[maritozzo]

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

[billy]

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]

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]

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]

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]

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]

@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]

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]

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

Anyway: these effects are only significant for low permeabilities.

[maritozzo]

I found the Navier-Stokes-Brinkman 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

[gschaider]

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

[maritozzo]

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?

[maritozzo]

Sorry G is the inverse of K tensor

[gschaider]

If you switch positions of U and G: Yes.

[maritozzo]

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]

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]

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]

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();

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

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);

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

[niklas]

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]

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