Thixo-elastoviscoplastic modeling of human blood

Human blood is a complex biofluid that consists of plasma, an aqueous solution, in which blood cells are suspended: the red blood cells (RBCs), the white blood cells (WBCs), and the platelets. Plasma possesses viscoelastic properties [1,2], but it is often approximated as a Newtonian fluid. Since the RBCs are the most abundant blood cells, they are the ones that contribute most to the rheological properties of human blood. Accurate estimation of these properties is crucial for disease detection and prevention [3]. The steady shear viscosity of blood is the most common material function that is utilized for this purpose because it is significantly affected by physiological parameters, mainly the hematocrit. Blood viscosity below physiological values (3– $4 mPa s$ at 37 °C) can be related to a reduced number of RBCs and, consequently, reduced transportation of oxygen and nutrients, which is the case in anemia disease. Since it also affects blood pressure, other diseases such as hypertension, diabetes, and obesity are connected to blood viscosity as well [4]. Because of their elastic membrane, the RBCs are flexible and highly deformable, possessing elastic properties. Moreover, at small rates of strain $( γ ˙)$⁠, they can adhere to each other, creating networks that are called rouleaux. However, at large values of $γ ˙$⁠, these networks break down, the structure is destroyed, and the RBCs are suspended independently in plasma [5,6]. The deformability of RBCs is important in the microcirculation, where their normal size can be smaller than the diameter of the capillaries. Because of the ability of the RBCs to form rouleaux structures, human blood behaves as an elastic solid under small stress. Thus, it also exhibits a yield stress. It should be noted that human blood is not an ideal yield stress (viscoplastic) fluid because its yield stress is not constant [7]. The solid-fluid transition is a continuous process, rather than an abrupt change, and depends on the history of deformation of the material. Accurate estimation of the yield stress is important for predicting blood clotting. Moreover, human blood is a thixotropic material whose microstructure changes in time [8] and, consequently, so do its rheological properties. This conclusion stems from the fact that aggregation and breakdown of rouleaux are time-dependent, and rather slow, processes in contrast to viscoelasticity, which is apparent at short time scales. In this way, there is a clear distinction between thixotropy and viscoelasticity. As a result, at the beginning of an experiment, the response of human blood is expected to be elastic, while at longer times, thixotropy is dominant and the material breaks down or rebuilds itself depending on its previous state [9–11].

Blood rheology has been studied for over a century both experimentally and theoretically with the development of constitutive models [12]. The well-known Poiseuille law was derived in 1840s to evaluate the pressure drop in the capillary flow by assuming that blood is a Newtonian fluid. The Newtonian fluid approximation has been used to simulate the blood flow in arteries [13,14], but this is a major simplification. During the 20th century, it was observed that human blood is a shear-thinning fluid [15–19]; its viscosity decreases when the rate of strain $( γ ˙)$ increases. It was also observed that the viscosity depends on the hematocrit [16,19,20]. Regarding the yield stress, it was initially estimated by extrapolation of the flow curve to zero $γ ˙$ [17,21], and models such as the Casson model were used for the rheological description of human blood. The yield stress was also evaluated experimentally [17,22,23] to be of the order of $O( 10 − 3 Pa)$ although deviations were significant because of different methodologies and variable physiological characteristics [22–25]. The viscoelastic nature of human blood was first observed during clotting [26] and later in oscillatory flow experiments in cylindrical tubes [27]. Since then, several experimental works have provided insight into the complex nature of human blood [25,28–35], while also introducing the concept of thixotropy. Clearly, human blood possesses non-Newtonian characteristics that cannot be ignored. Several models have been proposed to describe the rheology of human blood. Some of them are empirical formulas that correlate the shear viscosity with the concentration of particles based on Einstein's work [36]. However, these models only give a simplified estimation and do not account for flow dynamics. Moreover, they do not account for viscoelasticity. The shear-thinning behavior of blood has been treated with models that are commonly used in industrial applications, such as the power law, Carreau–Yasuda, and Cross models [36,37]. Some of these models can predict a yield stress value as well. The Heschel–Bulkley model is a common ideal viscoplastic model, while the Casson model has been used extensively in the case of human blood to estimate its yield stress and viscosity [17,21] and provide empirical expressions in terms of the hematocrit and fibrinogen concentration [38]. To predict the transient response of human blood, viscoelasticity must be considered. Generalized Oldroyd-B models have been used for this purpose [36,39], while the Giesekus model has been used to describe the viscoelasticity of blood plasma [2]. These models are also used as the basis for recently developed thixotropic models [40,41]. This is done by introducing a scalar parameter that accounts for the structure of the material. This structural parameter varies between [0, 1] and is governed by an evolution equation that resembles chemical kinetics. It includes rebuilding and destruction terms that describe the macroscopic response of the material to flow changes. Such equations are commonly used to model thixotropic materials, and numerous models have been developed with this methodology [35,42–49]. Detailed reviews of such models can be found in [11,50–52].

Regarding our previous TEVP models [53,54], they are both based on the elastoviscoplastic formalism by Saramito [55] with additional dependence of the material properties on the structure parameter, λ. Our first model [54] considers only the rate of strain dependence of λ, for both the rebuilding and the breakdown terms. Moreover, it encompasses the concept of kinematic hardening to account for back stresses upon flow reversal. The model consists of 15 parameters. Asymmetric behavior upon reversal of the flow has not been observed experimentally for human blood. For this reason, our second model [53], which is specific to blood, does not account for kinematic hardening. This is a significant simplification because of the use of four less parameters. Finally, in this model, the structure equation is purely stress controlled. A common feature of all TEVP models is the need for introducing several adjustable parameters (usually more than ten) to describe the different mechanisms in the micro-scale that induce the complex response of the material under different flow conditions. One of the goals of this work was to reduce the number of parameters to make the model as simple and computationally friendly as possible, without reducing its accuracy.

In this work, we introduce a structural parameter $λ$ and a corresponding evolution equation, which is based on experimental findings [5] and thermodynamic principles developed by Stephanou and Georgiou [56]. Regarding the constitutive equation for the stresses, it is based on the elastoviscoplastic (EVP) formalism by Saramito [55]. Additionally, we have replaced the viscoplastic max term with a continuous function of $λ$⁠. For simplicity, we consider that human blood plasma is a Newtonian fluid, and we treat it separately by introducing a constant solvent viscosity. The proposed model consists of nine adjustable parameters (TEVP 9), and it is an improvement of our previous thixo-elastoviscoplastic (TEVP) models [53,54]. The accuracy of the model is tested in numerous prototype flow conditions. This paper is organized as follows: in Sec. II, we discuss the formulation of the new model, the constitutive equation for the stresses in the tensorial form, and the evolution equation for the structural parameter. In Sec. III, we present the fitting procedure to experimental data on human blood samples and the resulting values of the adjustable parameters. We also fit our previous TEVP model for human blood (TEVP 11) [53] to the same experimental data. In Sec. IV, we present the predictions of the model in comparison to experiments and to TEVP 11 predictions. Several flow conditions are examined: steady shear, shear step-up/step-down, triangle, and oscillatory experiments. Predictions for uniaxial elongation are provided as well. We conclude and discuss ideas for future work in Sec. V.

There are three contributions in this equation. The first one is a restructuring term due to Brownian collisions of the cells and is manifested by the constant $k 1$⁠. The second one is also a rebuilding term that accounts for aggregation due to shear. As stated previously, the RBCs aggregate at small $γ ˙$ and form rouleaux. However, this phenomenon quickly diminishes as $γ ˙$ increases. This Gaussian bell-like function is a simple expression that describes this process, and $t R$ is a characteristic time for rouleaux formation. It is based on the experimental data [5] and analysis [63]. The function is presented in Fig. 1 for $t R=0.395 s$⁠, which is the value we obtained from the fitting procedure (Sec. III).

The last term in the equation is the breakdown term, which we expect to be dominant for moderate and large $γ ˙$⁠. This term is consistent with nonequilibrium thermodynamics following the formulation developed by Stephanou and Georgiou [56]. The double dot product $( ∇ _ u _ ) T: C _ _$ resembles the dissipation term $σ _ _: ∇ _ u _$ in the kinetic energy equation and accounts for the work done by viscoelastic stresses. $k 2$ and $m 2$ are additional adjustable parameters and $C _ _$ is the conformation tensor.

Finally, we introduce a constant solvent viscosity $η s=1.2× 10 − 3 Pa s$ to account for a Newtonian contribution by blood plasma. This value is typical for healthy human blood plasma [64]. Adding $η s+ η 0 , R B C$ results in the total zero-rate viscosity of nonaggregating blood, $η 0$⁠, considering the plasma contribution as well. This approximation is valid for flows inside almost all types of human vessels, except for capillaries and arterioles, where the shear rate increases dramatically, and the viscoelastic nature of the human blood plasma becomes important. Still, these extreme cases are out of the scope of our study. Overall, the model consists of nine adjustable parameters: $G 0, n g, ε P T T, η 0 , R B C, m 1, k 1, t R, k 2, and m 2$⁠. This is a considerable reduction compared to previous TEVP models that required the evaluation of 11–15 parameters. Consequently, this model is significantly simpler and computationally cheaper, while preserving or even improving the accuracy of previous models in predicting experimental results, as we will present in Sec. IV. In  Appendix C, we provide additional information regarding the effect of each parameter on the model predictions.

The parameters of the model are estimated via nonlinear regression on experimental data by Armstrong and Pincot [41]. In their work, five samples from healthy donors were prepared following standard procedures and protocols. The samples underwent numerous experiments: steady shear flow, shear rate step-up/step-down, hysteresis loops, and oscillatory flow. In this work, we will present the results of the fitting procedure using the experimental data from the sample of donor 5 for which the physiological parameters are given in Table I.

For each experiment, the constitutive equations for the stress components and $λ$ [Eqs. (8) and (10), respectively] are solved iteratively in time using the Runge–Kutta method (RK4) with adjustable time step, depending on the flow conditions. We also performed the calculations with the implicit Euler method, and the results were almost identical. We fit all experimental data for all flow conditions simultaneously to obtain a single set of parameters that best describes the rheology of human blood. The values of the TEVP 9 model parameters that resulted from the fitting procedure are given in Table II.

It is interesting to note that for $t R=0.395 s$⁠, we predict that the maximum of the Gaussian function (Fig. 1) is approximately at $γ ˙≈1.8 s − 1$⁠. This value is close to the value at which maximum rouleaux formation was observed experimentally $( γ ˙≈2.0 s − 1)$ [5]. However, we predict a more rapid decay for larger $γ ˙$ compared to the experiments. For a more accurate description of the aggregation process, we would need to introduce additional adjustable parameters.

To compare the new model with our previous formulation (TEVP 11), we perform the same procedure to evaluate the parameters of that model using the same experimental data. The resulting values are given in Table III:

In this section, we compare the predictions of the new model (TEVP 9) with experimental data and TEVP 11 predictions. The comparison involves the experimentally conducted flow conditions: steady shear, shear step-up/step-down, triangle, and oscillatory shear. Wherever we provide additional predictions, they always refer to TEVP 9.

First, we present the predictions in the steady shear flow. The equations are solved iteratively in time until steady state, which we assume is reached when the difference between the previous and current solution is below a threshold of $10 − 8$⁠. An initial condition of fully structured material $(λ=1)$ at rest $( σ i j=0)$ is imposed.

In Fig. 2, we present the predictions of the new model (red line) compared with the experimental data (black symbols) in a shear stress–shear rate plot. We also compare it with TEVP 11 (green dashed).

we obtain a similar result $σ y=4.14× 10 − 3 Pa$⁠. A more accurate comparison is obtained from the experiments of Picart et al. [23] and Chien et al. [19]. For $Ht=43.1%$⁠, the value of the yield stress in these works is approximately $σ y=2.2× 10 − 3 Pa$⁠, which is very close to our predicted result.

The structure parameter $λ$ is plotted with respect to $γ ˙$ in Fig. 3(i) for TEVP 9. At small shear rate, $λ$ is close to 1 (⁠ $λ=0.983$ for $γ ˙= 10 − 3 s − 1$⁠), and the material is almost structured. As $γ ˙$ increases, rouleaux breaks down and $λ$ decreases monotonically to 0.019 at $γ ˙= 10 3 s − 1$⁠.

Regarding normal stresses, no experimental information is provided. However, we present the predictions of the new model for the first normal stress difference $N 1= σ x x− σ y y$ in Fig. 3(ii). The Saramito model and all its variants predict zero $σ y y$ and $σ z z$⁠, and so does this model. Consequently, $N 1≡ σ x x$⁠, while the second normal stress difference is zero $( N 2= σ y y− σ z z=0)$⁠. We observe that the normal stress difference increases with $γ ˙$ as we would expect for a viscoelastic material. A normal stress plateau is also evident at low $γ ˙$⁠, which is two orders of magnitude smaller than its shear counterpart. Similar to the shear stress, an analytical expression for the normal stress can be derived as shown in  Appendix A. This small value does not contribute much to the von Mises criterion, and primarily the shear stress determines the yield stress. This response is expected when the yield strain $ε y= σ y G o$ is small [65], which is true in this case $( ε y ≈ 0.01)$⁠. The normal stress plateau is predicted by Saramito-type models [58] and has been reported before experimentally in various yield stress fluids [66,67]. Interestingly, $N 1$ increases faster than $σ x y$⁠, and at $γ ˙∼40 s − 1$⁠, it becomes greater than $σ x y$⁠. The significant normal stress value can be attributed mainly to the deformability and viscoelasticity of RBCs because the rouleaux structures are destroyed in this range of shear rates.

Regarding the new model (TEVP 9) in Figs. 2 and 3, we can clearly see that at approximately $γ ˙≈2 s − 1$⁠, there is a change in the slope of all curves. As we can see in Fig. 4, in this range of $γ ˙$ (between the gray vertical lines), the rebuilding term due to flow decays to zero very fast with $γ ˙$⁠, while the destruction term increases exponentially. Consequently, the structural parameter $λ$ decreases rapidly and the same holds for the viscosity of the material. The reduction in the viscosity implicitly affects the stress $( σ x y=η γ ˙)$ as we observe in Figs. 2 and 3(ii). Also, it is possible that the deviation of the flow curve from the experimental data is caused by the cancelation of these two rapidly changing terms, which is difficult to handle via the simultaneous fit of all parameters.

As we mentioned in Sec. III, we correctly predict the maximum of the aggregation term. However, without additional parameters, the term decays faster compared to the experiments [5]. We suspect that adding an extra adjustable parameter would mitigate this problem and eliminate the peculiar change of slope.

In shear step experiments, the sample is initially (for $t<0$⁠) subjected to a constant shear rate $γ ˙ i n$ for sufficient time to achieve steady state. At this steady state, the values of the stress components and $λ$ are used as initial conditions. Then, at time $t=0$⁠, $γ ˙$ changes abruptly to a new constant value, $γ ˙ f$⁠. In the step-up case, the final $γ ˙ f$ is larger, while in the step-down case, $γ ˙ f$ is smaller. This experiment is more appropriate for viscoelastic and thixotropic materials because the flow is time dependent.

In Fig. 5, we present the predictions of TEVP 9 (red line) and TEVP 11 (green dashed) models in comparison to experiments (black symbols) for three different step-up experiments. A stress overshoot is present in all cases, which is characteristic of the elastic and thixotropic nature of the material. The deformability of the rouleaux and the RBCs is responsible for this viscoelastic response. The ability of RBCs to aggregate or disaggregate contributes to the thixotropic response. Initially, the material resists this change and stores energy elastically. However, at longer times, the material cannot resist anymore, and the dissipative nature prevails. When $γ ˙ f$ is larger, the steady shear stress is higher $( σ x y=η γ ˙)$⁠, the overshoot occurs earlier, and the steady state is attained earlier. The response is faster at larger $γ ˙ f$ because the material is less structured and the relaxation time $( t r e l= η t/ G t)$ is smaller. Both models can accurately predict the evolution of $σ x y$ in time, but TEVP 9 performs slightly better in predicting the overshoot.

The prediction of the new model for the normal stress difference $N 1$ and the structural parameter $λ$ is presented in Fig. 6 for case $γ ˙ i n=0.1 s − 1$ to $γ ˙ f=20.0 s − 1$⁠. Similar results are obtained for the rest of the step-up cases. Regarding $N 1$⁠, the behavior is similar to that of the shear stress: an overshoot is present at approximately the same time $(t≈0.08 s )$⁠. Throughout the experiment, $N 1$ is smaller than $σ x y$⁠. In Fig. 6(ii), we can see that the structural parameter $λ$ decreases rapidly when we increase $γ ˙$ from $0.1$ to $20.0 s − 1$⁠, demonstrating the breakdown of rouleaux structures.

The predictions of the models for six step-down experiments are presented in Fig. 7. A characteristic of thixotropy is the presence of a stress undershoot in step-down experiments. As $γ ˙$ decreases to $γ ˙ f$⁠, the shear stress decreases. For pure viscoelastic materials, the decrease is exponentially monotonic until the new steady state. However, thixotropic materials partially rebuild their structure to a new equilibrium because $γ ˙$ has decreased. As a result, the material properties also change and particularly the viscosity increases as the material becomes more structured. The increase in the viscosity will lead to an increase in the shear stress $( σ x y=η γ ˙)$ at longer times. We can see that TEVP 9 is able to predict the experimental data accurately with a slight underestimation of the undershoot time. On the other hand, TEVP 11 fails to predict the stress undershoot in cases of small $γ ˙ f$ [Figs. 7(i) and 7(ii)] because the rebuilding terms are negligible at $γ ˙ f=0.1 s − 1$⁠. It is clear that the rebuilding term in TEVP 9 (which is responsible for the stress undershoot) is physically more accurate compared to the corresponding one in TEVP 11. In all cases, the new model predicts the overshoot (in step-up) or the undershoot (in step-down) earlier compared to the experiment. For a more accurate prediction, we would need to introduce multiple relaxation modes or a multi- $λ$ approach as in [43]. However, these methods would introduce additional adjustable parameters and would make the model unnecessarily complicated.

In Fig. 8, we present the predictions of TEVP 9 for $N 1$ and $λ$ for the case of $γ ˙ i n=300.0 s − 1$ to $γ ˙ f=20.0 s − 1$⁠. The evolution of the normal stress difference is similar to $σ x y$⁠. The undershoot occurs slightly later compared to the shear stress undershoot. Initially, $N 1$ is larger than $σ x y$ because $γ ˙$ is large. As $γ ˙$ decreases, $N 1$ decreases faster and eventually becomes smaller than $σ x y$ at steady state. In contrast to the step-up experiment, the structural parameter in the step-down experiment increases as the material rebuilds its structure under the imposition of smaller $γ ˙$⁠. It is worth noting that at steady state $λ=0.088$⁠, which is the same as the steady state value in Fig. 6(ii), which corresponds to the step-up $γ ˙ i n=0.1 s − 1$ to $γ ˙ f=20.0 s − 1$⁠. This means that at steady state, the value of $λ$ is the same if $γ ˙$ is the same, irrespectively of the previous state of the material. This result is expected because the breakdown of rouleaux is a reversible process.

In this new definition, $σ x y , s t$ is the shear stress value at the new steady state, $σ x y , m$ is the shear stress value at the undershoot, and $t s t$ and $t m$ are the corresponding time instants. Equation (15) is a normalized area in stress-time units, of the undershoot, which is consistent with our perception of thixotropy: the more $σ x y , s t$ increases because of rebuilding and the longer it takes to reach steady state, the more thixotropic the material is. When there is no undershoot: $σ x y , s t= σ x y , m$ and $t s t= t m$⁠, there is no rebuilding and $ξ=0$⁠. As the difference $σ x y , s t− σ x y , m$ (or $t s t− t m$⁠) increases, thixotropy becomes more important and $ξ$ increases. We evaluate $σ x y , s t$ at the point where $| d σ x y d t |<ε$⁠, with $ε= 10 − 4$⁠.

In Fig. 9, we present the value of $ξ$ for different values of $γ ˙$ for both models. In Fig. 9(i), we examine the effect of $γ ˙ i$ by keeping constant $γ ˙ f=0.1 s − 1$⁠. Thus, the value of $σ x y , s t$ is the same in these cases. We observe that as we increase $γ ˙ i$⁠, $ξ$ also increases. When the difference $γ ˙ i− γ ˙ f$ is large, the material undergoes a more abrupt change. As a result, the elastic response and, consequently, the stress undershoot are more prominent, i.e., $σ x y , m$ is smaller. For that reason, the thixotropic index is larger for larger $γ ˙ i$⁠. In Fig. 9(ii), we examine the effect of $γ ˙ f$ by keeping constant $γ ˙ i=20.0 s − 1$⁠. In this case, $ξ$ decreases with $γ ˙ f$⁠. For large $γ ˙ f$ the rebuilding term due to shear is negligible and the only contribution is due to Brownian motion. As a result, the thixotropic response is weakened compared to cases of small $γ ˙ f$⁠. Thus, we expect the thixotropic index to be large at small values of $γ ˙ f$⁠. In both cases, we note that the old model underestimates the importance of thixotropic rebuilding due to shear, while the new model performs much better.

The time evolution of $σ x y$ for four triangular experiments is shown in Fig. 10. We observe that as the parameters of the flow increase, i.e., $α$ and $t m$⁠, the shear stress increases because $γ ˙$ is ultimately higher. However, in all cases, the maximum value of $γ ˙$ is small (⁠ $γ ˙ m a x=1.2 s − 1$ in the case of $α=0.0448 s − 2, t m=26.8 s$⁠). Thus, the normal stresses are negligible. The new model performs well with minor underestimation of the shear stress. However, TEVP 11 clearly underestimates the shear stress in all cases.

The triangle experiment can provide insight into the viscoelastic and thixotropic character of the material under different flow conditions [29,34,36]. To observe this, we need to construct the graph of $σ x y− γ ˙$⁠, which visualizes better the hysteresis loop. This graph is presented in Fig. 11 for TEVP 9 and three different cases of a and $t m$⁠. Each graph can be divided into two curves: the ramp-up curve, where $γ ˙$ increases, and the ramp-down curve where $γ ˙$ decreases. The direction is denoted with arrows. The relative position of the ramp-up and ramp-down curves is an indicator of the rejuvenation and breakdown of the material. The steady flow curve is also presented with dashed lines for comparison. The first graph corresponds to one of the previous cases $(a=0.043 s − 2, t m=23.8 s )$⁠. We note that the up and down curves almost coincide with each other and with the steady state flow curve for most of the experiment. Consequently, the structure reaches its new equilibrium quickly whether it is breaking down (ramp-up) or rebuilding (ramp-down). However, at small $γ ˙$⁠, there is an overshoot in the ramp-up curve due to elasticity, which is not replicated in the ramp-down curve. It should be noted that the model slightly underpredicts ramp-up stress, which means that the ramp-up and ramp-down curves should be a bit more distant from each other. The response is different if we decrease $t m$ as we can see in Fig. 11(ii) for $a=0.043 s − 2, t m=5.0 s$⁠. In this case, the ramp-down curve is below the ramp-up curve, while the steady state flow curve lies between them. For a particular $γ ˙$⁠, during the ramp-up, the material breaks down, but the time scale for breakdown is large enough so that the material remains more structured compared to its equilibrium state. As a result, the viscosity and the shear stress are higher compared to the steady state. On the other hand, during ramp-down, the material is slowly rebuilding and remains more unstructured compared to its steady state. In both cases, the material needs more time to adapt to the flow conditions. Elasticity is also important since an overshoot is evident in the ramp-up curve. In the third case $(a=0.002 s − 2, t m=5.0 s )$⁠, the response is completely different. The ramp-down curve is higher than the ramp-up curve. This is a characteristic of dominant viscoelasticity [29] as the material accumulates stress close to its yielding point. These types of hysteresis loops have been observed experimentally [29,34] for different blood samples. We are not able to compare our results with these experiments because the samples had different physiological parameters. However, the order of magnitude of the stress is the same at similar values of $γ ˙$⁠, and we are confident that by adjusting our model to these samples, it can predict this response accurately.

Regarding the structural variable, the predictions of the model are similar to those of our previous formulation (TEVP 11) for which the corresponding figures can be found in [53]. As we expect, $λ$ decreases during the ramp-up part and increases during ramp-down.

In SAOS, the sample usually undergoes an amplitude sweep first to identify the LVE region. In this series of experiments, the frequency of oscillation is constant, while the amplitude is increased. In the LVE region, both $G ′$ and $G ′ ′$ are constant, independent of $γ 0$⁠. In the experiments of Armstrong and Pincot [41], the frequency of oscillation is $ω=4π rad s − 1$⁠, and the amplitude varies between $γ 0=0.10$ and $γ 0=29.86$⁠. In Fig. 12, we present the experimental values and the predictions of both models for $G ′$ and $G ′ ′$⁠. Experimentally, the two moduli are practically constant up to $γ 0=0.30$⁠, which we consider to be the limit of the LVE region. Beyond this value, the two moduli decrease with $γ 0$⁠, especially $G ′$⁠. For all strain amplitudes, $G ′ ′$ is larger than $G ′$⁠, which means that the viscous nature is dominant compared to the elastic one. The models can accurately predict the viscous response of human blood as the values of $G ′ ′$ are close to the experimental ones. Unfortunately, this is not true for $G ′$ because the models overpredict its value and cannot predict the LVE region. However, our new model performs much better compared to TEVP 11. Saramito-type models always predict $G ′> G ′ ′$ in the proximity of yielding, which is not desirable for human blood as we can observe from the experimental values. Clearly, we still lack a physical explanation of the yielding process that can be translated to an accurate mathematical expression.

A common deficiency of the Saramito model and its variants is that they predict $G ′= G o$ and $G ′ ′=0$ as $γ 0→0$⁠, in the absence of solvent viscosity. However, we would physically expect that a gel is characterized by viscous dissipation even in the fully structured solid state, which is expressed by the finite value of $G ′ ′$ in numerous experimental works [58,70–73]. In Fig. 13, we can observe that our new model can predict this response because of the introduction of the plasma (solvent) viscosity. The plateau of $G ′ ′≈0.015 Pa$ is reached smoothly because of the continuous yielding description through $λ$ [Eq. (8)]. This is much more realistic compared to the abrupt yielding that a max term would predict (see also Fig. 14 for the prediction of TEVP 11 at small $γ 0$⁠).

When we investigate elastoviscoplastic systems [54], it is common to exclude the viscous contribution of the solvent. However, doing so results in zero $G ′ ′$ for small $γ 0$ if the constitutive model contains the max term for viscoplasticity. On the contrary, our newly proposed model does not incorporate viscoplasticity through the max term, which has a major implication in the calculated $G ′ ′$⁠. In Fig. 14, we compare the behavior of the two models as $γ 0→0$⁠. In this figure, we do not account for the solvent viscosity both in TEVP 9 [Fig. 14(i)] and in TEVP 11 [Fig. 14(ii)]. We do this to compare the yielding predictions of the max term (found in TEVP 11) and the continuous yielding (found in TEVP 9). We also take into consideration the integration time step, $dt$⁠, of Eqs. (18) and (19) because it is important at this range of amplitude. We should note that we perform the integration using the trapezoidal rule. The arrows indicate the direction of decreasing integration time step. We see that for both models when we decrease $dt$⁠, $G ′ ′$ tends to zero. On the right, we see that our old model (TEVP 11) predicts $G ′ ′=0$ at approximately $γ 0≈0.015$ (provided that we decrease $dt$ very much). Therefore, it is clear that there is an abrupt change in $G ′ ′$ because of the max term that predicts sharp yielding of the material. Thus, this model suffers from the same deficiency as any Saramito-type model that contains the max term. However, on the left, we can see that our new model has a different response with decreasing $γ 0$⁠: $G ′ ′$ reaches a similar plateau as $γ 0→0$⁠, but its decrease is smooth. Moreover, it is important to state that $G ′ ′$ remains finite for very small values of $γ 0$⁠. We can observe that in the case of the smallest $dt$ (i.e., the most accurate prediction) where at $γ 0= 10 − 4$⁠, we have not reached the plateau $G ′ ′=0$ yet. Therefore, the response of the new model is much improved in predicting the viscous dissipation of a gel material because we describe the yielding process with a continuous function of the structure parameter $λ$⁠, rather than with a fixed value of yield stress and sudden yielding.

After the amplitude sweep test, the sample undergoes a frequency sweep test, where the amplitude of oscillation is constant, and the frequency varies. The frequency sweep is usually conducted at an amplitude that corresponds to the LVE region of the material, which has been identified from the amplitude sweep. However, in the experimental work of Armstrong and Pincot [41], the frequency sweep is performed at $γ 0=9.95$⁠, which is clearly in the non-LVE region (Fig. 12). The frequency range is $ω=0.5−40.0 rad s − 1$⁠. We plot the results for $G ′$ and $G ′ ′$ in Fig. 15. $G ′ ′$ is larger than $G ′$ through the whole range of frequencies, which means that human blood is mostly viscous. Again, the viscous nature is accurately predicted, but deviations are evident for $G ′$⁠. However, TEVP 9 is significantly more accurate than TEVP 11.

In  Appendix B, we provide additional predictions of TEVP 9 for the structure and stresses during a period of oscillation.

In the large amplitude oscillatory shear (LAOS) flow, the kinematics is similar to SAOS, but in this case, we are clearly in the non-LVE region. Consequently, the definitions of $G ′$ and $G ′ ′$ are not sufficient anymore because higher harmonics contribute significantly to the response of the material. Usually, the Lissajous–Bowditch diagrams are constructed, which are three-dimensional diagrams of the shear stress with respect to the strain and the rate of strain in the steady oscillatory state (alternance). The stress–strain plot is the elastic projection, and the stress-rate of the strain plot is the viscous projection.

The elastic projection is presented in Fig. 16 for three different frequency values at strain amplitudes $γ 0=0.5−100.0$⁠. The elastic projection is a circular loop in most experiments, indicating that the viscous nature is dominant compared to the elastic. The new model (red lines) is accurate for most experiments. Deviations arise for the highest frequency $(ω=5.0 rad s − 1)$ and for small strain amplitudes (⁠ $γ 0=0.5$ and $γ 0=1.0$⁠). For small $γ 0$⁠, the response is mostly elastic, and the model is not so accurate as we discussed in the SAOS test. For large $γ 0$⁠, the material is more viscous, and the model performs very well. The predictions of TEVP 11 (green lines) are similar and slightly better at large frequency and small amplitude.

The viscous projection is presented in Fig. 17 for the same range of frequency and amplitude. Again, TEVP 9 (red lines) slightly deviates from the experimental data at high frequency. The viscous projection loops are more linear, especially for large $γ 0$⁠, verifying the dominant viscous nature of human blood in this regime. Additionally, secondary loops are present in some cases (for example, at $ω=5.0 rad s − 1$ and $γ 0=1.0$ or at $ω=1.0 rad s − 1$ and $γ 0=10.0$⁠). These secondary loops have been reported in the past in viscoelastic materials [74] and are related to stress overshoots and nonlinear elasticity and thixotropy due to rebuilding and breakdown of the structure throughout a period of oscillation [44,70,75]. The predictions of TEVP 11 (green lines) are similar.

Having validated the accuracy of the model in shear flows, it is also important to test its predictions in extensional flows. Experimentally, little is known about the response of complex materials in extensional flows, and human blood is no exception [76]. Since we do not have experimental data to compare with, we present the predictions of the model for the values of the parameters that were evaluated from the regression on shear flow data (Table I).

For small $ε ˙$⁠, the viscosity decreases (extensional thinning). The rouleaux breaks down and the structure disintegrates leading to a reduction in the resistance to the flow. At large $ε ˙$⁠, the viscosity increases and reaches a local maximum, which indicates extensional hardening. In this region, the normal stress difference increases fast. The resistance to flow increases, possibly because the RBCs try to adapt to the flow. For large $ε ˙$⁠, the material is thinning again constantly, because the structure is destroyed, and the RBCs are aligned with the flow. The model is able to predict bounded extensional properties because of the PTT term that also accounts for thinning at very large $ε ˙$⁠. However, we do not have any indication as to whether these results are realistic specifically for human blood; related experiments are necessary.

The structural parameter as a function of $ε ˙$ is presented in Fig. 18(ii). The curve is similar to its shear flow counterpart. For small $ε ˙$⁠, $λ$ is close to unity and the material is structured. As $ε ˙$ increases, $λ$ decreases abruptly and approaches (but never reaches) zero. In this state, the material is unstructured, which is exactly what is expected.

The changes in the slope in Figs. 18(i) and 18(ii) are similar to those in Figs. 2 and 3. The interplay between the rebuilding and the destruction terms is responsible for this prediction as we discussed in Sec. IV A.

In realistic extensional flows, it is common that a steady state cannot be achieved. Consequently, it is important to be able to predict the transient aspect of the flow. In Fig. 19, we present the predictions of the model for startup uniaxial elongation for three different elongational rates. In Fig. 19(i), we plot the evolution of the first normal stress difference for different values of $ε ˙$⁠. We can see that for the largest $ε ˙$⁠, the stress difference is larger, which is expected. A stress overshoot is present similar to the shear flow case (Fig. 5). In all cases, we achieve a steady state, and the larger the extensional rate, the sooner the steady state. In Fig. 19(ii), the extensional viscosity is presented. Initially we observe only the Newtonian contribution, and the extensional viscosity of blood is three times the shear viscosity of plasma. The Trouton ratio is 3, which is the expected value for a Newtonian fluid in uniaxial elongation. As the elastic stresses evolve, we observe an overshoot in the extensional viscosity for all values of $ε ˙$⁠. For the lowest value of $ε ˙$⁠, the viscosity is larger. However, comparing the cases of $ε ˙=8.0 s − 1$ and $ε ˙=50.0 s − 1$⁠, we can see that in the latter, the viscosity is higher, because for $ε ˙=50.0 s − 1$⁠, we are in the region of extensional hardening [Fig. 18(i)]. Finally, the time evolution of $λ$ is plotted in Fig. 19(ii). Initially $λ=1$ because the material is at rest. As time goes by, the material is deformed and $λ$ decreases to its steady state value, which is higher for the smallest $ε ˙$⁠, as expected.

In this work, we developed a new constitutive model to describe the thixotropic, viscoelastic, and viscoplastic nature of human blood. For simplicity, we consider that blood plasma is Newtonian, so we add a solvent contribution to the total stress. The model was based on our previous formulations [53,54] and on the work of Stephanou and Georgiou for a thermodynamically consistent approach of thixotropy [56]. We have excluded the max term of the viscoplastic rate of deformation and replaced it with a continuous function of the structure of the material based on experimental reports [7]. This approach resembles the continuous process of yielding more accurately. The model consists of nine parameters, which is a major improvement compared to previous models, both in terms of computational requirements and simplicity.

The model was tested in numerous shear experiments and was able to accurately predict all major rheological characteristics. In steady shear, it predicts finite apparent yield stress despite not explicitly containing a yield stress parameter. The yield stress was of the same order of magnitude as previous experimental reports. In step shear experiments, the model predicts stress overshoots (in step-up) and undershoots (in step-down), which indicate viscoelastic and thixotropic response. In triangle experiments, the shear stress prediction was slightly underestimated, but the overall agreement was satisfactory. Moreover, we were able to reproduce hysteresis loops that have been observed experimentally [29]. In SAOS, $G ′ ′$ was accurately predicted, but $G ′$ slightly deviated from the experimental values. A significant improvement of the model is the ability to predict finite viscous dissipation for small values of the amplitude, i.e., close to the solid region, even in the absence of a solvent. The inability of the Saramito model and its variants to predict this response stems from the presence of the max term. By replacing the viscoplastic term with a continuous function, we were able to circumvent this deficiency and predict smooth variation of the dynamic moduli with $γ o$⁠. The elastic and viscous projections of the Lissajous–Bowditch diagram agree with experiments as well. Discrepancies were found for large frequency and small amplitude of oscillation. In all oscillatory experiments, the viscous response is dominant compared to the elastic one, and we also predict secondary loops in the viscous projection, an indication of interplay between nonlinear viscoelasticity and thixotropy. Additionally, we provide predictions of the model in uniaxial elongation. The model predicts finite normal stress difference and bounded elongational viscosity because of the PTT term. For large elongational rate of strain, $ε ˙$⁠, a region of extensional hardening exists where the elongational viscosity increases with $ε ˙$⁠. In any other case, the viscosity decreases with $ε ˙$⁠, i.e., the material is thinning. The transient response under these flow conditions is similar to the shear counterpart (shear step-up experiment). We compare this model to our previous TEVP formulation for human blood (TEVP 11) [53] and observe that the predictions of the present study are more accurate, especially in step-down, triangular and SAOS experiments.

In future works, we wish to use this model in realistic simulations of human blood in vessels. We have performed such simulations using TEVP 11 [77–79], and we are confident that we can extend that work by utilizing the present model and its improved predictability. Regarding the model itself, we should focus on two main goals in future formulations. First, we aim to use less adjustable parameters, in order to make the model simpler and efficient for large-scale simulations. It is important to remove unnecessary or peculiar parameters or replace them with others that have a clearer physical meaning. Second, we should focus on predicting experimental data more accurately, especially in the case of SAOS experiments where we observed significant deviations in the elastic response, i.e., in the $G ′$ predictions. This is inherent in all Saramito models, and we showed that the exclusion of the max term and the definite yield stress parameter improved the predictions of $G ′ ′$ at small $γ$ but did not completely alleviate the problem. Consequently, in order to predict the linear viscoelastic region correctly, we need to depart from this formulation and adopt a different approach to constitutive modeling. Additional experimental work and modeling from a microscale point of view can provide insight into the mechanisms of rouleaux formation and overall structural response of human blood under complex flow conditions.

The research work was supported by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the Project “MOFLOWMAT” (No. HFRI-FM17-2309). P.M. was supported by HFRI under the 3rd Call for HFRI PhD Fellowships (Fellowship No. 5854) and Y.D. was supported by H.F.R.I. under the Project “CARE” (No. HFRI-FM17-876). We are very thankful to the anonymous reviewers for their fruitful comments and suggestions that led to the improvement of the presentation of our work.

The authors have no conflicts to disclose.