CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Porosity and the energy equation

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

Like Tree1Likes
  • 1 Post By shock77

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 20, 2020, 11:39
Default Porosity and the energy equation
  #1
Senior Member
 
Join Date: Dec 2019
Posts: 215
Rep Power: 7
shock77 is on a distinguished road
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"
            }
        }
     }
}
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:


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)
    );
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
alainislas likes this.
shock77 is offline   Reply With Quote

Old   January 20, 2020, 18:07
Post
  #2
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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(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
Swagga5aur is offline   Reply With Quote

Old   January 24, 2020, 10:53
Default
  #3
Member
 
Join Date: Jun 2017
Posts: 73
Rep Power: 9
Friendly is on a distinguished road
Hi,


could you explain how you derived your change in enthalpie due to porosity?
Friendly is offline   Reply With Quote

Old   January 24, 2020, 11:23
Default
  #4
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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.
Swagga5aur is offline   Reply With Quote

Old   January 31, 2020, 10:00
Default
  #5
Senior Member
 
Join Date: Dec 2019
Posts: 215
Rep Power: 7
shock77 is on a distinguished road
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
shock77 is offline   Reply With Quote

Old   January 31, 2020, 16:15
Default
  #6
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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
Swagga5aur is offline   Reply With Quote

Old   February 4, 2020, 08:10
Default
  #7
Senior Member
 
Join Date: Dec 2019
Posts: 215
Rep Power: 7
shock77 is on a distinguished road
Thanks again for your reply.


Quote:
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.


Quote:
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
shock77 is offline   Reply With Quote

Old   March 7, 2020, 09:00
Default
  #8
New Member
 
Kunqing Jiang
Join Date: Mar 2020
Posts: 2
Rep Power: 0
Nikesilpers is on a distinguished road
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.
Nikesilpers is offline   Reply With Quote

Old   March 10, 2020, 10:15
Default
  #9
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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")];
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.
Swagga5aur is offline   Reply With Quote

Old   March 14, 2020, 10:33
Default
  #10
New Member
 
Kunqing Jiang
Join Date: Mar 2020
Posts: 2
Rep Power: 0
Nikesilpers is on a distinguished road
Dear Swagga5aur,

Thanks for your reply,it helps me a lot. I will have a try myself. Thank you very much indeed.
Nikesilpers is offline   Reply With Quote

Old   March 28, 2020, 16:59
Default
  #11
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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;
		}
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,.
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))
        )
    );
Let me know if there is any questions or issues.

Best regards,
Lasse
Swagga5aur is offline   Reply With Quote

Old   May 9, 2020, 16:52
Default
  #12
Senior Member
 
Lasse Brams Vinther
Join Date: Oct 2015
Posts: 118
Rep Power: 11
Swagga5aur is on a distinguished road
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;
    }
The following to createFields:
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))
        )
    );
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
Swagga5aur is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 14:55.