
[Sponsors] 
Darcy Forchheimer: trying to understand the code and how to run with it? 

LinkBack  Thread Tools  Display Modes 
February 12, 2015, 18:12 
Darcy Forchheimer: trying to understand the code and how to run with it?

#1 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
Hi all,
I'm facing strange behavior using fvOptions and simpleFoam. I mean by setting extrapolated values from experimental data I don't have same pressure losses on my simulation, but lower. The case is incompressible with MRF source for momentum. here is th fvOption: Code:
MRF1 { type MRFSource; active true; selectionMode cellZone; cellZone MRF; MRFSourceCoeffs { nonRotatingPatches (wallStat boccaglio AMI_ROT1 AMI_ROT2 AMI_ROT3);/ origin (0 0 0); axis (0 0 1); omega 300; } } ….... porosity2 { type explicitPorositySource; active true; //yes; selectionMode cellZone; cellZone filtro; explicitPorositySourceCoeffs { type DarcyForchheimer; DarcyForchheimerCoeffs { d d [0 2 0 0 0 0 0] (1790 1790 1790); f f [0 1 0 0 0 0 0] (1640 1640 1640); coordinateSystem { type cartesian; origin (0.05 0.11 0.15); coordinateRotation { type axesRotation; e1 (1 0 0); e2 (0 1 0); //e3 (0 0 1); } } } } } In order to set the D (linear?)& F(square?) I made some manipulation on experimental data. Experimental data consists of pressure losses (static) at different mean velocities. For each static pressure, I divided the value by thickness of filter (10mm), then by plotting the pressure/meters losses vs. velocity I was able to extrapolate a(=F) & b(=D) values by tendency line representation on the graph, forcing the curve passing by point (0,0). 1st question: is it correct to divided pressure losses by thickness of filter? I mean I followed a tutorial (not openfoam) and I think this division take into account the number of cell in order to interpolate pressure and velocity values, is it correct? Further by digging into the code I see on DarcyForchheimer.H: Code:
public porosityModel { private: // Private data // Darcy coeffient XYZ components (usersupplied) [1/m2] dimensionedVector dXYZ_; // Forchheimer coeffient XYZ components (usersupplied) [1/m] dimensionedVector fXYZ_; // Darcy coefficient  converted from dXYZ [1/m2] List<tensorField> D_; // Forchheimer coefficient  converted from fXYZ [1/m] List<tensorField> F_; // Name of density field word rhoName_; ….. 2nd question: how openfoam calculates rho as I specify only nu 1.5e5 m^2/s in transport properties? Further into DarcyForchheimer.C: Code:
void Foam::porosityModels::DarcyForchheimer::calcTranformModelData() { if (coordSys_.R().uniform()) { forAll (cellZoneIDs_, zoneI) { D_[zoneI].setSize(1); F_[zoneI].setSize(1); D_[zoneI][0] = tensor::zero; D_[zoneI][0].xx() = dXYZ_.value().x(); D_[zoneI][0].yy() = dXYZ_.value().y(); D_[zoneI][0].zz() = dXYZ_.value().z(); D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]); // leading 0.5 is from 1/2*rho F_[zoneI][0] = tensor::zero; F_[zoneI][0].xx() = 0.5*fXYZ_.value().x(); F_[zoneI][0].yy() = 0.5*fXYZ_.value().y(); F_[zoneI][0].zz() = 0.5*fXYZ_.value().z(); F_[zoneI][0] = coordSys_.R().transformTensor(F_[zoneI][0]); } } else { forAll(cellZoneIDs_, zoneI) { const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]]; D_[zoneI].setSize(cells.size()); F_[zoneI].setSize(cells.size()); forAll(cells, i) { D_[zoneI][i] = tensor::zero; D_[zoneI][i].xx() = dXYZ_.value().x(); D_[zoneI][i].yy() = dXYZ_.value().y(); D_[zoneI][i].zz() = dXYZ_.value().z(); // leading 0.5 is from 1/2*rho F_[zoneI][i] = tensor::zero; F_[zoneI][i].xx() = 0.5*fXYZ_.value().x(); F_[zoneI][i].yy() = 0.5*fXYZ_.value().y(); F_[zoneI][i].zz() = 0.5*fXYZ_.value().z(); } const coordinateRotation& R = coordSys_.R(mesh_, cells); D_[zoneI] = R.transformTensor(D_[zoneI], cells); F_[zoneI] = R.transformTensor(F_[zoneI], cells); } } 3rd question: can someone help me please? Thanks a lot! Bye 

February 12, 2015, 19:59 

#2 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
Hi,
I work a little on angleExplicit case on porousSimpleFoam folder. If in 0 folder I keep only U & p file, computation starts anyway. I'm a little in confusion about this topic, as I thought rho was calculated by perfect Gas equation, but removing T file, computation start and terminate with no errors. In any case, I discovered a mismatch between my BCs and the tutorial case: in the tutorial case a BC for the porosity walls is used and set to slip, but if I change it into wall with fixed uniform value (0 0 0) it works anyway...so what's the difference with it? thank you for any help. Bye 

February 13, 2015, 10:50 

#3 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
I test the DarcyForchHeimer on a test case based on measurments of pressure losses on a porous media using porousSimpleFoam and porosity dict:
https://www.dropbox.com/s/w27xlva4ftbw908/case.JPG?dl=0 it's 2D case of a channel of 21x2 meters with one porous zone of 1meter in the middle. As it's possibile to see, I calculated the D & F coefficient by plotting on a graph and drawing the tendecy line passing through (0,0). These values are used into porosity dict for the DarcyForchheimer model, but red lines show results of calculation, higly disagree with measruments, any suggestion? As it possible to see on the following picture, after the porous zone (fluid flows in x+ direction) velocity has a 1.06m/s value, while at inlet the BC is 0.88m/s: I was expecting a reduction on velocity after the porous media. https://www.dropbox.com/s/duxz3jeu4hjghjr/fdsf.png?dl=0 Can someone give me an explanation? Here is the case, with the mesh.unv https://www.dropbox.com/s/ixu1pdwevbpg88b/case.zip?dl=0 Thanks a lot 

February 13, 2015, 11:16 

#4 
Senior Member
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,409
Rep Power: 24 
Hi,
Why do you think that 5000 iterations is enough for the case to converge? 

February 13, 2015, 11:40 

#5 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
Hi,
thanks for reply. Here's a plot of the residual; https://www.dropbox.com/s/sr7kftxrvk...nning.png?dl=0 I see residual floating over very low values; honestly I haven't checked a physical quantity, as well I haven't set up a turbulent case, but a laminar one. In any case, I can't see how velocity should raise up its value after the porous media. Should it be as for the conservativness of the FV method? I mean as for the BC applied, mass conservation has to be kept, so higher pressure before porous media gives a "bump" to the fluid, but as for conservativness the calculation leads to have same mass flow between inlet and outlet and, as for constant section area, velocity can't decrease. 

February 13, 2015, 11:56 

#6 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
Hi,
here's a plot of the sum of inlet and outlet massflow. As you can see after few hundred iterations it comes to 0. https://www.dropbox.com/s/1r97qj99jj...sflow.png?dl=0 This should explain the local raise for velocity value. How can I the set up correctly the Darcy Forchheimer model? 

March 25, 2015, 10:22 

#7 
Senior Member
M. C.
Join Date: May 2013
Location: Italy
Posts: 122
Rep Power: 5 
Hi all,
I don't know why, but the DarcyForchHeimer model didn't work for me. I found this useful link: http://www.geocities.co.jp/SiliconVa...fvOptions.html and by changing model to FixedCoeff, I was able to solve my problem. Anyway, these are the steps I performed to test the model on a simple case: a straight 2D channel of constant section with a predefined thickness for the porous zone. 1  I divided pressure losses (experimental data) by the thickness of the porous zone: call this A 2  Using excel, I plotted A vs. velocity and ask to have the alfa e beta coefficient for the tendency line. 3  in the fvOptionDict I set these values for alfa & beta for the FixedCoeff model. Code:
porosity1 { type explicitPorositySource; active true; //yes; selectionMode cellZone; cellZone porous; explicitPorositySourceCoeffs { type fixedCoeff; fixedCoeffCoeffs { alpha alpha [0 0 1 0 0 0 0] (25000 25000 25000); //linear term beta beta [0 1 0 0 0 0 0] (10000 10000 10000); //squared term coordinateSystem { type cartesian; origin (1.005 0.25 0); rho 1.205; coordinateRotation { type axesRotation; e1 (1 0 0); e2 (0 1 0); //e3 (0 0 1); } } } } } Postprocessing the results I had the matching between the experimental and numerical results (within a small tolerance). Hope this can help. Regards 

Thread Tools  
Display Modes  

