We propose an enhanced model for the rheological characterization of human blood that accounts for thixotropy, viscoelasticity, and yield-stress. Blood plasma is assumed to act as a Newtonian solvent. We introduce a scalar variable, λ, to macroscopically describe the structure of blood. The temporal evolution of λ is governed by an equation that accounts for aggregation of red blood cells and breakdown of rouleaux structures. We introduce a Gaussian function that qualitatively describes experimental findings on rouleaux restructuring and the expression that was proposed by Stephanou and Georgiou for the breakdown term. The constitutive equation for stresses is based on the elastoviscoplastic formalism by Saramito. However, the max term of the viscoplastic deformation rate has been replaced by a continuous function of λ to account for smooth solid-fluid transition, following the experimental evidence. The continuous yielding description provides improved rheological predictions, especially in small amplitude oscillatory shear. The model predicts finite viscous dissipation at small amplitude oscillation, as we would expect from a gel material-like human blood. Overall, it has nine adjustable parameters that are fitted simultaneously to experimental data by nonlinear regression. The model can accurately predict numerous flow conditions: steady shear, step shear, hysteresis loops, and oscillatory shear. We compare this model (TEVP 9) to our previous formulation for human blood (TEVP 11), and we show that the predictions of the new model are more accurate, despite using fewer parameters. We provide additional predictions for uniaxial elongation, which include finite normal stress difference, extensional hardening at large values of the extensional rate, and extensional thinning at extremely large extensional rates.

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.

Our model is based on the elastoviscoplastic formulation by Saramito [55]. First, we define the total rate of strain tensor γ ˙ _ _,
γ ˙ _ _ = ( _ u _ ) + ( _ u _ ) T ,
(1)
where u _ is the velocity and _ is the nabla differential operator. The total strain rate is decomposed in an elastic and a viscoplastic part [57–59],
γ ˙ _ _ = γ ˙ _ _ e + γ ˙ _ _ v p .
(2)
The elastic part of the rate of strain is defined as
γ ˙ _ _ e = 1 G t σ _ _ .
(3)
The definition of the elastic part is similar to the Saramito model in accounting for Neo-Hookean behavior in the solid region. However, the elastic modulus G t depends on thixotropy (denoted with the subscript “t”) in the following way:
G t = G 0 λ n g .
(4)
In essence, the elastic modulus depends on the structure of the material that is characterized by the scalar variable λ, so that the material elasticity decreases, when it gets less structured. G 0 and n g are adjustable parameters. σ _ _ is the upper convected derivative of the extra stress tensor σ _ _, which is defined as
σ _ _ = σ _ _ t + u _ _ σ _ _ ( _ u _ ) T σ _ _ σ _ _ ( _ u _ ) .
(5)
The viscoplastic part of the rate of strain is defined as
γ _ ˙ _ v p = e ε P T T t r ( σ _ _ ) G t η t σ _ _ .
(6)
An exponential Phan–Thien Tanner (PTT) function is used to account for shear thinning and to bound the extensional viscosity at large extensional rates of strain. η t is the plastic viscosity that also depends on thixotropy, through the structure parameter λ, so that the material viscosity also decreases, when it gets less structured, but with a different dependence on λ,
η t = λ m 1 1 λ η 0 , R B C .
(7)
m 1 and η 0 , R B C are adjustable parameters. More specifically, η 0 , R B C is the zero-rate viscosity in the case of nonaggregating blood (in which no yield stress is present, and there exists a viscosity plateau as γ ˙ 0), excluding the plasma contribution. The major difference of this model compared to the Saramito variants and previous TEVP formulations is the replacement of the max term of the viscoplastic part of the deformation with a continuous function of the structure parameter λ. Instead of explicitly considering a yield stress value, we propose that the viscosity approaches infinity as λ 1, where the material is fully structured and behaves as a solid. This proposition is physically more appropriate because experiments have shown that there is no apparent yield stress [7], but its value depends on the structure and the history of the material. All in all, the tensorial stress constitutive equation is
1 G 0 λ n g σ _ _ + e ε P T T t r ( σ _ _ ) G t 1 λ η 0 , R B C λ m 1 σ _ _ = [ ( _ u _ ) + ( _ u _ ) T ] .
(8)
The constitutive equation combines the Saramito model and a PTT function, both of which are thermodynamically consistent [55,60]. The equation can also be written in terms of the conformation tensor, C _ _ = σ _ _ G t + I _ _, which accounts for the shape and orientation of the RBCs and the rouleaux structures and their deviation from equilibrium conditions,
C _ _ + e ε P T T t r ( C _ _ I _ _ ) 1 λ η 0 , R B C λ m 1 G t ( C _ _ I _ _ ) = 0 _ _ .
(9)
Regarding the structure of the material, the scalar variable λ accounts for microstructural changes, macroscopically. We note that the use of a scalar variable is a simplification. At small shear rates, the orientation of complex 3D rouleaux structures can lead to anisotropy. Anisotropy can also arise at large γ ˙ because of the diverse orientations of RBCs [61,62]. The significance of the phenomenon should be assessed in future studies, albeit the introduction of a tensorial structure variable would complicate the formulation immensely. The evolution equation of λ resembles a chemical kinetics equation that accounts for aggregation of RBCs and breakdown of rouleaux,
d λ d t = k 1 ( 1 λ ) B r o w n i a n a g g r e g a t i o n + | γ ˙ | e ( t R | γ ˙ | ) 2 ( 1 λ ) F l o w a g g r e g a t i o n k 2 [ ( _ u _ ) T : C _ _ ] λ m 2 D e s t r u c t i o n .
(10)

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

FIG. 1.

Representation of the flow aggregation function for t R = 0.395 s.

FIG. 1.

Representation of the flow aggregation function for t R = 0.395 s.

Close modal

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.

TABLE I.

Physiological parameters for sample of donor 5; see [41].

DonorHt (%)Fibrinogen ( g / dl )Total cholesterol ( mg / ml )Triglycerides ( mg / ml )HDL cholesterol ( mg / ml )LDL cholesterol ( mg / ml )
43.1 0.249 244 94 70 154 
DonorHt (%)Fibrinogen ( g / dl )Total cholesterol ( mg / ml )Triglycerides ( mg / ml )HDL cholesterol ( mg / ml )LDL cholesterol ( mg / ml )
43.1 0.249 244 94 70 154 

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.

TABLE II.

Values of the TEVP 9 model parameters for the sample of donor 5; see [41].

G0 (Pa)ng ε P T Tη0,RBC (Pa s)m1k1(s−1)tR(s)k2m2
0.147 0.281 0.011 0.028 0.668 0.370 0.395 628.7 3.961 
G0 (Pa)ng ε P T Tη0,RBC (Pa s)m1k1(s−1)tR(s)k2m2
0.147 0.281 0.011 0.028 0.668 0.370 0.395 628.7 3.961 

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:

TABLE III.

Values of the TEVP 11 model parameters for the sample of donor 5; see [41].

G (Pa)σy (Pa) ε P T Tη0 (Pa s)m1k1 (s−1) k 2 ( s n 1 1 )k3n1n2n3
0.211 0.003 0.001 0.015 0.505 0.010 0.688 2329.3 0.312 3.388 1.430 
G (Pa)σy (Pa) ε P T Tη0 (Pa s)m1k1 (s−1) k 2 ( s n 1 1 )k3n1n2n3
0.211 0.003 0.001 0.015 0.505 0.010 0.688 2329.3 0.312 3.388 1.430 

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

FIG. 2.

Steady shear flow predictions: shear stress as a function of shear rate. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

FIG. 2.

Steady shear flow predictions: shear stress as a function of shear rate. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

Close modal
The new model is accurate at large γ ˙ and slightly overpredicts the shear stress at small γ ˙. In other works, a contribution of independent RBCs has been introduced to account for viscoelasticity at large shear rates [41,52]. However, in our formulation, this term is not needed because the model is accurate in this region. Regarding the discrepancy at small γ ˙, we found that adding more weight on this data region did not improve but rather led to less accurate results in some transient rheometric flows. The model can also predict an apparent yield stress value under shear σ y 1.6 × 10 3 Pa despite not explicitly including a yield stress parameter. An analytical expression for the apparent yield stress can be determined under steady state in the limit γ ˙ 0 0 and ε P T T 0 (see  Appendix A),
σ y = η 0 G 0 k 1 k 2 .
(11)
This expression results in σ y 1.55 × 10 3 Pa, which matches the simulation result since ε P T T is small. Regarding TEVP 11, it is accurate for the whole range of γ ˙ but overpredicts the yield stress value ( 3.2 × 10 3 Pa ). The experiments do not show a tendency to reach a yield stress because they were conducted in a rate-controlled manner, in which the shear rate is imposed, and the stress is measured. In this way, the material will always flow since γ ˙ > 0, and the stress in the solid region cannot be measured. In order to measure the yield stress in a rate-controlled manner, one should initially impose a constant shear rate (ramp-up), followed by its gradual decrease to zero (ramp-down). The residual stress of the material as γ ˙ 0 indicates the yield stress. In the case of human blood, it is well documented that it is a yield-stress fluid. Several equations have been proposed for the estimation of the yield stress depending on the hematocrit and the concentration of fibrinogen, like the one by Morris et al. [22],
σ y = ( 0.091 + 0.47 H t + 0.22 c f 0.14 c f 2 + 0.48 H t c f ) 2 .
(12)
Here, σ y is the yield stress in dyn c m 2, H t is the volume fraction of the RBCs (hematocrit), and c f is the concentration of fibrinogen in g d L 1. Using the values of the hematocrit and fibrinogen from Table I, we calculate σ y = 4.38 × 10 3 Pa, which is close to the prediction of both models but clearly higher than the experimental data in Fig. 2. If we use the newer expression by Apostolidis and Beris [49],
σ y = { ( H t H t c ) 2 ( 0.5084 c f + 0.4517 ) 2 , H t H t c 0 , H t H t c ,
(13)
with the critical hematocrit,
H t c = { 0.3126 c f 2 0.468 c f + 0.1764 , H t H t c 0.0012 , H t H t c ,
(14)

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 H t = 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.

FIG. 3.

Dependence on the shear rate of (i) the structure parameter λ and (ii) the normal stress N 1 (TEVP 9).

FIG. 3.

Dependence on the shear rate of (i) the structure parameter λ and (ii) the normal stress N 1 (TEVP 9).

Close modal

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.

FIG. 4.

Dependence of the rebuilding and destruction terms, defined in Eq. (10), on γ ˙ and their effect on λ.

FIG. 4.

Dependence of the rebuilding and destruction terms, defined in Eq. (10), on γ ˙ and their effect on λ.

Close modal

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.

FIG. 5.

Shear step-up predictions: time evolution of shear stress. (i) γ ˙ i n = 0.1 s 1 to γ ˙ f = 5.0 s 1, (ii) γ ˙ i n = 0.1 s 1 to γ ˙ f = 10.0 s 1, and (iii) γ ˙ i n = 0.1 s 1 to γ ˙ f = 20.0 s 1. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

FIG. 5.

Shear step-up predictions: time evolution of shear stress. (i) γ ˙ i n = 0.1 s 1 to γ ˙ f = 5.0 s 1, (ii) γ ˙ i n = 0.1 s 1 to γ ˙ f = 10.0 s 1, and (iii) γ ˙ i n = 0.1 s 1 to γ ˙ f = 20.0 s 1. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

Close modal

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.

FIG. 6.

Shear step-up predictions: time evolution of (i) normal stress difference and (ii) structural parameter λ for γ ˙ i n = 0.1 s 1 to γ ˙ f = 20.0 s 1 (TEVP 9).

FIG. 6.

Shear step-up predictions: time evolution of (i) normal stress difference and (ii) structural parameter λ for γ ˙ i n = 0.1 s 1 to γ ˙ f = 20.0 s 1 (TEVP 9).

Close modal

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.

FIG. 7.

Shear step-down predictions: time evolution of shear stress. (i) γ ˙ i n = 5.0 s 1 to γ ˙ f = 0.1 s 1, (ii) γ ˙ i n = 10.0 s 1 to γ ˙ f = 0.1 s 1, (iii) γ ˙ i n = 20.0 s 1 to γ ˙ f = 0.1 s 1, (iv) γ ˙ i n = 300.0 s 1 to γ ˙ f = 5.0 s 1, (v) γ ˙ i n = 300.0 s 1 to γ ˙ f = 10.0 s 1, and (vi) γ ˙ i n = 300.0 s 1 to γ ˙ f = 20.0 s 1. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

FIG. 7.

Shear step-down predictions: time evolution of shear stress. (i) γ ˙ i n = 5.0 s 1 to γ ˙ f = 0.1 s 1, (ii) γ ˙ i n = 10.0 s 1 to γ ˙ f = 0.1 s 1, (iii) γ ˙ i n = 20.0 s 1 to γ ˙ f = 0.1 s 1, (iv) γ ˙ i n = 300.0 s 1 to γ ˙ f = 5.0 s 1, (v) γ ˙ i n = 300.0 s 1 to γ ˙ f = 10.0 s 1, and (vi) γ ˙ i n = 300.0 s 1 to γ ˙ f = 20.0 s 1. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

Close modal

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.

FIG. 8.

Shear step-down predictions: time evolution of (i) normal stress difference and (ii) structural parameter λ for γ ˙ i n = 300.0 s 1 to γ ˙ f = 20.0 s 1 (TEVP 9).

FIG. 8.

Shear step-down predictions: time evolution of (i) normal stress difference and (ii) structural parameter λ for γ ˙ i n = 300.0 s 1 to γ ˙ f = 20.0 s 1 (TEVP 9).

Close modal
To quantify the importance of thixotropy in step-down experiments, we introduce a thixotropic index that accounts for the importance of structural rebuilding [53,68,69],
ξ = σ x y , s t σ x y , m σ x y , s t t s t t m t s t .
(15)

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.

FIG. 9.

Thixotropic index predictions for the step-down experiment (i) against γ ˙ i with constant γ ˙ f = 0.1 s 1 and (ii) against γ ˙ f with constant γ ˙ i = 20.0 s 1.

FIG. 9.

Thixotropic index predictions for the step-down experiment (i) against γ ˙ i with constant γ ˙ f = 0.1 s 1 and (ii) against γ ˙ f with constant γ ˙ i = 20.0 s 1.

Close modal
Another common test for thixotropic materials is the hysteresis loop or triangular experiment. In this case, the material is initially at rest ( λ = 1, σ i j = 0), and at time t = 0, we increase γ ˙ linearly in time (ramp-up). At the half-time of the experiment, where γ ˙ is maximum, we decrease it with the same rate to γ ˙ = 0 (ramp-down). The expression for the time evolution of γ ˙ is
γ ˙ ( t ) = { α t , t t m a ( 2 t m t ) , t > t m ,
(16)
where α is the rate of increase of γ ˙ and t m is the half-time of the experiment.

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.

FIG. 10.

Triangular experiment predictions: time evolution of shear stress. (i) a = 0.0375 s 2 , t m = 20.0 s, (ii) a = 0.0399 s 2 , t m = 21.9 s, (iii) a = 0.043 s 2 , t m = 23.8 s, and (iv) a = 0.0448 s 2 , t m = 26.8 s. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

FIG. 10.

Triangular experiment predictions: time evolution of shear stress. (i) a = 0.0375 s 2 , t m = 20.0 s, (ii) a = 0.0399 s 2 , t m = 21.9 s, (iii) a = 0.043 s 2 , t m = 23.8 s, and (iv) a = 0.0448 s 2 , t m = 26.8 s. Comparison of TEVP 9 (continuous line) and TEVP 11 (dash line) with experiments [41] (symbols).

Close modal

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.

FIG. 11.

Triangular experiment predictions: hysteresis loops. (i) a = 0.043 s 2 , t m = 23.8 s, (ii) a = 0.043 s 2 , t m = 5.0 s, and (iii) a = 0.002 s 2 , t m = 5.0 s (TEVP 9).

FIG. 11.

Triangular experiment predictions: hysteresis loops. (i) a = 0.043 s 2 , t m = 23.8 s, (ii) a = 0.043 s 2 , t m = 5.0 s, and (iii) a = 0.002 s 2 , t m = 5.0 s (TEVP 9).

Close modal

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.

The small amplitude oscillatory shear (SAOS) experiment is useful for identifying the linear viscoelastic (LVE) region and for examining the dynamic behavior of a material in alternance (steady oscillating state). In oscillatory shear experiments, the material undergoes deformation by imposing either a strain (SAOStrain) or a stress (SAOStress) that varies sinusoidally in time. In their experiments, Armstrong and Pincot conducted strain-controlled tests [41]. In this case, the applied strain is given by
γ ( t ) = γ 0 sin ( ω t ) ,
(17)
where γ 0 is the amplitude and ω is the frequency of oscillation. The storage and loss moduli can be computed from the Fourier expansion of the measured shear stress,
G = ω π γ 0 t t + 2 π ω σ x y sin ( ω t ) d t ,
(18)
G = ω π γ 0 t t + 2 π ω σ x y cos ( ω t ) d t .
(19)

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.

FIG. 12.

Amplitude sweep at ω = 4 π rad s 1 for γ 0 = 0.10 29.86: values of G and G . Comparison of TEVP 9 (line) and TEVP 11 (dashed) with experiments [41].

FIG. 12.

Amplitude sweep at ω = 4 π rad s 1 for γ 0 = 0.10 29.86: values of G and G . Comparison of TEVP 9 (line) and TEVP 11 (dashed) with experiments [41].

Close modal

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

FIG. 13.

Predictions of TEVP 9 at ω = 4 π rad s 1 and small γ 0: values of G and G . Experimental values [41] are shown with symbols.

FIG. 13.

Predictions of TEVP 9 at ω = 4 π rad s 1 and small γ 0: values of G and G . Experimental values [41] are shown with symbols.

Close modal
FIG. 14.

Amplitude sweep at ω = 4 π rad s 1 for γ 0 = 10 4 29.86: values of G and G . (i) TEVP 9 and (ii) TEVP 11. Experimental values [41] are shown with symbols.

FIG. 14.

Amplitude sweep at ω = 4 π rad s 1 for γ 0 = 10 4 29.86: values of G and G . (i) TEVP 9 and (ii) TEVP 11. Experimental values [41] are shown with symbols.

Close modal

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, d t, 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 d t, 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 d t 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 d t (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.

FIG. 15.

Frequency sweep at γ 0 = 9.95 for ω = 0.5 40.0 rad s 1: values of G and G . Comparison of TEVP 9 (line) and TEVP 11 (dashed) with experiments [41].

FIG. 15.

Frequency sweep at γ 0 = 9.95 for ω = 0.5 40.0 rad s 1: values of G and G . Comparison of TEVP 9 (line) and TEVP 11 (dashed) with experiments [41].

Close modal

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.

FIG. 16.

Elastic projection: Stress-strain curves for different values of the frequency and the amplitude of oscillation. Black: experimental values [41], red: TEVP 9; green: TEVP 11.

FIG. 16.

Elastic projection: Stress-strain curves for different values of the frequency and the amplitude of oscillation. Black: experimental values [41], red: TEVP 9; green: TEVP 11.

Close modal

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.

FIG. 17.

Viscous projection: Stress-strain rate curves for different values of the frequency and the amplitude of oscillation. Black: experimental values [41], red: TEVP 9, and green: TEVP 11.

FIG. 17.

Viscous projection: Stress-strain rate curves for different values of the frequency and the amplitude of oscillation. Black: experimental values [41], red: TEVP 9, and green: TEVP 11.

Close modal

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

In Fig. 18(i), we present the first normal stress difference σ z z σ x x as a function of the extensional rate ε ˙. The model predicts an apparent extensional yield stress of approximately 2.7 × 10 3 Pa, which is 3 times the shear yield stress ( σ y = 1.6 × 10 3 Pa ). This value is exactly what the ideal viscoplastic theory predicts in uniaxial elongation. Varchanis et al. [65] showed that for an elastoviscoplastic material the value is the same when the yield strain ε y = σ y / G is smaller than 0.2. In our case, the yield strain is defined by using the value of G at λ = 1, i.e., G 0 = 0.147 Pa, so that ε y 0.01. Thus, the model is consistent with the predictions of elastoviscoplastic models. At large ε ˙, the stress difference increases abruptly, while at even larger ε ˙, the stress difference increases with a smaller slope. In Fig. 18(i), we also present the first extensional viscosity that is defined as
η 1 , e = σ z z σ x x ε ˙ .
(20)
FIG. 18.

Uniaxial elongation: (i) First normal stress difference and extensional viscosity and (ii) structural parameter λ as a function of the extensional rate. Model predictions for ε ˙ = 10 4 10 4 s 1 (TEVP 9).

FIG. 18.

Uniaxial elongation: (i) First normal stress difference and extensional viscosity and (ii) structural parameter λ as a function of the extensional rate. Model predictions for ε ˙ = 10 4 10 4 s 1 (TEVP 9).

Close modal

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.

FIG. 19.

Startup uniaxial elongation: (i) First normal stress difference, (ii) extensional viscosity, and (iii) structural parameter λ time evolution. Model predictions for ε ˙ = 1.0 s 1, ε ˙ = 8.0 s 1, and ε ˙ = 50.0 s 1 (TEVP 9).

FIG. 19.

Startup uniaxial elongation: (i) First normal stress difference, (ii) extensional viscosity, and (iii) structural parameter λ time evolution. Model predictions for ε ˙ = 1.0 s 1, ε ˙ = 8.0 s 1, and ε ˙ = 50.0 s 1 (TEVP 9).

Close modal

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.

Under steady simple shear, ( _ u _ ) T : C _ _ reads as
( _ u _ ) T : C _ _ = γ ˙ 0 C x y = γ ˙ 0 1 G 0 λ n g σ x y .
(A1)
Substituting Eq. (A1) in Eq. (10), the latter reduces to
0 = [ k 1 + | γ ˙ | e ( t R | γ ˙ | ) 2 ] ( 1 λ ) k 2 [ ( _ u _ ) T : C _ _ ] λ m 2 ( 1 λ ) = λ n g G 0 k 2 γ ˙ 0 σ x y λ m 2 k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 .
(A2)
Then, introducing Eq. (A2) into Eq. (8) yields
λ n g G 0 [ ( _ u _ ) T σ _ _ σ _ _ ( _ u _ ) ] + Y P T T λ n g + m 2 m 1 G 0 η 0 k 2 γ ˙ 0 σ x y k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 σ _ _ = [ ( _ u _ ) + ( _ u _ ) T ] ,
(A3)
where Y P T T = e ε P T T t r ( σ _ _ ) G t. For the shear component, we find that (given that the model predicts σ y y = 0)
λ n g G 0 ( γ ˙ 0 σ y y ) + Y P T T λ n g + m 2 m 1 G 0 η 0 k 2 γ ˙ 0 σ x y k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 σ x y σ x y = η 0 G 0 [ k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 ] k 2 ( λ n g + m 2 m 1 ) Y P T T .
(A4)
To obtain the total shear stress, we need to add in Eq. (A4) the Newtonian contribution from plasma ( σ x y , N = η s γ ˙ 0 ). In the limit γ ˙ 0 0 ( λ 1 ) and ε P T T 0 ( Y P T T 1 ),
σ x y , y = σ y = η 0 G 0 k 1 k 2 .
(A5)
Similarly, for the normal stress component,
λ n g G 0 ( 2 γ ˙ 0 σ x y ) + Y P T T λ n g m 2 m 1 G 0 η 0 k 2 γ ˙ 0 σ x y k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 σ x x = 0 σ x x = 2 η 0 [ k 1 + γ ˙ 0 e ( t R γ ˙ 0 ) 2 ] k 2 Y P T T λ m 2 m 1 .
(A6)
In the limit γ ˙ 0 0 ( λ 1 ) and ε P T T 0 ( Y P T T 1 ),
σ x x , y = 2 η 0 k 1 k 2 .
(A7)
Following the standard formalism for oscillatory flows, we can define an angle δ to quantify the ratio of the viscous and elastic response of the material, such that
t a n δ = G G .
(B1)

In Fig. 20, we present the angle δ in degrees as a function of the strain amplitude. We observe that for small γ o, the angle is almost zero, which indicates that the material is solid. Yielding occurs at approximately γ o 0.02, and the material transitions to a viscoelastic solid. As γ o increases further, the angle approaches 90 ° and the viscous nature prevails.

FIG. 20.

Angle δ as a function of the strain amplitude γ.

FIG. 20.

Angle δ as a function of the strain amplitude γ.

Close modal

For the amplitude sweep test at ω = 4 π rad s 1, we choose the limiting cases of γ 0 = 0.10 and γ 0 = 30.0 to examine the variation in the structure parameter and the stress components in a period of oscillation. For γ 0 = 0.10, TEVP 9 predicts δ = 35.0 ° and G > G . In Fig. 21, we present the periodic evolution of (i) λ and (ii) of the shear stress ( σ x y ) and normal stress difference ( N 1 ) components for γ 0 = 0.10. We can observe the phase shift of the shear stress with the strain γ (which is a sine function) and with the shear rate γ ˙ (which is a cosine function). The shear rate amplitude is small ( γ ˙ 0 = 0.20 s 1 ), and this observation supports the negligible normal stress. Concerning the structure parameter λ, it reaches maximum and minimum values when the combination of σ x y and γ ˙ is such that structural breakdown is minimum or maximum, respectively. The time at which this occurs depends on the phase lag between σ x y and γ ˙.

FIG. 21.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 0.10 and ω = 4 π rad s 1 (TEVP 9).

FIG. 21.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 0.10 and ω = 4 π rad s 1 (TEVP 9).

Close modal

For the large value of γ 0 = 30.0, the response is very different as we can see in Fig. 22. In this case, G < G and δ = 88.0 °. Thus, the material is mostly viscous, and σ x y is almost in phase with γ ˙ [Fig. 22(ii)]. However, since we are not in the LVE, the shear stress response is not purely sinusoidal but deviates from the standard cosine function. Moreover, γ ˙ 0 = 59.7 s 1, which means that N 1 is greater than τ x y (as discussed in Sec. IV A). Since σ x y and γ ˙ are almost in phase, when both are approximately zero (at t T / 4 and t 3 T / 4, where T = 2 π / ω is the period of oscillation), the structural parameter λ is maximum. Similarly, when σ x y and γ ˙ are maximum in absolute value (at t 0 and t T / 2), λ is minimum.

FIG. 22.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 30.0 and ω = 4 π rad s 1 (TEVP 9).

FIG. 22.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 30.0 and ω = 4 π rad s 1 (TEVP 9).

Close modal

For the frequency sweep at γ 0 = 9.95, the angle δ ranges from 81.3 ° (for ω = 40.0 rad s 1) to 88.9 ° (for ω = 0.5 rad s 1). Thus, the stress is mostly in phase with γ ˙ and is not purely sinusoidal because we are in the non-LVE region. γ ˙ 0 ranges from 0.79 to 63.3 s 1, and N 1 becomes important as ω increases. For the λ parameter, the evolution is similar to that of Fig. 22(i) throughout the frequency range. The periodic evolution of λ and the stress components for ω = 0.5 rad s 1 and ω = 40.0 rad s 1 are shown in Figs. 23 and 24, respectively.

FIG. 23.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 9.95 and ω = 0.5 rad s 1 (TEVP 9).

FIG. 23.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 9.95 and ω = 0.5 rad s 1 (TEVP 9).

Close modal
FIG. 24.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 9.95 and ω = 40.0 rad s 1 (TEVP 9).

FIG. 24.

Periodic evolution of (i) the structure parameter λ and (ii) the shear and normal stress difference for γ 0 = 9.95 and ω = 40.0 rad s 1 (TEVP 9).

Close modal

Here, we present the effect of each parameter on the predictions of the model. We keep the same values of the parameters from the fitting procedure, and we change only one parameter at a time by an order of magnitude (both higher and lower compared to the original value). We choose the case of a step-down experiment from γ ˙ = 300.0 to 20.0 s 1 to examine the evolution of the shear stress in time. The results are similar in all step-up and step-down cases.

In Fig. 25, we present the effect of the parameters of the elastic modulus, i.e., G 0, n g. We observe that increasing either one of them results in earlier undershoot. This result is expected because the elastic modulus G t = G 0 λ n g increases and the apparent relaxation time t r e l = η t G t decreases. We note that with increasing n g, the term λ n g decreases for the same λ because λ [ 0 , 1 ]. Moreover, the undershoot is more pronounced. To explain this, we need to compare the relaxation time with the thixotropic time, which is defined as t t h i x = 1 k 1. As it is argued in [80] when t r e l t t h i x 1, the thixotropic response is significant. In this case, the material is characterized by an initial fast elastic response followed by slow thixotropic rebuilding. Consequently, the undershoot is more evident when the relaxation time is small compared to the thixotropic time. These observations apply to the first normal stress difference N 1 as well. Finally, for larger G 0 or n g, the shear stress at steady state is larger [see Eq. (A4) of  Appendix A], while the first normal stress difference N 1 is smaller. In fact, this is the only difference in the predictions of σ x y and N 1. For all the other parameters, the effect is qualitatively the same on both the shear stress and the first normal stress difference. Hence, we will continue the discussion by referring only to σ x y.

FIG. 25.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) G 0 and (ii) n g.

FIG. 25.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) G 0 and (ii) n g.

Close modal

Next, we examine the effect of the viscosity parameters η 0 and m 1 in Fig. 26. Increasing η 0, the undershoot is observed later because the apparent relaxation time increases. For the largest η 0, the undershoot is barely evident, and the shear stress is larger. The effect of m 1 is opposite to that of η 0. With increasing m 1, the viscosity decreases because λ m 1 decreases.

FIG. 26.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) η 0 and (ii) m 1.

FIG. 26.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) η 0 and (ii) m 1.

Close modal

The effect of the kinetic parameters k 1 and k 2 is presented in Fig. 27. When the rebuilding constant k 1 increases, we observe larger σ x y overall [Fig. 27(i)] because λ (and consequently the viscosity) is larger [see also Eq. (A4) in  Appendix A for the steady state values]. Moreover, the apparent relaxation time is larger; the undershoot occurs later and it is less significant. The opposite response is predicted when the destruction constant k 2 increases [Fig. 27(ii)]. Larger values of k 2 lead to smaller viscosity and the shear stresses decrease.

FIG. 27.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) k 1 and (ii) k 2.

FIG. 27.

Predictions of TEVP 9 in a step-down experiment from 300.0 to 20.0 s 1. Effect of the model parameters: (i) k 1 and (ii) k 2.

Close modal

Finally, we examine the effect of the exponents m 2 and t R in Fig. 28. When m 2 increases, the destruction term becomes less dominant because λ m 2 decreases. Consequently, the viscosity is higher and the same holds for σ x y and the apparent relaxation time. Thus, the undershoot occurs later, and it is less pronounced [Fig. 28(i)]. The observations are the same for the first normal stress difference.

FIG. 28.

Predictions of TEVP 9 in a step-down experiment. Effect of the model parameters: (i) m 2: step-down from 300.0 to 20.0 s 1 and (ii) t R: step-down from 5.0 to 0.1 s 1.

FIG. 28.

Predictions of TEVP 9 in a step-down experiment. Effect of the model parameters: (i) m 2: step-down from 300.0 to 20.0 s 1 and (ii) t R: step-down from 5.0 to 0.1 s 1.

Close modal

To examine the effect of t R, we choose the step-down experiment from 5.0 to 0.1 s 1 [Fig. 28(ii)]. We do this to examine the range of γ ˙ where the flow rebuilding term is important. As we note in Sec. II, this term decreases fast with γ ˙, and at γ ˙ = 300.0 s 1, it is practically zero. When t R increases, the flow-rebuilding term becomes less significant and λ decreases. As a result, the viscosity σ x y and the relaxation time are smaller. As we can see in Fig. 28(ii), the undershoot occurs earlier.

We do not examine the effect of the ε P T T parameter because it does not contribute significantly to shear flows, except for large γ ˙. Its main contribution is to bound the stress in extensional flows and to introduce extensional thinning.

1.
Brust
,
M.
,
C.
Schaefer
,
R.
Doerr
,
L.
Pan
,
M.
Garcia
,
P. E.
Arratia
, and
C.
Wagner
, “
Rheology of human blood plasma: Viscoelastic versus Newtonian behavior
,”
Phys. Rev. Lett.
110
,
078305
(
2013
).
2.
Varchanis
,
S.
,
Y.
Dimakopoulos
,
C.
Wagner
, and
J.
Tsamopoulos
, “
How viscoelastic is human blood plasma?
,”
Soft Matter
14
,
4238
4251
(
2018
).
3.
Koenig
,
W.
, and
E.
Ernst
, “
The possible role of hemorheology in atherothrombogenesis
,”
Atherosclerosis
94
,
93
107
(
1992
).
4.
Letcher
,
R. L.
,
S.
Chien
,
T. G.
Pickering
,
J. E.
Sealey
, and
J. H.
Laragh
, “
Direct relationship between blood pressure and blood viscosity in normal and hypertensive subjects: Role of fibrinogen and concentration
,”
Am. J. Med.
70
,
1195
1202
(
1981
).
5.
Shiga
,
T.
,
K.
Imaizumi
,
N.
Harada
, and
M.
Sekiya
, “
Kinetics of rouleaux formation using TV image analyzer I.: Human erythrocytes
,”
Am. J. Physiol.
245
,
H252
H258
(
1983
).
6.
Chien
,
S.
, and
K.
Ming Jan
, “
Ultrastructural basis of the mechanism of rouleaux formation
,”
Microvasc. Res.
5
,
155
166
(
1973
).
7.
Balmforth
,
N. J.
,
I. A.
Frigaard
, and
G.
Ovarlez
, “
Yielding to stress: Recent developments in viscoplastic fluid mechanics
,”
Ann. Rev. Fluid Mech.
46
,
121
146
(
2014
).
8.
Kaliviotis
,
E.
,
I.
Ivanov
,
N.
Antonova
, and
M.
Yianneskis
, “
Erythrocyte aggregation at non-steady flow conditions: A comparison of characteristics measured with electrorheology and image analysis
,”
Clin. Hemorheol. Microcirc.
44
,
43
54
(
2010
).
9.
Barnes
,
H. A.
, “
Thixotropy—A review
,”
J. Non-Newton. Fluid Mech.
70
,
1
33
(
1997
).
10.
Mewis
,
J.
, and
N. J.
Wagner
, “
Thixotropy
,”
Adv. Colloid Interface Sci.
147–148
,
214
227
(
2009
).
11.
Larson
,
R. G.
, and
Y.
Wei
, “
A review of thixotropy and its rheological modeling
,”
J. Rheol.
63
,
477
501
(
2019
).
12.
Merrill
,
E. W.
, “
Rheology of blood
,”
Phys. Rev.
49
,
863
888
(
1969
).
13.
Olufsen
,
M. S.
,
C. S.
Peskin
,
W. Y.
Kim
,
E. M.
Pedersen
,
A.
Nadim
, and
J.
Larsen
, “
Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions
,”
Ann. Biomed. Eng.
28
,
1281
1299
(
2000
).
14.
Xu
,
D.
,
E.
Kaliviotis
,
A.
Munjiza
,
E.
Avital
,
C.
Ji
, and
J.
Williams
, “
Large scale simulation of red blood cell aggregation in shear flows
,”
J. Biomech.
46
,
1810
1817
(
2013
).
15.
Fåhræus
,
R.
, and
T.
Lindqvist
, “
The viscosity of the blood in narrow capillary tubes
,”
Am. J. Physiol. Legacy Content
96
,
562
568
(
1931
).
16.
Brundage
,
J. T.
, “
Blood and plasma viscosity determined by the method of concentric cylinders
,”
Am. J. Physiol. Legacy Content
110
,
659
665
(
1934
).
17.
Merrill
,
E. W.
,
E. R.
Gilliland
,
G.
Cokelet
,
H.
Shin
,
A.
Britten
, and
R. E.
Wells
, “
Rheology of blood and flow in the microcirculation
,”
J. Appl. Physiol.
18
,
255
260
(
1963
).
18.
Dintenfass
,
L.
, “
Viscosity of the packed red and white blood cells
,”
Exp. Mol. Pathol.
4
,
597
605
(
1965
).
19.
Chien
,
S.
,
S.
Usami
,
H. M.
Taylor
,
J. L.
Lundberg
, and
M. I.
Gregersen
, “
Effects of hematocrit and plasma proteins on human blood rheology at low shear rates
,”
J. Appl. Physiol.
21
,
81
87
(
1966
).
20.
Pirofsky
,
B.
, “
The determination of blood viscosity in man by a method based on Poiseuille’s law
,”
J. Clin. Invest.
32
,
292
298
(
1953
).
21.
Wayland
,
H.
, “
Rheology and the microcirculation
,”
Gastroenterology
52
,
342
355
(
1967
).
22.
Morris
,
C. L.
,
D. L.
Rucknagel
,
R.
Shukla
,
R. A.
Gruppo
,
C. M.
Smith
, and
P.
Blackshear
, “
Evaluation of the yield stress of normal blood as a function of fibrinogen concentration and hematocrit
,”
Microvasc. Res.
37
,
323
338
(
1989
).
23.
Picart
,
C.
,
J.-M.
Piau
,
H.
Galliard
, and
P.
Carpentier
, “
Human blood shear yield stress and its hematocrit dependence
,”
J. Rheol.
42
,
1
12
(
1998
).
24.
Merrill
,
E. W.
,
C. S.
Cheng
, and
G. A.
Pelletier
, “
Yield stress of normal human blood as a function of endogenous fibrinogen
,”
J. Appl. Physiol.
26
,
1
3
(
1969
).
25.
Huang
,
C. R.
,
H. Q.
Chen
,
W. D.
Pan
,
T.
Shih
,
D. S.
Kristol
, and
A. L.
Copley
, “
Effects of hematocrit on thixotropic properties of human blood
,”
Biorheology
24
,
803
810
(
1987
).
26.
Kaibara
,
M.
, and
E.
Fukada
, “
Non-Newtonian viscosity and dynamic viscoelasticity of blood during clotting
,”
Biorheology
6
,
73
84
(
1969
).
27.
Thurston
,
G. B.
, “
Viscoelasticity of human blood
,”
Biophys. J.
12
,
1205
1217
(
1972
).
28.
Thurston
,
G. B.
, “
Rheological parameters for the viscosity viscoelasticity and thixotropy of blood
,”
Biorheology
16
,
149
162
(
1979
).
29.
Bureau
,
M.
,
J. C.
Healy
,
D.
Bourgoin
, and
M.
Joly
, “
Rheological hysteresis of blood at low shear rate
,”
Biorheology
17
,
191
203
(
1980
).
30.
McMillan
,
D. E.
,
J.
Strigberger
, and
N. G.
Utterback
, “
Rapidly recovered transient flow resistance: A newly discovered property of blood
,”
Am. J. Physiol.
253
,
H919
H926
(
1987
).
31.
Vlastos
,
G.
,
D.
Lerche
,
B.
Koch
,
O.
Samba
, and
M.
Pohl
, “
The effect of parallel combined steady and oscillatory shear flows on blood and polymer solutions
,”
Rheol. Acta
36
,
160
172
(
1997
).
32.
Sousa
,
P. C.
,
J.
Carneiro
,
R.
Vaz
,
A.
Cerejo
,
F. T.
Pinho
,
M. A.
Alves
, and
M. S. N.
Oliveira
, “
Shear viscosity and nonlinear behavior of whole blood under large amplitude oscillatory shear
,”
Biorheology
50
,
269
282
(
2013
).
33.
Tomaiuolo
,
G.
,
A.
Carciati
,
S.
Caserta
, and
S.
Guido
, “
Blood linear viscoelasticity by small amplitude oscillatory flow
,”
Rheol. Acta
55
,
485
495
(
2016
).
34.
Huang
,
C.-R.
, and
S.-S.
Horng
, “
Viscoelastic-thixotropy of blood
,”
Clin. Hemorheol. Microcirc.
15
,
25
36
(
1995
).
35.
Horner
,
J. S.
,
M. J.
Armstrong
,
N. J.
Wagner
, and
A. N.
Beris
, “
Measurements of human blood viscoelasticity and thixotropy under steady and transient shear and constitutive modeling thereof
,”
J. Rheol.
63
,
799
813
(
2019
).
36.
Yilmaz
,
F.
, and
M.
Gundogdu
, “
A critical review on blood flow in large arteries; relevance to blood rheology, viscosity models, and physiologic conditions
,”
Korea-Australia Rheol. J.
20
,
197
211
(
2008
).
37.
Sousa
,
P. C.
,
F. T.
Pinho
,
M. A.
Alves
, and
M. S. N.
Oliveira
, “
A review of hemorheology: Measuring techniques and recent advances
,”
Korea Australia Rheol. J.
28
,
1
22
(
2016
).
38.
Apostolidis
,
A. J.
, and
A. N.
Beris
, “
Modeling of the blood rheology in steady-state shear flows
,”
J. Rheol.
58
,
607
633
(
2014
).
39.
Anand
,
M.
,
J.
Kwack
, and
A.
Masud
, “
A new generalized Oldroyd-B model for blood flow in complex geometries
,”
Int. J. Eng. Sci.
72
,
78
88
(
2013
).
40.
Armstrong
,
M.
, and
J.
Tussing
, “
A methodology for adding thixotropy to Oldroyd-8 family of viscoelastic models for characterization of human blood
,”
Phys. Fluids
32
,
094111
(
2020
).
41.
Armstrong
,
M.
, and
A.
Pincot
, “
Integration of thixotropy into Giesekus model for characterization of human blood
,”
AIP Adv.
11
,
035029
(
2021
).
42.
Quemada
,
D.
, and
R.
Droz
, “
Blood viscoelasticity and thixotropy from stress formation and relaxation measurements: A unified model
,”
Biorheology
20
,
635
651
(
1983
).
43.
Wei
,
Y.
,
M. J.
Solomon
, and
R. G.
Larson
, “
A multimode structural kinetics constitutive equation for the transient rheology of thixotropic elasto-viscoplastic fluids
,”
J. Rheol.
62
,
321
342
(
2018
).
44.
Armstrong
,
M. J.
,
A. N.
Beris
,
S. A.
Rogers
, and
N. J.
Wagner
, “
Dynamic shear rheology of a thixotropic suspension: Comparison of an improved structure-based model with large amplitude oscillatory shear experiments
,”
J. Rheol.
60
,
433
450
(
2016
).
45.
Wei
,
Y.
,
M. J.
Solomon
, and
R. G.
Larson
, “
Quantitative nonlinear thixotropic model with stretched exponential response in transient shear flows
,”
J. Rheol.
60
,
1301
1315
(
2016
).
46.
Dullaert
,
K.
, and
J.
Mewis
, “
A structural kinetics model for thixotropy
,”
J. Non-Newton. Fluid Mech.
139
,
21
30
(
2006
).
47.
Quemada
,
D.
, “
Rheological modelling of complex fluids: IV: Thixotropic and ‘thixoelastic’ behaviour.: Start-up and stress relaxation, creep tests and hysteresis cycles
,”
Eur. Phys. J. Appl. Phys.
5
,
191
207
(
1999
).
48.
Beris
,
A. N.
,
E.
Stiakakis
, and
D.
Vlassopoulos
, “
A thermodynamically consistent model for the thixotropic behavior of concentrated star polymer suspensions
,”
J. Non-Newton. Fluid Mech.
152
,
76
85
(
2008
).
49.
Apostolidis
,
A. J.
,
M. J.
Armstrong
, and
A. N.
Beris
, “
Modeling of human blood rheology in transient shear flows
,”
J. Rheol.
59
,
275
298
(
2015
).
50.
Armstrong
,
M.
,
J.
Horner
,
M.
Clark
,
M.
Deegan
,
T.
Hill
,
C.
Keith
, and
L.
Mooradian
, “
Evaluating rheological models for human blood using steady state, transient, and oscillatory shear predictions
,”
Rheol. Acta
57
,
705
728
(
2018
).
51.
Beris
,
A. N.
,
J. S.
Horner
,
S.
Jariwala
,
M. J.
Armstrong
, and
N. J.
Wagner
, “
Recent advances in blood rheology: A review
,”
Soft Matter
17
,
10591
10613
(
2021
).
52.
Armstrong
,
M.
,
A.
Pincot
,
S.
Jariwala
,
J.
Horner
,
N.
Wagner
, and
A.
Beris
, “
Tensorial formulations for improved thixotropic viscoelastic modeling of human blood
,”
J. Rheol.
66
,
327
347
(
2022
).
53.
Giannokostas
,
K.
,
P.
Moschopoulos
,
S.
Varchanis
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Advanced constitutive modeling of the thixotropic elasto-visco-plastic behavior of blood: Description of the model and rheological predictions
,”
Materials
13
,
4184
(
2020
).
54.
Varchanis
,
S.
,
G.
Makrigiorgos
,
P.
Moschopoulos
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Modeling the rheology of thixotropic elasto-visco-plastic materials
,”
J. Rheol.
63
,
609
639
(
2019
).
55.
Saramito
,
P.
, “
A new constitutive equation for elastoviscoplastic fluid flows
,”
J. Non-Newton. Fluid Mech.
145
,
1
14
(
2007
).
56.
Stephanou
,
P. S.
, and
G. G.
Georgiou
, “
A nonequilibrium thermodynamics perspective of thixotropy
,”
J. Chem. Phys.
149
,
244902
(
2018
).
57.
Dimitriou
,
C. J.
, and
G. H.
McKinley
, “
A comprehensive constitutive law for waxy crude oil: A thixotropic yield stress fluid
,”
Soft Matter
10
,
6619
6644
(
2014
).
58.
Fraggedakis
,
D.
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Yielding the yield stress analysis: A thorough comparison of recently proposed elasto-visco-plastic (EVP) fluid models
,”
J. Non-Newton. Fluid Mech.
236
,
104
122
(
2016
).
59.
Dimitriou
,
C. J.
,
R. H.
Ewoldt
, and
G. H.
McKinley
, “
Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress)
,”
J. Rheol.
57
,
27
70
(
2013
).
60.
Phan-Thien
,
N.
, “
A nonlinear network viscoelastic model
,”
J. Rheol.
22
,
259
283
(
1978
).
61.
Bambardekar
,
K.
,
J. A.
Dharmadhikari
,
A. K.
Dharmadhikari
,
T.
Yamada
,
T.
Kato
,
H.
Kono
,
Y.
Fujimura
,
S.
Sharma
, and
D.
Mathur
, “
Shape anisotropy induces rotations in optically trapped red blood cells
,”
J. Biomed. Opt.
15
,
041504
(
2010
).
62.
Steenbergen
,
W.
,
F.
de Mul
, and
R.
Kolkman
, “
Light-scattering properties of undiluted human blood subjected to simple shear
,”
J. Opt. Soc. Am. A
16
,
2959
2967
(
1999
).
63.
Moschopoulos
,
P.
,
E.
Chrysou
,
K.
Giannokostas
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “Modeling of the human blood rheology and simulation of its flow in elastic micro-vessels,” 91th SOR meeting, Raleigh, NC, October 2019.
64.
International Committee for Standardization in Haematology, “
Recommendation for a selected method for the measurement of plasma viscosity: International committee for standardization in haematology
,”
J. Clin. Pathol.
37,
1147
1152
(
1984
).
65.
Varchanis
,
S.
,
S. J.
Haward
,
C. C.
Hopkins
,
A.
Syrakos
,
A. Q.
Shen
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Transition between solid and liquid state of yield-stress fluids under purely extensional deformations
,”
Proc. Natl. Acad. Sci. U.S.A.
117
,
12611
12617
(
2020
).
66.
Thompson
,
R. L.
,
L. U. R.
Sica
, and
P. R.
de Souza Mendes
, “
The yield stress tensor
,”
J. Non-Newton. Fluid Mech.
261
,
211
219
(
2018
).
67.
de Cagny
,
H.
,
M.
Fazilati
,
M.
Habibi
,
M. M.
Denn
, and
D.
Bonn
, “
The yield normal stress
,”
J. Rheol.
63
,
285
290
(
2019
).
68.
Dimakopoulos
,
Y.
, and
K.
Giannokostas
, “
From biorheology to biofluid mechanics: Elucidating the behavior of biofluids in complex flows
,”
Sci. Talks
5
,
100139
(
2023
).
69.
Joly
,
M.
,
C.
Lacombe
, and
D.
Quemada
, “
Application of the transient flow rheology to the study of abnormal human bloods
,”
Biorheology
18
,
445
452
(
1981
).
70.
Min Kim
,
J.
,
A. P. R.
Eberle
,
A.
Kate Gurnon
,
L.
Porcar
, and
N. J.
Wagner
, “
The microstructure and rheology of a model, thixotropic nanoparticle gel under steady shear and large amplitude oscillatory shear (LAOS)
,”
J. Rheol.
58
,
1301
1328
(
2014
).
71.
Donley
,
G. J.
,
P. K.
Singh
,
A.
Shetty
,
S. A.
Rogers
, and
D. A.
Weitz
, “
Elucidating the G″ overshoot in soft materials with a yield transition via a time-resolved experimental strain decomposition
,”
Proc. Natl. Acad. Sci. U.S.A.
117
,
21945
21952
(
2020
).
72.
Kamani
,
K.
,
G. J.
Donley
, and
S. A.
Rogers
, “
Unification of the rheological physics of yield stress fluids
,”
Phys. Rev. Lett.
126
,
218002
(
2021
).
73.
Kamani
,
K. M.
,
G. J.
Donley
,
R.
Rao
,
A. M.
Grillet
,
C.
Roberts
,
A.
Shetty
, and
S. A.
Rogers
, “
Understanding the transient large amplitude oscillatory shear behavior of yield stress fluids
,”
J. Rheol.
67
,
331
352
(
2023
).
74.
Hyun
,
K.
,
M.
Wilhelm
,
C. O.
Klein
,
K. S.
Cho
,
J. G.
Nam
,
K. H.
Ahn
,
S. J.
Lee
,
R. H.
Ewoldt
, and
G. H.
McKinley
, “
A review of nonlinear oscillatory shear tests: Analysis and application of large amplitude oscillatory shear (LAOS)
,”
Prog. Polym. Sci.
36
,
1697
1753
(
2011
).
75.
Ewoldt
,
R. H.
, and
G. H.
McKinley
, “
On secondary loops in LAOS via self-intersection of Lissajous–Bowditch curves
,”
Rheol. Acta
49
,
213
219
(
2010
).
76.
Sousa
,
P. C.
,
R.
Vaz
,
A.
Cerejo
,
M. S. N.
Oliveira
,
M. A.
Alves
, and
F. T.
Pinho
, “
Rheological behavior of human blood in uniaxial extensional flow
,”
J. Rheol.
62
,
447
456
(
2018
).
77.
Giannokostas
,
K.
,
Y.
Dimakopoulos
,
A.
Anayiotos
, and
J.
Tsamopoulos
, “
Advanced constitutive modeling of the thixotropic elasto-visco-plastic behavior of blood: Steady-state blood flow in microtubes
,”
Materials
14
,
367
(
2021
).
78.
Giannokostas
,
K.
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Shear stress and intravascular pressure effects on vascular dynamics: Two-phase blood flow in elastic microvessels accounting for the passive stresses
,”
Biomech. Model. Mechanobiol.
21
,
1659
1684
(
2022
).
79.
Giannokostas
,
K.
,
D.
Photeinos
,
Y.
Dimakopoulos
, and
J.
Tsamopoulos
, “
Quantifying the non-Newtonian effects of pulsatile hemodynamics in tubes
,”
J. Non-Newton. Fluid Mech.
298
,
104673
(
2021
).
80.
Jamali
,
S.
, and
G. H.
McKinley
, “
The Mnemosyne number and the rheology of remembrance
,”
J. Rheol.
66
,
1027
1039
(
2022
).