CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Darcy-Forchheimer law for specifying Porous Zones (https://www.cfd-online.com/Forums/openfoam-solving/78705-darcy-forchheimer-law-specifying-porous-zones.html)

 Ger_US July 30, 2010 06:53

Darcy-Forchheimer law for specifying Porous Zones

Hi all,

I would like to use Darcy-Forchheimer law for specifying Porous Zones in Exhaust System application. But I am confused with coordinate system specification. Can anyone please explain what is e1 and e2?

and I have the values

a = 9.367
b = 1.029E7
alpha = 0.5 * a * density [kg/m^4]
beta = viscocity * b [kg/m^3s]

How can I calculate d and f parameters from the above data?

My guess: d= beta/viscocity [1/m^2]
f=alpha/density [1/m]

Am I correct?

if not, please tell me how to calculate d and f parameters.

thank you in advance

 VdG July 30, 2010 08:35

e1 and e2 are the vectors that are used to specify the porosity. In the porousZones file, you have to specify three components of f and d. The first component is in the direction of e1, the second in the direction of e2 and the third in the direction perpendicular to e1 and e2. An example can be found in tutorials/incompressible/porousSimpleFoam/angledDuctImplicit.

Furthermore,
d= beta/viscocity [1/m^2]
f=2*alpha/density [1/m]

 olesen August 3, 2010 04:30

Quote:
 Originally Posted by VdG (Post 269516) e1 and e2 are the vectors that are used to specify the porosity. In the porousZones file, you have to specify three components of f and d. The first component is in the direction of e1, the second in the direction of e2 and the third in the direction perpendicular to e1 and e2. An example can be found in tutorials/incompressible/porousSimpleFoam/angledDuctImplicit. Furthermore, d= beta/viscocity [1/m^2] f=2*alpha/density [1/m]
Note there is also useful little trick in darcy/forchheimer specification: since negative coefficients are physically meaningless, they can (mis)used to specify a multiplication factor for strongly anisotropic porous media. This can be quite convenient.

 Ger_US August 3, 2010 05:34

1 Attachment(s)
Hi,

thank you very much for your replies...here I am giving my data...can anyone please tell me e1 and e2 vectors are correct...I attached my porous part file (I gave not fully because it is confidential)

coordinateSystem
{

e1 (1 0 0);
e2 (0 0 1);
}

Darcy
{
d d [0 -2 0 0 0 0 0] (1.029e7 1 1);
f f [0 -1 0 0 0 0 0] (0 0 0);
}

Thanks

 VdG August 3, 2010 05:42

Hi Ger_US,

you have specified a Darcy resistance in x-direction, and virtually no resistance in y- and z-direction. If your intention is to model an anisotropic medium, this might be correct. Note that for isotropic media, all three components of d (and f) should have the same value.

 olesen August 3, 2010 06:44

Quote:
 Originally Posted by Ger_US (Post 269977) Hi, thank you very much for your replies...here I am giving my data...can anyone please tell me e1 and e2 vectors are correct...I attached my porous part file (I gave not fully because it is confidential) coordinateSystem { e1 (1 0 0); e2 (0 0 1); } Darcy { d d [0 -2 0 0 0 0 0] (1.029e7 1 1); f f [0 -1 0 0 0 0 0] (0 0 0); } Thanks
The e1/e2/e3 are the axes (x/y/z) of the local coordinate systems. There are some details here:
http://foam.sourceforge.net/doc/Doxy....html#_details

You may also wish to specify the 'origin' of your local coordinate system.
BTW: in you Darcy term, you have close to no resistance in the local y/z directions, but a fairly large one in the x-direction -- what type of porous media should this be?

 Ger_US August 3, 2010 09:47

Hi,

I think I am wrong...that is ceramic porous media...

this is information which I have

a = 9.367
b = 1.029E7
alphax = 0.5 * a * density [kg/m^4]
betax = viscocity * b [kg/m^3s]

d= beta/viscocity [1/m^2]
f=2*alpha/density [1/m]

alphay=1e6
betay=1e6

alphaz=1e6
betaz=1e6

viscocity=1.663e-5 kg/ms
density=0.61935 at Inlet

coordinateSystem
{

e1 (1 0 0);
e2 (0 0 1);
}

Darcy
{
d d [0 -2 0 0 0 0 0] (1.029e7 (6.01e10??) (6.01e10??));
f f [0 -1 0 0 0 0 0] (18.367 (?) (?));
}

please tell me if anything wrong because I am new to OpenFoam and CFD

I forgot to mention that it is anisotropic medium

Thanks

 VdG August 3, 2010 10:02

Hi,

from the description I assume that the medium has a certain resistance in the x-direction (d = 1.029e7, f = 18.367) and is impermeable in the other directions. My usual approach to this situation is to set the y- and z-resistance to a finite but much larger value than the x-resistance. Typically, I choose to set the y- and z-components of f to 1000 times the x-value. In your case, this would become:

f f [0 -1 0 0 0 0 0] (18.367 18e3 18e3 );

 Ger_US August 6, 2010 04:33

Hi,

as you suggested I tried without success. I don't why rho is not stable and bounding p.

Time = 92

smoothSolver: Solving for Ux, Initial residual = 0.00292406, Final residual = 2.10951e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00583214, Final residual = 0.000111279, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.013854, Final residual = 0.000225707, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0369572, Final residual = 0.000357835, No Iterations 1
GAMG: Solving for p, Initial residual = 0.165261, Final residual = 0.00597301, No Iterations 2
time step continuity errors : sum local = 2.06879, global = 0.061706, cumulative = 16.778
bounding p, min: -24708.8 max: 824884 average: 145678
rho max/min : 12.1778 0.0591955
smoothSolver: Solving for epsilon, Initial residual = 0.00295326, Final residual = 9.7231e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00285544, Final residual = 1.06368e-06, No Iterations 2
ExecutionTime = 354.46 s ClockTime = 355 s

Time = 93

smoothSolver: Solving for Ux, Initial residual = 0.00298496, Final residual = 2.15188e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.0047816, Final residual = 9.21071e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.0134613, Final residual = 0.000221754, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0368796, Final residual = 0.00035114, No Iterations 1
GAMG: Solving for p, Initial residual = 0.171999, Final residual = 0.00616598, No Iterations 2
time step continuity errors : sum local = 2.14248, global = 0.0565349, cumulative = 16.8345
bounding p, min: -28114.3 max: 848170 average: 145624
rho max/min : 12.769 0.0610522
smoothSolver: Solving for epsilon, Initial residual = 0.00294399, Final residual = 9.54477e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00287104, Final residual = 1.05367e-06, No Iterations 2
ExecutionTime = 358.29 s ClockTime = 359 s

Time = 94

smoothSolver: Solving for Ux, Initial residual = 0.00304401, Final residual = 2.19977e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00406858, Final residual = 7.70121e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00943227, Final residual = 0.000157461, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0366639, Final residual = 0.000342336, No Iterations 1
GAMG: Solving for p, Initial residual = 0.179248, Final residual = 0.00637306, No Iterations 2
time step continuity errors : sum local = 2.22181, global = 0.0517805, cumulative = 16.8863
bounding p, min: -29897.8 max: 870334 average: 145566
rho max/min : 14.1506 0.0643263
smoothSolver: Solving for epsilon, Initial residual = 0.00293563, Final residual = 9.37558e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00288938, Final residual = 1.04439e-06, No Iterations 2
ExecutionTime = 362.13 s ClockTime = 363 s

Time = 95

smoothSolver: Solving for Ux, Initial residual = 0.0030967, Final residual = 2.26804e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.0035207, Final residual = 6.69516e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00715987, Final residual = 0.000122852, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0362804, Final residual = 0.000331438, No Iterations 1
GAMG: Solving for p, Initial residual = 0.186293, Final residual = 0.00661338, No Iterations 2
time step continuity errors : sum local = 2.31546, global = 0.0470522, cumulative = 16.9333
bounding p, min: -32319.1 max: 893740 average: 145497
rho max/min : 15.0995 0.0689234
smoothSolver: Solving for epsilon, Initial residual = 0.00293385, Final residual = 9.21299e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00290919, Final residual = 1.03662e-06, No Iterations 2
ExecutionTime = 365.96 s ClockTime = 367 s

Time = 96

smoothSolver: Solving for Ux, Initial residual = 0.00314822, Final residual = 2.34295e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00319378, Final residual = 6.15854e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00594322, Final residual = 0.00010551, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0360695, Final residual = 0.000327196, No Iterations 1
GAMG: Solving for p, Initial residual = 0.194211, Final residual = 0.00690557, No Iterations 2
time step continuity errors : sum local = 2.42775, global = 0.0442164, cumulative = 16.9776
bounding p, min: -33417.8 max: 915435 average: 145415
rho max/min : 15.7051 0.0746939
smoothSolver: Solving for epsilon, Initial residual = 0.00293065, Final residual = 9.05813e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00293521, Final residual = 1.03033e-06, No Iterations 2
ExecutionTime = 369.81 s ClockTime = 370 s

Time = 97

smoothSolver: Solving for Ux, Initial residual = 0.00320052, Final residual = 2.42225e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00296225, Final residual = 5.77766e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00504761, Final residual = 8.91875e-05, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0357983, Final residual = 0.000325011, No Iterations 1
GAMG: Solving for p, Initial residual = 0.203505, Final residual = 0.00724398, No Iterations 2
time step continuity errors : sum local = 2.56053, global = 0.045573, cumulative = 17.0231
bounding p, min: -34279 max: 933095 average: 145318
rho max/min : 16.0785 0.0807459
smoothSolver: Solving for epsilon, Initial residual = 0.00291039, Final residual = 8.86735e-07, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.00298268, Final residual = 1.0239e-06, No Iterations 2
ExecutionTime = 373.66 s ClockTime = 374 s

Time = 98

smoothSolver: Solving for Ux, Initial residual = 0.00326285, Final residual = 2.48943e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.00276879, Final residual = 5.35671e-05, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.00437289, Final residual = 7.91479e-05, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0351593, Final residual = 0.000313212, No Iterations 1
[0] #0 Foam::error::printStack(Foam::Ostream&) in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #1 Foam::sigFpe::sigFpeHandler(int) in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #2 __restore_rt in "/lib64/tls/libc.so.6"
[0] #3 Foam::hPsiThermo<Foam::pureMixture<Foam::sutherlan dTransport<Foam::specieThermo<Foam::hConstThermo<F oam::perfectGas> > > > >::calculate() in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libbasicThermophysicalModels.so"
[0] #4 Foam::hPsiThermo<Foam::pureMixture<Foam::sutherlan dTransport<Foam::specieThermo<Foam::hConstThermo<F oam::perfectGas> > > > >::correct() in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libbasicThermophysicalModels.so"
[0] #5 main in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/applications/bin/linux64GccDPOpt/rhoPorousSimpleFoam"
[0] #6 __libc_start_main in "/lib64/tls/libc.so.6"
[0] #7 _start at ../sysdeps/x86_64/elf/start.S:116
[cn-12:16978] *** Process received signal ***
[cn-12:16978] Signal: Floating point exception (8)
[cn-12:16978] Signal code: (-6)
[cn-12:16978] Failing at address: 0x26c00004252
[cn-12:16978] [ 0] /lib64/tls/libc.so.6 [0x3e2842e2f0]
[cn-12:16978] [ 1] /lib64/tls/libc.so.6(gsignal+0x3d) [0x3e2842e25d]
[cn-12:16978] [ 2] /lib64/tls/libc.so.6 [0x3e2842e2f0]
[cn-12:16978] [ 3] /server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libbasicThermophysicalModels.so(_ZN4Foam10hPsiTher moINS_11pureMixtureINS_19sutherlandTransportINS_12 specieThermoINS_12hConstThermoINS_10perfectGasEEEE EEEEEE9calculateEv+0x4b7) [0x2a95597f77]
[cn-12:16978] [ 4] /server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libbasicThermophysicalModels.so(_ZN4Foam10hPsiTher moINS_11pureMixtureINS_19sutherlandTransportINS_12 specieThermoINS_12hConstThermoINS_10perfectGasEEEE EEEEEE7correctEv+0x33) [0x2a955a1303]
[cn-12:16978] [ 5] rhoPorousSimpleFoam [0x41cd9d]
[cn-12:16978] [ 6] /lib64/tls/libc.so.6(__libc_start_main+0xdb) [0x3e2841c3fb]
[cn-12:16978] [ 7] rhoPorousSimpleFoam [0x41ab29]
[cn-12:16978] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 16978 on node cn-12 exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------

Boundary conditions:

U:
In: type flowRateInletVelocity; flowRate 0.16448; value uniform (0 0 0);
wall: type fixedValue; value uniform (0 0 0);

P:
internalField uniform 134980;
In: type fixedValue; value uniform 152697;
Out: type fixedValue; value \$internalField;

T:
internalField uniform 882.8;
In: type fixedValue; value uniform 879.23;
Out: type fixedValue; value uniform 882.86;

K:
internalField uniform 64.3817;
In: type fixedValue; value uniform 64.3817;
wall:
type kqRWallFunction; value uniform 64.3817;

Epsilon:
internalField uniform 16975.07;
In: type fixedValue; value uniform 16975.07;
wall:
type kqRWallFunction; value uniform 16975.07;

Porous Zones:

Table_2
{
coordinateSystem
{
e1 (1 0 0);
e2 (0 0 1);
}

Darcy
{
d d [0 -2 0 0 0 0 0] (3.029e7 5.8275e10 5.8275e10);
f f [0 -1 0 0 0 0 0] (20.367 20367 20367);
}
}

Table_3
{
coordinateSystem
{
e1 (1 0 0);
e2 (0 0 1);
}

Darcy
{
d d [0 -2 0 0 0 0 0] (8.808e7 5.8275e10 5.8275e10);
f f [0 -1 0 0 0 0 0] (509.12 509120 509120);
}
}

But rho is stable when I used y and z coordinates of porouszones as 1 or 1000

I don't know why rho is not stable for the above Porous Zones parameters.

Can anyone please tell me where I am doing mistake

Thanks

 Chrisi1984 August 6, 2010 09:04

Hi,

check your discretization in fvSchemes and play a bit with the underrelaxation factors of p.

I realized that they should be really low (about 0.05) in exhaust systems.

And for a first calculation you should take a first order disrectization for the divSchemes (Gauss upwind).

Best regards

Chrisi

 rob3rt August 7, 2010 00:09

Hi Ger_US

I am wondering how do you get the a & b values? what are they?

I think your f value should be 9.367 instead of 18.367 ?

Kind Regards,
Robert

 Marc10 October 20, 2010 05:11

Quote:
 Originally Posted by rob3rt (Post 270624) Hi Ger_US I am wondering how do you get the a & b values? what are they? I think your f value should be 9.367 instead of 18.367 ? Kind Regards, Robert
Hi Ger_US!

I have the same question. What are the a & b values.
How can I calculate alpha and beta from a certain Permeability K?

 hesamgh March 11, 2013 05:09

hi
can anyone explain to me what are a,b which are use in equ.?

excuse me if i have mistake in my writing .

 maddalena May 3, 2013 08:20

references

a and b are the two coefficients in Deltap = a*v^2 + b*v, which defines the relation between pressure drop and velocity on the porous medium.
see here: https://www.sharcnet.ca/Software/Flu...e233.htm#36964
just as further reference for everybody interested on the subject:

 hesamgh May 12, 2013 03:56

where can i get these coefficient (a,b) for porous medium?

i have an example that have alpha ,beta which have different value in 3 aspects
would you please explain me where can i get them?

a = 9.367
b = 1.029E7
alphax = 0.5 * a * density [kg/m^4]
betax = viscocity * b [kg/m^3s]

d= beta/viscocity [1/m^2]
f=2*alpha/density [1/m]

alphay=1e6
betay=1e6

alphaz=1e6
betaz=1e6

Quote:
 Originally Posted by maddalena (Post 424959) a and b are the two coefficients in Deltap = a*v^2 + b*v, which defines the relation between pressure drop and velocity on the porous medium. see here: https://www.sharcnet.ca/Software/Flu...e233.htm#36964 just as further reference for everybody interested on the subject:

 maddalena May 13, 2013 02:30

Hi,
I really cannot understand what you miss on your example. Anyway:
Quote:
 Originally Posted by hesamgh (Post 426969) where can i get these coefficient (a,b) for porous medium?
a and b can be calculated with different techniques depending on your input. Simply scroll down the cited post and look for the one which best fit your need.

 hesamgh May 19, 2013 11:08

would you mind please tell me one of the way that i can reach to these 2 parameters (a,b)?
i have the velocity inlet and the pressure(zero gradient in of) ,are these values enough to calculate a,b?

i have confused for 2 months to calculate d & f coeff. in open foam ...

Quote:
 Originally Posted by maddalena (Post 427080) Hi, I really cannot understand what you miss on your example. Anyway: a and b can be calculated with different techniques depending on your input. Simply scroll down the cited post and look for the one which best fit your need. mad

 gooya_kabir August 30, 2013 09:30

same problem

did you solve your problem about how to get a and b ? If yes, please help me about that
Quote:
 Originally Posted by hesamgh (Post 428557) thanks again for your help maddalena would you mind please tell me one of the way that i can reach to these 2 parameters (a,b)? i have the velocity inlet and the pressure(zero gradient in of) ,are these values enough to calculate a,b? i have confused for 2 months to calculate d & f coeff. in open foam ...

 g.hamel November 11, 2014 05:50

Hello,

I reopen the post because the answers are bugging me.

if we have

Quote:
 alpha = 0.5 * a * density [kg/m^4] beta = viscocity * b [kg/m^3s]
and
Quote:
 d= beta/viscocity [1/m^2] f=2*alpha/density [1/m]
then d = b and f = a (!?), could somebody tell me if it is the case? Then why define alpha and beta?

I have read a few topics about this question but the definition of D and F remains unclear.