
[Sponsors] 
Viscoelastic Fluid Flows using OpenFOAM The solver viscoelasticFluidFoam 

LinkBack  Thread Tools  Display Modes 
May 16, 2009, 13:57 

#21 
New Member
Kerstin Heinen
Join Date: Mar 2009
Location: Ludwigshafen, Germany
Posts: 27
Rep Power: 8 
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 

May 17, 2009, 11:10 

#22 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 18, 2009, 18:47 

#23 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 19, 2009, 12:27 

#24 
New Member
Kerstin Heinen
Join Date: Mar 2009
Location: Ludwigshafen, Germany
Posts: 27
Rep Power: 8 
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.... 

May 19, 2009, 16:25 

#25 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
Hi,
Thank you for your reply. Do you mean D^d is equal to D? Then tau is a scalar? Thanks. Wendy 

May 19, 2009, 16:26 

#26 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
Sorry, I typed D : D


May 19, 2009, 16:32 

#27 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
I guess you are right, I agree with that tau= 2 eta( D,T) D. tau is a tensor


May 19, 2009, 21:06 

#28  
Member
MrFluent
Join Date: Mar 2009
Posts: 33
Rep Power: 8 
Quote:
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. 

May 19, 2009, 23:48 

#29 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 20, 2009, 01:44 

#30 
New Member
Kerstin Heinen
Join Date: Mar 2009
Location: Ludwigshafen, Germany
Posts: 27
Rep Power: 8 

May 20, 2009, 10:56 

#31 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 20, 2009, 11:01 

#32 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
forgot to attach it. Sorry.
Wendy 

May 20, 2009, 11:32 

#33 
Member

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 CarreauYassuda, for example: volScalarField etaEff = eta0 * Foam:ow( 1 + Foam:ow( 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 

May 20, 2009, 11:49 

#34 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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. 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 

May 20, 2009, 13:01 

#35 
Member

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 

May 20, 2009, 18:57 

#36 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 20, 2009, 21:08 

#37 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 20, 2009, 21:19 

#38 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
Hi, Jovani,
I found the PDF document of programmer's guide. Sorry for disturbing again. Wendy 

May 21, 2009, 22:49 
the way I changed the constitutive model to viscoplastic

#39 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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 

May 21, 2009, 22:50 

#40 
Member
xianghong wu
Join Date: Mar 2009
Posts: 57
Rep Power: 8 
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] 1e5; 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.224e8 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)); // ZenerHollomon 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; 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
VOF simulation of a viscoelastic fluid  sinah  OpenFOAM Running, Solving & CFD  10  November 25, 2010 12:02 
FREE SURFACE VISCOELASTIC FLOWS  Valdemir G. Ferreira  Main CFD Forum  6  December 18, 2009 07:14 
Viscoelastic flow modeling in OpenFOAM  vulda  OpenFOAM Running, Solving & CFD  1  March 17, 2008 08:32 
Polyflow & OpenFoam on Viscoelastic flow modeling  Sumeshen  Main CFD Forum  0  March 14, 2008 09:29 
Viscoelastic fluid codes  joel davison  Main CFD Forum  0  November 6, 2001 06:09 