Darcy-Forchheimer law for specifying Porous Zones

 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 30, 2010, 07:53 Darcy-Forchheimer law for specifying Porous Zones #1 New Member   usn Join Date: Jul 2010 Location: Deutschland Posts: 9 Rep Power: 13 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 FelicityNWU likes this.

 July 30, 2010, 09:35 #2 New Member   Vincent de Graaf Join Date: Jun 2010 Location: Germany Posts: 12 Rep Power: 13 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]

August 3, 2010, 05:30
#3
Senior Member

Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,372
Rep Power: 33
Quote:
 Originally Posted by VdG 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.

August 3, 2010, 06:34
#4
New Member

usn
Join Date: Jul 2010
Location: Deutschland
Posts: 9
Rep Power: 13
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
Attached Images
 ex-1.jpg (17.9 KB, 891 views)

 August 3, 2010, 06:42 #5 New Member   Vincent de Graaf Join Date: Jun 2010 Location: Germany Posts: 12 Rep Power: 13 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. Rojj likes this.

August 3, 2010, 07:44
#6
Senior Member

Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,372
Rep Power: 33
Quote:
 Originally Posted by Ger_US 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?

 August 3, 2010, 10:47 #7 New Member   usn Join Date: Jul 2010 Location: Deutschland Posts: 9 Rep Power: 13 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

 August 3, 2010, 11:02 #8 New Member   Vincent de Graaf Join Date: Jun 2010 Location: Germany Posts: 12 Rep Power: 13 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 );

 August 6, 2010, 05:33 #9 New Member   usn Join Date: Jul 2010 Location: Deutschland Posts: 9 Rep Power: 13 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:rintStack(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 > > > >::calculate() in "/server/appl/Programme/OpenFOAM/OpenFOAM-1.6/lib/linux64GccDPOpt/libbasicThermophysicalModels.so" [0] #4 Foam::hPsiThermo > > > >::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); Out: type zeroGradient; wall: type fixedValue; value uniform (0 0 0); P: internalField uniform 134980; In: type fixedValue; value uniform 152697; Out: type fixedValue; value \$internalField; wall: type zeroGradient; T: internalField uniform 882.8; In: type fixedValue; value uniform 879.23; Out: type fixedValue; value uniform 882.86; wall: type zeroGradient; K: internalField uniform 64.3817; In: type fixedValue; value uniform 64.3817; Out:type zeroGradient; wall: type kqRWallFunction; value uniform 64.3817; Epsilon: internalField uniform 16975.07; In: type fixedValue; value uniform 16975.07; Out:type zeroGradient; 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

 August 6, 2010, 10:04 #10 Senior Member   Join Date: Jan 2010 Location: Stuttgart Posts: 150 Rep Power: 13 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

 August 7, 2010, 01:09 #11 New Member   Robert Join Date: Apr 2010 Posts: 16 Rep Power: 13 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

October 20, 2010, 06:11
#12
New Member

Marc-Florian Uth
Join Date: Jan 2010
Posts: 10
Rep Power: 14
Quote:
 Originally Posted by rob3rt 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?

 March 11, 2013, 05:09 #13 New Member   hesam Join Date: Mar 2013 Posts: 13 Rep Power: 10 hi can anyone explain to me what are a,b which are use in equ.? excuse me if i have mistake in my writing .

 May 3, 2013, 09:20 references #14 Senior Member     maddalena Join Date: Mar 2009 Posts: 436 Rep Power: 20 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: esujby, manuc and sourav90 like this. Last edited by maddalena; May 3, 2013 at 10:49.

May 12, 2013, 04:56
#15
New Member

hesam
Join Date: Mar 2013
Posts: 13
Rep Power: 10
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 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:

May 13, 2013, 03:30
#16
Senior Member

Join Date: Mar 2009
Posts: 436
Rep Power: 20
Hi,
I really cannot understand what you miss on your example. Anyway:
Quote:
 Originally Posted by hesamgh 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.

May 19, 2013, 12:08
#17
New Member

hesam
Join Date: Mar 2013
Posts: 13
Rep Power: 10
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 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

August 30, 2013, 10:30
same problem
#18
Member

Reza
Join Date: Feb 2012
Posts: 67
Rep Power: 11
Quote:
 Originally Posted by hesamgh 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 ...

November 11, 2014, 05:50
#19
New Member

gaetan Hamel
Join Date: May 2014
Posts: 4
Rep Power: 12
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?

 November 11, 2014, 12:54 #20 New Member   Marco Join Date: Nov 2014 Location: Germany Posts: 14 Rep Power: 9 Yes, check out: http://www.tfd.chalmers.se/~hani/kur...ukurReport.pdf And therefore f = a and d = b sourav90 likes this.