Power-law viscosity model for blood
 

Hello friends,

Hope everything is ok.

I'm opening a new topic despite knowing that there are some other very similar. It happens that none clarifies me as I need...

So, I'm using solver simpleFoam to simulate blood in a 2D geometry as a newtonian and non-newtonian fluid and the only model that is giving me problems is power -law.

The power-law viscosity model has two parameters and, for blood, I found in the paper of Neofytou & Tsangaris (2006) the following values:

• k = 14.67×10−3 Pa.s^n
• n = 0.7755

I know that OpenFOAM uses kinematic viscosity and when I analysed the code of power-law model, noticed that the consistency index, k, as the same units of nu. So, I simply divided the value of k by rho (1040 kg/m3), and got a value of 1.411×10−5 m2/s for k.

In a first phase, I used this value in my simulation and it diverged rapidly.

In a second phase I have initialized all fields with potentialFoam and my simulation still doesn't converge.

I noticed, in some tutorials, that k has values in the order of 2500 and I have some questions:

1. Is anything escaping me? Am I supposed to obtain a value of k of this magnitude for blood with the data that I have available?
2. One thing I do not understand is the n exponent in the units of k. Can somebody explain me why is this exponent there? What does it mean?

I really need your help to proceed with my work. Sorry for my newbie question!

Cheers,

FS

People, can someone give me a help on this question, please?

It's really important to get my doubts solved!

FS

Greetings FS,

I quickly went to see your past posts and it seems I told you some months ago how important it is to increase the complexity in your test cases gradually, instead of jumping directly at the final problem ;)

Personally I never use the Power-law viscosity model, therefore I'm not familiar with it. If you can create and/or provide a simple example case, preferably with the expect analytical solution, for example something like the very first "cavity" tutorial explained in the User Guide, it should become pretty clear as to what is what ;).

Best regards,
Bruno

Hi Bruno!

First of all thank you for your answer! I remember well your advices! ;)

However, I think now the problem is slightly different and I think you could help me...

I'm using a 2D T-shaped geometry (2 inlets and 1 outlet) and studying, with simpleFoam solver:

1) 3 levels of mesh refinement;
2) different schemes for divergence and pressure gradient;
3) Newtonian, Power-law and Bird-Carreau viscosity models.


In a first phase, I studied the coarsest mesh: my simulations successfully converged when applying both type of schemes and the three viscosity models.

For the intermediate and refined mesh, I only can attain convergence when applying newtonian viscosity model. I tried to manipulate relaxation factors for Power-law and Bird-Carreau but without success... Residuals for Ux and P stabilize far away from 1e-05!
I'm using potentialFoam to initialize the fields in all my non-Newtonian simulations...

Attached, there are several files for the intermediate refined mesh using a Power-law viscosity model, as well as a checkMesh file and my residuals evolution...

Can you see them and give me some hints, please?

Do you think that leaving simpleFoam and use nonNewtonianIcoFoam could be a solution?

FS

P.S.: as I said in my first message, I don't understand very well why there's an exponent "n" in the units of the consistency index "k"... In http://en.wikipedia.org/wiki/Power-law_fluid itit is possible to see that "K is the flow consistency index (SI units Pa•sn)"... Also, I don't know very well how this influences the value to input in powerLawCoeffs (transportProperties file)...


Attachment 33997

Attachment 34002

Attachment 34004

Attachment 34005

Attachment 34006

Hi.

the things about k, nu and n. When I have a look into /src/transportModels/incompressible/viscosityModels/powerLaw/powerLaw.C I see this:

// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::powerLaw::calcNu() const
{
return max
(
nuMin_,
min
(
nuMax_,
k_*pow
(
max
(
dimensionedScalar("one", dimTime, 1.0)*strainRate(),
dimensionedScalar("VSMALL", dimless, VSMALL)
),
n_.value() - scalar(1.0)
)
)
);
}


k_: flow consistency index
n_: flow behaviour index
strain-rate: strainRate()
nu_: kin. viscosity

k_ has the dimension of nu_. , because the strainRate() [1/sec] is multyplied by 1sec. ( "dimensionedScalar("one", dimTime, 1.0)"). The equation is like this:

nu_ = k*{strain-Rate()*sec.}^(n-1)

The flow behaviour index n is dimensionless. So everything is in the proper dimension.

Is this the answer of your question??

kind regards
Florian

Hi,

-----------fvSchemes---------

Why you use this:

laplacianSchemes
{
default Gauss linear orthogonal;
}

snGradSchemes
{
default orthogonal;
}

instead of:

laplacianSchemes
{
default Gauss linear corrected;
}

snGradSchemes
{
default corrected;
}

????


Please post the complete case. Then we can see

Kind regards
Florian

Florian, thank you very much for your answer!!

Regarding the units of "k", I understood your explanation! Thank you very much! :)

So, I'm using "orthogonal" and not "corrected" because my mesh is orthogonal. I saw (I don't remember where) that in these cases we can use this definition... Do you suggest me to change this?

By the complete case, you mean my mesh and all the other folders ("0", "constant" and "system")?

Cheers,

FS

Hi,

Ah I see in your checkMesh-file:

Mesh non-orthogonality Max: 0 average: 0

o.k. then everything is fine with that.

yes with complete case I mean "0", "constant","system". Mesh is in PolyMesh. I don't need the blockMesh file.

I will give it a shot at my PC. Then I will see what's going wrong with this case.

kind regards
Florian

Hi 'fs.chemech',

A comment on your original problem. I also faced a similar issue and could not figure out why the solution was diverging for non-newtonian law. Initalise your flow field with a newtonian model and then switch non-newtonian and it wont diverge. This solved the problem for me. Hope it helps!

Florian,

my compressed case (zip) is much bigger (~800 KB) than the max size allowed by the Forum. How can I send it to you?




Nadish,

thank you for your suggestion. You change to non-Newtonian model while Newtonian simulation is running?

In the controlDict file you may set the option rumTimeModifiable to yes and then switch during runtime. But i guess the better option is to run the computation for say 0.1 sec, halt it, change the power law model and start from 0.1 sec. it will give a good initial value to the solution and i am pretty sure it wont diverge (given you have defined appropriate BC's)

Nadish,

I tried your suggestion... I know that with Newtonian model, my simulation converges in 290 iterations...

So, I ran it until the 260ª iteration, stoped it, changed viscosity model to powerLaw and adapted startTime in controlDict file and resumed simulation again...

Unfortunately, residuals Ux and p began to rise and stabilized again near 1e-04...

Do you have any idea why this happens?

Thank you for your time! ;)

Hi,

You just need to initialize the values using Newtonian model. I don't think running upto 260 iterations is wise. Just run a few... Say 20... Iterations with Newtonian and then switch

Nadish,

I've done it... Residuals continue to stabilize around 1e-04...

Hi,

do you have any experimental results for checking accuracy??? Pherhaps your results are good in case of experimental data.

Residuals are not good for checking accuracy. For example there can be some instationarities due to "fluctuations". (Keep in mind, laminar non-Newtonian flow is very similar to turbulent flow). Pherhaps you have transient effects in your simulation.

One hint:
If you want to have good results in case of residuals, then you can do a pseudo time-stepping. The only thing you have to do is to simulate your case unsteady (fvSchemes-ddtSchemes: Euler or backward). I would use backward for pseudo steady cases. You solve your case unsteady until there is no change in your residuals. Some people call this pseudo-transient continuation. This could help you a lot. Furthermore, with pseudo time-stepping, you do not need to init with newtonian flow.

Please let me know if this will work. In a lot of my cases this works really good. :)

kind regards
Florian

Hi again,

if the fluctuations in your residuals will not decay, then you have transient effects in your flow. In this case your problem is unsteady and not steady. To get good avarage data, a steady simulation is not the best way. A better way is to track your field data over time and make an time-avarage (similar to LES or DNS). This can be very tricky, because the data are correlated over time (Autocorrelation). You need several integral time scale, until your results are good due to statistically accuracy.

kind regards
Florian

Glad to have helped out :)

Hi Florian,

I don't have any experimental results to be based on...

Yes, one of my suspicions is also that I have transient effects affecting my results...
To solve my case as an unsteady case I have to change my solver to nonNewtonianIcoFoam isn't it? simpleFoam, doesn't run unsteady simulations, I think...

You said that residuals are not good for checking accuracy! What, in your opinion, should I be based on?

Cheers,

FS

Florian, can you explain me this a little better?
Sorry for my newbie questions...

Hi Nadish,

I wanted to say that I performed your suggestion. However, the simulation still does not converge ... Residuals continues to stabilize well above the desirable, around 1e-04!

Thank you very much for your help! ;)