Introduction

Traditional treatments have proven impotent as drug molecules act non-specifically by simply diffusing freely throughout the body, leading to undesirable side effects, and are detrimental in achieving the required doses for positive outcomes. On the other hand, medical applications of nanotechnology have proven proficient enough to deliver protagonistic clinical breakthroughs in the diagnosis and treatment of several diseases using drug-carrying particles (DCPs) (Irvine and Dane 2020). For example, very recently a lipid nanoparticle carrying messenger RNA has been approved by the US Food and Drug Administration for its first clinical trial aiming to treat the SARS-CoV-2 coronavirus (Cohen 2020). The use of DCPs provides promising and more effective alternative treatments, compared to traditional ones, to fight several diseases of the circulatory system, such as atherosclerosis (when the vascular wall thickens due to deposition of dead cells or cholesterol on the vascular wall), thrombosis (formation of a blood clot inside a vessel), and cerebrovascular amyloid angiopathy (when amyloid proteins build up on the walls of brain arteries) (Agyare and Kandimalla 2014).

The execution of in silico trials (the execution of clinical human trials in computers, i.e., on virtual patients) for prognostic purposes could promote the reduction of the time required to perform and the cost of design experiments, such as preclinical animal ones, that are a prerequisite for the approval of a new drug (Moradi Kashkooli et al. 2021). In addition, it could facilitate the tailor-design of each drug for optimum delivery to the infected areas (Stillman et al. 2020). More recently, in silico trials using sophisticated algorithms have emerged as a promising methodology with potential benefits for expanding computer modeling to drugs and other biomedical products (Viceconti et al. 2016).

While the blood vessel itself is seldom the target, almost all nanomedicines are injected into the microcirculation system and then make their way towards the targeted vascular area or tissue (e.g., the atherosclerotic plaque or the tumor). Several factors render such a task quite challenging. The transport of DCPs is affected by the local hemodynamics. In large vessels, the flow may be considered homogenous, and thus the transport of DCPs will strongly depend on the local velocity. On the other hand, the flow in microvessels, like capillaries, is inhomogeneous because of the migration of red blood cells (RBCs) towards the vessel centerline (Goldsmith et al. 1989) and the migration of DCPs towards the vascular wall (Tilles and Eckstein 1987; Eckstein et al. 1988). To achieve optimum-targeting DCPs, they must first successfully travel from the blood midstream to the RBC-free layer close to the vessel wall. As such, the physical properties of DCPs that result in optimal margination to the vessel wall must be considered. For example, a recent study (D’Apolito et al. 2015) has shown that the margination phenomenon is size- and shape-dependent as larger spherical particles and more elongated particles of the same size are more effectively marginated both in in vitro and in vivo experiments. Also, the margination effect is hematocrit-dependent since it occurs faster as the hematocrit increases (Fitzgibbon et al. 2015) (the hematocrit, \(H\), is the volume fraction of RBCs which, under physiological conditions, is between 42 and 47%). It would thus be highly desirable if an in silico study could investigate the transport of DCPs in vascular environments. This would allow for a faster design of suitable DCPs, reduce the need to perform numerous trials in animals, and alleviate all constraining factors associated with it, and/or minimize the error when conclusions of laboratory studies drawn on animals are extended to humans. Given the multitude of length and time scales involved, it is necessary to employ a multiscale modeling approach (Eckmann et al. 2020).

Although numerous coarse-grained approaches to simulate the flow behavior of particles in blood have been employed up to date, such as the immersed finite elements method (Tan et al. 2012), the lattice Boltzmann-immersed boundary method (Liu et al. 2018, 2019), and the dissipative particle dynamics method (Müller et al. 2016), these may be computationally expensive enough to forbid an elaborate parametrization of the complicated vascular transport of DCPs. In such approaches, the RBCs are modeled as flexible thin membranes enclosing a fluid, with the surface consisting of a triangular network, whereas DCPs are modeled as solid particles. Thus, numerical simulations, i.e., simulations wherein the continuity, momentum, and constitutive equations are solved numerically, usually using the finite elements method (FEM), are the only other tool available to address this issue. In such simulations, the rheological behavior of the fluid is described by a constitutive equation. Therefore, to be able to perform the aforementioned in silico trials, a reliable and accurate constitutive model is needed that properly and accurately addresses the rheology of nanofluids comprising of particles immersed in blood, which we shall henceforth refer to as nanoblood.

In most theoretical studies, the methodology employed is too simplistic. Decuzzi and co-workers (Decuzzi et al. 2005; Decuzzi and Ferrari 2008) analyzed the margination of a particle circulating in the bloodstream by considering the contribution of buoyancy, and of steric interactions between particles and the vascular wall. However, they considered blood to be a Newtonian fluid, i.e., with a constant viscosity. Hossain et al. (2013) attempted to predict the deposition of particles in a patient-specific arterial tree, but also considered blood to behave as a Newtonian fluid. The same approach was followed by others as well (Shaw et al. 2014; Elnaqeeb et al. 2019; Abdelsalam et al. 2020; Ahmed et al. 2020; Tripathi et al. 2021) who considered nanoblood to behave as a Newtonian fluid with a viscosity that depends on the particle’s volume fraction using the Brickman model. Gentile et al. (Gentile et al. 2008; Gentile and Decuzzi 2010) analyzed mathematically the longitudinal transport of particles in blood vessels by considering blood to be rheologically described as a Casson fluid and derived an expression for the effective longitudinal diffusion. Their analysis showed (Gentile et al. 2008) that rheology plays a decisive role in the transport of particles in the vascular network. More recently, Jafarzadeh et al. (2020) studied the accumulation of particles once they are injected into an abdominal aortic aneurysm for the duration of a cardiac cycle under the imposition of a pulsatile bloodstream. They considered blood to behave as a shear-thinning fluid (by using the phenomenological model proposed by Ameenuddin et al. (2019) according to which the shear viscosity depends on both the shear rate and the hematocrit) to study the hemodynamical flow in the vicinity of an abdominal aortic aneurysm. They noted that as the \(H\) and the size of the particle decrease the accumulation of the DCPs in the dilatation part of the artery is greater. Abdelsalam and co-workers have also used models that bear a shear-thinning behavior, despite being phenomenological, such as the Sutterby (Abdelsalam and Bhatti 2020; Bhatti et al. 2020) and Eyring-Powell (Abdelsalam et al. 2021) models. In all cases, the employed phenomenological expression refers to the viscosity of neat blood and not of the nanofluid as a whole. Dubey et al. (2020) performed computational fluid dynamics (CFD) simulations, using the FEM, of blood flow through a diseased artery with a stenosis followed by an aneurysm. They adopted two different models to address the rheological behavior of blood: they employed the Casson viscoplastic model (Casson 1959; Mitsoulis 2007) in the arterial core region, and the viscoelastic Sisko model in the peripheral (close to the arterial wall) region. This is necessary as due to the accumulation of RBCs in the core region, the fluid there behaves as a semi-solid; hence, the velocity of the fluid is constant (i.e., zero velocity gradient). They noted that the velocity at both the core and peripheral regions increases as the size of the particles decreases irrespective of whether nanoblood flows through a stenosis or through an aneurysm. A similar study was published by Shaw (2020) who employed the Herschel–Bulkley viscoplastic model (Herschel and Bulkley 1926; Mitsoulis 2007) in the arterial core region, whereas in the peripheral region nanoblood was considered to behave as a Newtonian fluid with a viscosity that depends not only on the particle volume fraction but also on the particle and the blood vessel diameters. Finally, Ponalagusamy and Priyadharshini (2019) studied the unsteady flow of the transport of magnetic nanoparticles through an artery with stenosis in the presence of a magnetic field by considering blood as a Herschel–Bulkley viscoplastic model and noted that when the yield stress and/or the particle concentration increases, then the wall shear stress and the flow resistance increase as well.

The majority of the abovementioned studies have selected to employ in their analysis simple rheological models. This necessity emanates from the need to have an analytical expression for the velocity profile within the vessel. Although considering blood to behave as a Newtonian fluid is an adequate approximation in arteries, it completely fails to capture the correct hemodynamics in capillaries (Sriram et al. 2014). Overall, all these works completely ignore the molecular characteristics of blood and their effect on its rheological behavior. This certainly raises significant limitations to our understanding of the rheological footprints of nanoblood. It should also be stressed that all the abovementioned studies have not derived a constitutive model for nanoblood, whose predictions should be compared against rheological data, but merely employed a previously proposed phenomenological model without first benchmarking its validity. The importance of correctly predicting, via mathematical modeling, the rheology of nanoblood is exemplified by the fact that the highest accumulation of DCPs is noted in regions of disturbed flow, e.g., branch points and curvatures in the vasculature (Gomez-Garcia et al. 2018), and at the shoulder regions of the atherosclerotic plaque (i.e., the junction between the plaque and the histologically healthy part of the vessel wall) (Peters et al. 2009); in both cases, the stresses there are minimum.

In this work, we will employ non-equilibrium thermodynamics (NET) to develop a sophisticated mathematical model addressing the rheological behavior of nanoblood. To accomplish this, we will develop the model using the generalized bracket formalism (Beris and Edwards 1994) of NET (Beris and Edwards 1994; Grmela and Öttinger 1997; Öttinger and Grmela 1997; Edwards et al. 2003; Öttinger 2005), by means of which several microstructured systems have been addressed up to date, such as, but not limited to, immiscible complex fluids (Edwards and Dressler 2003; Dressler and Edwards 2004; Dressler et al. 2008; Grmela et al. 2014; Mwasame et al. 2017), polymer melts and solutions (Beris and Edwards 1994; Grmela and Öttinger 1997; Öttinger and Grmela 1997; Öttinger 2005; Stephanou et al. 2009, 2016, 2020b), polymer nanocomposites (Rajabian et al. 2005; Eslami et al. 2010; Stephanou et al. 2014; Stephanou 2015), micellar systems (Germann et al. 2013; Stephanou et al. 2020a), blood (Tsimouri et al. 2018; Stephanou 2020; Stephanou and Tsimouri 2020), drilling fluids (Stephanou 2018), and thixotropic fluids (Stephanou and Georgiou 2018). The use of a NET formalism has the compelling advantage, over other approaches, that the constitutive model as a whole is guaranteed consistency with the laws of thermodynamics (as extended for beyond equilibrium systems). Furthermore, it provides a systematic way of developing and self-consistently coupling the various components (e.g., in our present work, RBCs and particles). Our ultimate goal in the present work is to provide an avenue to conduct in silico simulations for the testing of DCPs in virtual patients that will allow the tailor-design of DCPs for treating various diseases. This would ultimately reduce the size and duration of, as well as enable the design of more effective, human clinical trials, and lower both the development and market costs of new DCPs.

The paper is structured as follows: in “Model derivation”, the new model is introduced, whereas in “Results and discussion” we present the model predictions along with a comparison with available rheological experimental data for both spherical and non-spherical (rod-like) particles over a wide range of volume fractions (in the dilute regime) and shear rates. The paper concludes with “Conclusions,” where we elaborate on the significance of our work and future plans are highlighted and discussed.

Model derivation

The vector of state variables

Throughout this work, we consider an isothermal and incompressible flow, meaning that both the total mass density, ρ, and the entropy density (or temperature) are excluded from the vector of state variables. At first, to describe the rheology and microstructure of nanoblood, we follow previous works (Stephanou 2020, 2021; Stephanou and Tsimouri 2020) and consider RBCs as emulsions with a droplet morphology by using a constrained contravariant second-rank tensor, \(\overline{\mathbf{S} }\), such that det \(\overline{\mathbf{S} }\) is equal to the squared volume of a RBC. We also define the dimensionless tensor \(\mathbf{S}=\overline{\mathbf{S} }/{\left(\mathrm{det}\overline{\mathbf{S} }\right)}^{1/3}\) so that \(\mathrm{det}\mathbf{S}=1\). We consider RBCs as spheroids (i.e., ellipsoids of which the lengths of two axes are the same), as it is the simplest possible approximation to the actual shape of RBCs, i.e., biconcave disks. At equilibrium, the tensor \(\overline{\mathbf{S} }\) is equal to \(\mathrm{diag}[{a}^{2},{a}^{2},{c}^{2}]\) where \(a\) and \(c\) are the principal semi-axes of the RBC (see Fig. 1 of Stephanou and Tsimouri (2020)), whose determinant is \(\mathrm{det}{\overline{\mathbf{S}} }_{eq}={a}^{4}{c}^{2}={\left(\frac{3}{4\pi }{V}_{RBC}\right)}^{2}\); thus, \({\mathbf{S}}_{eq}={\overline{\mathbf{S}} }_{eq}/{\left(\mathrm{det}{\overline{\mathbf{S}} }_{eq}\right)}^{1/3}=\mathrm{diag}[{\left(a/c\right)}^\frac{2}{3},{\left(a/c\right)}^\frac{2}{3},{\left(c/a\right)}^\frac{4}{3},]\) (Stephanou 2020; Stephanou and Tsimouri 2020). If we then consider the typical values \(2a\) = 7.5 μm and \({V}_{RBC}=\) 90 μm3, we must consider \(2c=\) 3.06 μm (Stephanou and Tsimouri 2020). The mere use of the conformation tensor \(\mathbf{S}\) allows us to consider the deformability of RBCs but not their aggregation in normal blood, which leads to the exhibition of a yield stress at sufficiently low shear rates as a result of the network formed. We have recently proposed (Stephanou 2020) that to accomplish this we need to employ one additional scalar structural variable, \(\lambda\), to properly characterize the network formed by RBCs. When the network is fully formed, then \(\lambda\)=1, while in a completely broken state, i.e., when only single RBCs remain, it vanishes (Mewis and Wagner 2009; Stephanou and Georgiou 2018). Next, we consider particles as rigid spheroids; to describe the average particle orientation, we consider the orientation tensor \(\mathbf{a}\) defined (Stephanou et al. 2014; Stephanou 2015) as \(\mathbf{a}=\int \mathbf{n}\mathbf{n}\psi \left(\mathbf{n},\mathbf{r},t\right)dV\) where \(\mathbf{n}\) is the particle’s director vector and \(\psi \left(\mathbf{n},\mathbf{r},t\right)\) is the orientational distribution function for the vector \(\mathbf{n}\) at position \(\mathbf{r}\) and time \(t\). The particle volume fraction is given via \(\phi ={n}_{a}{v}_{a}\) where \({n}_{a}\) is the number density of particles, and \({v}_{a}=\left(\pi /6\right){d}^{2}l\) their volume, with \(l\) the length of the spheroid (for oblates their thickness, whereas for prolates their length) and \(d\) its diameter (see Fig. 1 of Stephanou (2018)). The orientation tensor is constrained to have a constant unit trace due to the rigidity of the particles considered. At equilibrium, the orientation distribution function of particles is isotropic; thus, \({\mathbf{a}}_{eq}=\mathbf{I}/3\), where \(\mathbf{I}\) is the unit tensor. Finally, we consider the momentum density \(\mathbf{M}\) as the hydrodynamic variable, so that overall the vector \(\mathbf{x}\) of state variables is expressed as \(\mathbf{x}=\left\{\mathbf{M},\mathbf{S},\mathbf{a},\lambda \right\}\).

Fig. 1
figure 1

Representative model predictions for the dimensionless steady-state a shear viscosity, b first normal stress coefficient, and c second normal stress coefficient, as a function of the dimensionless shear rate \(\mathrm{Ca}\). The results are shown for various values of the model parameters associated with the particles, while keeping \(H=40\%,\) \(\varepsilon \to \infty , \xi =0.01,\) and \({f}_{0}=1\)

Full size image

The resulting evolution equations

The complete derivation is presented in the Supplementary Information. The final evolution equations for each structural variable are as follows:

$$\frac{\partial \mathbf{S}}{\partial t}=-\mathbf{u}\bullet \nabla \mathbf{S}+\mathbf{S}\bullet \nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{T}\bullet \mathbf{S}+\frac{1}{3}\mathbf{S}\mathbf{u}\bullet \left({\mathbf{S}}^{-1}:\nabla \mathbf{S}\right)-\frac{\xi }{2}\left(\mathbf{S}\bullet \dot{{\varvec{\upgamma}}}+\dot{{\varvec{\upgamma}}}\bullet \mathbf{S}\right)-\frac{3\left(1-\lambda \right)}{{I}_{2}^{C}{\tau }_{S}}\left[-\left(\mathbf{S}\bullet {\mathbf{S}}_{eq}^{-1}\bullet \mathbf{S}-\frac{{\rm I}_{1}^{C}}{3}\mathbf{S}\right)-\frac{2{\rm I}_{2}^{C}}{3}\left({\mathbf{S}}_{eq}-\frac{{\rm I}_{2}^{C}}{3}\mathbf{S}\right)\right],$$
(1a)
$$\frac{\partial \mathbf{a}}{\partial t}=-\mathbf{u}\bullet \nabla \mathbf{a}+\mathbf{a}\bullet \nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{T}\bullet \mathbf{a}+\frac{\theta -1}{2}\left(\mathbf{a}\bullet \dot{{\varvec{\upgamma}}}+\dot{{\varvec{\upgamma}}}\bullet \mathbf{a}\right)-\theta \mathbf{a}\left(\mathbf{a}:\dot{{\varvec{\upgamma}}}\right)-\frac{{\Lambda }_{0}^{\mathrm{a}}}{{\tau }_{\mathrm{a}}}\left(\mathbf{a}-\frac{1}{3}\mathbf{I}\right),$$
(1b)
$$\frac{\partial \lambda }{\partial t}=-\mathbf{u}\bullet \nabla \lambda -\frac{1}{{\tau }_{\lambda }}\mathrm{ln}\lambda -\frac{2\lambda }{\mathrm{tr}\mathbf{S}}{\left(\nabla \mathbf{u}\right)}^{T}:\mathbf{S},$$
(1c)

where \(\nabla \mathbf{u}\) is the velocity gradient tensor (\({\mathbf{\rm X}}^{T}\) is the transpose of \(\mathbf{\rm X}\)), \(\dot{{\varvec{\upgamma}}}\equiv \nabla \mathbf{u}\)+\({\left(\nabla \mathbf{u}\right)}^{T}\) is the rate-of-strain tensor, \(\xi\) is the slip/non-affine parameter, \({I}_{i}^{C}\) is the ith invariant of tensor \(\mathbf{C}\), for which \({\mathbf{C}}_{\mathrm{eq}}=\mathbf{I}\), such that \(\mathbf{S}=\mathbf{C}\bullet {\mathbf{S}}_{\mathrm{eq}}\Rightarrow \mathbf{C}=\mathbf{S}\bullet {\mathbf{S}}_{\mathrm{eq}}^{-1}\) (see Stephanou and Tsimouri (2020)), \(\theta =\left({r}_{e}^{2}-1\right)/\left({r}_{e}^{2}+1\right)\) with \({r}_{e}=l/d\) the particle’s aspect ratio, \({\Lambda }_{0}^{\mathrm{a}}\) is a numerical constant, and \({\tau }_{S},{\tau }_{\mathrm{a}}\), and \({\tau }_{\lambda }\) are characteristic relaxation times. Here, we consider \({\tau }_{\mathrm{a}}={\tau }_{\mathrm{a}0}\phi\) as the characteristic particle rotational time, where \({\tau }_{\mathrm{a}0}\) is the particle rotational time at infinite dilution and \(\phi\) the particle volume fraction (Stephanou 2018), \({\tau }_{S}\) the characteristic relaxation time related to the deformability of RBCs (Stephanou 2020), and \({\tau }_{\lambda }\) the characteristic relaxation time for regeneration once the flow field has ceased. Note that, following Doi and Edwards (1988) and Larson (1999), we consider the particle rotational time to increase with the volume fraction since crowding leads to an increase of the time needed to complete one rotation. Also, we consider \({\tau }_{\lambda }=\varepsilon {\tau }_{S}\), where the parameter \(\varepsilon\) quantifies the relative importance between the regeneration/buildup and the flow-induced breakup of the RBC network formed at low shear rates or at stasis (Stephanou 2020). When \(\varepsilon \gg 1\), the time needed for regeneration is much larger than the one needed for breakup, and thus, the scalar structural variable vanishes, \(\lambda \to 0\), in which case a yield point is not predicted. The thermodynamic admissibility of the model is provided in the Supplementary Information. Finally, the stress tensor is given as

$${\varvec{\sigma}}=\left[1+f\left(\phi \right)\right]H\Gamma \left[2\left(1-\xi \right)\left({\rm I}_{1}^{C}\mathbf{C}-\mathbf{C}\bullet \mathbf{C}-\frac{2{\rm I}_{2}^{C}}{3}\mathbf{I}\right)-\frac{\lambda \mathrm{ln}\lambda }{\mathrm{tr}\mathbf{S}}\mathbf{S}\right]+{n}_{\mathrm{a}}{k}_{B}T\theta \left(3\mathbf{a}-\mathbf{I}\right)+{\eta }_{s}\left[1+f\left(\phi \right)\right]P\left(H,{\lambda }_{\eta }\right)\dot{{\varvec{\upgamma}}}+{\eta }_{s}\phi \left[\frac{B\left({r}_{e}\right)}{3}\left(\mathbf{a}\bullet \dot{{\varvec{\upgamma}}}+\dot{{\varvec{\upgamma}}}\bullet \mathbf{a}\right)+\frac{\overline{A }\left({r}_{e}\right)}{9}\mathbf{a}\left(\mathbf{a}:\dot{{\varvec{\upgamma}}}\right)\right].$$
(2a)

Here, \({\eta }_{s}\) is the viscosity of the Newtonian medium, e.g., plasma or other solvents, \({k}_{B}\) is Boltzmann’s constant, \(T\) is the absolute temperature, and \(\Gamma\) is the surface energy density, which can be shown to be given as \(\Gamma \left(\phi \right)=3{{f}_{0}\left(\phi \right)\eta }_{s}/{\tau }_{s}\)(Stephanou 2020; Stephanou and Tsimouri 2020). Note that the parameter \({f}_{0}\) has been considered to be a function of \(\phi\) as the surface energy density is expected to be affected by the particle concentration. A comparison against experimental data depicted in “Comparison with experimental data” will point out that \({f}_{0}\left(\phi \right)\), at least for suspension of rod-like nanoparticles in non-aggregating blood, is a decreasing function of \(\phi\). Given the uncertainty, in lieu of theoretical considerations, we will simply consider \({f}_{0}\left(\phi \right)\) as an adjustable parameter. Furthermore, following Stephanou (2020),

$$P\left(H,{\lambda }_{\eta }\right)={\left(1-\frac{H}{{H}_{m}}\right)}^{-T\left({\lambda }_{\eta }\right){H}_{m}}$$
(2b)

where \({H}_{m}\) is the maximum packing hematocrit, which is expected to have a value close to 0.72 (Stephanou 2020). Note that this expression as \(H\to 0\) boils down, as it should, to Taylor’s expression who showed that the viscosity of an emulsion, in the dilute regime, is given as \(P\left(H,{\lambda }_{\eta }\right)=1+T\left({\lambda }_{\eta }\right)H\), where \(T\left({\lambda }_{\eta }\right)=\left(5{\lambda }_{\eta }+2\right)/\left[2\left({\lambda }_{\eta }+1\right)\right]\) (Taylor 1932) and \({\lambda }_{\eta }\) is the ratio between the internal and external fluid viscosities. It is known that RBCs are very deformable since they are composed of a thin elastic membrane (lipid bilayer) enclosing the cytoplasm (a hemoglobin solution) (Yilmaz and Gundogdu 2008) which has a higher viscosity (equal to about 3–10 mPa s) than that of the surrounding blood plasma. Thus, the ratio between the internal (cytoplasm) and external (plasma) viscosities, \({\lambda }_{\eta }\), is between 2.5 and 8.3. As will be shown later, when we neglect the elasticity of RBCs (neglecting their deformability) under very low RBC concentrations, \(H\to 0\), and considering spherical rigid particles, Eq. (2a) boils down to the viscosity of a dilute suspension of spherical particles in a Newtonian solvent and generalizations thereof, i.e., \({\varvec{\sigma}}={\eta }_{s}\left[1+f\left(\phi \right)\right]\dot{{\varvec{\upgamma}}}\); as \(\phi \to 0\) this should boil down to Einstein’s expression (Einstein 1906, 1911; Woolard et al. 1928) by considering \(f\left(\phi \right)=\frac{5}{2} \phi\), which is accurate up until \(\phi \approx 0.02\), whereas Batchelor and Green (Batchelor and Green 1972a, b) extended Einstein’s equation to higher volume fractions (up until \(\phi \approx 0.1\)) by considering \(f\left(\phi \right)=\frac{5}{2} \phi +\overline{c}{\phi }^{2}\), with \(\overline{c }\approx 6.2\). Although other expressions could be employed (see Eq. (26c) of Stephanou et al. (2014)), these suffice for the present application since the volume fraction of DCPs in nanoblood is not expected to exceed 10%. Finally, the parameters \(\overline{A }\left({r}_{e}\right)\) and \(B\left({r}_{e}\right)\) are functions of the aspect ratio (see Eq. (2) in Stephanou (2018)), which vanish for spherical particles. Also, in the case of a dilute suspension of spheroidal particles then \(f\left(\phi \right)=C \phi\) where \(C\) is also a function of the aspect ratio (see Eq. (2) in Stephanou (2018)); in the case of spherical particles, \(C=5/2\). By noting that \({n}_{\mathrm{a}}{k}_{B}T=\phi \left({k}_{B}T/{v}_{\mathrm{a}}\right)=\phi X\left({\eta }_{S}/{\tau }_{\mathrm{a}0}\right)\) where \(X=6/\left(\pi {r}_{e}\right)\), then

$${\varvec{\sigma}}=3\left[1+f\left(\phi \right)\right]{f}_{0}H\left({\eta }_{s}/{\tau }_{s}\right)\left[2\left(1-\xi \right)\left({\rm I}_{1}^{C}\mathbf{C}-\mathbf{C}\bullet \mathbf{C}-\frac{2{\rm I}_{2}^{C}}{3}\mathbf{I}\right)-\frac{\lambda \mathrm{ln}\lambda }{\mathrm{tr}\mathbf{S}}\mathbf{S}\right]+\phi X\left({\eta }_{s}/{\tau }_{\mathrm{a}0}\right)\theta \left(3\mathbf{a}-\mathbf{I}\right)+{\eta }_{s}\left[1+f\left(\phi \right)\right]P\left(H,{\lambda }_{\eta }\right)\dot{{\varvec{\upgamma}}}+{\eta }_{s}\phi \left[\frac{B\left({r}_{e}\right)}{3}\left(\mathbf{a}\bullet \dot{{\varvec{\upgamma}}}+\dot{{\varvec{\upgamma}}}\bullet \mathbf{a}\right)+\frac{\overline{A }\left({r}_{e}\right)}{9}\mathbf{a}\left(\mathbf{a}:\dot{{\varvec{\upgamma}}}\right)\right]$$
(3)

or when defining the dimensionless stress via \(\stackrel{\sim }{{\varvec{\sigma}}}=\left({\eta }_{S}/{\tau }_{s}\right){\varvec{\sigma}}\):

$$\stackrel{\sim }{{\varvec{\sigma}}}={\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{Newtonian}}+{\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{RBCs}}+{\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{Particles}}$$
(4a)

where the contributions to the overall viscosity are as follows:

\({\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{Newtonian}}=\left[1+f\left(\phi \right)\right]P\left(H,{\lambda }_{\eta }\right)\stackrel{\sim }{\dot{{\varvec{\upgamma}}}}\), (4b)

$${\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{RBCs}}=3\left[1+f\left(\phi \right)\right]{f}_{0}H\left[2\left(1-\xi \right)\left({\rm I}_{1}^{C}\mathbf{C}-\mathbf{C}\bullet \mathbf{C}-\frac{2{\rm I}_{2}^{C}}{3}\mathbf{I}\right)-\frac{\lambda \mathrm{ln}\lambda }{\mathrm{tr}\mathbf{S}}\mathbf{S}\right],$$
(4c)
$${\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{Particles}}=\phi X\beta \theta \left(3\mathbf{a}-\mathbf{I}\right)+\phi \left[\frac{B\left({r}_{e}\right)}{3}\left(\mathbf{a}\bullet \stackrel{\sim }{\dot{{\varvec{\upgamma}}}}+\stackrel{\sim }{\dot{{\varvec{\upgamma}}}}\bullet \mathbf{a}\right)+\frac{\overline{A }\left({r}_{e}\right)}{9}\mathbf{a}\left(\mathbf{a}:\dot{{\varvec{\upgamma}}}\right)\right],$$
(4d)

where \(\beta ={\tau }_{s}/{\tau }_{\mathrm{a}0}\) and \(\stackrel{\sim }{\dot{{\varvec{\upgamma}}}}={\tau }_{s}\dot{{\varvec{\upgamma}}}\) is the dimensionless rate-of-strain tensor. The first contribution, \({\stackrel{\sim }{{\varvec{\sigma}}}}_{\mathrm{Newtonian}}\), refers to the Newtonian contribution to the stress tensor that depends on the concentration of both RBCs and particles, and the RBC viscosity ratio \({\lambda }_{\eta }\). The latter two contributions arise due to RBCs (Eq. (4b)) and particles (Eq. (4c)), respectively (Stephanou et al. 2014; Stephanou 2015, 2018, 2020; Stephanou and Tsimouri 2020).

In total, our model has 11 parameters. Four of these parameters characterize the particles:

  1. 1.

    The aspect ratio between the length and the diameter of the spheroid, \({r}_{e}=l/d\),

  2. 2.

    The particle volume fraction, \(\phi\),

  3. 3.

    The particle rotational time at infinite dilution, \({\tau }_{\mathrm{a}0}\), and

  4. 4.

    \({\Lambda }_{0}^{\mathrm{a}}\) which is a numerical constant.

If all the information concerning the nanoblood sample under study is provided, then the first two parameters should not be considered as adjustable ones. Further note that if the particles are spherical, then only the volume fraction is relevant as the orientation tensor is always equal to its equilibrium value (see “Reduction to special cases” and “Comparison with experimental data”) (Stephanou et al. 2014; Stephanou 2015). Seven additional parameters characterize the RBCs:

  1. 1)

    The slip/non-affine parameter, \(\xi\),

  2. 2)

    The characteristic relaxation time related to the deformability of RBCs, \({\tau }_{S}\),

  3. 3)

    The parameter \(\varepsilon\) that quantifies the relative importance between the regeneration/buildup and the flow-induced breakup of the RBC network,

  4. 4)

    The hematocrit, \(H\),

  5. 5)

    The ratio between the internal (cytoplasm) and external (solvent) fluid viscosities, \({\lambda }_{\eta }\),

  6. 6)

    The external (solvent) viscosity, \({\eta }_{S}\), and

  7. 7)

    \({f}_{0}\), an adjustable parameter.

Again, if all the information concerning the nanoblood sample under study is provided, then three parameters (\(H,{\lambda }_{\eta },{\eta }_{S}\)) are known constants. Note that in case the transient shear viscosity does not present a damping behavior, we may set \(\xi =0\) (Stephanou 2020; Stephanou and Tsimouri 2020); also, if we are interested in nanoblood that involves non-aggregating RBCs (e.g., when the solvent is not plasma), then \(\varepsilon \gg 1\). Overall, there are six adjustable parameters: \({\tau }_{\mathrm{a}0},{\Lambda }_{0}^{\mathrm{a}},\xi ,\varepsilon ,{\tau }_{S}\), and \({f}_{0}\).

Asymptotic behavior of the model in steady-state shear

The asymptotic values of the material functions, in both the zero- and infinite-shear-rate limits, can be obtained analytically. In particular, the zero-shear-rate values of the material functions when \(\varepsilon \to \infty\) are:

$$\frac{{\eta }_{0}}{{\eta }_{s}}=6{f}_{0}H{\left(1-\xi \right)}^{2}+\left[1+f\left(\phi \right)\right]P\left(H,{\lambda }_{\eta }\right)+\frac{2B\left({r}_{e}\right)}{3}\phi +X\phi {\theta }^{2}\frac{\phi }{{\Lambda }_{0}^{\mathrm{a}}},$$
(5a)
$$\frac{{\Psi }_{\mathrm{1,0}}}{{\eta }_{s}{\tau }_{s}}=12{f}_{0}H{\left(1-\xi \right)}^{2}+\frac{2X\phi {\theta }^{2}}{\beta }{\left(\frac{\phi }{{\Lambda }_{0}^{\mathrm{a}}}\right)}^{2},$$
(5b)
$$\frac{{-\Psi }_{\mathrm{2,0}}}{{\eta }_{s}{\tau }_{s}}=6{f}_{0}H{\xi \left(1-\xi \right)}^{2}+\frac{X\phi {\theta }^{2}\left(1-\theta \right)}{\beta }{\left(\frac{\phi }{{\Lambda }_{0}^{\mathrm{a}}}\right)}^{2}.$$
(5c)

On the other hand, at very large shear rates, the shear viscosity, for any value of \(\varepsilon\), approaches the infinite-shear-rate viscosity:

$$\frac{{\eta }_{\infty }}{{\eta }_{s}}=\left[1+f\left(\phi \right)\right]P\left(Ht,{\lambda }_{\eta }\right)+B\left({r}_{e}\right)\phi ,$$
(6)

whereas the two normal stress coefficients vanish.

Reduction to special cases

It is highly important to showcase how the model as derived here, which considers the base fluid to be aggregating blood and with non-spherical (spheroidal) particles, reduces to special cases:

  1. 1.

    Suspension of non-spherical particles in non-aggregating blood

    The case of non-aggregating blood, e.g., when RBCs are immersed in phosphate-buffered saline (PBS), can be obtained when \(\varepsilon \gg 1\) so that \(\lambda =0\) (Stephanou 2020).

  2. 2.

    Suspension of spherical particles in aggregating or non-aggregating blood

    In the case that we limit the analysis only to spherical particles then, \(\mathbf{a}=\left(1/3\right)\mathbf{I}\) and \({r}_{e}=1\), meaning that \(\overline{A }=B=0\) so that the stress tensor simplifies to,

    $${\varvec{\sigma}}=\left[1+f\left(\phi \right)\right]{\varvec{\sigma}}\left(\phi =0\right)\Rightarrow \eta =\left[1+f\left(\phi \right)\right]\eta \left(\phi =0\right),$$
    (7a)

    where \({\varvec{\sigma}}\left(\phi =0\right)\) is the stress expression of neat blood (Stephanou 2020):

    $${\varvec{\sigma}}\left(\phi =0\right)=3{f}_{0}H\left[2\left(1-\xi \right)\left({\rm I}_{1}^{C}\mathbf{C}-\mathbf{C}\bullet \mathbf{C}-\frac{2{\rm I}_{2}^{C}}{3}\mathbf{I}\right)-\frac{\lambda \mathrm{ln}\lambda }{\mathrm{tr}\mathbf{S}}\mathbf{S} \right]+{\eta }_{s}P\left(H,{\lambda }_{\eta }\right)\dot{{\varvec{\upgamma}}}.$$
    (7b)
  3. 3.

    Suspension of spheroidal particles in a Newtonian solvent

When the elasticity of the RBCs is neglected, forcing the conformation tensor \(\mathbf{S}\) to be always equal to its equilibrium value, then blood behaves as a Newtonian fluid; indeed, blood can be considered to be Newtonian in certain conditions, e.g., in large arterial blood vessels. Then, the stress tensor simplifies to

$${\varvec{\sigma}}={\eta }_{s}\left[1+f\left(\phi \right)\right]P\left(H,{\lambda }_{\eta }\right)\dot{{\varvec{\upgamma}}}+{\eta }_{s}\phi \left[\frac{B\left({r}_{e}\right)}{3}\left(\mathbf{a}\bullet \dot{{\varvec{\upgamma}}}+\dot{{\varvec{\upgamma}}}\bullet \mathbf{a}\right)+\frac{\overline{A }\left({r}_{e}\right)}{9}\mathbf{a}\left(\mathbf{a}:\dot{{\varvec{\upgamma}}}\right)\right]+X\left({\eta }_{S}/{\tau }_{\mathrm{a}0}\right)\phi \theta \left(3\mathbf{a}-\mathbf{I}\right),$$
(8)

which is identical to the one presented by Stephanou (2018). In this case, the zero- and infinite-shear-rate shear viscosities, \({\eta }_{0}\) and \({\eta }_{\infty }\), respectively, are given via (Stephanou 2018):

$$\frac{{\eta }_{0}}{{\eta }_{s}}=1+f\left(\phi \right)+\frac{2B\left({r}_{e}\right)}{3}\phi +\frac{X\phi {\theta }^{2}}{{\Lambda }_{0}^{\mathrm{a}}},$$
(9a)
$$\frac{{\eta }_{\infty }}{{\eta }_{s}}=1+f\left(\phi \right)+B\left({r}_{e}\right)\phi .$$
(9b)

Results and discussion

In this section, we present the predictions of the new model in the case of homogenous simple shear, described by the kinematics \(\mathbf{u}=\left(\dot{\gamma }y,\mathrm{0,0}\right)\), and how they compare with available experimental data. The material functions to analyze are as follows: (a) the shear viscosity \(\eta ={\sigma }_{yx}/\dot{\gamma }\) and the two normal stress coefficients \({\Psi }_{1}=\left({\sigma }_{xx}-{\sigma }_{yy}\right)/{\dot{\gamma }}^{2} \, \mathrm{and} \, {\Psi }_{2}=\left({\sigma }_{yy}-{\sigma }_{zz}\right)/{\dot{\gamma }}^{2}\). The results have been obtained by numerically solving the constitutive model (Eq. (1)), using MATLAB. In particular, as we are interested here only in steady-state simple shear, the system of ordinary differential equations boils down to a system of algebraic equations that is solved using the fsolve command of MATLAB. Once the solutions are obtained, then the stress tensor is calculated through Eq. (2a).

Material functions in steady-state simple shear flow

In this section, we will keep the following parameters, associated with RBCs, as constants: \({\lambda }_{\eta }=5\), meaning \(T({\lambda }_{\eta })\approx 2.17\), and \({H}_{m}=72\%\)(Stephanou 2020). We first consider the case of non-aggregating blood (\(\varepsilon \to \infty\), in which case a yield point is not predicted), while we later consider the implications of considering finite values of \(\varepsilon\).

In Fig. 1, we depict the predictions of the dimensionless steady-state material functions as a function of the dimensionless shear rate, \(\mathrm{Ca}={\tau }_{s}\dot{\gamma }\), for various values of the model parameters associated with the particles, while keeping fixed the ones associated with the RBCs: \(H=40\%,\) \(\varepsilon \to \infty , \xi =0.01,\) and \({f}_{0}=1\) (note that here we omit the volume fraction dependency of \({f}_{0}\) for simplicity). We note that the shear viscosity (Fig. 1a) is seen to slightly shear thicken before shear thinning, leading to the exhibition of a maximum, which is seen to intensify as the particle volume fraction is increased. Also, by increasing \(\phi\), both the zero- and infinite-shear-rate viscosities increase following Eqs. (5a) and (6). On the other hand, when the particle aspect ratio is decreased the shear-thickening behavior is seen to vanish without affecting the zero- and infinite-shear-rate viscosities. We should note here that the neat blood predictions do not predict any shear thickening (see Fig. 4 of Stephanou and Tsimouri (2020)) indicating that the shear thickening observed is particle-induced. Overall, these predictions, except for the slight shear-thickening behavior, are in line with available rheological measurements for non-aggregating blood, as will be more carefully documented in the next section. In the case of the first normal stress coefficient (Fig. 1b), we note that the zero-shear-rate value is noted to increase when either \(\beta\) or \({\Lambda }_{0}^{\mathrm{a}}\) decrease, in line with Eq. (5b), whereas the predictions in the range \(1\le \mathrm{Ca}\le 10\) are noted to be insensitive to the parameter values. Also, at large shear rates, a power-law behavior is predicted: when \({r}_{e}=20\) we find that \({\Psi }_{1}\propto {\mathrm{Ca}}^{-1.65},\) whereas when \({r}_{e}=10\) we note that \({\Psi }_{1}\propto {\mathrm{Ca}}^{-2}\). Finally, the second normal stress coefficient (Fig. 1c) is noted to be completely independent of the values of all model parameters varied here, exhibiting a power-law behavior of the form \({-\Psi }_{2}\propto {\mathrm{Ca}}^{-2},\) except for the zero-shear-rate value that is noted to increase when either of \(\beta\) or \({\Lambda }_{0}^{\mathrm{a}}\) decrease, in line with Eq. (5c). A remarkable feature is the appearance of a maximum at about \(\mathrm{Ca}\approx 1\); this is attributed to the RBCs as it can also be noted in the case of neat non-aggregating blood (see Fig. 5b of Stephanou (2020)).

When selecting a finite value for the \(\varepsilon\) parameter, the scalar structural variable, \(\lambda\), decreases as the shear rate increases, as noted in Fig. 2. When the parameter \(\varepsilon\) increases, \(\lambda\) is noted to decrease faster with the shear rate, which is the expected behavior given the definition of \(\varepsilon\). Note that as \(\varepsilon \to 0\) then the structural variable will remain equal to unity irrespective of the applied shear rate, whereas when \(\varepsilon \gg 1\) then the structural vanishes even if a very small shear rate is applied, meaning that a yield-stress point is not predicted. Finally, when increasing the slip parameter to \(\xi =0.01\) we note that the prediction remains invariant at small shear rates and reaches a constant value at large shear rates; this is because after a critical shear rate is reached the RBCs begin to tumble without deforming, which ceases the further destruction of the network.

Fig. 2
figure 2

Representative model predictions for the scalar structural variable, \(\lambda\), as a function of the dimensionless shear rate for various values of \(\xi\) and \(\varepsilon\)

Full size image

When selecting a finite value for the \(\varepsilon\) parameter, a yield point is predicted for all material functions (see panel a in Figs. 3, 4, and 5) irrespective of the value of the remaining parameters. The yield shear stress, \({\upsigma }_{y}\), and yield first, \({N}_{1,y}\), and negative second, \({-N}_{2,y}\), normal stress differences are all seen to be decreasing when \(\xi\) and \(H\) increase and \({f}_{0}\) decreases. The value of the \(\varepsilon\) parameter is, however, completely irrelevant at large shear rates, as expected. In line with Eq. (6), the infinite-shear-rate shear viscosity depends only on the \(H\) when the particle volume fraction and the viscosity ratio \({\lambda }_{\eta }\) are constant. On the other hand, we note that both the first (Fig. 4a) and negative second (Fig. 5a) normal stress differences at large shear rates shift downwards when \(\xi\) increases or \({f}_{0}\) decreases, whereas both curves shift upwards when the \(H\) increases, which is the expected outcome. As mentioned by Stephanou (2020), despite the plethora of shear stress rheological measurements of blood, no data for the normal stresses are available in the literature. The same statement holds also to nanoblood for which only very scarce rheological measurements are available, such as the ones reported by Antonova et al. (2014) (see Figs. 6, 7, 8, and 9).

Fig. 3
figure 3

Representative model predictions for the steady-state dimensionless a shear stress and b shear viscosity as a function of the dimensionless shear rate. The results are shown for various values of the model parameters associated with the RBCs, while keeping \(\phi =1\%,{r}_{e}=10, \beta =1,\) and \({\Lambda }_{0}^{\mathrm{a}}={10}^{-2}\)

Full size image
Fig. 4
figure 4

Representative model predictions for the steady-state dimensionless a first normal stress difference and b first normal stress coefficient as a function of the dimensionless shear rate. Same parameter values as in Fig. 3

Full size image
Fig. 5
figure 5

Representative model predictions for the steady-state dimensionless negative a second normal stress difference and b second normal stress coefficient as a function of the dimensionless shear rate. Same parameter values as in Fig. 3

Full size image
Fig. 6
figure 6

Comparison of the model predictions for the shear viscosity with the experimental rheological data of neat non-aggregating blood with H = 40% (Antonova et al. 2014). The dotted line depicts the Newtonian contribution \(P\left(H,{\lambda }_{\eta }\right)\dot{\gamma }\), where \(P\left(H,{\lambda }_{\eta }\right)\) is given from Eq. (2b), at large shear-rates

Full size image
Fig. 7
figure 7

Comparison of the model predictions for a the shear viscosity of a suspension of spherical nanoparticles in saline (a Newtonian fluid) and b the shear viscosity of suspensions of spherical nanoparticles in non-aggregating blood at various concentrations with the experimental rheological data of Antonova et al. (2014)

Full size image
Fig. 8
figure 8

Comparison of the model predictions for the shear viscosity of a suspension of rod-like nanoparticles in saline (a Newtonian fluid) with the experimental rheological data of Antonova et al. (2014). The dotted lines depict the saline viscosity value \({\eta }_{s}=0.739 \mathrm{mPa s}\)

Full size image
Fig. 9
figure 9

Comparison of the model predictions for the shear viscosity of a suspension of rod-like nanoparticles in non-aggregating blood (blue line) with the experimental rheological data of Antonova et al. (2014). Note that the neat blood prediction, as presented in Fig. 6, is also presented (red line)

Full size image

Comparison with experimental data

Rheological data for neat non-aggregating blood

Antonova et al. (2014) measured the viscosity of a RBC-PBS solution, with a hematocrit value adjusted to 40%, that they use as their base fluid for preparing their nanoblood suspensions (see Figs. 7b and 9). Given that the viscosity of PBS (external fluid) at 37 °C (where all rheological measurements were made) is \({\eta }_{s}=0.73 \, \mathrm{mPas}\) (Chien et al. 1966) and that the cytoplasmic fluid that resides inside the RBC has a viscosity in the range of 3–10 mPa s (Cokelet and Meiselman 1968), then the ratio between the internal (cytoplasm) and external (plasma) viscosities is between 4.1 and 13.7; we will here consider the average value \({\lambda }_{\eta }=8.9\), meaning that \(T({\lambda }_{\eta })\approx 2.3\). Furthermore, since the measurements were made on a RBC-PBS solution, i.e., using as a base fluid non-aggregating blood, we consider \(\varepsilon \to \infty\). Thus, only three more parameters remain to be selected:\({f}_{0}, \xi\), and \({\tau }_{s}\).We compare the model prediction with these experimental data in Fig. 6, considering \({f}_{0}=3.8, \xi =0, \, \mathrm{and} \, {\tau }_{s}=0.5 \mathrm{s}\), and note that the model compares adequately well with the experimental data.

Rheological data for a suspension of spherical nanoparticles in non-aggregating blood

We first consider the viscosity of a suspension of core–shell-type star polymers in a saline solution whose interior is composed of hyperbranched-polystyrene-bearing arms of poly(acrylic acid); these star polymers resemble spherical nanoparticles with a (hydrodynamic) radius equal to 14 nm (Antonova et al. 2014). As seen in Fig. 7a, these solutions behave as Newtonian fluids at various concentrations. The average viscosity of the saline solution is equal to \({\eta }_{s}=0.739\,\mathrm{mPas}\). Since the volume fraction for each Newtonian suspension is not provided by Antonova et al., (2014), we use the shear viscosity data and the Batchelor and Green (1972a, b) expression, \({\varvec{\sigma}}={\eta }_{s}\left(1+\frac{5}{2} \phi +\overline{c}{\phi }^{2}\right)\dot{{\varvec{\upgamma}}}\) with \(\overline{c }\)=6.2, in order to obtain an estimate of the volume fraction for each solution; we find \(\phi \left(0.02\frac{\mathrm{mg}}{\mathrm{ml}}\right)=0.0565, \phi \left(0.1\frac{\mathrm{mg}}{\mathrm{ml}}\right)=0.0594\), and \(\phi \left(0.2\frac{\mathrm{mg}}{\mathrm{ml}}\right)=0.082\). Then, and given that the volume fraction for each nanoblood suspension is small (same volume fractions are assumed as those obtained in Fig. 7a), we consider that the parameters \(\xi ,\varepsilon ,{\tau }_{S}\), and \({f}_{0}\) remain equal to those selected in the case of neat blood (Fig. 6). Thus, the model predictions are obtained by simply rescaling the pure blood predictions, as dictated by Eq. (7a), by a factor \(1+f\left(\phi \right)\). The comparison against the shear viscosity data of the suspension of spherical nanoparticles in a Newtonian fluid is provided in Fig. 7a which showcases that the selected values for the volume fractions are adequately accurate. Then, we proceed to predict, without additional adjustable parameters, the rheological predictions in the case of suspensions of spherical nanoparticles in non-aggregating blood at various concentrations. This comparison, shown in Fig. 7b, is again noted to be adequately good.

Rheological data for a suspension of rod-like nanoparticles in non-aggregating blood

Antonova et al. (2014) measured the viscosity of linear  poly(acrylic acid) (PAA) chains, that exhibit a rod-like conformation, immersed in a saline solution for various molecular weights and two concentrations (0.02 and 0.1 mg/mL). Unfortunately, as no information is given concerning the aspect ratio or the volume fraction for either sample, we will obtain these values by comparison against the rheological data. Antonova et al. (2014) noted that in their experiments, carried out with rather dilute PAA solutions, no measurable effect of the polymer molecular weight on the viscosity was registered. By simply fitting their experimental data against our model, we obtain the following values for the volume fractions \(\phi \left(0.02\frac{\mathrm{mg}}{\mathrm{ml}}\right)=0.005\) and \(\phi \left(0.1\frac{\mathrm{mg}}{\mathrm{ml}}\right)=0.025\), and the following values of all remaining parameters: aspect ratio \({r}_{e}=120, {\Lambda }_{0}^{\mathrm{a}}={10}^{-5}\), and \({\tau }_{\mathrm{a}0}=0.5\, \mathrm{s}\). The comparison is presented in Fig. 8 where we note a satisfactory comparison against the experimental data. Note that the non-spherical shape of these particles induces a shear-thinning behavior (cf. Figure 7a where spherical particles are used).

Finally, in Fig. 9, we compare against the experimental rheological data of Antonova et al. (2014) with \(H=0.4\) for a suspension of rod-like (PAA) nanoparticles in non-aggregating blood, with the solvent being PBS (\(\eta_s=0.73\,\mathrm{mPa}s\)), when the particle concentration is 0.1 mg/mL (\(\phi =0.025\)). As mentioned by Antonova et al. (2014), the rheological measurements are the same irrespective of the molecular weight of PAA, as also noted in Fig. 8; however, the nanoblood rheological measurements are lower than the neat blood ones. Since PAA particles are immersed in PBS (i.e., the base fluid is non-aggregating blood), we consider \(\varepsilon \to \infty\). Again, as the volume fraction is small, we employ the same values for the parameters \(\xi\) and \({\tau }_{S}\) as in the case of neat blood (Fig. 6), i.e., \(\xi =0 \,\mathrm{and}\, {\tau }_{s}=0.5\, \mathrm{s}\), but consider \({f}_{0}\left(\phi =0.025\right)=2\). The fact that this parameter has to be reduced is a direct consequence of the viscosity reduction of the nanoblood sample relative to the neat blood one. Finally, we consider the same values for the parameters \({r}_{e},{\Lambda }_{0}^{\mathrm{a}}\), and \({\tau }_{\mathrm{a}0}\) as in the case of the rod-like suspensions compared in Fig. 8, i.e., \({r}_{e}=120, {\Lambda }_{0}^{\mathrm{a}}={10}^{-5}\), and \({\tau }_{\mathrm{a}0}=0.5\, \mathrm{s}\). As noted in Fig. 9, the new model is seen to compare favorably with the experimental data of Antonova et al. (2014).

Conclusions

Given the growing interest in the use of drug-carrying particles to fight numerous diseases, it becomes absolutely necessary to have in our arsenal a rheological model that adequately addresses the rheology of nanofluids comprising of particles immersed in blood, what we refer to in this article as nanoblood. The availability of such a model would, undoubtedly, offer the opportunity to execute in silico human trials to bolster the reduction of both the time needed to perform and the cost of design experiments, and further the tailor-design of each drug for optimum delivery to the infected areas. Such a model should be properly derived to accommodate the coupling between particles and blood, not to violate physical laws, and be properly parameterized against available rheological data.

We have herein introduced such a constitutive model derived via the use of non-equilibrium thermodynamics. Thus, the state variables have been properly coupled with the flow field, while concurrently guaranteeing the model's thermodynamic admissibility. To the best of the author’s knowledge, this is the very first rheological model derived for nanoblood using such an approach. In deriving the model, we have considered the actual morphologies of both the RBCs, modeled as deformed droplets with a constant volume that are able to aggregate, and particles, modeled as rigid spheroids. We have thoroughly established its capacity to favorably compare against available suspension rheological data (Antonova et al. 2014) with both spherical (Fig. 7b) and rod-like nanoparticles (Fig. 9) immersed in non-aggregating blood (RBC suspension in PBS), given especially the fact that the parameters (except for \({f}_{0}\) in Fig. 9) are carried over from previous comparisons: the comparisons in Fig. 7b and Fig. 9 are purely predictive ones as the parameters are kept the same to the ones used in Figs. 6 and 7a, for Fig. 7b, and Fig. 8, for Fig. 9, respectively. In addition, we have showcased that the model is noted to predict a positive first normal stress difference (Figs. 1 and 4), and a negative second normal stress difference (Figs. 1 and 5). Recent experimental evidence (Kim et al. 2019) suggests that the migration of particles is a direct result of the appearance of normal stress differences. Alas, if optimization of the tailor-design of DCPs is to be sought, it is extremely paramount to properly predict normal stresses in vascular environments; however, it has been ignored by the experimental rheological community altogether, as no such data can be found in the literature. This surely leaves a gap in our effort to develop a constitutive model able to accurately address the rheological behavior of both neat blood and nanoblood. We, therefore, urge the experimental rheological community to be more actively interested in studying the rheological behavior of nanoblood comprising particles immersed in normal (aggregating) blood.

We do admit that the current version of our model bears certain limitations. At first, in our previous works (Stephanou 2020; Stephanou and Tsimouri 2020) the model parameters were considered constants. We could use carefully executed experiments (Fischer 2007) wherein single RBCs are suspended in fluids with viscosities ranging from 12.9 to 109 mPa s. This would allow the specification of the \(\xi\) parameter as a function of the viscosity ratio \({\lambda }_{\eta }\). Also, since the yield (shear) stress can be expressed (Apostolidis and Beris 2014) as a function of the \(H\) and the concentration of fibrinogen, \({c}_{f}\), a key blood protein responsible for the aggregation of RBCs, we could specify the \(\varepsilon\) parameter as a function of \(H\) and \({c}_{f}\). These generalizations would reduce the number of adjustable parameters. Secondly, its validation has been made using rheological data (Antonova et al. 2014) that refer to non-aggregating blood. This stands as a huge void in our efforts to showcase the model capacity to address the rheological behavior of nanoblood in actual vascular environments. However, as mentioned above, this void stems from the unavailability of such data in the literature. Thirdly, in our present derivation we have neglected coupling terms between the structural variables for mere simplicity. Given the fact that nanoparticles are noted to interact with blood cells (Vu et al. 2020), and that blood proteins interact with nanoparticles by forming an absorbed layer on the nanoparticles that leads to nanoparticle aggregation (De Paoli Lacerda et al. 2010), such an amendment is particularly important. The use of NET will allow us to properly couple the structural variables while, at the same time, guaranteeing the thermodynamic admissibility of the revised version, as we have done previously for polymer nanocomposites with remarkable success (Stephanou et al. 2014; Stephanou 2015).

In the future, we aim to employ our model in CFD simulations in vascular environments to provide answers as to how the tailor-design of DCPs should take place. This could include the use of numerical techniques, such as the FEM or finite-volume method (FVM) solvers. Commercial software such as Fluent (Ansys Inc.) or FreeFem +  + , an open-source FEM solver used by Dubey et al. (2020), could also be employed in performing the aforementioned CFD simulations. This will enable us to understand what specific factors result in the optimal margination of DCPs to the vessel wall. As we have exemplified prior, when one targets atherosclerotic plaques, the highest accumulation of DCPs is noted at the shoulder regions of the plaque (Peters et al. 2009) where the stresses are minimum. This result should be verified via CFD simulations. Overall, the use of our model in CFD simulations would allow for the execution of in silico human trials of several treatments to facilitate a faster and more economic design of suitable DCPs.