NonNewtonian Flow: is it implemented right?
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 nonnewtonian fluid from a cylindroconical 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 24% error) or the massflux at the orifice (13% 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 nonnewtonian 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 CrossPowerLaw 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 CrossPowerLaw the fluid is flowing far to quick. For example, in reality it takes the different mixtures between 3756 seconds, in the corresponding simulations between 1520 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 CrossPowerLawCoefficients 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 CrossPowerLaw model? I would really appreciate your help and I am looking forward to your answers My best regards Alex 
Addition
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.cfdonline.com/Forums/ope...1vs15a.html and http://www.cfdonline.com/Forums/ope...ateerror.html Best regards again... Alex 
Quote:
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. liddriven cavity flow)! I hope this will help you. Cheers, Primoz. 
Quote:
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.cfdonline.com/Forums/ope...1vs15a.html it is always difficult to use a 3D definition of shear rate and apply it on a 1D Model like Crosspowerlaw 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 Crosspowerlaw 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 nonnewtonian simulations and perhaps even validated the results? Thanks for your answers and time Best regards Alex 
Hi Alex,
My interest is nonNewtonian 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,
thank you very much for your reply. I already implemented this definition and it also worked for me. Thanks for caring and sharing ;) Alex 
Quote:
nice to see that my post regarding correct shear rate (or better said 2nd invariant of symmetrical rateofdeformation tensor) did help you! Enjoy, Primoz 
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 rateofdeformation tensor was corrected.

Quote:
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: 7374 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 nonNewtonian models I would strongly suggest that you check the source code to determine what it is solving for. Cheers, CB 
Incidentally OpenFOAM 1.7.0 still has the original incorrect shear rate definition.
Regards, Lachlan 
Does anyone know whether this error has been corrected in OpenFoam, please? I have just downloaded v 2.1.
Thanks, MJ. 
Yes it has been corrected MJ. See the following file
$FOAM_SRC/transportModels/incompressible/viscosityModels/viscosityModel/viscosityModel.C Chris 
Hi,
Thank you very much for such a quick reply. Best regards, MJ. 
All times are GMT 4. The time now is 15:27. 