- Email: [email protected]
A new generalized Oldroyd-B model for blood flow in complex geometries
International Journal of Engineering Science 72 (2013) 78–88
Contents lists available at SciVerse ScienceDirect
International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci
A new generalized Oldroyd-B model for blood flow in complex geometries M. Anand a,⇑, J. Kwack b, A. Masud b a b
Department of Chemical Engineering, Indian Institute of Technology Hyderabad, Yeddumailaram, Andhra Pradesh 502205, India Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
a r t i c l e
i n f o
Article history: Received 21 May 2013 Received in revised form 18 June 2013 Accepted 21 June 2013 Available online 3 August 2013 Keywords: Blood Oldroyd-B model Thermodynamic framework Variational multiscale method Convergence
a b s t r a c t A new generalization of the Oldroyd-B model is developed to describe the flow of blood. The model is developed within a thermodynamic framework which recognizes that a viscoelastic fluid can remain stress free in multiple configurations. The new model is an improvement over an earlier model developed within the same framework to describe the response characteristics of blood. It captures the shear-thinning and deformationdependent viscoelastic behavior of blood just like the previous model. More importantly, unlike the previous model, it does not have the shortcoming of an abrupt transition of the material properties at low shear rates: instead, it allows for a smooth variation of the rate of dissipation, and therefore viscosity, over the entire range of physically feasible shear-rates. This feature is very attractive for developing high-fidelity numerical methods for application to the complex geometries that are typically encountered in the human vasculature. Convergence of the numerical method in simple geometries shows its superior properties as compared to the earlier model: this demonstration of model performance is a precursor to its use in 3D geometries. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Human blood is a suspension of formed elements (erythrocytes, leukocytes, and platelets, that together constitute 45% by volume) in an aqueous solution (plasma). It exhibits complex rheological behavior namely shear-thinning by orders of magnitude (Chien, Usami, Taylor, Lundberg, & Gregersen, 1966), and deformation-dependent viscoelastic behaviour (recorded as shear-rate dependent relaxation times (McMillan, Utterback, Nasrinasrabadi, & Lee, 1986) or normal stress differences (Copley & King, 1976)). An accurate model that captures this range of behavior in a single framework can be invaluable in generating precise knowledge of flow patterns and stress distributions in the heart and the vasculature. Such a model can be a tool for patient-specific simulations wherein the geometry information is gleaned from either angiograms or intra-vascular-ultrasound images. Patient-specific simulations are an important tool in tracking the progress of cardiovascular pathologies and identifying candidates for further medical intervention. Examples of such cardiovascular pathologies include aneurysms and atherosclerotic plaques: these pathologies are influenced in their development by the flow patterns of blood (e.g. aneurysms (Inzoli, Boschetti, Zappa, Longo, & Fumero, 1993; Raghavan, Vorp, Federle, Makaroun, & Webster, 2000), and atherosclerosis (Zarins et al., 1983)). Let us consider the case of an (intra-cranial saccular) aneurysm: the dilation of the artery wall due to mechanical weakness leads to the formation of a balloon-shaped aneurysm with a local circulation (slow flow with low shear rates) within it. Calculating the stresses imposed by blood flow on the walls of the aneurysm is a key measure in determining how much further it will stretch, and a constitutive model for blood used in this simulation ⇑ Corresponding author. Tel.: +91 40 23016090; fax: +91 40 23016032. E-mail addresses: [email protected] (M. Anand), [email protected] (J. Kwack), [email protected] (A. Masud). 0020-7225/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2013.06.009
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
79
must be accurate at the low shear rates that are encountered in this geometry. Further, the dilatation occurs adjacent to the mainstream flow of the cerebral artery where shear rates are high. The constitutive model has to capture stress distributions in the high shear rate regimes of the mainstream cerebral flow as well. This example illustrates the need for a constitutive model for blood that makes accurate predictions of flow over a wide range of shear rates: from extremely low values (at near stagnant conditions) to high values. A systematic thermodynamic framework for the development of rate-type fluid models was presented by Rajagopal and Srinivasa (2000). This thermodynamic framework results from the melding together of two central ideas. The first idea is that a body’s ‘‘natural configuration’’ (the configuration in the absence of external stimuli) evolves as the body undergoes a process, and that it accesses different sets of these configurations depending on the class of processes it is subject to. The stress field of the body in the current configuration is determined with respect to its deformation from the natural configurations for the particular process. The second idea is that the form of the stress tensor can be determined uniquely by specifying the forms- in terms of the kinematic variables and temperature - for the stored energy function and the rate of dissipation, and then maximising the rate of dissipation subject to relevant constraints. The constraints being that (i) the rate of dissipation is equal to the difference between the stress power and the rate of change of the stored energy, and (ii) the material is incompressible. Using such a framework, Anand and Rajagopal (2004) developed a model (which we refer to from now on as the Anand–Rajagopal model) that was intended to describe the flow characteristics of blood. The Anand–Rajagopal model was one among a few other shear-thinning viscoelastic models, for describing blood flow that were developed within a few years of each other. These other models include those by Sun and DeKee (2001), and by Owens (2006). The model in Sun and DeKee (2001) uses a shear rate-dependent structural parameter combined with a standard frame-invariant convected derivative to incorporate the shear-rate dependent viscosity and relaxation times. The model in Owens (2006) also uses a shear rate-dependent structural parameter to model the viscosity and relaxation time of blood, but the parameter is derived using theories in polymer network formation. Owens considered the erythrocyte as analagous to a monomer, and the erythrocyte aggregates (rouleaux) in blood as analogous to a ‘k’-mer (a network of ‘k’ cross-linked monomer units). The difference between these other models and the Anand–Rajagopal model lies in the framework adopted by the latter: they used the microstructural evolution of erythrocyte aggregates (rouleaux) and the elasticity of the erythrocyte to motivate forms for the stored energy function and rate of dissipation characterizing blood. The Anand–Rajagopal model was a generalization of the classical Oldroyd-B model that captured the rheological behavior of blood over a wide range of shear rates (0.06 s1 to 650 s1 ). However, the model equations, by themselves, presented the problem of infinite viscosity being recorded at zero-shear rate. In order to correct this anomaly, a low-shear rate criterion that transitioned between two forms of dissipation (and apparent viscosity forms), was included to ensure a finite zero-shear viscosity and to capture viscosity variation in the low shear-rate range (0 s1 to 0.06 s1 ). The Anand–Rajagopal model has recently been used to simulate blood flow in a straight tube with a stenosis, as well as in a 90° curved tube. Bodnar, Rajagopal, and Sequeira (2011) used a finite-volume discretization along with an explicit Runge– Kutta time marching scheme: they thereby obtained numerical results at three flow rates (Reynolds number less than 100) that illustrated the dynamics and flow patterns generated. They found that in a straight tube with stenosis the relative difference of the axial velocities, between that predicted by the Anand–Rajagopal model and that by the Newtonian model, increased as the flow rate decreased. They also found that in a 90° curved tube the axial velocity contour is pushed away from the central axis, and there is a slowing down of velocity towards the outer radius along with a corresponding speed-up towards the inner radius. The secondary vortices that are set up are slower than the corresponding ones in a Newtonian fluid due to the higher viscosity of the model. Further, they found that there is a flattening of the velocity at the central core, both in the stenosed pipe and the curved pipe, due to the shear-thinning nature of the model. This was a reiteration of the results that their group had obtained earlier in Bodnar, Sequeira, and Prosi (2011) and Bodnar, Sequeira, and Pirkl (2009) using an Oldroyd-B model with a Cross model-type shear-thinning viscosity. However, they were not concerned with the convergence rate of the simulations using the model. We performed numerical tests and found difficulties in the numerical implementation of the Anand–Rajagopal model (along with the low shear rate criterion) – specifically non-convergence of the residual of the mass conservation, momentum balance, and constitutive equations when used in a bent tube geometry. Our investigations have shown that it is non-trivial to obtain an appropriate low-shear criterion1 that will be acceptable from a modeling standpoint, and will also satisfy the requirements of the particular numerical scheme. In this paper, we therefore introduce a new model which has a finite zeroshear viscosity and allows for a smooth variation of viscosity at all possible shear-rates from zero-upwards: we corroborate the model with existing data, and generate parameters to fit this data. We then show numerical convergence plots for a simplified version of the model in two sample geometres (straight cylindrical tube, and bent tube) to demonstrate the potential of the model in the numerical scheme as well as to set the stage for simulations in more complex geometries. The rest of the paper is organized as follows: The notation used to describe the kinematics of a viscoelastic fluid with multiple stress-free states is given in Section 2. We motivate the need for, and propose, the new model in Section 3. The new model is subject to a corroboration procedure in Section 4 to get a set of constants that will best fit experimental data. Numerical results for sample geometries are documented in Section 5 to demonstrate the promise of the new model for 3D simulations. The results are discussed in Section 6.
1
Note that we do not say that such a criterion does not exist.
80
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
2. Preliminaries Consider the following schematic (Fig. 1) of the deformation of a viscoelastic fluid body with a single relaxation mechanism, and capable of instantaneous elastic response. In Fig. 1, jR ðBÞ and jt ðBÞ denote the reference and the current configuration (at time t) of the body B, respectively. jpðtÞ ðBÞ denotes the stress-free configuration that is reached by instantaneously unloading the body from its current configuration. The motion of the body is defined through a one-to-one mapping between each point in the reference configuration (X 2 jR ðBÞ) and a corresponding point in the current configuration at time instant t (x 2 jt ðBÞ):
x ¼ vjR ðX; tÞ:
ð1Þ
We assume the motion to be sufficiently smooth and invertible. We will also, for convenience, suppress B in the notation
jt ðBÞ etc. The deformation gradient FjR , and the left Cauchy-Green stretch tensor BjR are defined as
FjR ¼
@ vjR ; @XjR
and BjR ¼ FjR FTjR :
ð2Þ
The tensor FjpðtÞ denotes the deformation gradient between the instantaneous stress-free (natural) configuration and the current configuration, provided the deformation is homogeneous. Further, the left Cauchy-Green stretch tensor associated with the instantaneous elastic response from the stress-free configuration jpðtÞ is then given by:
BjpðtÞ ¼ FjpðtÞ FTjpðtÞ :
ð3Þ
The mapping G is defined by:
G ¼ FjR !jpðtÞ ¼ F1 jpðtÞ FjR :
ð4Þ
The velocity gradients L and LjpðtÞ are defined by:
_ 1 L ¼ F_ jR F1 jR ; andLjpðtÞ ¼ GG :
ð5Þ
where the dot represents the material time derivative. The symmetric parts of L and LjpðtÞ are obtained in the usual manner, and are represented by the symbols D and DjpðtÞ , respectively. r The upper convected Oldroyd derivative of BjpðtÞ ; BjpðtÞ , is defined as: r
BjpðtÞ ¼ B_ jpðtÞ LBjpðtÞ BjpðtÞ LT :
ð6Þ
Additional constraints, that are enforced during the afore-mentioned maximization procedure for the rate of dissipation, are those due to incompressibility (of blood); these imply:
trðDÞ ¼ 0
ð7Þ
trðDjpðtÞ Þ ¼ 0:
ð8Þ
and
Fig. 1. Configurations accessed by a viscoelastic fluid body with instantaneous elasticity.
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
81
3. Constitutive model for blood Recent numerical tests have shown that the low shear-rate criterion given in Anand and Rajagopal (2004) is not suitable for numerical implementation in complex geometries. This is because it does not allow for convergence of the residuals of the numerical equations. The abrupt transition (using a Heaviside function) in the model terms for relaxation time and shear-thinning index (the parameters K and n, respectively, in the model equations: see Anand & Rajagopal, 2004) has been pin-pointed as the reason for the non-convergence. This abrupt transition results in numerical oscillations in the internal fields described by these terms. Because of the fully coupled nature of the system of equations, these internal ripples affect the computed fields and lead to the divergence of the solution from the actual solution. For these reasons, we tried to develop an alternate low-shear rate criterion which would both ensure a smooth transition in relaxation time and shear-thinning index, through the cut-off shear-rate of 0.06 s1 , as well as ensure a finite apparent viscosity at zero shear rate. To get a smooth transition in the model equations of Anand and Rajagopal (2004), K and n must be functions of Bjp ðtÞ . Imposing constraints from both a modeling standpoint and a numerical standpoint resulted in a procedure to specify K(Bjp ðtÞ ) and n(Bjp ðtÞ ). However, we found that it is not trivial to identify appropriate forms for K and n: we tested power law forms among others (which are natural choices for shear-thinning behavior) but were unable to satisfy the constraints in the procedure. We therefore decided not to seek further for acceptable functional forms for the relaxation time and shearthinning index, and instead proceeded to construct a new model which does not pose such numerical difficulties requiring additional workarounds. The forms for the stored energy, and rate of dissipation, proposed for the new model are:
W¼
l
ðIB 3Þ; 2 n ¼ g Djp ðtÞ Bjp ðtÞ Djp ðtÞ þ g1 D D:
ð9Þ ð10Þ
We choose:
g ¼ g^ ðBjp ðtÞ Þ ¼ a trBjp ðtÞ
m
;
m > 0:
ð11Þ
We point out that the rate of dissipation contains a power-law term, and a viscous term (akin to a Newtonian fluid), whereas the stored energy form is that of the neo-Hookean type. Srinivasa (2000) used a pure power-law term in the rate of dissipation expression along with a neo-Hookean type stored energy function, and developed a Maxwell-type model from these forms. The model proposed here uses a different method compared to Srinivasa (2000) to incorporate the shear-thinning viscosity: namely, a power-law dependence in the rate of ^ ðBjp ðtÞ Þ, instead of a power-law dependence on the term Djp ðtÞ BjpðtÞ Djp ðtÞ . dissipation coefficient g ¼ g The model equations derived using the constrained maximisation procedure in Rajagopal and Srinivasa (2000) are:
T ¼ p1 þ S;
ð12Þ
S ¼ lBjp ðtÞ þ g1 D; r 2l Bjp ðtÞ ¼ Bjp ðtÞ k1 ;
ð13Þ
g
k¼
3 trB1 jp ðtÞ
:
ð14Þ ð15Þ
4. Model corroboration In this section we corroborate the model with experimental data for both apparent viscosity, as well as a shear flow loop hysteresis, so as to come up with a set of model parameters that matches both sets of experimental data. Ensuring that the model fits both sets of data (apparent viscosity, and shear flow loop hysteresis) allows us to identify a unique set of parameters for the model. 4.1. Apparent viscosity We consider the flow between two concentric cylinders separated by a small gap, with the inner cylinder being held at rest. The flow between the cylinders is assumed to be of the form (we use some of the notation that is used in Anand & Rajagopal (2004) to enable easy comparison):
v ¼ uðrÞ^eh ¼ rxðrÞ^eh ;
ð16Þ
which automatically satisfies the balance of mass. Substituting the constitutive Eqs. (12)–(15) in the equations for balance of linear momentum, and assuming that the stress field depends only on ‘r’, we get:
82
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
gk þ g1
T rh ¼
2
du u : dr r
ð17Þ
In a concentric-cylinder rheometer, it is possible to approximate the shear-rate term in the following manner:
du u dx ¼ c_ ¼ R0 : dr r dr
ð18Þ
The expression for apparent viscosity (gapp ) for this model is thus:
gapp ¼
gk þ g1 2
ð19Þ
;
where g is a function to be determined using Eq. (11). k is determined using the incompressibility condition (detBjp ðtÞ = 1) as:
1 k¼ 1=3 : g2 c_ 2 1 þ 4l2
ð20Þ
Noting that:
2c_ 2 g2 1 trBjp ðtÞ ¼ 3 þ 1=3 ; 4l2 g2 c_ 2 1 þ 4l2
ð21Þ
We obtain the following non-linear equation for g to be solved at each value of c_ :
2 3m 2g2 c_ 2 1 6 7 g ¼ a4 3 þ 1=3 5 : 4l2 g2 c_ 2 1 þ 4l2
ð22Þ
We fix the four unknown parameters (a; m; g1 ; l) by minimising the error between predicted viscosity (Eq. 19), and measured viscosity at a set of shear-rates (data reported in Yeleswarapu (1996)). The minimization is done using the MATLAB function ‘fminsearch’, and a residual (sum of least square of error for both apparent viscosity data and shear loop hysteresis data) of 0.0052 is reported. Since gk ¼ 0 as c_ ! 1, we set g1 ¼ 2g1 ¼ 0:01 Pa s. There are multiple sets of parameters that fit the data equally well, and one set that we find is a = 953.655 Pa s, m = 7.7779, l = 0.0283 Pa, g1 = 0.01 Pa s; this set predicts a zero-shear viscosity of g0 = 0.0978 Pa s. The data fit for these parameters is given in Fig. 2: we see that the proposed model fits the data as well as the Anand–Rajagopal model (which uses the parameters: K = 1.2056 s1 , n = 0.7525, l = 0.0227 Pa, g1 = 0.01 Pa s). We note that the relaxation time associated with the evolution of Bjp ðtÞ is given by:
trel ¼
g aðtrBjp ðtÞ Þm ¼ : 2l 2l
ð23Þ
A zero-shear relaxation time of 3.278 s is obtained for the proposed model, and this compares favorably with the value of 3.022 s obtained for the Anand–Rajagopal model. −1
10
Apparent Viscosity (Pa⋅ s)
Experimental data Anand−Rajagopal model Proposed model
−2
10
−3
10
−2
10
−1
10
0
10
1
10
2
10
3
10
Shear rate (sec−1)
Fig. 2. Apparent viscosity vs c_ for proposed model: comparison with experimental data in Yeleswarapu (1996).
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
83
4.2. Shear loop hysteresis Experimental data is available for the shear stress at various instants during, and residual shear stress (hysteresis) at the end of, the application of a triangular loop of shear rate in Bureau, Healy, Bourgoin, and Joly (1980) (see also Bureau, Healy, Bourgoin, & Joly, 1978). In order to corroborate the proposed model with this data, we first derive the equations governing it’s time variation of shear stress in a time-dependent Couette flow. We then compare the model predictions with the experimental data for a triangular loop of shear rate. The applied shear rate varies in the following (triangular) manner:
c_ ðtÞ ¼ c€c t; t 6 t0 ; c_ ðtÞ ¼ 2c€c t0 c€c t; t0 < t 6 2t0 :
ð24Þ ð25Þ
Here, for the data in Bureau et al. (1980), we have c€c ¼ 0:0185s2 , and t0 = 6.5 s. We seek a time-dependent Couette flow solution for the velocity field:
v ¼ c_ ðtÞy^ex :
ð26Þ
Substituting the above velocity field in the constitutive model, and assuming that the stress is a function of only ‘y’, we obtain:
dBxy 2l ¼ B þ c_ Byy ; dt g xy dBxx 2l ¼ ðB kÞ þ 2c_ Bxy ; dt g xx dByy 2l ¼ ðB kÞ; dt g yy
ð27Þ ð28Þ ð29Þ
where:
Sxy ¼ lBxy þ
g1 c_
ð30Þ
; 2 g ¼ g^ ðBjp ðtÞ Þ ¼ aðBxx þ 2Byy Þm ; k¼
ð31Þ
3Byy ðBxx Byy B2xy Þ B2yy þ 2Bxx Byy B2xy
ð32Þ
and the other stress components are given by Bzz ¼ Byy ; Bxz ¼ Byz ¼ 0. The stress components at initial time (t = 0) are calculated by setting the initial elastic stretch to be unity (Bjp ðtÞ ¼ 1), and the initial shear rate to be zero (c_ ¼ 0s1 ). Hence, Bxy ðt ¼ 0Þ ¼ 0; Bxx ðt ¼ 0Þ ¼ 1, and Byy ðt ¼ 0Þ ¼ 1. We use a simple ExplicitEuler scheme to solve this ODE-IVP, and obtain the stress components at each instant, t, in 0 6 t 6 2t 0 . Although Runge–Kutta methods of higher accuracy are available to solve the ODE-IVP, we have settled for the simplest method which yields a solution. We only need to ensure that the results are grid independent to have confidence in the results from the scheme. We determined Dt = 0.01 s to be the smallest time step for which the results are reproducible, and this time step is sufficient for our purpose of comparison with the data. −3
8
x 10
7
Shear Stress (Pa)
6 5 4 3 2 Experimental data Anand−Rajagopal model Proposed model
1 0
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
−1
Shear rate (sec )
Fig. 3. Shear stress variation with time during, and hysteresis at the end of, a triangular loop of shear rate. Data from Bureau et al. (1980).
84
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
The shear stress, Sxy , predicted by the model is compared with the experimental data, and the parameters that give the best match with the relevant data in Bureau et al. (1980) (as well as for the apparent viscosity data) are reported. The parameters are thus the same as those given earlier for the apparent viscosity data, namely a = 953.655 Pa s, m = 7.7779, l = 0.0283 Pa, g1 = 0.01 Pa s. The predictions from both models are compared with the data in Fig. 3: we see that the match between the proposed model and data is excellent, while that between the Anand–Rajagopal model and data is not as good. We hasten to point out that the less than exemplary data match with this one set of data does not in any way detract from the capability of the Anand–Rajagopal model in its ability to model blood. The reason for the poor fit is that the parameters used for the Anand–Rajagopal model were obtained by minimizing the error between model predictions with both apparent viscosity data (used in subSection 4.1), and the data for oscillatory flow in a pipe given in Thurston (1975) (which is different from the data used in this subsection). Therefore, to get a better fit using the Anand–Rajagopal model, we have to rerun the optimizations to get best fit with the two sets of data given in this paper. However, we will not undertake this exercise, as the object of this paper is to show the improved efficacy of the proposed model in high-fidelity simulations in complex geometries. Obtaining new parameters for the Anand–Rajagopal model will still not improve it’s performance in such simulations, as we shall show in the next section, because of the abrupt transition of properties at low shear rates which cannot be avoided. We thus use the same parameters for each model used in Figs. 2 and 3 to demonstrate the numerical results as well. 5. Stabilized three-field formulation and numerical results A three-field multiscale finite element method for incompressible viscoelastic fluids was developed by Kwack and Masud (2010). This method is comprised of velocity, pressure and conformation-tensor as the unknown fields. The method finds roots in the Variational Multiscale framework proposed by Hughes (1995) and Hughes, Feijoo, Mazzei, and Quincy (1998), Masud and Khurram (2006), Masud and Franca (2008) and Masud and Calderer (2009), along with SUPG (Streamlined-Upwind Petrov Galerkin) type ideas applied to the constitutive model Brooks and Hughes (1982). High Weissenberg number flows were investigated in Kwack and Masud (2010) and the computed results were compared with four-field formulations wherein rate-of-deformation tensor was the fourth field. These competing methods include the EVSS (Elastic-Viscous Split Stress) (Rajagopalan, Armstrong, & Brown, 1990), DEVSS methods (Discrete Elastic-Viscous Split Stress) (Guenette & Fortin, 1995), and the four-field GLS (Galerkinn Least-squares) (Coronado, Arora, Behr, & Pasuali, 2006) methods. It is important to note that the three-field method (Kwack & Masud, 2010) leads to ten degrees of freedom per node in 3D calculations, while the four-field formulations result in sixteen degrees of freedom per node. Consequently, not only the bandwidth for the three-field method is smaller as compared to the four field methods, it also leads to a substantially small assembled system for the same spatial discretization. The enhanced stability properties of the three-field formulation (Kwack & Masud, 2010) help resolve the flow with the same accuracy as the competing methods, however at substantially reduced computational costs. Both qualitatively and quantitatively, the three-field method shows excellent agreement with the published results on different benchmark problems. In several cases considered, higher Weissenberg numbers were reached than were reported previously in the literature. We have now extended the ideas presented in Kwack and Masud (2010) to full VMS (Variational Multiscale) based three-field formulation that results in enhanced numerical stability in 3D simulations: Details of this formulation will be presented in a forthcoming paper. By setting k = 1, the weak form of the finite element formulation can be written as follows:
qðw; v t Þ þ qðw; v gradv Þ þ ðgradw; 2gðvÞÞ ðdiv w; pÞ þ ðgradw; lBjp ðtÞ Þ þ ðq; divv Þ þ ðw; B;t;jp ðtÞ Þ T
þ ðw; v gradðBjp ðtÞ ÞÞ ðw; Bjp ðtÞ gradv Þ ðw; ðgradv Þ Bjp ðtÞ Þ þ ðw; fB1 ðBjp ðtÞ Þ½Bjp ðtÞ 1Þ þ ðSTB1 þ STB2 ; srÞ þ ðSTBCR ; sCR sÞ ¼ ðw; fÞ þ ðw; hÞCh :
ð33Þ
In the above equation, ‘’ symbolizes direct multiplication, w denotes the test function, v denotes the velocity, w is the test function for the stretch tensor field Bjp ðtÞ and is a second order tensor, s is the residual of Eq. (27) and is also a second order tensor, while fB1 ; fB2 ; STB1 ; STB2 , and STBCR are given by the following expressions:
fB1 ¼ 2
l ðtrBjp ðtÞ Þm ; a
ð34Þ
fB2 ¼ k ¼ 1;
ð35Þ
STB1 ¼ qðgradv w þ ðdivv Þw þ v gradwÞ þ gðgradðdiv wÞ þ DwÞ þ gradq; STB2 ¼ gradBjp ðtÞ : w þ Bjp ðtÞ : ðgradw þ gradðwT ÞÞ þ w ðdiv Bjp ðtÞ Þ þ ðdiv Bjp ðtÞ Þ w ; T STBCR ¼ v gradw þ ðdivv Þw þ gradv w þ w ðgradv Þ w : Bjp ðtÞ 1 GB1 fB1 w;
ð36Þ
GB1
dfB1 l ¼ ¼ 2 mðtrBjp ðtÞ Þm1 1 dBjp ðtÞ a
ð37Þ ð38Þ ð39Þ
and s; sCR are stabilization tensors. The function fB1 is the major difference between the proposed model and the constitutive model in Anand and Rajagopal (2004).2 2
The function fB1 for the constitutive model in Anand and Rajagopal (2004) is given by: fB1 ¼ 2KðtrBjp ðtÞ 3Þn .
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
85
Fig. 4. Geometries used in the numerical simulations.
In order to compare the numerical convergence of the model considered in Anand and Rajagopal (2004) and the proposed model, two geometric configurations in 3D are employed. The first test case is a straight tube (Fig. 4(a)), and the second test case is an arbitrary bent tube with a radius r = 3.5 mm as shown in Fig. 4(b). Inflow velocity is based on mean value of blood flow in the human carotid artery U 0 (=0.2 m/ s) and it serves as the reference inflow velocity. The inflow boundary has a parabolic velocity profile, and the problem is run with gradually increasing the mean inflow velocity. For both constitutive models, we used the same computational mesh and employed identical increments in the inflow mean velocity. The convergence of the residual of the momentum-balance equation, mass balance equation, and the viscoelastic constitutive equations during Newton–Raphson iterations for the straight tube is presented in Fig. 5. Fig. 5(a) shows the very first Newton–Raphson iteration loop for the Anand–Rajagopal model (Anand & Rajagopal, 2004). It can be seen that although at the beginning of the loop there is a reduction in the L2 -norm of the residual, it stagnates and does not converge to the desired pre-defined tolerance in subsequent iterations. Consequently, numerical simulations stop at this point. The new constitutive model presented in this paper shows optimal convergence trends during the nonlinear iterations. Several representative steps during the course of the simulation are shown in Fig. 5(b). The very first Newton–Raphson loop converges in only seven iterations, and this trend in convergence is maintained as the inflow velocity is gradually increased up to five times the inflow mean velocity U 0 . The super-linear convergence of the current model is attributed to the smoother transition in the internal constitutive variables at the low shear rates. Fig. 6 shows convergence of the residual for the bent tube geometry. Unlike the straight tube case presented above, the Anand–Rajagopal model diverges during the very first Newton–Raphson loop. On the other hand, the current constitutive model shows super-linear convergence all the way up to an inflow mean velocity of U in = 3.29U 0 . In order to understand the variation of the nonlinear relaxation time for different inflow velocities, the mean relaxation time (t rel ) at an inflow surface is computed as given below:
R
trel ¼ SRinflow
t rel dA
Sinflow
dA
:
ð40Þ
Recall that t rel is the non-linear relaxation time of the proposed model as given in Eq. (23). The Weissenberg number (Wi) and the Reynolds number corresponding to the flow are calculated using the above values of the mean relaxation time and the mean inflow velocity (U 0 ) in Eqs. (41) and 42, respectively.
t rel U 0 2r 2qU 0 r Re ¼ : g1 =2
Wi ¼
ð41Þ ð42Þ
Fig. 7(a) shows the relation of the mean relaxation time versus the mean inflow velocity: this was computed for simulations in the straight tube given in Fig. 4(a). As the magnitude of the inflow velocity increases, the mean relaxation time rapidly decreases from 0.0811 s to 0.00291 s. Fig. 7(b) shows variations of Wi and Re when mean inflow velocity increases from 0.1U 0 to 5U 0 . As inflow velocity increases, Wi increases slightly from 0.22 to 0.395, while Re increases linearly from 28 to 1,400.
86
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
Fig. 5. Convergence of the residual for the straight tube geometry.
Fig. 6. Convergence of the residual for the bent tube geometry.
Fig. 7. (Proposed Model) Variation of mean relaxation time (t rel ), Wi, & Re with mean inflow velocity (U 0 ) in a straight tube.
The proposed model has a propensity to reduce the relaxation time for increasing inflow velocity. Therefore, the Weissenberg number stays within a narrow range as shown in Fig. 7(b) even when the inflow velocity increases 50 times. As the velocity increases linearly, the viscous stress also increases rapidly, but the viscoelastic stress stays bounded. Consequently,
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
87
the contribution of the viscoelastic stress to the total stress decreases rapidly as the flow velocity is increased. This trend agrees with the observation from hemorheology that viscoelastic fluid response is dominant in the regions of low velocity (equivalently shear rate), while it is negligible in the flow domain with high velocity. 6. Discussion We have motivated the need for, and developed, a new model to describe the rheological characteristics of blood (namely shear-thinning and deformation-dependent viscoelasticity) in both steady and unsteady flows, respectively. Further, we have corroborated the model with a set of data for both types of flow: apparent viscosity data under steady flow conditions, and shear stress & hysteresis data under time-dependent shear flow conditions. This exercise has yielded a set of parameters that can be used to characterize blood flow in bent tube geometry as well as in more complex 3D geometries likely to be encountered in the human vasculature. The parameter set is especially relevant to simulate the flow in a patient-specific aneurysm wherein there is both a wide range of shear rates as well as pulsatile flow conditions. With this objective in mind, we have simulated the flow in a straight tube and a bent tube and shown the superior convergence properties of the model. By way of comparison, the model in Anand and Rajagopal (2004) displays much poorer convergence characteristics in the straight tube, and diverges in the bent tube. Further, the simulations in a straight tube show that the viscoelastic properties of the proposed model dominate at low mean velocities but are negligible at high mean velocities: this is further corroboration of experimental data in blood rheology. We opine that researchers intending to use a model in extensive numerical simulations of blood flow will find the proposed model more amenable to study than the model in Anand and Rajagopal (2004): this holds particularly true for very low Reynolds number flows, and start-up flows. We will take up the features of the proposed model’s flow in a subsequent paper. Acknowledgement MA thanks Dr. C. S. Sastry for useful discussions. JK was partially funded by a grant from the NCSA, and this support is acknowledged. We thank Professor K. R. Rajagopal for discussions on the presentation of the results. References Anand, M., & Rajagopal, K. R. (2004). A shear-thinning viscoelastic fluid model for describing the flow of blood. International Journal of Cardiovascular Medicine and Science, 4(2), 59–68. Bodnar, T., Sequeira, A., & Pirkl, L. (2009). Numerical simulations of blood flow in a stenosed vessel under different flow rates using a generalized Oldroyd-B model. In Numerical Analysis and Applied Mathematics (vol. 2, pp. 645–648) Melville, New York. Bodnar, T., Rajagopal, K. R., & Sequeira, A. (2011). Simulation of the three-dimensional flow of blood using a shear-thinning viscoelastic fluid model. The Mathematical Modelling of Natural Phenomena, 6(5), 1–24. Bodnar, T., Sequeira, A., & Prosi, M. (2011). On the shear-thinning and viscoelastic effects of blood flow under various flow rates. Applied Mathematics and Computation, 217(11), 5055–5067. Brooks, A. N., & Hughes, T. J. R. (1982). Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 32, 199–259. Bureau, M., Healy, J., Bourgoin, D., & Joly, M. (1978). Etude expérimentale in vitro du comportment rhéologique du sang et régime transitoire à faible vitesse de cisaillement. Rheologica Acta, 17(6), 612–625. Bureau, M., Healy, J., Bourgoin, D., & Joly, M. (1980). Rheological hysteresis of blood at low shear rate. Biorheology, 17(1–2), 191–203. Chien, S., Usami, S., Taylor, H. M., Lundberg, J. L., & Gregersen, M. I. (1966). Effect of hematocrit and plasma proteins on human blood rheology at low shear rates. Journal of Applied Physiology, 21(1), 81–87. Copley, A. L., & King, R. G. (1976). On the viscoelasticity of anticoagulated whole human blood in steady shear as tested by rheogoniometric measurements of normal forces. Biorheology, 12, 5–10. Coronado, O. M., Arora, D., Behr, M., & Pasuali, M. (2006). Four-field Galerkin/least-squares formulation for viscoelastic fluids. Journal of Non-Newtonian Fluid Mechanics, 140, 132–144. Guenette, R., & Fortin, M. (1995). A new mixed finite element method for computing viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 60, 27–52. Hughes, T. J. R. (1995). Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127, 387–401. Hughes, T. J. R., Feijoo, G. R., Mazzei, L., & Quincy, J. B. (1998). The variational multiscale method- a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 166, 3–24. Inzoli, F., Boschetti, F., Zappa, M., Longo, T., & Fumero, R. (1993). Biomechanical factors in abdominal aortic aneurysm rupture. European Journal of Vascular and Endovascular Surgery, 7(6), 667–674. Kwack, J., & Masud, A. (2010). A three-field formulation for incompressible viscoelastic fluids. International Journal of Engineering Science, 48, 1413–1432. Masud, A., & Calderer, R. (2009). A variational multiscale stabilized formulation for the incompressible Navier–Stokes equations. Computational Mechanics, 44, 145–160. Masud, A., & Franca, L. P. (2008). A hierarchical multiscale framework for problems with multiscale source terms. Computer Methods in Applied Mechanics and Engineering, 197, 2692–2700. Masud, A., & Khurram, R. (2006). A multiscale finite element method for the incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 195, 1750–1777. McMillan, D. E., Utterback, N. G., Nasrinasrabadi, M., & Lee, M. M. (1986). An instrument to evaluate the time dependent flow properties of blood at moderate shear rates. Biorheology, 23, 63–74. Owens, R. G. (2006). A new microstructure-based constitutive model for human blood. Journal of Non-Newtonian Fluid Mechanics, 140(1–3), 57–70. Raghavan, M. L., Vorp, D. A., Federle, M. A., Makaroun, M. S., & Webster, M. W. (2000). Wall stress distribution on three dimensionally reconstructed models of human abdominal aortic aneurysm. Journal of Vascular and Endovascular Surgery, 31(4), 760–769. Rajagopalan, D., Armstrong, R. C., & Brown, R. A. (1990). Finite element methods for calculation of steady, viscoelastic flow using constitutive equations with a Newtonian viscosity. Journal of Non-Newtonian Fluid Mechanics, 36, 159–192. Rajagopal, K. R., & Srinivasa, A. R. (2000). A thermodynamic framework for rate-type fluids. Journal of Non-Newtonian Fluid Mechanics, 88(3), 207–227.
88
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
Srinivasa, A. R. (2000). Flow characteristics of a multiconfigurational shear-thinning viscoelastic fluid with reference to the orthogonal rheometer. Theoretical and Computational Fluid Dynamics, 13(5), 305–317. Sun, N., & DeKee, D. (2001). Simple shear, hysteresis, and yield stress in biofluids. The Canadian Journal of Chemical Engineering, 79, 36–41. Thurston, G. B. (1975). Elastic effects in pulsatile blood flow. Microvascular Research, 9, 145–157. Yeleswarapu, K. K. (1996). Evaluation of Continuum Models for characterizing the constitutive behavior of blood, Ph.D. Thesis, University of Pittsburgh, Pittsburgh, PA. Zarins, C. K., Giddens, D. P., Bharadvaj, B. K., Sottiurai, V. S., Mabon, R. F., & Glagov, S. (1983). Carotid bifurcation atherosclerosis. Quantitative correlation of plaque localization with flow velocity profiles and wall shear stress. Circulation Research, 53(4), 502–514.
Contents lists available at SciVerse ScienceDirect
International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci
A new generalized Oldroyd-B model for blood flow in complex geometries M. Anand a,⇑, J. Kwack b, A. Masud b a b
Department of Chemical Engineering, Indian Institute of Technology Hyderabad, Yeddumailaram, Andhra Pradesh 502205, India Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
a r t i c l e
i n f o
Article history: Received 21 May 2013 Received in revised form 18 June 2013 Accepted 21 June 2013 Available online 3 August 2013 Keywords: Blood Oldroyd-B model Thermodynamic framework Variational multiscale method Convergence
a b s t r a c t A new generalization of the Oldroyd-B model is developed to describe the flow of blood. The model is developed within a thermodynamic framework which recognizes that a viscoelastic fluid can remain stress free in multiple configurations. The new model is an improvement over an earlier model developed within the same framework to describe the response characteristics of blood. It captures the shear-thinning and deformationdependent viscoelastic behavior of blood just like the previous model. More importantly, unlike the previous model, it does not have the shortcoming of an abrupt transition of the material properties at low shear rates: instead, it allows for a smooth variation of the rate of dissipation, and therefore viscosity, over the entire range of physically feasible shear-rates. This feature is very attractive for developing high-fidelity numerical methods for application to the complex geometries that are typically encountered in the human vasculature. Convergence of the numerical method in simple geometries shows its superior properties as compared to the earlier model: this demonstration of model performance is a precursor to its use in 3D geometries. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Human blood is a suspension of formed elements (erythrocytes, leukocytes, and platelets, that together constitute 45% by volume) in an aqueous solution (plasma). It exhibits complex rheological behavior namely shear-thinning by orders of magnitude (Chien, Usami, Taylor, Lundberg, & Gregersen, 1966), and deformation-dependent viscoelastic behaviour (recorded as shear-rate dependent relaxation times (McMillan, Utterback, Nasrinasrabadi, & Lee, 1986) or normal stress differences (Copley & King, 1976)). An accurate model that captures this range of behavior in a single framework can be invaluable in generating precise knowledge of flow patterns and stress distributions in the heart and the vasculature. Such a model can be a tool for patient-specific simulations wherein the geometry information is gleaned from either angiograms or intra-vascular-ultrasound images. Patient-specific simulations are an important tool in tracking the progress of cardiovascular pathologies and identifying candidates for further medical intervention. Examples of such cardiovascular pathologies include aneurysms and atherosclerotic plaques: these pathologies are influenced in their development by the flow patterns of blood (e.g. aneurysms (Inzoli, Boschetti, Zappa, Longo, & Fumero, 1993; Raghavan, Vorp, Federle, Makaroun, & Webster, 2000), and atherosclerosis (Zarins et al., 1983)). Let us consider the case of an (intra-cranial saccular) aneurysm: the dilation of the artery wall due to mechanical weakness leads to the formation of a balloon-shaped aneurysm with a local circulation (slow flow with low shear rates) within it. Calculating the stresses imposed by blood flow on the walls of the aneurysm is a key measure in determining how much further it will stretch, and a constitutive model for blood used in this simulation ⇑ Corresponding author. Tel.: +91 40 23016090; fax: +91 40 23016032. E-mail addresses: [email protected] (M. Anand), [email protected] (J. Kwack), [email protected] (A. Masud). 0020-7225/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijengsci.2013.06.009
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
79
must be accurate at the low shear rates that are encountered in this geometry. Further, the dilatation occurs adjacent to the mainstream flow of the cerebral artery where shear rates are high. The constitutive model has to capture stress distributions in the high shear rate regimes of the mainstream cerebral flow as well. This example illustrates the need for a constitutive model for blood that makes accurate predictions of flow over a wide range of shear rates: from extremely low values (at near stagnant conditions) to high values. A systematic thermodynamic framework for the development of rate-type fluid models was presented by Rajagopal and Srinivasa (2000). This thermodynamic framework results from the melding together of two central ideas. The first idea is that a body’s ‘‘natural configuration’’ (the configuration in the absence of external stimuli) evolves as the body undergoes a process, and that it accesses different sets of these configurations depending on the class of processes it is subject to. The stress field of the body in the current configuration is determined with respect to its deformation from the natural configurations for the particular process. The second idea is that the form of the stress tensor can be determined uniquely by specifying the forms- in terms of the kinematic variables and temperature - for the stored energy function and the rate of dissipation, and then maximising the rate of dissipation subject to relevant constraints. The constraints being that (i) the rate of dissipation is equal to the difference between the stress power and the rate of change of the stored energy, and (ii) the material is incompressible. Using such a framework, Anand and Rajagopal (2004) developed a model (which we refer to from now on as the Anand–Rajagopal model) that was intended to describe the flow characteristics of blood. The Anand–Rajagopal model was one among a few other shear-thinning viscoelastic models, for describing blood flow that were developed within a few years of each other. These other models include those by Sun and DeKee (2001), and by Owens (2006). The model in Sun and DeKee (2001) uses a shear rate-dependent structural parameter combined with a standard frame-invariant convected derivative to incorporate the shear-rate dependent viscosity and relaxation times. The model in Owens (2006) also uses a shear rate-dependent structural parameter to model the viscosity and relaxation time of blood, but the parameter is derived using theories in polymer network formation. Owens considered the erythrocyte as analagous to a monomer, and the erythrocyte aggregates (rouleaux) in blood as analogous to a ‘k’-mer (a network of ‘k’ cross-linked monomer units). The difference between these other models and the Anand–Rajagopal model lies in the framework adopted by the latter: they used the microstructural evolution of erythrocyte aggregates (rouleaux) and the elasticity of the erythrocyte to motivate forms for the stored energy function and rate of dissipation characterizing blood. The Anand–Rajagopal model was a generalization of the classical Oldroyd-B model that captured the rheological behavior of blood over a wide range of shear rates (0.06 s1 to 650 s1 ). However, the model equations, by themselves, presented the problem of infinite viscosity being recorded at zero-shear rate. In order to correct this anomaly, a low-shear rate criterion that transitioned between two forms of dissipation (and apparent viscosity forms), was included to ensure a finite zero-shear viscosity and to capture viscosity variation in the low shear-rate range (0 s1 to 0.06 s1 ). The Anand–Rajagopal model has recently been used to simulate blood flow in a straight tube with a stenosis, as well as in a 90° curved tube. Bodnar, Rajagopal, and Sequeira (2011) used a finite-volume discretization along with an explicit Runge– Kutta time marching scheme: they thereby obtained numerical results at three flow rates (Reynolds number less than 100) that illustrated the dynamics and flow patterns generated. They found that in a straight tube with stenosis the relative difference of the axial velocities, between that predicted by the Anand–Rajagopal model and that by the Newtonian model, increased as the flow rate decreased. They also found that in a 90° curved tube the axial velocity contour is pushed away from the central axis, and there is a slowing down of velocity towards the outer radius along with a corresponding speed-up towards the inner radius. The secondary vortices that are set up are slower than the corresponding ones in a Newtonian fluid due to the higher viscosity of the model. Further, they found that there is a flattening of the velocity at the central core, both in the stenosed pipe and the curved pipe, due to the shear-thinning nature of the model. This was a reiteration of the results that their group had obtained earlier in Bodnar, Sequeira, and Prosi (2011) and Bodnar, Sequeira, and Pirkl (2009) using an Oldroyd-B model with a Cross model-type shear-thinning viscosity. However, they were not concerned with the convergence rate of the simulations using the model. We performed numerical tests and found difficulties in the numerical implementation of the Anand–Rajagopal model (along with the low shear rate criterion) – specifically non-convergence of the residual of the mass conservation, momentum balance, and constitutive equations when used in a bent tube geometry. Our investigations have shown that it is non-trivial to obtain an appropriate low-shear criterion1 that will be acceptable from a modeling standpoint, and will also satisfy the requirements of the particular numerical scheme. In this paper, we therefore introduce a new model which has a finite zeroshear viscosity and allows for a smooth variation of viscosity at all possible shear-rates from zero-upwards: we corroborate the model with existing data, and generate parameters to fit this data. We then show numerical convergence plots for a simplified version of the model in two sample geometres (straight cylindrical tube, and bent tube) to demonstrate the potential of the model in the numerical scheme as well as to set the stage for simulations in more complex geometries. The rest of the paper is organized as follows: The notation used to describe the kinematics of a viscoelastic fluid with multiple stress-free states is given in Section 2. We motivate the need for, and propose, the new model in Section 3. The new model is subject to a corroboration procedure in Section 4 to get a set of constants that will best fit experimental data. Numerical results for sample geometries are documented in Section 5 to demonstrate the promise of the new model for 3D simulations. The results are discussed in Section 6.
1
Note that we do not say that such a criterion does not exist.
80
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
2. Preliminaries Consider the following schematic (Fig. 1) of the deformation of a viscoelastic fluid body with a single relaxation mechanism, and capable of instantaneous elastic response. In Fig. 1, jR ðBÞ and jt ðBÞ denote the reference and the current configuration (at time t) of the body B, respectively. jpðtÞ ðBÞ denotes the stress-free configuration that is reached by instantaneously unloading the body from its current configuration. The motion of the body is defined through a one-to-one mapping between each point in the reference configuration (X 2 jR ðBÞ) and a corresponding point in the current configuration at time instant t (x 2 jt ðBÞ):
x ¼ vjR ðX; tÞ:
ð1Þ
We assume the motion to be sufficiently smooth and invertible. We will also, for convenience, suppress B in the notation
jt ðBÞ etc. The deformation gradient FjR , and the left Cauchy-Green stretch tensor BjR are defined as
FjR ¼
@ vjR ; @XjR
and BjR ¼ FjR FTjR :
ð2Þ
The tensor FjpðtÞ denotes the deformation gradient between the instantaneous stress-free (natural) configuration and the current configuration, provided the deformation is homogeneous. Further, the left Cauchy-Green stretch tensor associated with the instantaneous elastic response from the stress-free configuration jpðtÞ is then given by:
BjpðtÞ ¼ FjpðtÞ FTjpðtÞ :
ð3Þ
The mapping G is defined by:
G ¼ FjR !jpðtÞ ¼ F1 jpðtÞ FjR :
ð4Þ
The velocity gradients L and LjpðtÞ are defined by:
_ 1 L ¼ F_ jR F1 jR ; andLjpðtÞ ¼ GG :
ð5Þ
where the dot represents the material time derivative. The symmetric parts of L and LjpðtÞ are obtained in the usual manner, and are represented by the symbols D and DjpðtÞ , respectively. r The upper convected Oldroyd derivative of BjpðtÞ ; BjpðtÞ , is defined as: r
BjpðtÞ ¼ B_ jpðtÞ LBjpðtÞ BjpðtÞ LT :
ð6Þ
Additional constraints, that are enforced during the afore-mentioned maximization procedure for the rate of dissipation, are those due to incompressibility (of blood); these imply:
trðDÞ ¼ 0
ð7Þ
trðDjpðtÞ Þ ¼ 0:
ð8Þ
and
Fig. 1. Configurations accessed by a viscoelastic fluid body with instantaneous elasticity.
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
81
3. Constitutive model for blood Recent numerical tests have shown that the low shear-rate criterion given in Anand and Rajagopal (2004) is not suitable for numerical implementation in complex geometries. This is because it does not allow for convergence of the residuals of the numerical equations. The abrupt transition (using a Heaviside function) in the model terms for relaxation time and shear-thinning index (the parameters K and n, respectively, in the model equations: see Anand & Rajagopal, 2004) has been pin-pointed as the reason for the non-convergence. This abrupt transition results in numerical oscillations in the internal fields described by these terms. Because of the fully coupled nature of the system of equations, these internal ripples affect the computed fields and lead to the divergence of the solution from the actual solution. For these reasons, we tried to develop an alternate low-shear rate criterion which would both ensure a smooth transition in relaxation time and shear-thinning index, through the cut-off shear-rate of 0.06 s1 , as well as ensure a finite apparent viscosity at zero shear rate. To get a smooth transition in the model equations of Anand and Rajagopal (2004), K and n must be functions of Bjp ðtÞ . Imposing constraints from both a modeling standpoint and a numerical standpoint resulted in a procedure to specify K(Bjp ðtÞ ) and n(Bjp ðtÞ ). However, we found that it is not trivial to identify appropriate forms for K and n: we tested power law forms among others (which are natural choices for shear-thinning behavior) but were unable to satisfy the constraints in the procedure. We therefore decided not to seek further for acceptable functional forms for the relaxation time and shearthinning index, and instead proceeded to construct a new model which does not pose such numerical difficulties requiring additional workarounds. The forms for the stored energy, and rate of dissipation, proposed for the new model are:
W¼
l
ðIB 3Þ; 2 n ¼ g Djp ðtÞ Bjp ðtÞ Djp ðtÞ þ g1 D D:
ð9Þ ð10Þ
We choose:
g ¼ g^ ðBjp ðtÞ Þ ¼ a trBjp ðtÞ
m
;
m > 0:
ð11Þ
We point out that the rate of dissipation contains a power-law term, and a viscous term (akin to a Newtonian fluid), whereas the stored energy form is that of the neo-Hookean type. Srinivasa (2000) used a pure power-law term in the rate of dissipation expression along with a neo-Hookean type stored energy function, and developed a Maxwell-type model from these forms. The model proposed here uses a different method compared to Srinivasa (2000) to incorporate the shear-thinning viscosity: namely, a power-law dependence in the rate of ^ ðBjp ðtÞ Þ, instead of a power-law dependence on the term Djp ðtÞ BjpðtÞ Djp ðtÞ . dissipation coefficient g ¼ g The model equations derived using the constrained maximisation procedure in Rajagopal and Srinivasa (2000) are:
T ¼ p1 þ S;
ð12Þ
S ¼ lBjp ðtÞ þ g1 D; r 2l Bjp ðtÞ ¼ Bjp ðtÞ k1 ;
ð13Þ
g
k¼
3 trB1 jp ðtÞ
:
ð14Þ ð15Þ
4. Model corroboration In this section we corroborate the model with experimental data for both apparent viscosity, as well as a shear flow loop hysteresis, so as to come up with a set of model parameters that matches both sets of experimental data. Ensuring that the model fits both sets of data (apparent viscosity, and shear flow loop hysteresis) allows us to identify a unique set of parameters for the model. 4.1. Apparent viscosity We consider the flow between two concentric cylinders separated by a small gap, with the inner cylinder being held at rest. The flow between the cylinders is assumed to be of the form (we use some of the notation that is used in Anand & Rajagopal (2004) to enable easy comparison):
v ¼ uðrÞ^eh ¼ rxðrÞ^eh ;
ð16Þ
which automatically satisfies the balance of mass. Substituting the constitutive Eqs. (12)–(15) in the equations for balance of linear momentum, and assuming that the stress field depends only on ‘r’, we get:
82
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
gk þ g1
T rh ¼
2
du u : dr r
ð17Þ
In a concentric-cylinder rheometer, it is possible to approximate the shear-rate term in the following manner:
du u dx ¼ c_ ¼ R0 : dr r dr
ð18Þ
The expression for apparent viscosity (gapp ) for this model is thus:
gapp ¼
gk þ g1 2
ð19Þ
;
where g is a function to be determined using Eq. (11). k is determined using the incompressibility condition (detBjp ðtÞ = 1) as:
1 k¼ 1=3 : g2 c_ 2 1 þ 4l2
ð20Þ
Noting that:
2c_ 2 g2 1 trBjp ðtÞ ¼ 3 þ 1=3 ; 4l2 g2 c_ 2 1 þ 4l2
ð21Þ
We obtain the following non-linear equation for g to be solved at each value of c_ :
2 3m 2g2 c_ 2 1 6 7 g ¼ a4 3 þ 1=3 5 : 4l2 g2 c_ 2 1 þ 4l2
ð22Þ
We fix the four unknown parameters (a; m; g1 ; l) by minimising the error between predicted viscosity (Eq. 19), and measured viscosity at a set of shear-rates (data reported in Yeleswarapu (1996)). The minimization is done using the MATLAB function ‘fminsearch’, and a residual (sum of least square of error for both apparent viscosity data and shear loop hysteresis data) of 0.0052 is reported. Since gk ¼ 0 as c_ ! 1, we set g1 ¼ 2g1 ¼ 0:01 Pa s. There are multiple sets of parameters that fit the data equally well, and one set that we find is a = 953.655 Pa s, m = 7.7779, l = 0.0283 Pa, g1 = 0.01 Pa s; this set predicts a zero-shear viscosity of g0 = 0.0978 Pa s. The data fit for these parameters is given in Fig. 2: we see that the proposed model fits the data as well as the Anand–Rajagopal model (which uses the parameters: K = 1.2056 s1 , n = 0.7525, l = 0.0227 Pa, g1 = 0.01 Pa s). We note that the relaxation time associated with the evolution of Bjp ðtÞ is given by:
trel ¼
g aðtrBjp ðtÞ Þm ¼ : 2l 2l
ð23Þ
A zero-shear relaxation time of 3.278 s is obtained for the proposed model, and this compares favorably with the value of 3.022 s obtained for the Anand–Rajagopal model. −1
10
Apparent Viscosity (Pa⋅ s)
Experimental data Anand−Rajagopal model Proposed model
−2
10
−3
10
−2
10
−1
10
0
10
1
10
2
10
3
10
Shear rate (sec−1)
Fig. 2. Apparent viscosity vs c_ for proposed model: comparison with experimental data in Yeleswarapu (1996).
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
83
4.2. Shear loop hysteresis Experimental data is available for the shear stress at various instants during, and residual shear stress (hysteresis) at the end of, the application of a triangular loop of shear rate in Bureau, Healy, Bourgoin, and Joly (1980) (see also Bureau, Healy, Bourgoin, & Joly, 1978). In order to corroborate the proposed model with this data, we first derive the equations governing it’s time variation of shear stress in a time-dependent Couette flow. We then compare the model predictions with the experimental data for a triangular loop of shear rate. The applied shear rate varies in the following (triangular) manner:
c_ ðtÞ ¼ c€c t; t 6 t0 ; c_ ðtÞ ¼ 2c€c t0 c€c t; t0 < t 6 2t0 :
ð24Þ ð25Þ
Here, for the data in Bureau et al. (1980), we have c€c ¼ 0:0185s2 , and t0 = 6.5 s. We seek a time-dependent Couette flow solution for the velocity field:
v ¼ c_ ðtÞy^ex :
ð26Þ
Substituting the above velocity field in the constitutive model, and assuming that the stress is a function of only ‘y’, we obtain:
dBxy 2l ¼ B þ c_ Byy ; dt g xy dBxx 2l ¼ ðB kÞ þ 2c_ Bxy ; dt g xx dByy 2l ¼ ðB kÞ; dt g yy
ð27Þ ð28Þ ð29Þ
where:
Sxy ¼ lBxy þ
g1 c_
ð30Þ
; 2 g ¼ g^ ðBjp ðtÞ Þ ¼ aðBxx þ 2Byy Þm ; k¼
ð31Þ
3Byy ðBxx Byy B2xy Þ B2yy þ 2Bxx Byy B2xy
ð32Þ
and the other stress components are given by Bzz ¼ Byy ; Bxz ¼ Byz ¼ 0. The stress components at initial time (t = 0) are calculated by setting the initial elastic stretch to be unity (Bjp ðtÞ ¼ 1), and the initial shear rate to be zero (c_ ¼ 0s1 ). Hence, Bxy ðt ¼ 0Þ ¼ 0; Bxx ðt ¼ 0Þ ¼ 1, and Byy ðt ¼ 0Þ ¼ 1. We use a simple ExplicitEuler scheme to solve this ODE-IVP, and obtain the stress components at each instant, t, in 0 6 t 6 2t 0 . Although Runge–Kutta methods of higher accuracy are available to solve the ODE-IVP, we have settled for the simplest method which yields a solution. We only need to ensure that the results are grid independent to have confidence in the results from the scheme. We determined Dt = 0.01 s to be the smallest time step for which the results are reproducible, and this time step is sufficient for our purpose of comparison with the data. −3
8
x 10
7
Shear Stress (Pa)
6 5 4 3 2 Experimental data Anand−Rajagopal model Proposed model
1 0
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
−1
Shear rate (sec )
Fig. 3. Shear stress variation with time during, and hysteresis at the end of, a triangular loop of shear rate. Data from Bureau et al. (1980).
84
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
The shear stress, Sxy , predicted by the model is compared with the experimental data, and the parameters that give the best match with the relevant data in Bureau et al. (1980) (as well as for the apparent viscosity data) are reported. The parameters are thus the same as those given earlier for the apparent viscosity data, namely a = 953.655 Pa s, m = 7.7779, l = 0.0283 Pa, g1 = 0.01 Pa s. The predictions from both models are compared with the data in Fig. 3: we see that the match between the proposed model and data is excellent, while that between the Anand–Rajagopal model and data is not as good. We hasten to point out that the less than exemplary data match with this one set of data does not in any way detract from the capability of the Anand–Rajagopal model in its ability to model blood. The reason for the poor fit is that the parameters used for the Anand–Rajagopal model were obtained by minimizing the error between model predictions with both apparent viscosity data (used in subSection 4.1), and the data for oscillatory flow in a pipe given in Thurston (1975) (which is different from the data used in this subsection). Therefore, to get a better fit using the Anand–Rajagopal model, we have to rerun the optimizations to get best fit with the two sets of data given in this paper. However, we will not undertake this exercise, as the object of this paper is to show the improved efficacy of the proposed model in high-fidelity simulations in complex geometries. Obtaining new parameters for the Anand–Rajagopal model will still not improve it’s performance in such simulations, as we shall show in the next section, because of the abrupt transition of properties at low shear rates which cannot be avoided. We thus use the same parameters for each model used in Figs. 2 and 3 to demonstrate the numerical results as well. 5. Stabilized three-field formulation and numerical results A three-field multiscale finite element method for incompressible viscoelastic fluids was developed by Kwack and Masud (2010). This method is comprised of velocity, pressure and conformation-tensor as the unknown fields. The method finds roots in the Variational Multiscale framework proposed by Hughes (1995) and Hughes, Feijoo, Mazzei, and Quincy (1998), Masud and Khurram (2006), Masud and Franca (2008) and Masud and Calderer (2009), along with SUPG (Streamlined-Upwind Petrov Galerkin) type ideas applied to the constitutive model Brooks and Hughes (1982). High Weissenberg number flows were investigated in Kwack and Masud (2010) and the computed results were compared with four-field formulations wherein rate-of-deformation tensor was the fourth field. These competing methods include the EVSS (Elastic-Viscous Split Stress) (Rajagopalan, Armstrong, & Brown, 1990), DEVSS methods (Discrete Elastic-Viscous Split Stress) (Guenette & Fortin, 1995), and the four-field GLS (Galerkinn Least-squares) (Coronado, Arora, Behr, & Pasuali, 2006) methods. It is important to note that the three-field method (Kwack & Masud, 2010) leads to ten degrees of freedom per node in 3D calculations, while the four-field formulations result in sixteen degrees of freedom per node. Consequently, not only the bandwidth for the three-field method is smaller as compared to the four field methods, it also leads to a substantially small assembled system for the same spatial discretization. The enhanced stability properties of the three-field formulation (Kwack & Masud, 2010) help resolve the flow with the same accuracy as the competing methods, however at substantially reduced computational costs. Both qualitatively and quantitatively, the three-field method shows excellent agreement with the published results on different benchmark problems. In several cases considered, higher Weissenberg numbers were reached than were reported previously in the literature. We have now extended the ideas presented in Kwack and Masud (2010) to full VMS (Variational Multiscale) based three-field formulation that results in enhanced numerical stability in 3D simulations: Details of this formulation will be presented in a forthcoming paper. By setting k = 1, the weak form of the finite element formulation can be written as follows:
qðw; v t Þ þ qðw; v gradv Þ þ ðgradw; 2gðvÞÞ ðdiv w; pÞ þ ðgradw; lBjp ðtÞ Þ þ ðq; divv Þ þ ðw; B;t;jp ðtÞ Þ T
þ ðw; v gradðBjp ðtÞ ÞÞ ðw; Bjp ðtÞ gradv Þ ðw; ðgradv Þ Bjp ðtÞ Þ þ ðw; fB1 ðBjp ðtÞ Þ½Bjp ðtÞ 1Þ þ ðSTB1 þ STB2 ; srÞ þ ðSTBCR ; sCR sÞ ¼ ðw; fÞ þ ðw; hÞCh :
ð33Þ
In the above equation, ‘’ symbolizes direct multiplication, w denotes the test function, v denotes the velocity, w is the test function for the stretch tensor field Bjp ðtÞ and is a second order tensor, s is the residual of Eq. (27) and is also a second order tensor, while fB1 ; fB2 ; STB1 ; STB2 , and STBCR are given by the following expressions:
fB1 ¼ 2
l ðtrBjp ðtÞ Þm ; a
ð34Þ
fB2 ¼ k ¼ 1;
ð35Þ
STB1 ¼ qðgradv w þ ðdivv Þw þ v gradwÞ þ gðgradðdiv wÞ þ DwÞ þ gradq; STB2 ¼ gradBjp ðtÞ : w þ Bjp ðtÞ : ðgradw þ gradðwT ÞÞ þ w ðdiv Bjp ðtÞ Þ þ ðdiv Bjp ðtÞ Þ w ; T STBCR ¼ v gradw þ ðdivv Þw þ gradv w þ w ðgradv Þ w : Bjp ðtÞ 1 GB1 fB1 w;
ð36Þ
GB1
dfB1 l ¼ ¼ 2 mðtrBjp ðtÞ Þm1 1 dBjp ðtÞ a
ð37Þ ð38Þ ð39Þ
and s; sCR are stabilization tensors. The function fB1 is the major difference between the proposed model and the constitutive model in Anand and Rajagopal (2004).2 2
The function fB1 for the constitutive model in Anand and Rajagopal (2004) is given by: fB1 ¼ 2KðtrBjp ðtÞ 3Þn .
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
85
Fig. 4. Geometries used in the numerical simulations.
In order to compare the numerical convergence of the model considered in Anand and Rajagopal (2004) and the proposed model, two geometric configurations in 3D are employed. The first test case is a straight tube (Fig. 4(a)), and the second test case is an arbitrary bent tube with a radius r = 3.5 mm as shown in Fig. 4(b). Inflow velocity is based on mean value of blood flow in the human carotid artery U 0 (=0.2 m/ s) and it serves as the reference inflow velocity. The inflow boundary has a parabolic velocity profile, and the problem is run with gradually increasing the mean inflow velocity. For both constitutive models, we used the same computational mesh and employed identical increments in the inflow mean velocity. The convergence of the residual of the momentum-balance equation, mass balance equation, and the viscoelastic constitutive equations during Newton–Raphson iterations for the straight tube is presented in Fig. 5. Fig. 5(a) shows the very first Newton–Raphson iteration loop for the Anand–Rajagopal model (Anand & Rajagopal, 2004). It can be seen that although at the beginning of the loop there is a reduction in the L2 -norm of the residual, it stagnates and does not converge to the desired pre-defined tolerance in subsequent iterations. Consequently, numerical simulations stop at this point. The new constitutive model presented in this paper shows optimal convergence trends during the nonlinear iterations. Several representative steps during the course of the simulation are shown in Fig. 5(b). The very first Newton–Raphson loop converges in only seven iterations, and this trend in convergence is maintained as the inflow velocity is gradually increased up to five times the inflow mean velocity U 0 . The super-linear convergence of the current model is attributed to the smoother transition in the internal constitutive variables at the low shear rates. Fig. 6 shows convergence of the residual for the bent tube geometry. Unlike the straight tube case presented above, the Anand–Rajagopal model diverges during the very first Newton–Raphson loop. On the other hand, the current constitutive model shows super-linear convergence all the way up to an inflow mean velocity of U in = 3.29U 0 . In order to understand the variation of the nonlinear relaxation time for different inflow velocities, the mean relaxation time (t rel ) at an inflow surface is computed as given below:
R
trel ¼ SRinflow
t rel dA
Sinflow
dA
:
ð40Þ
Recall that t rel is the non-linear relaxation time of the proposed model as given in Eq. (23). The Weissenberg number (Wi) and the Reynolds number corresponding to the flow are calculated using the above values of the mean relaxation time and the mean inflow velocity (U 0 ) in Eqs. (41) and 42, respectively.
t rel U 0 2r 2qU 0 r Re ¼ : g1 =2
Wi ¼
ð41Þ ð42Þ
Fig. 7(a) shows the relation of the mean relaxation time versus the mean inflow velocity: this was computed for simulations in the straight tube given in Fig. 4(a). As the magnitude of the inflow velocity increases, the mean relaxation time rapidly decreases from 0.0811 s to 0.00291 s. Fig. 7(b) shows variations of Wi and Re when mean inflow velocity increases from 0.1U 0 to 5U 0 . As inflow velocity increases, Wi increases slightly from 0.22 to 0.395, while Re increases linearly from 28 to 1,400.
86
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
Fig. 5. Convergence of the residual for the straight tube geometry.
Fig. 6. Convergence of the residual for the bent tube geometry.
Fig. 7. (Proposed Model) Variation of mean relaxation time (t rel ), Wi, & Re with mean inflow velocity (U 0 ) in a straight tube.
The proposed model has a propensity to reduce the relaxation time for increasing inflow velocity. Therefore, the Weissenberg number stays within a narrow range as shown in Fig. 7(b) even when the inflow velocity increases 50 times. As the velocity increases linearly, the viscous stress also increases rapidly, but the viscoelastic stress stays bounded. Consequently,
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
87
the contribution of the viscoelastic stress to the total stress decreases rapidly as the flow velocity is increased. This trend agrees with the observation from hemorheology that viscoelastic fluid response is dominant in the regions of low velocity (equivalently shear rate), while it is negligible in the flow domain with high velocity. 6. Discussion We have motivated the need for, and developed, a new model to describe the rheological characteristics of blood (namely shear-thinning and deformation-dependent viscoelasticity) in both steady and unsteady flows, respectively. Further, we have corroborated the model with a set of data for both types of flow: apparent viscosity data under steady flow conditions, and shear stress & hysteresis data under time-dependent shear flow conditions. This exercise has yielded a set of parameters that can be used to characterize blood flow in bent tube geometry as well as in more complex 3D geometries likely to be encountered in the human vasculature. The parameter set is especially relevant to simulate the flow in a patient-specific aneurysm wherein there is both a wide range of shear rates as well as pulsatile flow conditions. With this objective in mind, we have simulated the flow in a straight tube and a bent tube and shown the superior convergence properties of the model. By way of comparison, the model in Anand and Rajagopal (2004) displays much poorer convergence characteristics in the straight tube, and diverges in the bent tube. Further, the simulations in a straight tube show that the viscoelastic properties of the proposed model dominate at low mean velocities but are negligible at high mean velocities: this is further corroboration of experimental data in blood rheology. We opine that researchers intending to use a model in extensive numerical simulations of blood flow will find the proposed model more amenable to study than the model in Anand and Rajagopal (2004): this holds particularly true for very low Reynolds number flows, and start-up flows. We will take up the features of the proposed model’s flow in a subsequent paper. Acknowledgement MA thanks Dr. C. S. Sastry for useful discussions. JK was partially funded by a grant from the NCSA, and this support is acknowledged. We thank Professor K. R. Rajagopal for discussions on the presentation of the results. References Anand, M., & Rajagopal, K. R. (2004). A shear-thinning viscoelastic fluid model for describing the flow of blood. International Journal of Cardiovascular Medicine and Science, 4(2), 59–68. Bodnar, T., Sequeira, A., & Pirkl, L. (2009). Numerical simulations of blood flow in a stenosed vessel under different flow rates using a generalized Oldroyd-B model. In Numerical Analysis and Applied Mathematics (vol. 2, pp. 645–648) Melville, New York. Bodnar, T., Rajagopal, K. R., & Sequeira, A. (2011). Simulation of the three-dimensional flow of blood using a shear-thinning viscoelastic fluid model. The Mathematical Modelling of Natural Phenomena, 6(5), 1–24. Bodnar, T., Sequeira, A., & Prosi, M. (2011). On the shear-thinning and viscoelastic effects of blood flow under various flow rates. Applied Mathematics and Computation, 217(11), 5055–5067. Brooks, A. N., & Hughes, T. J. R. (1982). Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 32, 199–259. Bureau, M., Healy, J., Bourgoin, D., & Joly, M. (1978). Etude expérimentale in vitro du comportment rhéologique du sang et régime transitoire à faible vitesse de cisaillement. Rheologica Acta, 17(6), 612–625. Bureau, M., Healy, J., Bourgoin, D., & Joly, M. (1980). Rheological hysteresis of blood at low shear rate. Biorheology, 17(1–2), 191–203. Chien, S., Usami, S., Taylor, H. M., Lundberg, J. L., & Gregersen, M. I. (1966). Effect of hematocrit and plasma proteins on human blood rheology at low shear rates. Journal of Applied Physiology, 21(1), 81–87. Copley, A. L., & King, R. G. (1976). On the viscoelasticity of anticoagulated whole human blood in steady shear as tested by rheogoniometric measurements of normal forces. Biorheology, 12, 5–10. Coronado, O. M., Arora, D., Behr, M., & Pasuali, M. (2006). Four-field Galerkin/least-squares formulation for viscoelastic fluids. Journal of Non-Newtonian Fluid Mechanics, 140, 132–144. Guenette, R., & Fortin, M. (1995). A new mixed finite element method for computing viscoelastic flows. Journal of Non-Newtonian Fluid Mechanics, 60, 27–52. Hughes, T. J. R. (1995). Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127, 387–401. Hughes, T. J. R., Feijoo, G. R., Mazzei, L., & Quincy, J. B. (1998). The variational multiscale method- a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 166, 3–24. Inzoli, F., Boschetti, F., Zappa, M., Longo, T., & Fumero, R. (1993). Biomechanical factors in abdominal aortic aneurysm rupture. European Journal of Vascular and Endovascular Surgery, 7(6), 667–674. Kwack, J., & Masud, A. (2010). A three-field formulation for incompressible viscoelastic fluids. International Journal of Engineering Science, 48, 1413–1432. Masud, A., & Calderer, R. (2009). A variational multiscale stabilized formulation for the incompressible Navier–Stokes equations. Computational Mechanics, 44, 145–160. Masud, A., & Franca, L. P. (2008). A hierarchical multiscale framework for problems with multiscale source terms. Computer Methods in Applied Mechanics and Engineering, 197, 2692–2700. Masud, A., & Khurram, R. (2006). A multiscale finite element method for the incompressible Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 195, 1750–1777. McMillan, D. E., Utterback, N. G., Nasrinasrabadi, M., & Lee, M. M. (1986). An instrument to evaluate the time dependent flow properties of blood at moderate shear rates. Biorheology, 23, 63–74. Owens, R. G. (2006). A new microstructure-based constitutive model for human blood. Journal of Non-Newtonian Fluid Mechanics, 140(1–3), 57–70. Raghavan, M. L., Vorp, D. A., Federle, M. A., Makaroun, M. S., & Webster, M. W. (2000). Wall stress distribution on three dimensionally reconstructed models of human abdominal aortic aneurysm. Journal of Vascular and Endovascular Surgery, 31(4), 760–769. Rajagopalan, D., Armstrong, R. C., & Brown, R. A. (1990). Finite element methods for calculation of steady, viscoelastic flow using constitutive equations with a Newtonian viscosity. Journal of Non-Newtonian Fluid Mechanics, 36, 159–192. Rajagopal, K. R., & Srinivasa, A. R. (2000). A thermodynamic framework for rate-type fluids. Journal of Non-Newtonian Fluid Mechanics, 88(3), 207–227.
88
M. Anand et al. / International Journal of Engineering Science 72 (2013) 78–88
Srinivasa, A. R. (2000). Flow characteristics of a multiconfigurational shear-thinning viscoelastic fluid with reference to the orthogonal rheometer. Theoretical and Computational Fluid Dynamics, 13(5), 305–317. Sun, N., & DeKee, D. (2001). Simple shear, hysteresis, and yield stress in biofluids. The Canadian Journal of Chemical Engineering, 79, 36–41. Thurston, G. B. (1975). Elastic effects in pulsatile blood flow. Microvascular Research, 9, 145–157. Yeleswarapu, K. K. (1996). Evaluation of Continuum Models for characterizing the constitutive behavior of blood, Ph.D. Thesis, University of Pittsburgh, Pittsburgh, PA. Zarins, C. K., Giddens, D. P., Bharadvaj, B. K., Sottiurai, V. S., Mabon, R. F., & Glagov, S. (1983). Carotid bifurcation atherosclerosis. Quantitative correlation of plaque localization with flow velocity profiles and wall shear stress. Circulation Research, 53(4), 502–514.