porousSimpleFoam - comparison to Darcy's law
 

Dear community,

I want to solve a the flow through a porous medium with constant permeability. For the beginning I wanted to start with an easy test case and try to reproduce the pressure drop given by Darcy's law.
Therefore I choose a rectangular geometry with given inflow velocity (0.01 m/s) in x-direction and in the middle a porous medium with constant thickness (1 cm) and constant permeability (1e-10 m^2). The boundary conditions at the walls are chosen as symmetric in y-direction and as empty in z-direction.
As fluid I use water (rho=1000 kg/m^3, mu = 1 mPas).
For this setup (and the assumption of an one-dimensional flow), Darcy's law yields a pressure drop delta p of 1000Pa.

If I use the porousSimpleFoam solver from openFoam, the results are quite different depending on the number of discretization cells for the porous medium (the files you find attached).

I have the feeling that in the openFoam solution one discretization cell is missing. If the discretization cells are equistant, the pressure drop over all cells are the same. Thus Darcy's law yields for the pressure drop in one discretization cell 1000 Pa / N, where N is the number of cells. After some computations for different number of cells, I think that openFoam is just getting the pressure drop of N-1 cells thus, the result from openFoam is (N-1)/N*1000 Pa.
Let me illustrate this with a few numbers from my computations:

https://picload.org/image/rowcriaa/p...comparison.png

Since I am just started with openFoam, I am not sure, if I did something wrong. Therefore I attached my files. Maybe someone can tell me, if they are correct.
If my results would be correct, is there some bug in openFoam for this example? Is openFoam neglecting or forgetting about one cell? Is this in some way related to some interpolation scheme?

If some more data are needed, please tell me.

Thank you in advance for your answers and your help.

Regards

Sebastian

Dear community,

if nobody knows a solution or reason for the above problem, can I consider it as bug in openFoam?

Maybe someone can check, if my geometry is correct. Otherwise I might expect something completely different from the solution.

Some more information: I am using version v1606+ on a virtual machine on windows. Is something wrong with this version?

Thank you in advance for your answers.

Regards

Sebastian

Porosity issues with Openfoam v5
 

Hi Sebastian,

I found the same strange behaviour for porosities as you have. I'm using simpleFoam, but also here, the following things happen:

- with Darcy coefficients d and f calculated as already often stated here in the forum, I get nonsense (much too high) pressure drop. The pressure is nonuniform over the porous zone with spots etc and the convergence is horrible. Speed behaves accordingly. Pressure drop is grid dependent (by a factor of 2!! Dense grid 60Pa, coarse grid 120Pa). Calculations done with Openfoam 5.
- same case setup calculated with OF2.3.x: nice pressure values, no spots, good convergence, minimal grid dependence.

Did you get your problem resolved? Is there a special setting for OF5?

Regards,
-Thomas

Dear Thomas,

to answer your last questions first:
I have not resolved the problem and I do not know about some new features in OF5.
I can not give you an answer, if there is a difference between the solver in OF5 and OF2.3, maybe someone else can give us some more insight here.

Last year, I wrote to the developers to report this bug, however they said that there is no capacity in the moment to work on this issue. They recommended to use some different discretization scheme (pointCellsLeastSquares grad) but this had no influence on the solution, at least in my case. Maybe you can also have a try.
Therefore I would expect any improvement or fix in measurable time.

But maybe a short question to your case.
In my experience usually the pressure drop is too low because (as reported) openFoam seems to miss one cell. Therefore I am wondering that you are reporting a too high pressure drop. Maybe you can give me some information on the setup and geometry you are using.

Regards

Sebastian

Hello Sebastian,

I assumed to have a pressure drop curve following dp=4*v^2+13*v with a porous zone of thickness 94mm. I assumed a kinematic viscosity of 1.e-5 and a density of 1.0, and by comparing the coefficients of my pressure drop curve and Darcys law, I found my values to be f=85.11 and d=13829787.2. With these values put into fvOptions file and using simpleFoam I get the following results on the same grid:

v[m/s] dpCalc[Pa] dpOF23x[Pa] dpOF5[Pa]
4 116 116 8.35e3
8 360 359 1.64e4
12 732 731 2.72e4

The results for OF5 are close to nonsense. I cannot exclude an user error but as I cannot find a difference in the documentation between OF2.3.x and OF5 how porosity is implemented I'd say that porosity in OF5 and its variants is useless.

If you'd like to look at my cases, I'd pack them up and pass them to you.

Regards,
-Thomas

One more thing I've found: in the FvOptions file, the setting for d is (13829787.2 -1000 -1000) in order to block the directions perpendicular to the x direction. If I change the values from -1000 to 0, the presure drop is much lower, for 8m/s velocity for example down to 816Pa, which is still way off from what is expected.

Hello Thomas,

I am not an expert in computing the coefficients since I am mostly using a given permeability value for the porous medium and then d is directly clear and my flow is usually that slow that I am not using f. However, in https://www.cfdsupport.com/OpenFOAM-...t/node232.html the computation is different to https://www.cfd-online.com/Forums/op...tml#post664224.

Maybe this is the difference between the versions of openFOAM?

I always set all three values for d to the same value since I assume that I have an isotropic (or orthotropic) porous medium which has the same resistance in all three directions. Depending on your geometry, I don't know if it makes sense to try to block directions with the coefficients. Usually one assumes that the flow in the porous medium should be approximately assigned to the surface normal direction.
Going to the underlying equations a negative sign for d means that you have instead of a velocity sink a velocity source and therefore this might give different results since you are accelerating your flow (compare the definition of the source https://www.cfdsupport.com/OpenFOAM-...t/node232.html). You compared also the velocity fields of the two openFoam versions?

Kind regards

Sebastian

Hi guys,
Have you found a solotion of the mentioned problem?
I have few additional questions.
1) Why there is a porousSimpleFoam solver if the porous region can be implemented with fvOptions using any solver?
2) How I should modify codeAddSup in order to allow D and F variable coefficients? in particular, i do not fully understand the coordinateSystem set with regards to the source term: https://openfoamwiki.net/index.php/DarcyForchheimer
thanks in advance

Hey, sorry for my late reply.
Concerning the last reply I got from the developers, there is no fix of the mentioned problem. Also using a different scheme does not work and using a high number of cells for discretizing the porous medium does not work for me.
Concerning your questions:
1) Why there are two ways to use a porous region, I can not tell you. I have not experienced any differences in the use of the possibilities. But maybe there are some applications where this is the case. Maybe someone else can comment on this.
2) What do you mean with codeAddSup here? You might get varying coefficients when you define a field for the coefficients and use some setFields utility. However, I am not sure if you then also have to modify your solver to recognize those fields. If you have only a finite number of different coefficients, you might define multiple porous zones.
The coordinate system is basically defining your resistance tensor. Lets say the diagonal of your tensor is given as d = (d1 d2 d3).
Then using e1 = (1 0 0) and e2 = (0 1 0) yields e3 = (0 0 1) (cross product) and d1 will be the resistance value in the direction of e1, d2 in e2 and d3 in e3. If you use a different coordinate system, the directions will be changed. If for example, e1 = (0 1 0), e2 = (1 0 0), the third one is e3 = (0 0 -1) and the resistance in direction of e1 is d1 and so on. However, looking in the global coordinate system, in the first case e1 is the x-direction and in the second case e1 is the y-direction. This changes your results.

You can make a easy test, if you only put d1>0 and d2=d3=0. Then in the case of one-dimensional Darcy, changing the local coordinate system of the tensor will influence your results. If e1 is not (1 0 0), then the flow might experience no resistance at all.


I hope I could at least partly answer your questions.


Regards


Sebastian

Hi sepp.zell,

I run the same case on simpleFoam using fvOptions with a explicitPorositySource on the porous_media1 cell zone. Here are the results:

dp=951,6392 (10 cell in the porous_media1)
dp=995,0351 (100 cell in the porous_media1)
dp=999,268232 (1000 cell in the porous_media1) <-taken from paraview

These are the pressure difference from the center of the first and last cells of the porous zone. If you consider the pressure difference from face to face the results doesn't change significantly.
I don't know, it doesn't seem to me a matter of mesh refinement since the difference from 10 to 1000 is quite large for this slow improvement on the results.

P.S.: I set walls_Y and walls_Z as empty so the case is now trully 1D. I also add a residual control for the SIMPLE. All the cases converged to 1e-6.