
[Sponsors] 
October 19, 2005, 16:19 
Does nayone knows if OpenFoam

#1 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
Does nayone knows if OpenFoam is able to simulate a porous material.
Please answer 

October 19, 2005, 23:29 
I think it depends on what do

#2 
Senior Member
Billy
Join Date: Mar 2009
Posts: 167
Rep Power: 9 
I think it depends on what do you want to simulate.
OpenFOAM is flexible and can simulate many things. If there is no application that suits your needs, you can always build your own. 

October 20, 2005, 05:27 
You simply have to add the Dar

#3 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
You simply have to add the Darcyterm to the velocityequations (in other words: write a solver for this)
But there is no outofthebox solver for this (and there is always the problem how you define which parts of the geometry are porous and which are not)
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 20, 2005, 12:09 
so substantially you are telli

#4 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
so substantially you are telling me that i have to connect two different solvers, one when i'm out of the box and another one when i'm in the box, have you ever done something like this?


October 20, 2005, 12:22 
Sorry. With "outofthebox" I

#5 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
Sorry. With "outofthebox" I meant "a solver available in the OFdistribution". You'll have to write a solver yourself.
One approach would be to simply extend an existing solver by adding Darcy as a sourceterm. In the Darcyterm there is a permeability/resistivity (whatever formulation you prefer). By using a field for that and specifying appropriate values for certain regions you can define porous/nonporouszones. That approach has some problems at the porous/nonporousinterfaces.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 20, 2005, 12:40 
You're telling me to treat it

#6 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
You're telling me to treat it as a boundary condition?
I have to make the flux pass through a zone which i define as a porous zone and i have to evaluate the consequent pressure loss Sorry if i'm a little bit boring but i have to do this as a thesis job, so what kind of approach do you think would be best? 

October 24, 2005, 07:14 
@boundary condition: No. You j

#7 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
@boundary condition: No. You just implement an additional sourceterm. I just wanted to warn you, that there might be problems at the porous/nonporousinterface
Is your whole simulation domain a porous body with a homogenous permeability or is your setup of the type: inlet/NaviesStokesFlow/porous body (Darcy)/NSFlow/outlet? My above comments were always meant for the second case.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 24, 2005, 15:22 
I'm in the second condition,wi

#8 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
I'm in the second condition,with inlet, NS flow, porous body, NS flow, and outlet.
Do you know a different approach to avoid the problems at the porous non porous inteface? 

October 24, 2005, 18:31 
No. Sorry. But I'll tell you w

#9 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
No. Sorry. But I'll tell you when I find out.
Anyway: these effects are only significant for low permeabilities.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 25, 2005, 06:48 
I found the NavierStokesBrin

#10 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
I found the NavierStokesBrinkman equation
which is U*grad(U)nu*laplacian(U)nu*U/K=(1/rho)*grad(p) where nu is the kinematic viscosity and K represents the permeability tensor. So I should simply add the term nu*U/K I looked at simpleFoam solver and i found this tmp<fvvectormatrix> UEqn ( fvm::div(phi, U) + turbulence>divR(U) ); UEqn().relax(); solve(UEqn() == fvc::grad(p)); I have some questions 1)Where is the division by rho regarding grad(p), because i foundend in the programmer guide that R represents nu(eff)*gradU and nu is the kinematic visvosity. 2)If the equation is right i should simply define the K tensor in the porous media domain the change the equation to this form tmp<fvvectormatrix> UEqn ( fvm::div(phi, U) + turbulence>divR(U)+(nu/K)*U ); UEqn().relax(); solve(UEqn() == fvc::grad(p)); Please tell me if i've made some syntax error 

October 25, 2005, 07:12 
This looks good. I would sugge

#11 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
This looks good. I would suggest using fvm::Sp or fvm::SuSp for an implicit treatment of the source term.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 25, 2005, 07:21 
You mean like this
tmp

#12 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
You mean like this
tmp<fvvectormatrix> UEqn ( fvm::div(phi, U) + turbulence>divR(U)+nu*fvm::SuSp(U,G) ); UEqn().relax(); solve(UEqn() == fvc::grad(p)); P.S. ; I'dont have to worry about the lack of density? 

October 25, 2005, 07:23 
Sorry G is the inverse of K te

#13 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
Sorry G is the inverse of K tensor


October 25, 2005, 07:28 
If you switch positions of U a

#14 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
If you switch positions of U and G: Yes.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

October 25, 2005, 11:57 
Sorry Bernhard but from my que

#15 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
Sorry Bernhard but from my questions i think you understood i'm very new in using OpenFoam, so i need some more hints, i know i'm boring, i'm trying to define the field of the permeability tensor and i really don't know how to begin.
I thaught to convert setGammaDambreak utility to my case but it only defines the initial field and it works on an already existing field it doesn't define a new one. In createFields.H i also observed that it creates the fields reading from existing files. Futhermore i think I should define the field as a volTensorField so that it refers to cell centers or i have to define it as a pointfield Please answer 

October 25, 2005, 13:31 
Density: the incompressible co

#16 
Super Moderator
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 17 
Density: the incompressible codes tend to solve for U/rho, p/rho etc. (see the dimensions in the fields) with density set to 1.


October 28, 2005, 06:59 
Hi Muzio!
When developing a

#17 
Assistant Moderator
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,926
Rep Power: 41 
Hi Muzio!
When developing a new solver you should be prepared to edit the files for the initial conditions. Copy an existing fieldfile (p for example) and edit it: 1. set the correct dimensions 2. set correct boundary conditions (zeroGradient should be alright for the permeability). The you can use damBreakderived utility on that (or you might have a look at the setFieldsutility which is new in 1.2 and might just do what you're looking for). @tensorField: do you have directed permeabilities? if not a scalar would be sufficient.
__________________
Note: I don't use "Friend"feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request 

November 7, 2005, 08:31 
Hi Bernhard I wrote the solver

#18 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
Hi Bernhard I wrote the solver and i used it but i had very bad results (negative pressure for example).
I found that when i write nu*fvm::SuSp(U,G) this only discretizes G in a implicit or explicit way depending on the sign of U, and that is not correct, so i discretized it as an explicit term This is the .C file I used #include "fvCFD.H" #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H" #include "incompressible/turbulenceModel/turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" # include "createFields.H" # include "created.H" # include "createG.H" # include "createNu.H" # include "initContinuityErrs.H" //mesh.clearPrimitives(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; # include "readSIMPLEControls.H" p.storePrevIter(); // Pressurevelocity SIMPLE corrector { // Momentum predictor tmp<fvvectormatrix> UEqn ( fvm::div(phi, U) + (1.0+2.5*(1.0d))*turbulence>divR(U)+nu*G*U ); UEqn().relax(); solve(UEqn() == fvc::grad(p)); p.boundaryField().updateCoeffs(); volScalarField AU = UEqn().A(); U = UEqn().H()/AU; UEqn.clear(); phi = fvc::interpolate(U) & mesh.Sf(); adjustPhi(phi, U, p); // Nonorthogonal pressure corrector loop for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::laplacian(1.0/AU, p) == fvc::div(phi) ); fvScalarMatrix::reference pRef = pEqn.setReference(pRefCell, pRefValue); pEqn.solve(); pEqn.unsetReference(pRef); if (nonOrth == nNonOrthCorr) { phi = pEqn.flux(); } } # include "continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Momentum corrector U = fvc::grad(p)/AU; U.correctBoundaryConditions(); } turbulence>correct(); runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s\n\n" << endl; } Info<< "End\n" << endl; return(0); } // ********** PLEASE HELP ME find out where is my error 

November 7, 2005, 11:41 
what error?
first though, n

#19 
Super Moderator
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 21 
what error?
first though, nu*fvm::SuSp(U,G) is not the same as fvm::SuSp(nu*G, U) which is what you should use, second: why cant the pressure be negative? 

November 7, 2005, 12:20 
I'm analysing the flow of an i

#20 
Member
Muzio Grilli
Join Date: Mar 2009
Posts: 36
Rep Power: 9 
I'm analysing the flow of an incompressible fluid through a porous media inside and outside it, the therm which i have to introduce is ((nu)/K)*U
writing fvm::SuSp(nu*G,U)where G=(K)^1 i discretize only U and not the product ((nu)*G)*U 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
nonisotropic porous material  gmmh  FLUENT  0  September 4, 2007 06:38 
Error while setting up porous material.  Kiddo  CFX  1  October 10, 2005 10:42 
porous material  ioana  CFX  2  March 10, 2005 08:52 
model for porous material  sleepinglily  CFX  3  October 19, 2004 10:45 
Material in Porous media  Rajab Rajab  FLUENT  3  July 4, 2003 13:42 