CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (

maritozzo October 19, 2005 16:19

Does nayone knows if OpenFoam
Does nayone knows if OpenFoam is able to simulate a porous material.
Please answer

billy October 19, 2005 23:29

I think it depends on what do
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.

gschaider October 20, 2005 05:27

You simply have to add the Dar
You simply have to add the Darcy-term to the velocity-equations (in other words: write a solver for this)

But there is no out-of-the-box solver for this (and there is always the problem how you define which parts of the geometry are porous and which are not)

maritozzo October 20, 2005 12:09

so substantially you are telli
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?

gschaider October 20, 2005 12:22

Sorry. With "out-of-the-box" I
Sorry. With "out-of-the-box" I meant "a solver available in the OF-distribution". You'll have to write a solver yourself.

One approach would be to simply extend an existing solver by adding Darcy as a source-term. In the Darcy-term 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/non-porous-zones.

That approach has some problems at the porous/non-porous-interfaces.

maritozzo October 20, 2005 12:40

You're telling me to treat it
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?

gschaider October 24, 2005 07:14

@boundary condition: No. You j
@boundary condition: No. You just implement an additional source-term. I just wanted to warn you, that there might be problems at the porous/non-porous-interface

Is your whole simulation domain a porous body with a homogenous permeability or is your setup of the type: inlet/Navies-Stokes-Flow/porous body (Darcy)/NS-Flow/outlet? My above comments were always meant for the second case.

maritozzo October 24, 2005 15:22

I'm in the second condition,wi
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?

gschaider October 24, 2005 18:31

No. Sorry. But I'll tell you w
No. Sorry. But I'll tell you when I find out.

Anyway: these effects are only significant for low permeabilities.

maritozzo October 25, 2005 06:48

I found the Navier-Stokes-Brin
I found the Navier-Stokes-Brinkman equation
which is


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)


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


solve(UEqn() == -fvc::grad(p));

Please tell me if i've made some syntax error

gschaider October 25, 2005 07:12

This looks good. I would sugge
This looks good. I would suggest using fvm::Sp or fvm::SuSp for an implicit treatment of the source term.

maritozzo October 25, 2005 07:21

You mean like this tmp
You mean like this

tmp<fvvectormatrix> UEqn
fvm::div(phi, U)
+ turbulence->divR(U)+nu*fvm::SuSp(U,G)


solve(UEqn() == -fvc::grad(p));

P.S. ; I'dont have to worry about the lack of density?

maritozzo October 25, 2005 07:23

Sorry G is the inverse of K te
Sorry G is the inverse of K tensor

gschaider October 25, 2005 07:28

If you switch positions of U a
If you switch positions of U and G: Yes.

maritozzo October 25, 2005 11:57

Sorry Bernhard but from my que
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

mattijs October 25, 2005 13:31

Density: the incompressible co
Density: the incompressible codes tend to solve for U/rho, p/rho etc. (see the dimensions in the fields) with density set to 1.

gschaider October 28, 2005 06:59

Hi Muzio! When developing a
Hi Muzio!

When developing a new solver you should be prepared to edit the files for the initial conditions.

Copy an existing field-file (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 damBreak-derived utility on that (or you might have a look at the setFields-utility 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.

maritozzo November 7, 2005 08:31

Hi Bernhard I wrote the solver
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

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"


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

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

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

# include "readSIMPLEControls.H"


// Pressure-velocity SIMPLE corrector
// Momentum predictor

tmp<fvvectormatrix> UEqn
fvm::div(phi, U)
+ (1.0+2.5*(1.0-d))*turbulence->divR(U)+nu*G*U


solve(UEqn() == -fvc::grad(p));

volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);

// Non-orthogonal 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);

if (nonOrth == nNonOrthCorr)
phi -= pEqn.flux();

# include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector

// Momentum corrector
U -= fvc::grad(p)/AU;



Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;

Info<< "End\n" << endl;


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

PLEASE HELP ME find out where is my error

niklas November 7, 2005 11:41

what error? first though, n
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?

maritozzo November 7, 2005 12:20

I'm analysing the flow of an i
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
fvm::SuSp(nu*G,U)where G=(K)^-1

i discretize only U and not the product -((nu)*G)*U

All times are GMT -4. The time now is 10:36.