Viscoelastic Fluid Flows using OpenFOAM The solver viscoelasticFluidFoam
 

Hello Foamers,

My name is Jovani and I want announce that will soon be available a solver for treatment of viscoelastic fluid in OpenFOAM. I was working in the development of a viscoelasticFluidFoam solver (with the great orientation of the Dr. Hrvoje Jasak, thank you very much Hrvoje!!) to simulate viscoelastic flows.

The follow models are implemented: Maxwell, Oldroyd-B, Giesekus, PTT linear and exponential,
FENE-P, FENE-CR and DCPP. There are others models already implemented but not tested (Pom-Pom, SXPP, DXPP, WM ...).

The metodology used to obtain solutions in high We numbers was the DEVSS. The test geometry was mainly the 4:1 planar contraction. The solver are able to simulate multimode cases with n modes.

For the interested: A presentation that will soon be available on line to give more informations.

Jovani

Brilliant work, Jovani! Do
 

Brilliant work, Jovani!

Do you think I could take a look at the sources? This is something we were looking to implement, but it looks like you've saved us the trouble!

Hello Menon, Yes, the solve
 

Hello Menon,

Yes, the solver will be available in the -dev version of OpenFOAM, but you need wait untill the solver is put into the source. Don't be worried, this will be done soon.

Jovani

Hi Jovani, Great work! This
 

Hi Jovani,

Great work! This will speed up research in viscoelastic simulations a lot. 8-)
Really great!
Greetings
Kerstin

A piece of propaganda: Jovani
 

A piece of propaganda: Jovani did this work under a Masters project at Universidade Federal Rio Grande do Sul, Brazil and we are trying to convince him to stay on and do a PhD along similar lines.

Can I ask the audience to kindly help me in this.

Hrv

Hi Jovani, great job! And t
 

Hi Jovani,

great job! And the fact that it was made under a master's project is simply terrific! Great.

And... You should just listen to Hrvoje! Go on with CFD and do a PhD. So we can wait for more good to come in the field of viscoelastic flow simulations! http://www.cfd-online.com/OpenFOAM_D...part/happy.gif

best
Holger

Thanks everyone!! Thanks Hr
 

Thanks everyone!!

Thanks Hrvoje for the propaganda!

Marschall, I'm thinking to a PhD. There are a lot of studies to be done in this area.
As there isn't much in this area into commercial softwares, I think that if we do any efforts, with the help of Hrvoje of course, we can create something top of line.

Jovani

Hello Jovani, Thank you for
 

Hello Jovani,

Thank you for your good work. When will the next solver be released?

Ning

Hi Jovani, that all looks v
 

Hi Jovani,

that all looks very promising.

However, could you give some details about the stability of your solver:

- up to which We (or De)-Number stable solutions can be obtained in your 4to1 planar geometry?

- any tests on realistic 3-D geometries and if, what is the performance in terms of CPU-time?

cheers Max

Hello All, Ning, the solver
 

Hello All,

Ning, the solver will released in this month or the maximum in April. I believe I will not have time to test some models (white-Metzner and derivations, PP, XPP and Feta-PTT) but their will be into the solver too. Is necessary validate these models.


Max,

About the De numbers: I make any test to explore this question, but only for some models (PTT and Giesekus). But to give to you a number I simulate a case that I simulate was a LDPE flow with De=250 using Giesekus 8-mode model, this was't the limit for this case, I simply don't tried higher values. A complete study to this maximum Deborah values need to be done. Other, the DEVSS methodology sugest an arbitrary value for the the term added to moment equation, generally is used the value of etaP (value that is used in the solver), but we know that the better value for this term is related with each model and your parameters, this need more study too and I think with it's possible obtain any improvements. The interpolation shemes used also afect it.

About 3-D simulations: I do a test with a realistic case that was the flow into a capillary with a contraction. The mesh was about 150000 volumes, considering 2 simetry planes. As is make the transient simulation the Courant number control the time step, but the simulation time is acceptable for me. I used multigrid and parallelization tecniques for it.

file:///home/jovani/Desktop/PresentationSolver/pictures/magtautri2.png

file:///home/jovani/Desktop/PresentationSolver/pictures/meshtri.png


Jovani

Hello All, Ning, the solver
 

Hello All,

Ning, the solver will released in this month or the maximum in April. I believe I will not have time to test some models (white-Metzner and derivations, PP, XPP and Feta-PTT) but their will be into the solver too. Is necessary validate these models.


Max,

About the De numbers: I make any test to explore this question, but only for some models (PTT and Giesekus). But to give to you a number I simulate a case that I simulate was a LDPE flow with De=250 using Giesekus 8-mode model, this was't the limit for this case, I simply don't tried higher values. A complete study to this maximum Deborah values need to be done. Other, the DEVSS methodology sugest an arbitrary value for the the term added to moment equation, generally is used the value of etaP (value that is used in the solver), but we know that the better value for this term is related with each model and your parameters, this need more study too and I think with it's possible obtain any improvements. The interpolation shemes used also afect it.

About 3-D simulations: I do a test with a realistic case that was the flow into a capillary with a contraction. The mesh was about 150000 volumes, considering 2 simetry planes. As is make the transient simulation the Courant number control the time step, but the simulation time is acceptable for me. I used multigrid and parallelization tecniques for it.

file:///home/jovani/Desktop/PresentationSolver/pictures/magtautri2.png

file:///home/jovani/Desktop/PresentationSolver/pictures/meshtri.png


Jovani

Ops, Sorry!!! An error occu
 

Ops, Sorry!!!

An error occurred, I sent twice the same message and the pictures do not is showed. I believe the pictures appears now.

U magnitude:


http://img24.imageshack.us/img24/174/magutri2.png][IMG]http://img24.imageshack.us/img24/174/magutri2.th.png[/IMG][/URL]


Stress magnitude:

http://img6.imageshack.us/img6/6859/magtautri2.png][IMG]http://img6.imageshack.us/img6/6859/magtautri2.th.png[/IMG][/URL]

Mesh:

http://img22.imageshack.us/img22/9122/meshtri.png][IMG]http://img22.imageshack.us/img22/9122/meshtri.th.png[/IMG][/URL]

Jovani

Hello Jovani

I am interested in knowing if the solver is released yet.

I also would like to know, before switching to OpenFOAM, if it possible to model two phase flows, one of the phase having specific viscoelastic properties.

Thank you

Mike

Hi Michael,
In principle nothing is impossible in OpenFOAM. ;-)
I guess handling the viscoelastic stresses correctly in the region, where you don't have a sharp interface is the most tricky part...I mean to make it really physical.
But how do you intend to do it, if you don't switch to OpenFOAM? Everything I can imagine to do as alternative (commercial code or own code) will be much more difficult than in OpenFOAM...
Greetings
Kerstin

Hello Mike,

Kerstin is rigth about OpenFOAM :)

The solver was not release yet! I am finishing the tests Hrvoje ask me to do and I need talk with him about it in this days.
About free surface viscoelastic flow: Yes is possible in OpenFOAM, I am simulate rod-climming and die swell, but this work are just beginning, few tests for the methodology was done untill now!!

I advert you when the solver is released.

Best,
Jovani

viscoplastic constitutive model
 

Hi,
Jovani,

I am planing to use OpenFOAM to simulate aluminum extrusion which is viscoplastic constitutive model. Will you give the viscoplastic constitutive model in the future? I think maybe for you it is a very easy work to do. But I just don't know how to do it and where I can get started to do it by myself. Could you please give me some advice?
Attatchment is the constitutive model. Thank you.

Wendy.

Hi Wendy,

What do you mean mathematically with D^d in the .doc?
D is a tensor and what is meant with the power of d???

At the moment I don't see any viscoelastic part in your equation. Viscoelastic properties are related to the left hand side (time derivative) of a constitutive equation. You should have a look on the theory of linear viscoelasticity, where material behaviour and equations are explained with spring's and dash-pots...(keyword: Maxwell fluids... ) you only have an explicit equation for stress tensor as a function of deformation gradient tensor with non constant viscosity...

but please explain the D^d further, bacause your material law looks interesting...
Kerstin

Hi, thank you.

It is a viscoplastic model , not viscoelastic.
D is rate of deformation tensor, the superscript d represents deviator.
the equation for D is as attathment.

Thank you very much.

Wendy

Hi, I have another constitutive model , maybe it is easier to realize. please see the attathment.
Thank you .
Have a nice day.

Wendy

No pressure... but Jovani has got his examination on Monday - just registering support in public :)

Hrv

Hi,
Oh sorry, I didn't read correctly... viscoplastic makes sense for this equation...
I am still confused. In pricipal your equations don't look very difficult because explicit.
But your .doc files are not self explaining... your epsilon with "overline and dot" is a scalar elongation rate? Calculated e.g. from trace of D tensor? or is this something else?

and in the first equation, still I don't get the meaning of "d"... "D" is clear.
\tau is a tensor... and D is a tensor...but what is deviator? I am really interested in understanding this approach. Looks for me in principal as a extension of newtonian stress tensor law with non constant viscosity. And viscostity is a function of temperature and shear rate...
Can you give a literature reference for the first equation?
Kerstin

Hi,

they are two different model. d means deviator, I don't know exactly what does it mean also. I will check later.
the first is coming from the thesis of Proefschirift: finite element simulation of the aluminum extrusion process.
the author uses it according to: L.A. Lalli and A. J. DeArdo. Experimental assessment of stucture and property prediction s during hot working.
and R.N.Wright , G.G.Lea, and F.F.Kraft. Constitutive equations and flow stress characterization concepts for aluminum extrusion.

the second is coming from the software Superforge.
thank you

Wendy

Hi,
I did the details for the first model. please see the attachment .
Please check if it is correct.
as for the second model I will do it right after.

Thank you.

Wendy

Hi,
ah ok. D^d is dyadic product... in OpenFOAM nomenclature it is D && D

But still your equation is then not very universal. D : D is a scalar quantity. So your stress tau is also scalar, because you multiply scalar viscosity with D^d...
If you like to include the new stress _tensor_ in the momentum balance, you need a tensor equation for your material law. tau= 2 eta(D,T) D^d is a scalar equation. I would guess you mean tensor tau = 2 eta(shearrate, temperature) D and D is also tensor...but in principle you have only explicit equations and you need to update the calculation of viscosity in each iteration step with actual values of shearrate and temperature.

define the shearrate e.g. like

volScalarField shearrate
(
IOobject
(
"shearrate",
runTime.timeName(),
mesh,
IOobject::NO_READ
),
pow(2.0/3.0* (Dtensor && Dtensor), 0.5)
);

define Z and eta similar and update calculation of Z and eta in the time loop....

Hi,
Thank you for your reply.
Do you mean D^d is equal to D:D? Then tau is a scalar?
Thanks.

Wendy

Sorry, I typed D : D

I guess you are right, I agree with that tau= 2 eta( D,T) D. tau is a tensor

is the gamma dot that you have written is shear rate.
you wrote gamma dot = sqrt( (2 /3) D : D)

I think if your gamma dot is shear rate than the above expression is wrong as shear rate should be :
shear rate = sqrt( 2 D : D ).

I mean this is how i used it for power law fluid in my solver.

I don't know if it is wrong or not, I copied it from a PhD thesis of B.J.E. van Rens.
Thank you.

Wendy

D^d = D : D is what _you_ typed in your word .doc ... I have never seen this definition for D^d before...

Sorry, I did not mean that.
I only has not finish D^d equation because I was not sure how to do it.
Now please see the attachment for the equation of D^d

Thank you very much.

Wendy

forgot to attach it. Sorry.

Wendy

Hello Andy,

This seems confused. What you need to have is the correct formulation for each term of your model (mainly for D^d)!! I think you get it where this model was introduced!!

The gamma dot can be defined in different ways, the more used is shearRate = sqrt( 2 D : D ), but see what is correct in the definition of your model.

I believe tau=2eta(D,T)D, where D is a tensor obviously (tau is a tensor!!!!) or D^d is the deviatoric part of D, but still is a tensor!!

Then after you have the correct formulation, follow icoFoam solver:

1) create the temperature T (volScalarField) like the pressure p in icoFoam;

Info<< "Reading field T\n" << endl;
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

2) Read all the properties for your model like for nu:

dimensionedScalar nu
(
transportProperties.lookup("nu")
);

3) Put your T equation into the time loop, before or after UEqn?? follow your algorithm, if before:

fvScalarMatrix TEqn
(
the equation
);


4) Update you viscosity, using your model like for Carreau-Yassuda, for example:

volScalarField etaEff = eta0 * Foam::pow( 1 + Foam::pow( K* sqrt(2.0)*mag(symm(fvc::grad(U))),a), (m - 1)/a );

here sqrt(2.0)*mag(symm(fvc::grad(U))) = sqrt( 2 D : D ), if is sqrt( 2/3 D : D ), no problem, do this modification!!


5) Use etaEff into your moment equation

6) Enjoy

Jovani

Hi, Jovani,

Thank you very much. I now understand that D^d is a deviative part of D. it is a tensor, tau is a tensor, yes. :o
I will study your code right now. I am really a beginner. I will ask you when I have questions about it.
Thank you very much, I really appreciate that. :)

Wendy

Hello Wendy,

Sorry for "Andy" in my last post!!

Another suggestions for you: You can obtain the deviatoric part of a tensor directally in the OpenFoam, for example, if you have a tensor T the function dev(T) give to you the deviatoric of this tensor, see page 21 and 24 of the programer's guide!

Have a good job,

Jovani

Hi, Jovani,

It is really a good suggestion for me.
And I know you are talking to me. so never mind. :)

Thank you very much.
Have a nice day.

Wendy

Hi, Jovani,
When you mentioned programmer's guide, do you mean the webpage? http://www.nabla.co.uk/
or is there another document?
I can not see page number on the programmer's guide
Thank you.

Wendy

Hi, Jovani,

I found the PDF document of programmer's guide.
Sorry for disturbing again.

Wendy

the way I changed the constitutive model to viscoplastic
 

Hi, Jovani,

Thank you for your help (and Kerstin also ) , I have done the modification based on icoFoam slover to solve viscoplastic problem.
the attachment is the steps I did.
would you have some time to see if it is right or not? thanks again.
Have a nice weekend.

Wendy

Sorry , I can not upload it.
so just paste here

1 ) set up viscoPlasticProperties file

FoamFile
{
version 2.0;
format ascii;
class dictionary;
object viscoPlasticProperties;
}
// modified based on viscoelasticProperties
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


//InnerIterations 1; //Solver Algorithm control


//nuS nuS [0 2 -1 0 0 0 0] 1e-5;


rho rho [1 -3 0 0 0 0 0] 2700;


//numberOfModes 3;
R [0 0 0 -1 0 0 0] 8.321 // J mol^-1 K^-1
A [0 0 -1 0 0 0 0]1.445e+12
n 4.750
beta [-1 1 -2 0 0 0 0] 3.224e-8
deltaH 1.772e+5



// PTTExp on; //PTTExp on, off;
// PTTLin off;//;on; //PTTLin on, off;4.750




// ************************************************** *********************** //


2) read the properties in createFieldPlastic.H

Info<< "Reading viscoPlasticProperties\n" << endl;
IOdictionary viscoPlasticProperties
(
IOobject
(
"viscoPlasticProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);

dimensionedScalar rho ( viscoPlasticProperties.lookup("rho") );

dimensionedScalar R ( viscoPlasticProperties.lookup("R") );

dimensionedScalar A ( viscoPlasticProperties.lookup("A") );

dimensionedScalar beta ( viscoPlasticProperties.lookup("beta") );


Scalar deltaH ( readScalar(viscoPlasticProperties.lookup("deltaH") )
);




scalar n (
readScalar(viscoPlasticProperties.lookup("n"))
);


//effective shear rate
volScalarField etaEff
(
IOobject
(
"etaEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE //why
),
sqrt(2.0/3)*mag(symm(fvc::grad(u))); //
);



volScalarField Z
(
IOobject
(
"Z",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
etaEff*exp(deltaH/(R*T)); // Zener-Hollomon parameter
);

// viscosity
volScalarField eta
(
IOobject
(
"eta",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
asinh(pow(Z/A,1.0/n))/(3*beta*etaEff); // viscosity
);


//dynamic viscosity
volScalarField eta_d
(
IOobject
(
"eta_d",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
eta/rho; //dynamic viscosity
);



3) my_icoFoam_viscoPlastic.C




Info<< "\nStarting time loop\n" << endl;


for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;


# include "readPISOControls.H"
# include "CourantNo.H"


// here add viscosity calculation
etaEff=sqrt(2.0/3)*mag(symm(fvc::grad(u)));
Z=etaEff*exp(deltaH/(R*T));
eta=asinh(pow(Z/A,1.0/n))/(3*beta*etaEff);
eta_d=eta/rho;

// end of adding


fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(eta_d, U) // right?? change nu to eta_d
);


...........................


// and add T field calculation
fvScalarMatrix TEqn
(
fvm::ddt(T)
+fvm::div(phi,T)
-fvm::laplacian(DT,T)
);
TEqn.solve();
runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;