
[Sponsors] 
January 20, 2020, 10:39 
Porosity and the energy equation

#1 
Senior Member
Join Date: Dec 2019
Posts: 215
Rep Power: 7 
Hi,
I am trying to model a porous media using the darcy forchheimer model. I am using Openfoam 4 and the solver sonicFoam. For this purpose I have added a fvOptions dict: Code:
porosity_wall { type explicitPorositySource; active true; explicitPorositySourceCoeffs { type DarcyForchheimer; selectionMode cellZone; cellZone porosity; // Specify the name of the cellZone DarcyForchheimerCoeffs { // Negative coeffs are multiplied by largest positive coeff, // taking the magnitude, e.g. for 1000, coeff = 1e7*1000 = 1e10 d [0 2 0 0 0 0 0] (1e10 1e10 1e10); f [0 1 0 0 0 0 0] (0 0 0); coordinateSystem // Cartesian coordinates for the cellZone { x (1 0 0); y (0 1 0); #includeEtc "caseDicts/general/coordinateSystem/cartesianXY" } } } } But the problem is, that the temperature seems not to be affected by the porosity. I have looked up EEqn.H for sonicFoam: Code:
fvScalarMatrix EEqn ( fvm::ddt(rho, e) + fvm::div(phi, e) + fvc::ddt(rho, K) + fvc::div(phi, K) + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")  fvm::laplacian(turbulence>alphaEff(), e) == fvOptions(rho, e) ); Does anyone know whether there is a source term in the energy equation? Kind regards, shock77 

January 20, 2020, 17:07 

#2 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
Hello shock77,
The fvOptions in EEqn, regards the potential sources, constraints or corrections you could add to your energy equation, however, the thermal part of the porous media is not described in openFoam. You can do this your self by developing the fvOptions library required, or implement it directly in the solver (this is what I have done previously). Lastly you could do a coded source that does it all, however it requires you develop the entire sink term yourself from a empirical formula or similarly. I have done something alike the last part as well, just a source term that depended on the geometrical location of the cell in the domain, volume of the cell, velocity in the cell and so on. I have attached an example of the code I have used, where a1 through a5 are just numerical constants from a fit. Code:
energySource { type scalarCodedSource; //scalarSemiImplicitSource active true; name sourceTime; scalarCodedSourceCoeffs //scalarSemiImplicitSourceCoeffs S(x) = Su + Sp*x // q in [W]; or in [W/m³] if you use specific mode { selectionMode cellZone; cellZone porousity1; fields (h); fieldNames (h); name sourceTime; codeInclude #{ #}; codeCorrect #{ // Pout<< "**codeCorrect**" << endl; #}; codeAddSup #{ // const Time& time = mesh().time(); const scalarField& V = mesh_.V(); const vectorField& C = mesh_.C(); const volVectorField& U = mesh().lookupObject<volVectorField>("U"); const volScalarField z = U.mesh().C() & vector(0,0,1); scalarField& hSource = eqn.source(); forAll(C, i) { hSource[i] = a1*((a2 a3*exp(a4*mag(a5z[i])))*V[i]); } // Pout << "***codeAddSup***" << endl; #}; codeSetValue #{ // Pout<< "**codeSetValue**" << endl; #}; // Dummy entry. Make dependent on above to trigger recompilation code #{ $codeInclude $codeCorrect $codeAddSup $codeSetValue #}; } sourceTimeCoeffs { // Dummy entry } } Lasse 

January 24, 2020, 09:53 

#3 
Member
Join Date: Jun 2017
Posts: 73
Rep Power: 8 
Hi,
could you explain how you derived your change in enthalpie due to porosity? 

January 24, 2020, 10:23 

#4 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
My approach is somewhat based upon the following source: https://www.slideshare.net/Cypiii/op...trainingv51en
With the heat transfer coefficient being based upon empirical formulas. I believe I found applicable formulas in Chemical reactor analysis and design by Froment, Bischoff. And then only executing the catalyst code if in the porous zone. Let me know if you have any additional questions. 

January 31, 2020, 09:00 

#5 
Senior Member
Join Date: Dec 2019
Posts: 215
Rep Power: 7 
Dear Swagga5aur,
thank you very much for your reply and sry for my very late response. So what you do is treating the heat transfer within the porous media when the gas and the solid have different temperatures right? What I am wondering about is, that I cant find any treatment of the temperature change due to the increased friction inside a porous material. It seems to be always negleted. Maybe I am mistaken and its not important. What I also cant find is how turbulence models are affected by porosity in OpenFOAM. Do you know how it works? Kind regards, shock77 

January 31, 2020, 15:15 

#6 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
Yes, heat transfer requires a difference in temperatures to exist.
Do you suggest that due to the lower cross sectional flow area of the fluid a frictional heat generation occurs? To my understanding frictional heat generation are usually neglectable compared to the convective heat transfer. I have mainly dealt with low velocities, however, perhaps some the of the following links may help. http://scientiairanica.sharif.edu/ar...f69f2ffd37.pdf https://www.cfdonline.com/Wiki/V2f_models In general I would consider the resistance of the porous media outweighs the influence of the turbulent flow. https://vbn.aau.dk/ws/portalfiles/po..._equations.pdf Check page 242 for a description of the critical Reynolds number. Hope this helps. Best regards Lasse 

February 4, 2020, 07:10 

#7  
Senior Member
Join Date: Dec 2019
Posts: 215
Rep Power: 7 
Thanks again for your reply.
Quote:
Quote:
It is a difficult topic, at least for me. I am lacking experience here. But thank you very much for your kind help! It helped me a lot. Kind regards, shock77 

March 7, 2020, 08:00 

#8 
New Member
Kunqing Jiang
Join Date: Mar 2020
Posts: 2
Rep Power: 0 
Dear Swagga5aur and shock77
I have met the same problem,and I benefited a lot from your discussion.But it's hard for me to do this "And then only executing the catalyst code if in the porous zone".Because I haven't learned openfoam long. Could you please tall me how to do that ,thanks. 

March 10, 2020, 09:15 

#9 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
Hello Kunqing,
I currently don't have the time to make this for you, however, just a brief description of what I would do. I would add cells from a porosity zone as the following. Code:
const labelList& cells = mesh.cellZones() [mesh.cellZones().findZoneID("porosity")]; Hope it helps. 

March 14, 2020, 09:33 

#10 
New Member
Kunqing Jiang
Join Date: Mar 2020
Posts: 2
Rep Power: 0 
Dear Swagga5aur,
Thanks for your reply,it helps me a lot. I will have a try myself. Thank you very much indeed. 

March 28, 2020, 15:59 

#11 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
To anyone of interest here is an example of how to make a zone specific Y solver where you can remove for example the reaction term.
Its in no way improved for speed just briefly tested. Code:
tmp<fv::convectionScheme<scalar> > mvConvection ( fv::convectionScheme<scalar>::New ( mesh, fields, phi, mesh.divScheme("div(phi,Yi)") ) ); { reaction.correct(); Qdot = reaction.Qdot(); volScalarField Yt(0.0*Y[0]); const labelList& cellZ = mesh.cellZones() [mesh.cellZones().findZoneID("catalyst")]; const label zoneID = mesh.cellZones().findZoneID("catalyst"); if (zoneID != 1) { forAll(cellZ, i) { zF=1; } forAll(Y, i) { if (Y[i].name() != inertSpecie) { volScalarField& Yi = Y[i]; fvScalarMatrix YiEqn ( mvConvection>fvmDiv(phi, Yi)  fvm::laplacian(turb.muEff(), Yi) == reaction.R(Yi) + fvOptions(rho, Yi) ); YiEqn.relax(); YiEqn.relax(mesh.equationRelaxationFactor("Yi")); fvOptions.constrain(YiEqn); YiEqn.solve(mesh.solver("Yi")); fvOptions.correct(Yi); Yi.max(0.0); Yt += Yi; } } } forAll(Y, i) { if (Y[i].name() != inertSpecie) { volScalarField& Yi = Y[i]; fvScalarMatrix YiEqn ( mvConvection>fvmDiv(phi, Yi)  fvm::laplacian(turb.muEff(), Yi) == fvOptions(rho, Yi) ); YiEqn.relax(); YiEqn.relax(mesh.equationRelaxationFactor("Yi")); fvOptions.constrain(YiEqn); YiEqn.solve(mesh.solver("Yi")); fvOptions.correct(Yi); Yi.max(0.0); Yt += Yi; } } Y[inertIndex] = scalar(1)  Yt; Y[inertIndex].max(0.0); } Just remove the following code to make it compile able. Code:
forAll(cellZ, i) { zF=1; } Note this is made for a custom multiRegion solver but the main idea is the same,. Code:
PtrList<volScalarField> zFReactor(reactorRegions.size()); Code:
Info<< " Adding to zFReactor\n" << endl; zFReactor.set ( i, new volScalarField ( IOobject ( "zF", runTime.timeName(), reactorRegions[i], IOobject::NO_READ, IOobject::AUTO_WRITE ), reactorRegions[i], dimensionedScalar("zF", dimless, scalar(0.0)) ) ); Best regards, Lasse 

May 9, 2020, 15:52 

#12 
Senior Member
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 112
Rep Power: 10 
Hello all,
Turned out this doesn't work, however, the following does but is less elegant regarding skipping unnecessary steps. Added this to setRegionFields: Code:
word cellSetName = "catalyst"; label zoneID = mesh.cellZones().findZoneID(cellSetName); if (zoneID == 1) { FatalErrorIn("yourFunctionName") << "Cannot find cellZone " << cellSetName << endl << "Valid cellZones are " << mesh.cellZones().names() << exit(FatalError); } const labelList& cells = mesh.cellZones()[zoneID]; Info << "Cells in cellzone " << cellSetName << ":" << endl; forAll(cells, i) { const label cell = cells[i]; zF[cell]=1; } Code:
Info<< " Adding to zF\n" << endl; zF.set ( i, new volScalarField ( IOobject ( "zF", runTime.timeName(), reactorRegions[i], IOobject::NO_READ, IOobject::AUTO_WRITE ), Regions[i], dimensionedScalar("zF", dimless, scalar(0.0)) ) ); Best regards, Lasse 

Thread Tools  Search this Thread 
Display Modes  

