[shock77]

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:

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"
            }
        }
     }
}

It works fine concerning the pressure drop. The darcy forchheimer modell is an additional source term in the momentum equation.


But the problem is, that the temperature seems not to be affected by the porosity. I have looked up EEqn.H for sonicFoam:

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)
    );

It seems the solver considers fvOptions in the energy equation too. Unfortunately I cannot find how or whether there is an additional source term in the energy equation aswell.


Does anyone know whether there is a source term in the energy equation?




Kind regards,
shock77

[Swagga5aur]

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.

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(a5-z[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
    }
}

Best regards,

Lasse

[Friendly]

Hi,


could you explain how you derived your change in enthalpie due to porosity?

[Swagga5aur]

My approach is somewhat based upon the following source: https://www.slideshare.net/Cypiii/op...training-v51en

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.

[shock77]

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

[Swagga5aur]

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 neglect-able 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.cfd-online.com/Wiki/V2-f_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

[shock77]

Thanks again for your reply.

Do you suggest that due to the lower cross sectional flow area of the fluid a frictional heat generation occurs?

Yes excactly. I have a high velocity flow and I am not sure whether the friction plays a role or not. I have read lots of papers and every author neglects the additional friction. Maybe you are right and its not very important.

In general I would consider the resistance of the porous media outweighs the influence of the turbulent flow.

I have read both. Either you neglect the turbulence, whenever you have a very fine porosity or you have to implement it into your turbulence model.


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

[Nikesilpers]

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.

[Swagga5aur]

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.

    const labelList& cells = mesh.cellZones()
        [mesh.cellZones().findZoneID("porosity")];

And then add a for loop or while loop for the cells where the normal y equation is solved and then everywhere else the y equation is solved with removed reactions such that no reactions occur, but the mixing of the specie still does.

Hope it helps.

[Nikesilpers]

Dear Swagga5aur,

Thanks for your reply,it helps me a lot. I will have a try myself. Thank you very much indeed.

[Swagga5aur]

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.

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.

 		forAll(cellZ, i)
		{
		   zF=1;
		}

If interested for a simple indicator of the functionality, leave this code and add the following to createFields directory to make a field zF being 1 if its the specified zone or 0 if outside it.

Note this is made for a custom multiRegion solver but the main idea is the same,.

PtrList<volScalarField> zFReactor(reactorRegions.size());

    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))
        )
    );

Let me know if there is any questions or issues.

Best regards,
Lasse

[Swagga5aur]

Hello all,
Turned out this doesn't work, however, the following does but is less elegant regarding skipping unnecessary steps.

Added this to setRegionFields:

    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;
    }

The following to createFields:

   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))
        )
    );

and then just uses as a multiplier to disable fx. reaction.R(Yi) call in YEqn to remove reactions in cells outside of the catalyst cellzone.

Best regards,
Lasse