Advanced Constitutive Modeling of the Thixotropic Elasto-Visco-Plastic Behavior of Blood: Steady-State Blood Flow in Microtubes

The present work focuses on the in-silico investigation of the steady-state blood flow in straight microtubes, incorporating advanced constitutive modeling for human blood and blood plasma. The blood constitutive model accounts for the interplay between thixotropy and elasto-visco-plasticity via a scalar variable that describes the level of the local blood structure at any instance. The constitutive model is enhanced by the non-Newtonian modeling of the plasma phase, which features bulk viscoelasticity. Incorporating microcirculation phenomena such as the cell-free layer (CFL) formation or the Fåhraeus and the Fåhraeus-Lindqvist effects is an indispensable part of the blood flow investigation. The coupling between them and the momentum balance is achieved through correlations based on experimental observations. Notably, we propose a new simplified form for the dependence of the apparent viscosity on the hematocrit that predicts the CFL thickness correctly. Our investigation focuses on the impact of the microtube diameter and the pressure-gradient on velocity profiles, normal and shear viscoelastic stresses, and thixotropic properties. We demonstrate the microstructural configuration of blood in steady-state conditions, revealing that blood is highly aggregated in narrow tubes, promoting a flat velocity profile. Additionally, the proper accounting of the CFL thickness shows that for narrow microtubes, the reduction of discharged hematocrit is significant, which in some cases is up to 70%. At high pressure-gradients, the plasmatic proteins in both regions are extended in the flow direction, developing large axial normal stresses, which are more significant in the core region. We also provide normal stress predictions at both the blood/plasma interface (INS) and the tube wall (WNS), which are difficult to measure experimentally. Both decrease with the tube radius; however, they exhibit significant differences in magnitude and type of variation. INS varies linearly from 4.5 to 2 Pa, while WNS exhibits an exponential decrease taking values from 50 mPa to zero.


Introduction
Understanding blood flow is of high theoretical and practical importance as it is directly associated with the pathophysiology and the development of diseases such as endotheliitis, microangiopathy, COVID-19 [1,2] in the microvasculature of any human being. Consequently, a consistent dynamic model is necessary for assessing the hemodynamic resistance and its regulation in the microcirculation [3]. We have already developed an integrated constitutive model of blood rheology in our recent work [4]. Blood has a pronounced non-Newtonian character, primarily attributed to aggregation, disaggregation, deformation, orientation, and migration of the erythrocytes [5][6][7][8][9]. Our model encom-Previous mathematical models [26,43] treated both fluid phases as generalized Newtonian fluids, which is a quite simplified approach for blood flows in microtubes. Das et al. [44] proposed a two-phase Casson model to describe the in vitro velocity profiles of blood flow in low flow rates, while Sriram et al. [27] adopted the well-known Quemada model [45] for the determination of the CFL thickness under various rheological conditions. More recently, Moyers-Gonzalez and Owens [46] used principles of the kinetic theory to derive a nonhomogeneous hemorheological model, which was applied for the evaluation of the CFL thickness using empirical laws for the reduction of the hematocrit within the microtube. The model was adopted for the first time in multidimensional flows by Dimakopoulos et al. [47]. Another approach is through multi-particle flow models [48][49][50][51][52], which can also give accurate predictions of local variation of the CFL width, with the major drawback being the high computational cost. Additionally, Qi et al. [53] and Narshimhan et al. [54] based on a coarse-grained theory predicted migration effects and the erythrocytes concentration along with the CFL layer in Couette flow and pressure-driven flow in tubes. This model takes into consideration the wall-induced hydrodynamic lift and the cell-cell interactions, predicting the Fähraeus-Lindqvist effect.
In this work, we combine both hemorheological features in a comprehensive mathematical formulation for the two-phase blood flow in microtubes. We invoke a thixotropic elasto-visco-plastic (TEVP) constitutive model to accurately predict the stresses in the blood in the RBC-rich central region. This model is described in detail in our previous work [4] and fitted on the steady and transient experimental data of McMillan et al. [55]. Following experimental evidence, the plasma phase is treated as a viscoelastic fluid, represented by the linear form of the well-known Phan-Tien-Tanner (or PTT) model. Through our advanced modelling, we try to address some of the open questions in microcirculation regarding the various migration mechanisms and the complex microstructure of blood. Particularly, we investigate the interplay between aggregation and the formation of the cell-free layer and the effect of aggregability on velocity profiles in microtubes with cross-sections comparable to the RBC diameter. Moreover, we examine the intensity of the migration and how it is affected by the diameter of the tube or the imposed pressure gradient. Additionally, we examine the distribution of the relaxation times of blood within the tube along with the aggregation size, which varies with the available space. Finally, we predict the normal and shear viscoelastic stresses in the flow field and how the rheological conditions affect the thixotropic microstructural state of blood.
Our contribution is divided into four sections: In Section 2, we briefly describe the problem formulation, the hemodynamic constraints, and the Thixotropic Elasto-Visco-Plastic (TEVP) constitutive modeling of the whole blood and the viscoelastic modeling of the plasma phase. Section 3 presents the validation of our model with experimental data. A thorough parametric analysis follows in Section 4, which describes the effect of the rheological parameters on blood behavior extensively.

Problem Formulation
We consider the transient blood flow in two separate regions (i = 1,2), one for the RBC rich core (i = 1), placed at the center of the microtube, and a second one for the plasma layer, placed adjacent to the microtube's wall (i = 2) ( Figure 1). The momentum balance for the two-phase blood/plasma flow of density ρ (i) in a microtube of radius R is then written as where U  Due to the different composition of each layer, we assume a different constitutive model for each phase. For the blood core, its modeling is based on the formulation reported in our tensorial TEVP form [4] and especially for the blood constitutive modeling. In this recent work, we invoked a consistent and validated model for TEVP materials proposed by Varchanis et al. [56], which features a coupling of the tensorial constitutive model by Saramito [57] for EVP materials with thixotropy. Another form of the proposed TEVP formulation, without thixotropy, was recently used to evaluate the transition between the solid and liquid state of elasto-viscoplastic fluids under extensional flow [58]. In addition, we assume a viscoelastic behavior for the proteinic plasma phase (CFL), which is mostly determined by fibrinogen concentration. The starting point is the decomposition of the extra stress tensor τ = (i) to solvent and viscoelastic terms as Regarding the modeling of the RBC-rich region, our model accounts for the viscoelastic contribution of the RBCs and the viscoelasticity of blood plasma since it has been fitted to high-shear rate plateau [4]. Consequently, there is no explicit Newtonian contribution ( τ = n (1) = 0) as long as the solvent (plasma) is included in the viscoelastic term (τ = ve (1) ).
Additionally, the presence of plasma proteins [59,60] dominate the stress contribution in the pure plasma phase (CFL), and hence the Newtonian contribution (water) is assumed to be negligible there as well τ = n (2) = 0.

Whole Blood Constitutive Modeling
The total rate of deformation tensor is decomposed into an elastic contribution D The elastic term of the total deformation rate tensor D = e (1) accounts for memory effects and can be written as where the upper-convected time derivative is (1) , while G denotes the elastic modulus of blood. The effect of plasticity is is the magnitude of the stress tensor, η t is the plastic viscosity, and τ y is the blood yield stress (see [4]). The l-PTT function is defined as where ε PTT is a fitted parameter, introducing shear-thinning along with bounding of the extensional viscosity. Given that the total rate of deformation tensor is equal to The time evolution of the structure (or rouleaux) parameter is then given by Wei et al. [61] as dλ dt = (k 1 + k 2 ϕ n 1 )(1 − λ) − k 3 ϕ n 2 λ n 3 (8) The first term in the RHS of Equation (8) represents the rebuilding of rouleaux, while the second one refers to the disintegration process. We base the dependency of the level of structure explicitly on the stresses, via the parameter ϕ [61] Yielding occurs when ϕ > 0, and this is the von Mises criterion. The plastic viscosity of blood η t (λ) is a thixotropy-dependent variable via the structural variable λ η t (λ) = η 0 λ m 1 (10) where η 0 , m 1 are fitting parameters. Similarly, we can define the relaxation time of blood, χ, as the ratio of the plastic viscosity to the shear modulus as To extract realistic values for the set of the parameters of the constitutive model, we adopt a non-linear regression procedure [62] on experimental data for steady-state and transient experiments. The steady-state experiment is the simple shear flow one providing the shear stress response as a function of the imposed shear rate. Our model is fitted on the experimental data of McMillan et al. [55] for a hematocrit equal to 45%, and the set of the adjustable parameters (eleven in number) are reported in Table 1, along with the nomenclature. Alternatively, a modified form of the Casson model [63] is commonly used to simulate blood behavior under steady-conditions. It encompasses the property of the yield stress along with the shear-thinning description of blood and is given by where τ y and µ are the yield stress and the viscosity, respectively, while they both are fitted to the steady state experimental data of McMillan et al. [55] and tabulated in Table 2. A generalized Newtonian constitutive model, such as the Casson model [63], needs only steady state experiments so as for the adjustable parameters to be determined. A plethora of steady shear experiments with different hematocrits reported in the literature, enables the extraction of a mathematical correlation between core and discharged hematocrit.

Plasma Constitutive Modeling
Although in most investigations, blood plasma has been considered a Newtonian fluid, recent shear and extension dominated flow experiments prepared by Brust et al. [59] revealed that blood plasma features bulk viscoelasticity. They reported phenomena that have been widely studied, both experimentally and theoretically, in polymeric solutions [64][65][66], clearly indicating that human blood plasma has viscoelastic properties. In contrast to Newtonian fluids, blood plasma exhibited a shear-thinning viscosity. Despite the evidence of the existence of blood plasma viscoelasticity, only a few investigations take this into account, such as the work of Varchanis et al. [60], who provided a complete data set of the rheometric material functions of plasma rheology in simple shear and elongational flows. Finally, combining their results with previous in vivo measurements, they additionally found that the viscoelasticity of human blood plasma must not be ignored when examining the flow of whole blood in micro-tube, such as arterioles or capillaries. Consequently, we adopt viscoelastic constitutive modelling via the Phan-Thien-Tanner (or PTT) model [67][68][69]. The model is an extension of the Maxwell model accompanied with the upper convected derivative to include a function dependent upon tr(τ = (2) ve ), the trace of the where λ pl is the relaxation time of plasma and η pl is the viscosity of plasma, which are both adjustable parameters while ε PTT,pl is the PTT parameter which is responsible for bounding the extensional viscosity. Fitting the l-PTT formula to the data reported in [60], we extract the parameters of the model illustrated in Table 3. The one-dimensional forms of the above equations are presented in detail in Appendix A.

Hemodynamical Constraints
The model presented above is not complete, and its equations cannot be solved unless the location of the core/plasma interface or the core radius δ, can be determined. Hence, this location must be computed along with all the other unknowns, making this a moving boundary problem. To calculate it, we need to estimate the impact of the hemorheological parameters on the velocity and stress fields by introducing hemodynamical constraints, which quantify major microcirculation effects [70]. As such, we refer to the Fähraeus [71] and Fähraeus-Lindqvist [17] phenomena. Both are related to the increase of the CFL thickness as the tube radius decreases, attributed to the increasing tendency of the RBCs to migrate towards the center of the vessel [72]. Unfortunately, neither effect can be predicted using first principles, so we have to rely on experiments and develop suitable correlations. The Fähraeus effect is the reduction of the bulk hematocrit, whereas the Fähraeus-Lindqvist effect accounts for the decrease in the apparent viscosity (η app ).
Before presenting these correlations, we need to clarify the different existing definitions of the hematocrit, which represents the volume fraction of RBCs in whole blood expressed as a percentage. Namely, these are the discharged hematocrit H d , the tube hematocrit, H t , and the core hematocrit, H c . The discharged hematocrit is the velocity-weighted average of the local volume fraction of the erythrocytes given by (14) where H(r) is the local volume fraction of RBCs across a section of the tube. Discharged hematocrit is a measure of its bulk value accounting for the fact that the local volume fraction at a specific radial position is carried by the local axial velocity, which also varies radially. Since the higher RBC concentration is near the tube center due to cell aggregation and it is carried by the higher local velocity, H d is larger than H t , which is defined next. The tube hematocrit, H t represents the average hematocrit within the vessel, assuming that the velocity has a plug-flow (radially constant) profile. Hence, H t is given by In general and in this work, given that the RBC migration away from the wall leads to a considerable decrease of the local viscosity, which cannot be described due to the absence of any realistic differential model, we assume that H(r) follows a simple distribution which assigns a constant value to the hematocrit in the core region denoted by H c and a zero value in the CFL. If we substitute Equation (16) into Equation (15), we readily get a linear relationship between H t and H c given by where δ is the core radius shown in Figure 1. Similarly, H d can be related to H c using a simple mass balance of RBCs in any cross-section The difference between the three types of hematocrit diminishes, i.e. H t ≈ H c ≈ H d as δ/R → 1 . Tube and core hematocrits are forms that have been introduced exclusively in arteriolar blood flows.
Next, we rely on extensive experimental observations to quantify macroscopic hemodynamical properties [73,74] and develop the necessary correlations as follows. The Fähraeus-Lindqvist effect has been examined in detail by Pries and Secomb [73]. To this end, a series of experiments has been undertaken for different hemodynamical conditions (i.e., vessel diameter, discharged hematocrit, etc.). The non-Newtonian blood flow has been correlated with the well-known Hagen-Poiseuille law for laminar flow of Newtonian fluids in a tube. To match the blood flow data given in [73], the relationship of the flow-rate with the pressure-drop for different tube diameters is revised as follows Q = JD 4 128 η rel η pl (19) where Q and J are the volume flow rate and pressure gradient, respectively, and η rel is the relative apparent viscosity (= η app /η pl ). For a tube of a given diameter, Q, J and η pl are measured, and Equation (19) is used to determine η rel . Then, we invoke the Fähraeus-Lindqvist phenomenon, i.e., the strong dependence of η rel on the tube diameter and the discharged hematocrit (H d ), by the empirical law derived from the experiments conducted in glass tubes by Pries and Secomb [73]. This is a three-part relationship, which is relatively complicated, and hence, we have replaced it with our simplified form (the so called GDAT model) given by a single equation after a non-linear fitting on experimental data [73], which is provided by where D is the tube diameter in µm. A more detailed description of this relationship and comparison with the original form of Pries and Secomb [73] (Equations (A11)-(A13)) is presented in Appendix B. Equation (18) determines H d . Measurements of the tube hematocrit are typically reported relative to the discharged hematocrit, i.e., the ratio H t /H d is measured. This ratio is related to the Fähraeus effect and was measured experimentally by Pries et al. [74] for different values of the discharged hematocrit and a wide range of tube diameters. The relationship that describes this ratio is a function of H d and R (in µm) and is given by Equation (21) determines H t and combing its value with the value of H d and Equation (17) yields the CFL thickness. We discretize the above system of equations (Equations (1)-(21)) using central finite differences [75] in space and a backward Euler scheme for their time integration. Since we assume that the location of the interface (δ) is part of the solution, we deal with a moving boundary problem and hence we need to introduce a linear spatial transformation to solve the system in a fixed domain. The one-dimensional forms of the entire equation set, along with the boundary conditions, are presented in detail in Appendix A.
Although we present the steady-state values, our model is inherently transient, and thus all the predictions must be obtained through a transient solution. For all the examined cases, we assume an initially unperturbed and unyielded state (λ |t=0 = 1), while the stress and the velocity field are equal to zero (τ ve|t=0 = 0, U z|t=0 = 0), as it is described in detail in Appendix C.

Validation
Initially, we need to validate the predictions of our model against relevant flows reported in the literature. It must be pointed out that throughout the validation procedure, the examined flow conditions do not exactly match those of the corresponding experiments regarding the hematocrit value. As described in Section 2 and ref. [4], the parametrization of the TEVP model in this study corresponds to a rheological description of blood with a bulk hematocrit equal to 45%, and thus, the simulations are implemented by using a constant core hematocrit equal to this value. Most of the investigations assume a constant H d rather than a constant H c . In the following Figures, blue solid curves represent the steady-state predictions of our TEVP model, while the symbols are data obtained by experimental studies.
Firstly, our predictions for velocity profiles are compared with experimental micro-PIV measurements of human blood flow in glass tubes. We select experimental data where flow conditions are close enough to ours regarding the core hematocrit H c . Apparently, this is a challenging task since the literature lacks systematic works that present complete hemorheological data and hemodynamical measurements. To this end, we compare our steady-state results with those reported in the work of Bugliarello and Sevilla [76] for two cross-sections with different hemorheological conditions. Results under steady conditions are shown in Figure 2a,b for a discharged hematocrit equal to 40% in arterioles of 20 µm and 35 µm radius, respectively, while our simulation is under a constant core hematocrit value equal to 45%. Additionally, the experiments are conducted for an imposed mean axial velocity equal to U mean = 13 mm/s and U mean = 3.8 mm/s for 20 µm and 35 µm radius, respectively. Despite the somewhat different values of core hematocrit, our model is in excellent agreement with the experimental observations for both radii. We also provide the predictions of the Casson model [63] for the aforementioned experiments. We can observe that the predictions of the Casson model are not quantitatively good, especially for the case of R = 35 µm, where the predicted velocity profile is more plug-like than that observed in the experiments. Additionally, the predictions of the inelastic Casson model are not close enough to the experimental observation regarding the CFL thickness. This is more obvious in Figure 2a,b where the CFL thickness is highlighted with light red, yellow and grey for the Casson, TEVP and experimental value, respectively. Indicatively, for the flow field in a microtube with a radius equal to 20 µm the predictions for the CFL thickness are 5.2 µm, 3.3 µm of the Casson and the TEVP models, while the corresponding values for the case of 35 µm radius, they are 5.25 µm and 3.85 µm, respectively. Our model is in excellent agreement with experimental observations of the CFL as Bugliarello and Sevilla [76] reported in their work that the CFL thickness is 3.3 µm and 3.36 µm for 20 µm and 35 µm tube radius, respectively. On the contrary, the predictions of the CFL thickness by the Casson model are not as good. In Table 4 we provide the predictions of our model for CFL thickness (w), flow rate (Q), Interfacial Normal (INS) and Shear Stresses (ISS), Wall Normal Stress (WNS) as well as the Wall Shear Stress (WSS) regarding the simulations with 20 µm and 35 µm radius, respectively. In Figure 2c, we compare our numerical predictions for the dimensionless thickness of the cell-free peripheral layer with those experimentally observed [5,77,78] for a discharged hematocrit equal to 45%, and a pressure gradient J = 10 4 Pa/m. The experimental points refer to tube diameters ranging from 30.8 µm and 132.3 µm with a pseudo-shear-rate varying between 5.38 s −1 and 15.91s −1 , respectively. Our model captures very accurately the observed decrease in the relative CFL thickness (w/R), when the microtube radius increases, depicting an overall deviation of about 2%. This discrepancy is quite insignificant if one considers that our simulation encompasses a constant H c , while the experiments employ a constant H d . In the same Figure we provide the predictions of a Newtonian model [26], the model of Casson [63], and the model proposed by Moyers-Gonzalez and Owens [46] for the normalized CFL thickness. It is evident that our model provides more accurate results compared to those predicted by the other models. Although both the Newtonian and Casson models predict a similar behavior with that presented by our model, they both deviate considerably from the experimental values throughout, and the model of Moyers-Gonzalez and Owens [46] does not follow the trend of the experiments for radii below 30 µm. All these are absolutely normal since Newtonian and Casson models have only one and two adjustable parameters for the description of one or two out of five major mechanisms of the rheological response of blood. On the other hand the Owens' model has eleven parameters as our TEVP model but it does not account for blood plasticity making it inaccurate at low shear-rates. Definitely, the coupled appearance of several viscous, elastic, plastic and thixotropic (RBC's aggregation and disaggregation) phenomena and the underlying mechanisms necessitates the introduction of an adequate number of adjustable parameters to accurate represent their contribution in a thixotropic elasto-visco-plastic (TEVP) model and consequently the study of hemodynamics. Sevilla [76], (b) velocity profile for R = 35 µm under the hemorheological conditions reported in Bugliarello and Sevilla [76], (c) Normalized thickness of CFL from various experimental investigations [5,77,78] accompanied with the predictions of a Newtonian model [26], the model of Casson [63] and Moyers-Gonzalez and Owens [46], (d) Total flow rate for the flow conditions reported in Bugliarello and Sevilla [76] for R = 20 µm, (e) Comparison of wall shear stress predictions with respect to pseudo shear rate . γ between TEVP and experimental observations from [79]. Comparison of relative wall shear stress WSS rel for various pseudo shear rates . γ between TEVP, Casson models and the experimental data from [80] for (f) R = 15 µm and (g) R = 50 µm. Figure 2d illustrates the experimentally measured average volumetric flow rate of blood with respect to the imposed pressure gradient for 40% of discharged hematocrit. The experimental data are provided by the work of Bugliarello and Sevilla [76], and the only available data which were close enough to our rheological conditions were those for a microtube with a radius equal to 20 µm. The predicted flow-rate demonstrates a quite good agreement with the experimental data with a small overestimation at low pressure gradients. This deviation is mainly attributed to the presence of the migration effects, which may provoke a considerable deviation in the hematocrit of the experiment to that imposed by our simulation. In Figure 2e, we present the prediction of the Wall Shear Stress (WSS) as a function of the imposed pseudo shear-rate . γ (U aver /R) along with the experimental observations of the same quantity reported in [79]. We observe that our predictions demonstrate a similar behavior with that reported by Merill et al. [79] for a hematocrit equal to 40% and R = 100 µm. In any steady, axisymmetric flow, the shear stress is given by τ rz = Jr 2 , where r is the radial position and J is the pressure gradient. Consequently, the shear stress on the wall demonstrates a linear dependence on both the pressure gradient and pseudo shear-rate. In experiments we observe a slight nonlinearity which is attributed to either the failure of the experimental techniques since WSS it difficult to be measured accurately or the loss of axisymmetry (e.g., sedimentation of RBCs). Figure 2f,g show the predictions of the relative wall shear stress (WSS rel ) of the TEVP and Casson models [63] along with the experimental measurements of Yang [80] for horse blood flow with a hematocrit equal to 40%, in tubes with R equal to 15 µm and 50 µm, respectively. According to Yang [80], the WSS rel is defined as the wall shear stress obtained by the experiment WSS exp divided by the wall shear stress exerted on the tube wall by the medium alone, i.e., blood in the absence of the RBCs WSS plasma . The latter quantity can be derived from the Poiseuille law as WSS plasma = 2 η pl U max /R. In Figure 2f,g, the TEVP and the Casson parametrizations correspond to human blood with a hematocrit equal to 45% [49]. Considering the different values of hematocrit and the difference in the blood type, the predicted WSS rel of our TEVP model demonstrates a reasonably good agreement with the experimental data for both radii. The larger deviation between TEVP predictions and experimental data at low shear-rates can be attributed to the limitations of Equations (20) and (21). Both expressions are valid only in the high shear-rate limit. On the contrary, the Casson model always overestimates the TEVP prediction for the WSS rel .

Numerical Results
This section presents the predictions of our one-dimensional steady-state simulations accounting for a continuous RBC-rich phase at the center of the tube and a plasma layer placed adjacent to the microtube wall. Arterioles are typically varying from 10 µm to 80 µm of radius, with a mean velocity of 2 mm/s. We perform a thorough parametric analysis under a wide range of microtube radii from 10 µm to 250 µm. The flow is driven by either a pressure gradient varying from 10 2 Pa/m to 10 5 Pa/m or a mean axial velocity U mean ranging between 0.1 mm/s and 5 mm/s. The mean axial velocity is related to the total flow rate Q by Q = U mean π R 2 . This model needs data from startup shear, cessation, and steady rheometry (or LAOS etc.), which can hardly be found in the open literature. Once the model parameters have been evaluated for a value of hematocrit (in this work is 45% given by McMillan's experiments [55]), our hemodynamic simulations (e.g., flows in microcirculation) should consistently follow this restriction. Since blood flows only in the core region of a vessel or a glass tube, the core hematocrit (neither the tube nor the discharge hematocrit) should be equal to 45% in all cases. Typically, blood flow investigations in which microcirculatory phenomena are included impose a constant discharged hematocrit (H d ) rather than its core counterpart (H c ).

Comparison with an Inelastic Model
It is very useful to compare the results of our model with that of the Casson inelastic model. In the following Figures, we present the predictions for both narrow and wide microtube with cross-section equal to 10 µm and 40 µm, respectively. The imposed pressure drop is J = 10 4 Pa/m. We can observe that the predicted velocity profiles for TEVP and Casson models exhibit considerable discrepancies in the velocity profile (Figure 3a,b), interfacial axial velocities (Figure 3c), and wall shear stress (Figure 3d). Regarding the profile of the axial velocity, we can see that both models capture the non-parabolic profile of the blood flow-a behavior which is consistent with experimental observations in vitro [33]. In narrow microtubes, both TEVP and Casson models predict a blunted profile (Figure 3a), which gradually shifts towards parabolic as the radius of the tube increases ( Figure 3b). The TEVP model predicts a smaller CFL thickness compared to that predicted by Casson, which in turn affects the maximum velocity prediction ( Figure 3b). As we can see from Figure 3c,d, the interfacial velocity U int and the WSS predicted respectively by the TEVP model are lower than those predicted by the inelastic model for R = 40 µm.

Effect of Proteinic Elasticity in the Plasma Layer
Since the non-elastic model of plasma is used extensively in blood flow modeling, we compare its predictions against those of our l-PTT model to estimate the error involved in neglecting the elastic contribution of plasma proteins. Although plasma exhibits much smaller elasticity than blood, as indicated by the shorter relaxation time of plasma λ pl (see Tables 1 and 3), it can induce significant differences in the stress field response. The latter is mainly affected by high shear rates or extensional phenomena. These phenomena are absent in one-dimensional flows but are substantial in more complex flows such as hemodynamics in arteriole branching, saccular aneurysms, and arterial bifurcations [42,81]. The scope of this work is limited to examining the crucial features of blood in one-dimensional flows offering a consistent model for more complex flows and not to probe extensional flows. Figure 4 illustrates the normal and shear viscoelastic stress profiles for abnormal hemodynamical conditions characterized by high pressure-gradients J. The prediction of normal and shear stress in the plasma phase is much higher than that of a Newtonian model, as depicted in Figure 4a.

Effect of Radius
Cell migration is an essential part of the mechanism behind both the development of the CFL formation and the Fåhraeus effect, as it has been stated in the previous sections. These phenomena are mainly governed by excluded volume effects [82,83] and cell-cell interactions or collisions [84,85], which are enhanced by aggregation [86,87]. Erythrocytes appear in two structural forms of individual cells and aggregated cells, the distribution of which is strongly related to the stress field applied. The dynamic equilibrium shifts toward more individual cells when the applied shear-rate increases, which affects the overall configuration of the velocity field, inducing a distinguishable difference between narrow and wide tubes. Figure 5 presents the profile of the axial velocity U z along the radial position r of the tube, for radius ranging from 10-80 µm for an applied pressure gradient equal to J = 10 4 Pa/m. Under these conditions, the dynamics of the system in narrow microtubes reveals that the velocity is more plug-like compared to wider ones, which is enhanced by the fact that in small tubes, aggregation is promoted [86]. The discontinuity in the shear rate between the fluid in the core region and that in the peripheral layer indicates the presence of the CFL. As the radius increases, the applied shear rates increase too and lead to a progressive rouleaux breakdown into individual cells promoting a more parabolic profile, as illustrated in Figure 5. In large tubes, the migration effects are less intense, and hence the CFL thickness is significantly smaller than that predicted for the narrower microtubes. When the applied shear stress is higher than the yield stress, the rouleaux network is broken, and the blood is free to flow like a liquid. This can also be evaluated by the structure parameter λ, which constitutes an indicator for the instantaneous state of blood within the tube.  Figure 6 demonstrates this thixotropic variable justifying the previous assertion. In general, blood is predicted to be in a fully structured state near the center of the tube, the extent of which is highly affected by the applied shear-rates. For very narrow tubes, such as those of R = 10 µm where the cross-section is comparable to RBC diameter, we observe two distinct zones with a sharp transition between them. The first zone is the RBC-rich region characterized by a fully structured state, and the second one is the plasma phase with λ = 0. It is evident that in this case, aggregation in the core is quite intense as blood does not demonstrate any change from its initial state and remains fully structured at steadystate conditions. The higher the imposed pseudo-shear-rate, the narrower is the plug flow region where blood depicts a fully structured form. At the center of the tube, λ is always equal to unity since the shear-rate is zero there. On the contrary, as the distance from the center increases, the shear-rate attains higher values. Consequently, the breakdown term dominates, enforcing an abrupt decrease in blood aggregates. This behavior is apparent in microtubes with wider cross-sections where the aggregation of RBCs is relatively weak. In comparison with the same blood flowing in smaller microtubes, λ attains lower values.
An indicator that demonstrates the deviation of the blunted velocity profile from a parabolic one is the parameter β. It characterizes the bluntness of the velocity profile in the core by correlating the average viscosities in the two phases given by where η b and η p are the mean viscosity of the core region and the plasma phase, respectively.
The mean shear viscosity of blood, η b , is determined as δ o τ ve,rz . γ rz rdr/ δ 0 rdr. Particularly, a value of β close to zero indicates that the velocity profile is nearly plug, while the pure parabolic profile is indicated by a β equal to unity, i.e., viscosities in the core and the plasma layer are the same and equal to the bulk viscosity. Figure 7 reveals that the bluntness parameter β increases as the diameter of the tube increases; in other words, the deviation from the parabolic profile decreases as the tube diameter increases. In the same Figure, we also illustrate the predictions of a two-phase model with a Newtonian representation for both blood and plasma. As it is expected, our predictions depict a considerable deviation from the pure Newtonian modeling, demonstrating a blunter profile for values of radius below 100 µm. The two curves converge above 100 µm of radius, predicting an almost equal bluntness parameter as the velocity profile reaches the parabolic form. Notably, our model predicts that the bluntness of the velocity profile is increased when β drops from 0.9 to 0.3 as well as the tube radius reduces from 120 µm to 10 µm. Further, for 100 µm to 150 µm in radius, the parameter β is increased from 0.9 to 0.96. For further increase of the radius, β asymptotically reaches unity. Thus, the velocity profile becomes more parabolic when the tube diameter is increased.    (Figure 8) depicts the same pattern for all examined cases. As the shear rate gradually increases along the radial position, τ ve,zz progressively increases from zero to a maximum value, which is strongly dependent on the local shear rates. The pick of normal stress component is observed at the blood/plasma interface, followed by an abrupt decrease within the plasma phase. Within the CF Layer, the normal stress is finite but relatively insignificant compared to that in the core region for the imposed pressure gradient. A considerable contribution of normal stress in the plasma phase should occur under extremely high shear rates. However, the viscoelastic behavior must not be underestimated. The only observable differentiation in normal stress predictions for different R is the maximum value at the blood/plasma interface as well as the phase change location, i.e., the CFL thickness w. The maximum value for each case at the blood/plasma interface is 21 mPa, 112 mPa, 525 mPa, 914 mPa for 10, 20, 50, and 80 µm in radius, respectively. In Figure 9, we present the spatial variation of the shear component of the viscoelastic stress tensor τ ve,rz for the same rheological conditions. The latter refers to the total shear stress applied to the system, as we assumed a negligible contribution from the solvent, highlighting a linear distribution along the radius of the microtube. Similarly to the normal stress prediction, the magnitude of τ ve,rz increases with r and R as a consequence of the appearance of higher values of shear rates. By comparing the maxima of normal and shear components, we observe that τ ve,zz is lower than τ ve,rz up to a tube radius equal to 20 µm. On the contrary, for cross-sections higher than 50 µm the normal stress contributes significantly to the total stress and overcomes the contribution of the shear component. Most of the blood constitutive modelling investigations do not present the normal stress prediction, and hence we are not able to make a comparison with other studies. Varchanis et al. [56] in their work reported a significant contribution of normal stress in simple shear tests and compared their findings with those predicted by the ML-IKH model [61], which was found to have similar behavior. The presence of normal stresses in suspensions is attributed to the intense interaction between the particles, whereas in plasma to protein stretching. Similar arguments are presented by Mall-Gleissle et al. [88] for suspensions with viscoelastic matrix fluids. The condition of whether blood is yielded or unyielded is defined through the von Mises criterion via the parameter ϕ. Figure 10a presents the spatial variation of ϕ along the radial position r for microtubes of different radii. If the stress components present in the system are large enough to satisfy the von Mises yielding criterion, the quantity ϕ acquires non-zero values indicating that the yield stress has been exceeded and blood is fluidized. On the contrary, a zero value of ϕ indicates unyielded blood, like this depicted in the case with R = 10 µm. Near the center of this tube, the stresses are insignificant, resulting in an unyielded region, the size of which depends on the radius of the tube for the same imposed pressure gradient. Clearly, ϕ = 0 in the CFL. Figure 10b presents the relaxation time χ = η t (λ) G variation along the radial position r. It is obvious that it follows similar dynamics with that depicted by the structure parameter due to our assumption that plastic viscosity depends on the instantaneous state of blood, and thus it is potentially a thixotropic property. have considered a plastic viscosity that depends on the instantaneous state of blood. Thus, as the radius increases, the microstructure of blood is more disintegrated.
(a) (b)  As we have already mentioned, the instantaneous state of the blood is defined through the parameter λ. In the current formulation, we have assumed that our model includes a stress-controlled structural parameter in the sense that the thixotropic behavior of blood is controlled by the applied stresses via the von Mises criterion. Figure 11a depicts the steadystate values of the mean structural parameter λ in the core region of the tube, as a function of its radius. We observe a continuous deconstruction of the RBCs aggregates as the radius of the tube increases due to higher shear and extensional stresses (see Figures 8 and 9). Across the whole range of the examined radii, λ experiences a reduction from λ = 1 for R = 10 µm to λ = 0.1 for R = 240 µm, with a higher decrease up to 100 µm, while beyond this point, the average structural parameter continuously seems to approach an asymptote. Figure 11b shows the variation of the fully structured fraction of the core region, where λ is equal to unity. In particular, we use the quantity r f , which stands for the percentage of the structured region out of the radius of the microtube. Also, it is used for quantifying the extension of the RBC aggregation. Blood aggregability exhibits an almost continuous reduction as blood flows in larger tubes. Indicatively, in a microtube of radius equal to R = 10 µm, we observe that fully structured material extends to 30% of R. r f follows an abrupt reduction when R = 80 µm for which its corresponding value is 6%. As the radius increases further to 240 µm, r f asymptotes to 5% approximately.
A comparison of our model results against plug flow predictions can be made by invoking the results of Gupta et al. [89]. They reported the experimental observations of the velocity profiles in microtubes. In particular, they measured the region for which the velocity profile follows a plug pattern. Figure 11c reports the experimental data of Gupta et al. [89] along with our predictions for the normalized plug flow radius r c as a function of the radius of the microtube. We observe an excellent agreement, from 40 µm up to 180 µm of radius with an overall deviation of about 4.5%. For tubes with a radius less than 40 µm, a non-monotonic behavior can be observed. This is mainly attributed to the variation of CFL with the radius of the tube, the width of which varies in a similar manner as r c . At this point, it is necessary to underline the difference of r c and r f . The former is a feature of the velocity profile and highly correlated with the thickness of the CF Layer, while the latter represents the region where the RBCs are structured in aggregated forms. In Figure 11d, we observe the prediction of the steady blood mean relaxation time χ = η t /G as a function of the microtube radius R. It exhibits a continuous reduction as the radius of the tube increases. From R = 10 µm to R = 100 µm the relaxation time experiences a steep decrease from χ = 31.4 ms to χ = 13.7 ms, while for wider radii, the mean value approaches an asymptote. This behavior comes mainly from the fact that we have considered a plastic viscosity that depends on the instantaneous state of blood. Thus, as the radius increases, the microstructure of blood is more disintegrated. One of the most important quantities in blood flow studies is the Wall Shear Stress (WSS), which is the total shear stress exerted on the microtube's wall. A proper calculation of WSS has an exceptional role, especially when blood flows in vivo because it is the stress applied on the internal Endothelial Cells (EC) surface. Vascular operations such as biochemical reactions are considerably affected by the WSS as it has been proved to be directly associated with Nitric Oxide (NO) production [43] and calcium activation in Smooth Muscle Cells (SMC) [90] by triggering the biochemical reactions that take place in vascular beds, leading to the regulation of vascular tone [91]. Although in one-dimensional blood flows, the prediction of WSS is not complex, we offer a consistent model for a proper prediction of WSS in more complicated flows. Figure 12 shows the distribution of ISS, WSS, INS, and WNS for different cross-sections of the arteriole for an imposed mean axial velocity equal to U mean = 1 mm/s. As it is expected, the wall shear stress demonstrates a non-linear dependence on arteriole radius. As we have neglected the solvent contribution in the plasma phase, the predicted WSS is related to the pure viscoelastic contribution of the proteinic phase. For a fixed mean velocity within the tube, the ISS and WSS demonstrate a gradual decrease until an asymptotic behavior is achieved. Our model predicts a significant contribution of normal stress for both interfacial and wall locations, which would be crucial for more complex blood flows such as arterial bifurcation or saccular aneurysms where the extensional phenomena are intense enough to promote the development of τ ve,zz . Both shear and normal stress components can be expressed as a function of the tube radius R in µm through a non-linear relationship of the form s n (R) = a n + b n R + c n R 2 (23) with n = 1, 2, 3, 4 refer to ISS, WSS, INS, and WNS respectively. The corresponding coefficients are presented in Table 5.    Figure 14 presents the maximum velocity U max and interfacial velocity U int as a function of the microtube radius. As blood flows under a constant pressure gradient in larger tubes, the axial velocity attains higher values, and hence the predicted maximum value demonstrates a continuous increase (Figure 14a). However, interfacial velocity does not exhibit a monotonic behavior, as illustrated in Figure 14b. This is attributed to the reduction of the aggregation effects, which enhance the transition of the velocity profiles from plug to parabolic ones. This behavior can also be justified by the predictions of velocity profiles for the two-phase blood flow in narrow tubes [27,46].
Migration effects in microcirculation are more clearly visible through the variation of the CFL thickness prediction and the evaluation of the discharged hematocrit. Figure 13 demonstrates the steady-state values of w, H d and η rel as a function of the tube radius R for J = 10 4 Pa/m and H c = 0.45. Apparently, the tendency of the erythrocytes to migrate towards the center of the tube implies an interrelation between the hemodynamical properties such as the CFL thickness, the discharged hematocrit, and the apparent viscosity. As the tube radius decreases from 250 µm down to 20 µm, the apparent viscosity (η rel ) drops due to the CFL formation next to the wall, leading to a decrease in discharged hematocrit. As the location of the CFL interface is coincident with the region of the highest shear rate within the flow, the presence of such a layer can significantly reduce the apparent viscosity. Irrespectively of the pressure drop used in the current simulation, both quantities are affected considerably by the radius R. Figure 13a demonstrates a nonmonotonic behavior of the computed cell-free layer width across the wide range of the examined cross-sections. CFL is determined by a balance between lateral migration and mass diffusion caused by cell-cell interactions. To this end, as the tube diameter decreases below a specific value, red cell migration becomes restricted due to strong interaction between RBCs in the core. This is the reason why we observe the non-monotonic behavior for the microtubes of radius below 20 µm. Although our continuum model agrees very well with the empirical observations of Pries et al. [92], it may cease to be valid once the tube radius becomes comparable to the diameter of the individual RBCs. To validate our prediction, we invoke the results of a related study in which a different constitutive model for the description of blood rheology has been used. In that work, Moyers-Gonzalez and Owens [46] conducted blood flow simulations regarding the CFL thickness under various hemodynamical conditions. Unfortunately, we are not able to completely compare our model results with those reported in their work because they simulated blood flows with different values of the discharged hematocrit and not with its core counterpart. As we impose a constant value of the core hematocrit in the current work, the only comparable CFL thickness is that predicted for a microtube with a radius of about 27 µm. At this radius, the discharged hematocrit predicted by our model is almost equal to 0.2 which is one of the examined values of the discharged hematocrits in [46]. Our model predicts w = 5.74 µm, while Moyers-Gonzalez and Owens reported a value equal to 5.68 µm, i.e., a discrepancy of only 1.04%. As the radius increases from 20 µm to 250 µm the CFL thickness follows a monotonic decrease as it is expected. When a radius lies in this range, migration phenomena become inappreciable leading to a more parabolic velocity profile and hence the apparent viscosity resembles that given by the Poiseuille law. These phenomena seem to be quite weak above a radius of 180 µm where an asymptote is approached. For a more detailed description of the interplay between the apparent viscosity and the CFL formation, the reader is referred to Appendix A. Determining the CFL thickness is of significant importance not only for a two phase blood flow simulation but also for other processes. To this end we provide a mathematical expression for the evaluation of w in µm as a function of the microtube radius R given by w(R) = w 1 + w 2 R 1 + w 3 R + w 4 R 2 (24) where w k with k = 1, 4 are adjustable parameters, which are presented in Table 6. The radius R is in µm.  Figure 13b describes the effective migration through the display of the H d as a function of the radius of the microtube, corresponding to a fixed value of H c . It is obvious that a constant H c does not imply a constnant H d for small radii. As it is expected, the increased CFL thickness in a narrower microtube indicates intensification of the migration of the erythrocytes, yielding a considerably lower value of the discharged hematocrit corresponding to a core hematocrit equal to 45%. Indicatively, from 10 µm up to 67 µm of radius, we observe a remarkable increase of H d from an extremely low value of about 0.135 to 0.44, respectively. Consequently, it is safe to claim that above a 70 µm radius, the migration effects are relatively insignificant when blood flows in vitro. function of the microtube radius. As blood flows under a constant pressure gradient in larger tubes, the axial velocity attains higher values, and hence the predicted maximum value demonstrates a continuous increase (Figure 13a). However, interfacial velocity does not exhibit a monotonic behavior, as illustrated in Figure 13b. This is attributed to the reduction of the aggregation effects, which enhance the transition of the velocity profiles from plug to parabolic ones. This behavior can also be justified by the predictions of velocity profiles for the two-phase blood flow in narrow tubes [27,46].

Effect of Pressure Gradient
Another significant model parameter is the flow rate or the applied pressure-gradient because it impacts the velocity profile [93] and, consequently, the instantaneous state of blood microstructure. In the following Figures, we demonstrate the effect of the imposed pressure gradient J on the velocity profile, the viscoelastic stress distribution, and the steady structure parameter λ. The simulation refers to a microtube with a radius of R = 20 µm while the pressure gradient ranges from 10 2 Pa/m to 10 5 Pa/m with a constant core hematocrit equal to H c = 0.45. Figure 15 demonstrates the steady-state profile of the axial velocity U z along the radial position r, when blood flows under the aforementioned rheological conditions. As it is expected, the imposed pressure gradient affects the bluntness of the velocity profile, indicating a plug-like flow for and a more parabolic-like behavior for J = 10 5 Pa/m. However, the variation is limited to profile skewness and the magnitude of the maximum velocity. Particularly, as the bluntness of the profile increases, the axial velocity field acquires higher values. However, the CFL thickness remains nearly constant, as it is implied by the experimental observations of Pries et al. [73], who argued that migration effects are not significantly affected by the applied shear rates. To further elucidate the impact of the pressure gradient on the rheological behavior of blood, we present the steady-state values of the structure parameter distribution along the radial position of the tube r ( Figure 16). For all imposed pressure gradients, the blood is initially at rest with λ |t=0 = 1. Beginning from the same state, we distinguish four different responses of the blood regarding its final microstructural configuration. When J = 10 2 Pa/m and J = 10 3 Pa/m, blood does not exhibit any change from its initial state, meaning that the stress has not exceeded the yield-stress value to disintegrate the rouleaux. Thus, in these cases, we have two distinct areas, a core region with λ = 1, and a plasma phase with λ = 0. Further increase of pressure gradient causes a partial deconstruction of rouleaux. Particularly, for a pressure gradient equal to J = 10 4 Pa/m, λ demonstrates a monotonic decrease from 1 to 0.48. It is obvious that the total stress near the center of the tube does not surpass the blood yield stress, and thus, blood remains in an unstructured state. Interestingly, with an imposition of a pressure gradient of one order of magnitude higher than the previous one, the thixotropic parameter λ experiences a steep decrease caused by the higher stress applied. In the vicinity of the blood/plasma interface, where the shear rates are high, the structure parameter approaches an asymptote at a low value. The size of the region where λ has an almost constant value is associated with the imposed pressure-gradient. The higher the J, the wider the region where λ is maintained at a constant low value. This is more obvious in the case with J = 10 5 Pa/m, where the microstructure of blood demonstrates a significant collapse. Here, blood has a constant value of about λ = 0.13 from 8 µm to 13.6 µm. In any case, blood never becomes fully   Figure 17 shows the normal viscoelastic stress distributions for different pressure gradients. An increase in J from 10 2 Pa/m to 10 5 Pa/m causes an increase in stress magnitude in both core and plasma regions. Regarding the RBC-rich central region, τ ve,zz demonstrates a continuous non-linear increase as the distance from the center is increased too. In the plasma phase, we can observe that the developed stress is quite insignificant, but when blood flows at higher velocities, as those depicted by the case with J = 10 5 Pa/m, our model predicts an observable contribution of normal stress in the plasma phase. In this case, the rheological behavior is quite reasonable as we do not expect a considerable viscoelastic contribution from plasma. However, further increase in J yields a considerable normal stress distribution along the plasma layer, which is comparable to that predicted for the RBCs. This observation is in excellent agreement with the findings in the work of Varchanis et al. [60]. They predicted a pronounced normal stress, caused by the extension of plasma proteins, especially in high shear rates. This extra elastic contribution to the rheological response of whole blood may have a significant impact on the red blood cell deformation and interaction when flowing in microtubes. The effect of the pressure gradient on the viscoelastic shear stress distribution is presented in Figure 18. As expected, τ ve,rz implies a linear dependence on the radial position r along the tube as it is the total stress, the magnitude of which is strongly associated with the imposed pressure gradient. One of the most important quantities in blood flow studies is the Wall Shear Stress (WSS), which is the total shear stress exerted on the microtube's wall. A proper calculation of WSS has an exceptional role, especially when blood flows in vivo because it is the stress applied on the internal Endothelial Cells (EC) surface. Vascular operations such as biochemical reactions are considerably affected by the WSS as it has been proved to be directly associated with Nitric Oxide (NO) production [43] and calcium activation in Smooth Muscle Cells (SMC) [90] by triggering the biochemical reactions that take place in vascular beds, leading to the regulation of vascular tone [91]. Although in onedimensional blood flows, the prediction of WSS is not complex, we offer a consistent model for a proper prediction of WSS in more complicated flows. Figure 12 shows the distribution of ISS, WSS, INS, and WNS for different cross-sections of the arteriole for an imposed mean axial velocity equal to U = 1 mm/s. As it is expected, the wall shear  Figure 19a presents the ϕ parameter for different pressure gradients. It reveals that the pressure gradient has a significant impact on the state of blood regarding its yielded or unyielded regions. Since plasma does not demonstrate plasticity, ϕ has a dual role indicating both the unyielded blood and the absence of plasticity in the annulus region. For the lowest imposed J, the predicted behavior implies that ϕ = 0 throughout the tube, i.e., unyielded blood. Increasing J to J = 10 3 Pa/m, results in partial fluidization with an unyielded region up to 6 µm of the tube and a yielded region from this point up to the phase change location δ. On the contrary, the remaining two cases depict the total fluidization of blood. Although, the transition from the unyielded J = 10 2 Pa/m to completely yielded blood J = 10 4 Pa/m is quite steep, under intermediate conditions we would observe situations with partial fluidization. For the cases with J = 10 3 Pa/m to J = 10 5 Pa/m, the region with zero ϕ is the plasma layer which is totally yielded demonstrating zero plasticity.
The developed stresses are quite sensitive to pressure gradient imposition, as it is depicted in Figure 20a,b, which demonstrates the viscoelastic stresses at the blood/plasma interface and on the microtube wall. We can observe the monotonic variation of stress when the pressure gradient is ranging between 10 Pa/m to 10 5 Pa/m for a microtube with a radius equal to R = 20 µm and a constant core hematocrit equal to H c = 0.45. From Figure 20a, we can see that ISS is always lower than the WSS. In Figure 20b, INS attains higher values than WNS does, as it is expected. The viscoelastic contribution of blood is more significant than that of pure plasma, and therefore, we cannot observe any excess of WNS.    Figure 21 demonstrates the prediction for the interfacial velocity U int and the maximum attainable velocity at the center of the tube U max under the aforementioned rheological conditions. For low-pressure gradients up to 10 3 Pa/m, the velocity profiles are almost plug and thus the U max and the U int acquire the same values. As the imposed pressure gradient increases, we observe that the two curves begin to deviate from each other up to J = 10 5 Pa/m where the velocity profile tends to be parabolic, the predictions for the U max and the U int are 5.49 mm/s and 3.33 mm/s, respectively. In Figure 22a,b, we illustrate the mean microstructural configuration λ and the plug flow radius r c , respectively. The flow conditions correspond to a pressure gradient range of 10-10 5 Pa/m, for microtube with a radius of R = 20 µm with a constant core hematocrit equal to H c = 0.45. Regarding Figure 22a, we can observe a quite expectable dependence of the mean structural variable on the pressure gradient J. For extremely low-pressure gradients, and hence low flow rates, the blood structure presents no change from the initial fully structured state up to a critical value of about J = 10 3 Pa/m. From this point onward, the microstructure of blood starts to disintegrate, and λ experiences an abrupt and continuous reduction from 1 to 0.25. From Figure 22b, we can observe that the prediction of the normalized plug-flow radius r c as a function of the imposed pressure gradient is quantitatively similar to that presented for the mean value of λ (Figure 22a). Particularly, we observe two distinct responses of blood regarding the normalized region for which the blood velocity presents a plug profile. Initially, for low-pressure gradients, the velocity profile is plug throughout the core region, and hence the r c is constant and about 0.68. As the pressure gradient increases from J = 10 3 Pa/m to J = 10 5 Pa/m the velocity profile gradually obtains a more parabolic pattern, and thus the plug flow region is reduced dramatically to r c = 0.32, which corresponds to the extreme pressure gradient value of our simulations.

Conclusions
This work presents steady-state predictions of blood flow in microtubes incorporating a two-layer fluid model consisting of an outer annulus filled with plasma and an inner RBCrich core. Two constitutive formulations are employed, a thixotropic elasto-visco-plastic (TEVP) model for the blood core and a viscoelastic (linear-PTT) for the plasma, offering a holistic formulation for the accurate two-phase rheological behavior of blood/plasma in microtubes. The blood model was parameterized by Giannokostas et al. [4], using the hemorheological data of donors with 45% core hematocrit [55]. This is the first time that a TEVP model is validated in non-rheometric flows or flows with spatial variation of the velocity and the stress fields. The model is accompanied with CFL existence via empirical expressions of Pries and Secomb [92], and we have proposed a new simplified form (Equation (20)) compared to that provided by Pries et al. [74] for the description of the apparent viscosity. We present the predictions of the new analytical relationship of the CFL thickness with the radius of the tube. Comparison with experiments enables us to claim that our formulation constitutes a consistent model capable of successfully addressing complex blood flow behavior in narrow microtubes. In particular, we compared our predictions with experimental data of velocity profiles [76], cell-free layer thickness [5,77,78], as well as wall shear-stress measurements [79], obtaining an excellent agreement. We offered a thorough parametric analysis examining the impact of the microtube diameter as well as the pressure gradient on velocity profiles, normal and shear viscoelastic stresses, and the structural state of blood. Among the findings, we can summarize those that are of greater importance. Elasticity is found to have a significant impact on both RBCs and plasma phases. The normal stress field would be important in flows with a high flow rate or where extensional phenomena are quite intense. Additionally, elasticity provoked considerable discrepancies between TEVP and Casson predictions regarding the flow rate and the interfacial velocity. It was shown that the radius of the tube and the pressure gradient have a significant impact on the configuration of the internal structure of blood. For narrow microtubes, the aggregation was intense, the blood state was almost fully structured, leading to a blunted velocity profile. In addition, the relaxation time of the blood was also affected by the radius of the tube depicting higher values for narrower microtubes. As the location of the CFL is coincident with the region of the highest shear rate within the flow, the presence of such a layer was found to reduce the apparent viscosity significantly. Our model encompassed migration through a hemodynamical mathematical description, which was verified by experimental evidence and hence was able to predict the variation in local hematocrits accurately. For a small microtube, the model predictions implied a considerable reduction in discharged hematocrit due to migration effects, which in some cases is up to 70% and was proved to be invariant with the pressure gradient. The pressure gradient had a significant impact on the velocity profile and the mean microstructure of blood. We show a transition from totally unyielded to fully yielded blood with intermediate states of partial fluidization. The size of each region is linked with the imposed pressure gradient. The reduction of the blood structure λ within a reasonable range of pressure gradient was almost 80%.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The reduced form of the z-component of Equation (1) in cylindrical coordinates is given by where the first equation refers to the blood-core region, while the second one to the blood plasma in the annulus. U z (1) and U z (2) are the velocities of the core region and plasma in the z-direction, respectively, τ ve,rz (1) and τ ve,rz (2) are the shear components of the viscoelastic stress tensor in core and plasma phases, respectively, ∂P ∂z is the pressure drop in the z-direction, and δ is the location of the blood/plasma interface. ρ (1) and ρ (2) are the densities of the blood and plasma equal to 1060 Kg/m 3 and 1025 Kg/m 3 , respectively. It must be pointed out that a two-phase representation is mandatory not only for the distinct rheological description of each phase but also for the prediction of the interface position and, consequently, the computation of the cell-free layer (CFL) thickness w = R − δ. We also impose two boundary conditions for the velocity, a zero gradient at the center of the tube and a no-slip condition on the wall given by At the blood/plasma interface, we impose continuity conditions for both stress and velocity as τ ve,rz As we have a moving boundary problem, due to the unknown location of the blood/plasma phase, the equations are solved in a transformed computational domain with fixed boundaries, and the solution is transformed back to the physical domain. The reduced one-dimensional forms of the constitutive law for the shear τ ve,rz (1) and normal τ ve,zz (1) components of the stress tensor in the core are given by

Appendix B
The original form of the empirical description of the in vitro apparent viscosity experiments has been given by Secomb & Pries [94]. They assembled results from several studies, and they developed algebraic relationships to describe the dependence of the relative apparent viscosity η rel (the ratio of the apparent viscosity to the suspending medium viscosity) on tube diameter and discharged hematocrit as where D is the diameter of the tube in µm, and C a dimensionless function of D. The apparent viscosity for a typical discharged hematocrit in humans H d = 0.45, η 0.45 , is given by Equation (A11). The variable C (Equation (A12)) describes the type of dependence of the relative viscosity on the hematocrit, which is approximately linear for diameters up to about 8 µm, but shows a highly nonlinear dependence at large diameters. Finally, the expression of relative apparent viscosity η rel is given by Equation (A13) based on the aforementioned quantities. This set of equations is much more complex compared to our proposed expression for the relative apparent viscosity given by Equation (20). Additionally, the dependence of parameter C on D 12 is quite unreliable, and prone to large errors. The prediction of the relative apparent viscosity with respect to the tube diameter for different values of discharge hematocrit, according to our equation, is illustrated in Figure A1a along with the corresponding experimental data given by Pries and Secomb [3]. The trend of the apparent viscosity to decrease with reducing the tube diameter continues down to about 7 µm. At smaller diameters, the relative apparent viscosity rises rapidly as the diameter approaches a critical minimum diameter, which is about 3 µm. Although highly deformable, RBCs are subject to constraints of constant volume and almost constant surface area. These constraints prevent passage of intact cells through tubes narrower than this critical diameter. As diameter increases, the relative viscosity attains higher values which tend to asymptotically reach the relative viscosity given by Poiseuille's law for a Newtonian fluid. Figure A1b illustrates the predictions of cell-free layer thickness in µm with respect to the tube diameter and for different imposed core hematocrits. As our analysis is valid only for a certain value of core hematocrit that of 45%, we provide the predictive capability of Equation (20) using the Casson model (Equation (12), Table 2) [63] for hematocrit values between 10% to 60%. The predictions are compatible with experimental indications exhibiting a continuous decrease as the diameter of the tube increases. At high cross-sections, the CFL thickness asymptotically reaches a small but finite constant value, depending on the imposed core hematocrit. Figure A1. (a) The predictions of our proposed expression (Equation (20)) for relative apparent viscosity against experimental data reported by Pries and Secomb [73]. (b) Cell-free layer variation with respect to the tube radius and different values of core hematocrit based on Equation (20) for two different constitutive models applied in the core-annular flow problem described in Section 2.

Appendix C
Our model is inherently transient, in the sense that there exists an infinitude of steady states, which depend in a continuous manner on the initial conditions and state. Even if it is just the steady-state that is sought, one cannot simply discard the time derivatives; the steady-state must be obtained by a time integration with appropriate initial conditions. The initial conditions imposed for the solution of the governing equations are those corresponding to initially unyielded and unperturbed blood. To this end, we assume a zero stress field ( τ ve,rz (i) t=0 = τ ve,zz (i) t=0 = 0) and velocity field ( U z (i) t=0 = 0), while the structure variable is assumed to be unity ( λ| t=0 = 1). Figure A2. Time evolution of (a) axial velocity U z , and (b) structure parameter λ from fully structured state to steady-state.
The evolution of the profile of the axial velocity U z is presented in Figure A2a for five-time instances. During the transient simulation, the velocity field is gradually formed from zero state to a blunted profile at a steady-state, as it is indicated by the respective symbols in the graph. As time increases, the velocity attains a higher flow rate with larger bluntness and lower CFL thickness. Figure A2b illustrates the evolution of the structure parameter λ for a radius equal to R = 50 µm and a pressure drop equal to J = 10 4 Pa/m to achieve steady-state conditions. The initial state for λ is λ| t=0 = 1 so that the blood is not affected by any prior history effect. The stresses develop rapidly, leading to a continuous breakdown of blood rouleaux and consequently promote a decrease of λ. Gradually, more of the material acquires values below unity until the steady-state is achieved, and a small region near the center of the tube still retains its original state of the fully structured configuration. The size of this region is clearly dependent on the imposed flow conditions and especially on the magnitude of the pressure gradient and, consequently, on the applied shear-rates. The higher the developed shear-rates, the narrower is the region of a fully structured state. Apparently, as shear-rate attains higher values, the bulk configuration of λ acquires lower values in the sense that the blood is getting softer, and the microstructural destruction is higher. Irrespective of the intensity of shear rate, λ never reaches 0, but asymptotes to zero for extreme values of the shear-rate. In other words, the blood never becomes fully unstructured. Additionally, the fully structured region is always present, no matter how shear-rate increases.