[alexey7783]

Hi forum members

Although I already did some simulations with OpenFOAM I would call me a real newbie, especially as my programming skills are still under development. So I hope that someone of you can help me:

for my PhD thesis (which is experimental) my professor wants me to implement some numerical simulations. It comprises the outflow of a non-newtonian fluid from a cylindro-conical vessel. I did some first experiments with newtonian fluids (water and a water/glycerol mixture at 22°C) and simulated them with OpenFOAM. The results are really really promising, as for example the outflow times are nearly the same (only 2-4% error) or the massflux at the orifice (1-3% error). You see, I was really happy at that point

But now the point came, where I had to do the "real" work, which means try all of that with a non-newtonian fluid. I chose a very simple one, namely the glycerol/water mixture mentioned above to which I added some carboxymethylcellulose (CMC). I recorded the flow curves with our rheometer in the range from shear stresses starting from approx. 0 to 500 1/s, because these are the shear rates I extracted from the newtonian simulations. Than I chose the Cross-Power-Law because the curve fits showed excellent agreement and (which is the more important reason) is already implemented in the interFoam solver (as the fluid is flowing out, air is coming in from above, that's why I need a twophase solver).

When I simulate the outflow with the parameters I computed for the Cross-Power-Law the fluid is flowing far to quick. For example, in reality it takes the different mixtures between 37-56 seconds, in the corresponding simulations between 15-20 seconds. That is really a huge deviation. And of course, the mass fluxes at the orifice are therefore also too high. Do you think this is the fault of the solver? I tried to do the following first: I changed the Cross-Power-Law-Coefficients in that way, that newtonian flow with the viscosity of the pure glycerol/water mixture was present (m=0, n and nuInf can be whatever) and Lo and behold! the same results as in the pure newtonian case. This means for me, that the solver is working right, but I had an idea:

I than tried to analyse the structure of the solver with my basic c++ knowledge and I understand it the following way: the solver is calculating the viscosity value of each cell by using the velocities that are currently present. Than he stores it as viscosity value of this cell and than he basically solves the NSE acting like the viscosity of this cell is constant (newtonian behaviour so to say) in this time step (the one calculated at the beginning). But isn't it more correct to calculate the viscosity iteratively, because when the velocity changes, the calculated nu should also?
I say so, because this for me could explain the faster outflow.

My question to you with more experience is: do you think it is worth spending so much time (for me it will take a lot of time, due to my lack of programming knowledge till now) on implementing the calcNu() in the iteration, or are the wrong results due to something else, for example the wrong shear rate range in which I determined my fit parameters for the Cross-Power-Law model?

I would really appreciate your help and I am looking forward to your answers

My best regards

Alex

[alexey7783]

Hi it's me again,

I searched the forum once again (but spending some more time) and I saw, that there were already some discussions on the calculation of the strainrate(). Perhaps that is the solution to the problem?!? What do you think?

http://www.cfd-online.com/Forums/ope...1-vs-15-a.html

and

http://www.cfd-online.com/Forums/ope...ate-error.html

Best regards again...

Alex

[ternik]

Originally Posted by alexey7783

Hi it's me again,

I searched the forum once again (but spending some more time) and I saw, that there were already some discussions on the calculation of the strainrate(). Perhaps that is the solution to the problem?!? What do you think?

http://www.cfd-online.com/Forums/ope...1-vs-15-a.html

and

http://www.cfd-online.com/Forums/ope...ate-error.html

Best regards again...

Alex

Hi Alex,
I think that correct implementation of shear rate (IInd invariant of symmetrical rate of deformation tensor) should be:

return sqrt(2.0)*mag(symm(fvc::grad(U_)));

I have done it like this and tested it - results are in agreement with other authors for "standard test (benchmark) cases" (e.g. lid-driven cavity flow)!

I hope this will help you.

Cheers,
Primoz.

[alexey7783]

Originally Posted by ternik

Hi Alex,
I think that correct implementation of shear rate (IInd invariant of symmetrical rate of deformation tensor) should be:

return sqrt(2.0)*mag(symm(fvc::grad(U_)));

I have done it like this and tested it - results are in agreement with other authors for "standard test (benchmark) cases" (e.g. lid-driven cavity flow)!

I hope this will help you.

Cheers,
Primoz.

Hi Primoz,

thank you for your reply. I read a lot yesterday and the more I read the more doubts I get. As Kerstin wrote in

http://www.cfd-online.com/Forums/ope...1-vs-15-a.html

it is always difficult to use a 3D definition of shear rate and apply it on a 1D Model like Cross-power-law whose parameter were obtained in a rheometer. Especially as a 1D rheometerflow doesn't comprise elongation parts but the 3D definition of shear rate does. One has to bear in mind, that elongation viscosity for newtonian fluids is three times higher than shear viscosity. Or do I understand something wrong?

Yesterday I tried the following: I computed the gradient tensor field L of the velocity field L=fvc::grad(U) and than I had a look at the entries L11, L22 and L33 which correspond to dU/dx, dV/dy and dW/dz. I thought for finite volumes the sum of these quantities has to be zero?!? Because of the continuity equation for incompressible fluids. But for all of my simulations (regardless which) this doesn't apply. Why is this so?

My newest idea is to define a new shear rate, which only uses shear velocities and excludes elongation parts.

Any tipps how to do that? I thought of something similar to the current definition, except that I substract all elongation parts. Why is it a "good" idea? As my CMC solution is flowing out to quickly the shear rates have to be too high (high shear rate = low calculated viscosity with Cross-power-law as CMC is shear thinning). When I exclude the elongation parts (which are quite high in the conical part) I will automaticly reduce my shear rate. But this sounds a bit like a "dirty approach" to me.

What is your opinion? Has somebody some experience with non-newtonian simulations and perhaps even validated the results?

Thanks for your answers and time

Best regards

Alex

[Locky3827]

Hi Alex,
My interest is non-Newtonian laminar flows in complex geometries. As a validation I calculated pipe flow of a powerlaw fluid for which there is an analytical solution. I used a fluid with shear stress=0.5*(shear rate)^0.5 at a superficial velocity of 0.2 m/s in a 100mm diameter pipe. The analytical pressure gradient is 0.0894 kPa/m. With OpenFOAM 1.6.x "out of the box" I got 0.105 kPa/m. With the shear rate altered as in the post by ternik above I got the correct value (0.09 kPa/m reading off the graph).
Hope this helps,
Lachlan

[alexey7783]

Hi Lachlan,

thank you very much for your reply. I already implemented this definition and it also worked for me.

Thanks for caring and sharing

Alex

[ternik]

Originally Posted by Locky3827

Hi Alex,
My interest is non-Newtonian laminar flows in complex geometries. As a validation I calculated pipe flow of a powerlaw fluid for which there is an analytical solution. I used a fluid with shear stress=0.5*(shear rate)^0.5 at a superficial velocity of 0.2 m/s in a 100mm diameter pipe. The analytical pressure gradient is 0.0894 kPa/m. With OpenFOAM 1.6.x "out of the box" I got 0.105 kPa/m. With the shear rate altered as in the post by ternik above I got the correct value (0.09 kPa/m reading off the graph).
Hope this helps,
Lachlan

Hi Lachlan,

nice to see that my post regarding correct shear rate (or better said 2nd invariant of symmetrical rate-of-deformation tensor) did help you!

Enjoy,
Primoz

[Locky3827]

Yes your post was most useful. I'm intending to try out the other rheological models when I get a chance as well. I've already tried the Herschel Bulkley model and that was OK also after the 2nd invariant of symmetrical rate-of-deformation tensor was corrected.

[chrisb]

Originally Posted by ternik

Hi Alex,
I think that correct implementation of shear rate (IInd invariant of symmetrical rate of deformation tensor) should be:

return sqrt(2.0)*mag(symm(fvc::grad(U_)));

I have done it like this and tested it - results are in agreement with other authors for "standard test (benchmark) cases" (e.g. lid-driven cavity flow)!

I hope this will help you.

Cheers,
Primoz.

I agree with this definition. A good resource for those wanting to understand why is:
@book{fung1993biomechanics,
title={{Biomechanics: mechanical properties of living tissues}},
author={Fung, Y.C.},
year={1993},
publisher={Springer}
}
pp: 73-74

Using this You arrive with the definition in "OpenFOAM" as:
return sqrt(2.0)*mag(symm(fvc::grad(U_)));

Note: This is only different to the implemented version by a factor of sqrt(2.0) therefore it may be accounted for easily in a model.
If you are using non-Newtonian models I would strongly suggest that you check the source code to determine what it is solving for.

Cheers,
CB

[Locky3827]

Incidentally OpenFOAM 1.7.0 still has the original incorrect shear rate definition.
Regards,
Lachlan

[mjfarnley]

Does anyone know whether this error has been corrected in OpenFoam, please? I have just downloaded v 2.1.

Thanks,

MJ.

[chrisb]

Yes it has been corrected MJ. See the following file

$FOAM_SRC/transportModels/incompressible/viscosityModels/viscosityModel/viscosityModel.C

Chris

[mjfarnley]

Hi,

Thank you very much for such a quick reply.

Best regards,

MJ.