Multiphase nonNewtonian viscosity of blood
Dear users,
I am working on multiphase blood simulations in patientspecific aortas using ansys FLUENT (version 19.1). I have succesfully implemented a twophase model on different geometries, but am struggling with the interpretation of some results.
I model the blood as being nonNewtonian using a modified CarreauYasuda model with two phases: plasma + RBC (as a granular phase). The UDF uses the shear rate of the RBCs and RBC volume fraction to calculate the RBC viscosity and this is then implemented in FLUENT under the specific settings for the RBCs in the granular phase. The bulk viscosity is assumed to be zero. The mixture viscosity is then (when postprocessing) calculated according to: mu_mix = mu_rbc*VOF_rbc + mu_plasma*VOF_plasma. The plasma viscosity is assumed to be constant.
To compare and verify my results, I did a single phase nonNewtonian simulation with the same UDF, but using a fixed VOF of 0.45 (therefore, the viscosity is only dependent on the strain rate magnitude.
When I compare the mixture viscosities for the single and multiphase model, I encounter lower viscosity magnitude values for the multiphase model. I have trouble finding a reason why this is the case, since (for as far as I know), the strain rate magnitude is only dependent on the velocity of the particles. There are no differences in the velocity contour plots, so I would expect the viscosity contours to be the same as well (at least during steady state). I listed the differences between the simulations below:
1. Single phase has overall plasma viscosity of 0.0037 Pa s, whereas the twophase has a plasma viscosity of 0.001 Pa s and the RBC viscosity is determined by the UDF.
2. The single phase model has an overall density of 1080 kg/m3 and the two phase model uses 1025 kg/m3 for the plasma + 1125 kg/m3 for the RBCs.
3. Obviously, the twophase model has the granular eulereuler multiphase approach enabled. Here, RBCs are assumed to be the granular phase with a diameter of 8e6 m. Granular bulk viscosity is assumed to be zero, only drag is enabled (Gidaspow model).
All other convergence settings are similar, and the results presented below are from converged steady state solutions.
Also, the UDF for 1) the single phase nonNewtonian model:
DEFINE_PROPERTY(cy_viscosity_singlephase, c, t)
{
float mu_plasma = 0.001;
float eps_RBC = 0.45;
float lambda = 0.110;
float n = (0.8092 * pow(eps_RBC,3))  (0.8246 * pow(eps_RBC,2))  (0.3503 * eps_RBC) + 1;
float m = (122.28 * pow(eps_RBC,3))  (51.213 *pow(eps_RBC,3)) + (16.305 * eps_RBC) + 1;
float exp = (n1)/2;
double gamma;
double mu;
gamma = C_STRAIN_RATE_MAG(c, t);
mu = mu_plasma*m*pow((1+ pow((lambda*gamma),2)),exp);
return mu;
}
And the UDF for the twophase model:
#include "udf.h"
DEFINE_PROPERTY(granular_viscosity_rbc , cell , thread_rbc)
{
int phase_domain_index , ID=9;
Thread *mixthread;
Thread *subthread;
Domain *mixture_domain;
mixture_domain = Get_Domain(1);
mixthread = Lookup_Thread(mixture_domain ,ID);
/* predefine variables */
real muplasma=0.001;
real murbc;
real mumix;
real eps_rbc;
real lambda=0.110;
real m;
real n;
real sr;
Domain *plasma_domain;
Domain *rbc_domain;
/* get domain numbers of plasma and rbc */
plasma_domain = Get_Domain(2);
rbc_domain = Get_Domain(3);
/* loop over all threads */
sub_thread_loop(subthread , mixthread , phase_domain_index)
{
/* viscosity in plasma thread is muplasma */
if ( subthread == Lookup_Thread(plasma_domain ,ID) )
{
muplasma = C_MU_L(cell , subthread);
}
else if ( subthread == Lookup_Thread(rbc_domain ,ID) )
{
sr = C_STRAIN_RATE_MAG(cell ,subthread); /* shear rate from rbc, could be the
wrong one */
eps_rbc = C_VOF(cell , subthread);
}
}
n = 0.8092*pow(eps_rbc ,3.)  0.8246*pow(eps_rbc ,2.)  0.3503*eps_rbc + 1.;
m = 122.28*pow(eps_rbc ,3.)  51.213*pow(eps_rbc ,2.) + 16.305*eps_rbc + 1.;
mumix = m* pow( (1.+ pow((lambda*sr) ,2.) ),((n1.)/2.) );
murbc = (mumix*muplasma (1.eps_rbc)*muplasma)/eps_rbc;
/*if (cell %50000==0)
{
Message("muplasma = %f\n", muplasma);
Message("strainrate = %f\n",sr);
Message("vof = %f\n",eps_rbc);
Message("n = %f\n",n);
Message("m = %f\n",m);
Message("mumix %f\n",mumix);
Message("murbc %f\n",murbc);
Message(" %f\n",);
}*/
return murbc;
}
Furthermore, I investigated the effect of taking the strain rate magnitude for the two phase model from the plasma phase / from the plasma + rbc phase. However, this gave no different viscosity magntide if deducted from the plasma phase, and if it was deducted from the plasma + rbc phase, the mixture viscosity was even lower (due to the shear thinning effect).
Therefore, I was wondering: what exactly are all the parameters effecting the shear rate magnitude (and this the mixture viscosity) in a multiphase model? Since everything I read / research leads to the conclusion that it's only a function of the velocity field, whereas those are similar for the different simulations.
I included images of the different Velocity Fields + RBC viscosities + mixture viscosities.
Thanks in advance!
Best Answer

Rob UKForum Coordinator Posts: 8,371
The granular phase will have some effect on the strain rate due to the drag etc of the solids, However I'm not sure how you've fixed the volume fraction in the single phase model?
You're correct that blood is a slurry, however, the nonNewtonian model is designed to account for this. Double check what plasma viscosity is as I suspect it's Newtonian. Where did you get the CY properties?
Answers
Hello,
I am not an expert in rheology and might not be able to answer all your questions. I just want to bring your attention into two points.
1) when we estimate the viscosity of a mixture, we use mass fraction instead of volume fraction. Not sure if this is applicable in your case or not. This will not make huge difference as your densities are close.
2) strain rate will have different values at different locations. It is function of the velocity gradient, not just velocity. For the bulk fluid, it might be the same but near the walls and interfaces, if you have any, the values are different. Check the velocities at the wall, you might find an answer to your question.
Hi!
Thanks for the quick response :)
Regarding the velocity gradient: do you mean that the velocity gradient can also be different at the plasmaRBC interface? In that case, that might be where strain rate differences occur (i think)?
if you are using VOF model, no it will not as it has only velocity of the mixture. I mean the cells near the wall (boundary layer) will have higher velocity gradient than the bulk fluid. plot the velocity vectors and zoom on the wall. You might find the answer as the velocity profile at the wall should be different.
Okay thanks :) I'm using EulerEuler btw.
Good, I was into recommending mixture or EulerEuler and solve for interfacial area.
Nice. What do you mean to solve for interfacial area (sorry to ask so much questions)?
No problem. This is another conservation equation that you can enable in the eulereuler model. you will find and option to select it when you choose Eulerian model. If your secondary phase does not keep constant shape, enable this option
The granular phase will have some effect on the strain rate due to the drag etc of the solids, However I'm not sure how you've fixed the volume fraction in the single phase model?
You're correct that blood is a slurry, however, the nonNewtonian model is designed to account for this. Double check what plasma viscosity is as I suspect it's Newtonian. Where did you get the CY properties?
Dear Rob,
Thank you for your answer. For the single phase model, I used one value in my UDF for the strain rate, which is epsilon_RBC = 0.45 (meaning it's just a fixed parameter).
Yes the plasma viscosity is Newtonian, which is why I model the RBC viscosity to be nonNewtonian. Then the total viscosity is calculated using: mu_mix = mu_rbc*VOF_rbc + mu_plasma*VOF_plasma.
The CY properties are deducted from the following paper:
Jonghwun Jung, Robert W. Lyczkowski, Chandrakant B. Panchal, and Ahmed Hassanein. Mul
tiphase hemodynamic simulation of pulsatile flow in a coronary artery.
Journal of Biomechanics, 39(11):2064 – 2073, 2006. ISSN 00219290. doi: https://doi.org/10.1016/j.jbiomech.2005.06.023
URL http://www.sciencedirect.com/science/article/pii/S0021929005002769
Wow I figured out why the single phase nonnewtonian has much higher values since there's a typo in the udf.. this explains a lot, nevertheless, your answers have been very useful in the thinking process so thanks a lot!
@michelletramper
Thanks, that's a more recent paper than
Gijsen F.J.H., van de Vosse F.N., Janssen J.D.; 1999; The Influence of the nonNewtonian properties of blood on the flow in large arteries: steady flow in a carotid bifurcation model; Jou of Biomechanics, Vol 32, pp601
Which I used to model this sort of application a "few" years ago. "The Heartbeat of Pulmonary Modelling" publication date Mar 1, 2003 publication description FLUENT Inc (now ANSYS Inc) publication description Fluent NEWS article showing the transfer of MRI data via stl into a CFD model of the pulmonary artery. Fluent NEWS, spring 2003, pp1617
Okay thanks! I will look into your paper then as well, might pick up some useful things :)
The work flow will still be valid, but the preprocessing, solver & hard ware have moved on a fair bit. Also remember that back then 1M cells was a big model and I did a lot of manual work on the mesh. The youth of today don't appreciate all of that! 😉