CFD Online Discussion Forums - Adding porosity to UEqn.H in porousInterFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Adding porosity to UEqn.H in porousInterFoam (https://www.cfd-online.com/Forums/openfoam-programming-development/140725-adding-porosity-ueqn-h-porousinterfoam.html)

Adding porosity to UEqn.H in porousInterFoam

Hi all,

I would like to modify UEqn.H in porousInterFoam to incorporate porosity in the solver.

I defined porosity by linking it to a zone called "porous" in blockMeshDict and set the field using setFieldsDict and zoneToCell.

I have defined porosity and porosityf in createFields.H as follows:
/////////////////createFields.H////////////////////////////////
...
// Read porosity from /0
volScalarField porosity
(
IOobject
(
"porosity",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

// Converting porosity to surfaceScalarField for operation compatability in UEqn.
surfaceScalarField porosityf
(
IOobject
(
"porosityf",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(porosity)
);
...
/////////////////createFields.H////////////////////////////////

and now in UEqn.H I update the values every time step:

/////////////////UEqn.H////////////////////////////////
// Calculate and cache mu for the porous media
volScalarField mu(twoPhaseProperties.mu());

volScalarField rho_porosity
(
"rho_porosity",
rho/porosity
);

surfaceScalarField rhoPhi_porosityf
(
"rhoPhi_porosityf",
rhoPhi/(porosityf*porosityf)
);

surfaceScalarField muEff_porosityf
(
"muEff_porosityf",
muEff/porosityf
);

fvVectorMatrix UEqn
(
//pZones.ddt(rho, U)
fvm::ddt(rho_porosity, U)
+ fvm::div(rhoPhi_porosityf, U)
- fvm::laplacian(muEff_porosityf, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf()))
==
fvOptions(rho, U)
);
...
/////////////////UEqn.H////////////////////////////////

This compiles fine and the solver runs but diverges shortly into the simulation as I am running a case identical (with the exception of the additional porosity terms) to:

http://www.cfdandco.com/resources/JR...sInter_v11.pdf

where the liquid is adjacent to the interface at t=0.

Without the porosity terms (unmodified solver), this case runs fine. with the porosity terms (1 in free flow and 0.32 in porous zone) the solution diverges.

My attempts to debug have led me to some good information that I do not know what to do with : )

When I make the porosity field uniform (both the volScalarField and the surfaceScalarField for porosity) with a value of 1 (free flow), the solution still diverges. This case is essentially dividing the modified terms by 1 and should have no effect on the solution, which is strange that it still diverges.

something/1 = something, but the divergence here seems like something/1 = something_else. my example is crude and probably doesn't represent whats actually happening.

Clearly there is something happening as a result of my modifications that I am too novice to realize.

I have experienced success in eliminating diverging temperature fields in an unrelated case by streamlining the arguments of operators such as fvm::div(rhoCpPhi, T) by shifting to rhoCp*fvm::div(phi, T). This is also something I'd like to implement here with UEqn.H

for example:

fvVectorMatrix UEqn
(
//pZones.ddt(rho, U)
fvm::ddt(rho, U) / porosity
+ fvm::div(rhoPhi, U) /(porosityf*porosityf)
- fvm::laplacian(muEff, U) / porosityf
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf()))
==
fvOptions(rho, U)
);

But i get a compile error when I try to compile porousInterFoam with the only modification being "fvm::ddt(rho, U) / porosity":

UEqn.H:33:28: error: no match for ‘operator/’ in ‘Foam::fvm::ddt(const volScalarField&, const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::Vector<double>, Foam::volScalarField = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]((*(const Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>*)(& U))) / porosity’

I would appreciate any insight regarding operations within fvVectorMatrix and the effects of adding arguments to something like fvm::div() vs. moving them outside of the divergence operator. Also I am curious as to why I get the compile error when operating on this term fvm::div()*something/something_else and why I don't when I define the field first and then pass it to the fvm::div() operator.

My over arching concern is naturally the divergence of the solution, but also, why dividing by a uniform field of 1 causes this divergence.

Thanks in advance, James

PS I'm new to the forums here and have benefited immensely from simply reading related problems! Thanks all

Hi James,

You might be better off replicating existing implementations of porosity models. For example, in foam 2.3 there is an extra term (fvOptions(U)) added onto the momentum equation in many of the solvers like this (from pimpleFoam):

Code:

tmp<fvVectorMatrix> UEqn

(

    fvm::ddt(U)

  + fvm::div(phi, U)

  + turbulence->divDevReff(U)

 ==

    fvOptions(U)

);

This lets you apply momentum sinks at runtime via the fvOptions library. An example of this would be put in a file '/system/fvOptions' like this:

Code:

porosity1

{

    type            explicitPorositySource;

    active          yes;

    selectionMode   cellZone;

    cellZone        stator;



    explicitPorositySourceCoeffs

    {

        type            DarcyForchheimer;



        DarcyForchheimerCoeffs

        {

            d   d [0 -2 0 0 0 0 0] (1e5 1000 1000);

            f   f [0 -1 0 0 0 0 0] (0 0 0);



            coordinateSystem

            {

                type    cartesian;

                origin  (0 0 0);

                coordinateRotation

                {

                    type    axesRotation;

                    e1  (1 0 0);

                    e2  (0 1 0);

                }

            }

        }

    }

}

You might have to play with the selection mode in order to set the correct cells to porosity 'on'.

I hope that helps you some!
Cheers,
Kyle

Thanks Kyle,

This is some great info! I looked into some models that incorporate porosity and they seem to modify the transient term only. Im set on modifying UEqn.H as the terms i need to modify are in this header. See the attached photo. Pi is already a feature in porousInterFoam but porosity is not. I would think, just as I have built custom functions and equations for something like the heat equation, I could just modify the terms in UEqn.H but the solver doesnt seem to like the modifications, even if it involves dividing by 1.

Ill keep looking and continue to explore fvoptions as well.

Thanks again!

Any more suggestions?

James

Hallo James,

You might want to look into the waves2Foam package. It comes with a porous solver, which looks a lot like the set of equations sketched above. One of the tutorials is even identical with the one you are trying (though other dimensions), but it comes with experimental data and a post-processing script for the validation.

Kind regards,

Niels

Thanks Niels, I'll take a look at that asap.

James

Quote:

Originally Posted by jameswilson620 (Post 507104)

Hi all,

I would like to modify UEqn.H in porousInterFoam to incorporate porosity in the solver.

I defined porosity by linking it to a zone called "porous" in blockMeshDict and set the field using setFieldsDict and zoneToCell.

I have defined porosity and porosityf in createFields.H as follows:
/////////////////createFields.H////////////////////////////////
...
// Read porosity from /0
volScalarField porosity
(
IOobject
(
"porosity",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

// Converting porosity to surfaceScalarField for operation compatability in UEqn.
surfaceScalarField porosityf
(
IOobject
(
"porosityf",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(porosity)
);
...
/////////////////createFields.H////////////////////////////////

and now in UEqn.H I update the values every time step:

/////////////////UEqn.H////////////////////////////////
// Calculate and cache mu for the porous media
volScalarField mu(twoPhaseProperties.mu());

volScalarField rho_porosity
(
"rho_porosity",
rho/porosity
 );

surfaceScalarField rhoPhi_porosityf
(
"rhoPhi_porosityf",
rhoPhi/(porosityf*porosityf)
 );

surfaceScalarField muEff_porosityf
(
"muEff_porosityf",
muEff/porosityf
 );

fvVectorMatrix UEqn
(
//pZones.ddt(rho, U)
fvm::ddt(rho_porosity, U)
+ fvm::div(rhoPhi_porosityf, U)
- fvm::laplacian(muEff_porosityf, U)
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf()))
==
fvOptions(rho, U)
);
...
/////////////////UEqn.H////////////////////////////////

This compiles fine and the solver runs but diverges shortly into the simulation as I am running a case identical (with the exception of the additional porosity terms) to:

http://www.cfdandco.com/resources/JR...sInter_v11.pdf

where the liquid is adjacent to the interface at t=0.

Without the porosity terms (unmodified solver), this case runs fine. with the porosity terms (1 in free flow and 0.32 in porous zone) the solution diverges.

My attempts to debug have led me to some good information that I do not know what to do with : )

When I make the porosity field uniform (both the volScalarField and the surfaceScalarField for porosity) with a value of 1 (free flow), the solution still diverges. This case is essentially dividing the modified terms by 1 and should have no effect on the solution, which is strange that it still diverges.

something/1 = something, but the divergence here seems like something/1 = something_else. my example is crude and probably doesn't represent whats actually happening.

Clearly there is something happening as a result of my modifications that I am too novice to realize.

I have experienced success in eliminating diverging temperature fields in an unrelated case by streamlining the arguments of operators such as fvm::div(rhoCpPhi, T) by shifting to rhoCp*fvm::div(phi, T). This is also something I'd like to implement here with UEqn.H

for example:

fvVectorMatrix UEqn
(
//pZones.ddt(rho, U)
fvm::ddt(rho, U) / porosity
+ fvm::div(rhoPhi, U) /(porosityf*porosityf)
- fvm::laplacian(muEff, U) / porosityf
- (fvc::grad(U) & fvc::grad(muEff))
//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf()))
==
fvOptions(rho, U)
);

But i get a compile error when I try to compile porousInterFoam with the only modification being "fvm::ddt(rho, U) / porosity":

UEqn.H:33:28: error: no match for ‘operator/’ in ‘Foam::fvm::ddt(const volScalarField&, const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::Vector<double>, Foam::volScalarField = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]((*(const Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>*)(& U))) / porosity’

I would appreciate any insight regarding operations within fvVectorMatrix and the effects of adding arguments to something like fvm::div() vs. moving them outside of the divergence operator. Also I am curious as to why I get the compile error when operating on this term fvm::div()*something/something_else and why I don't when I define the field first and then pass it to the fvm::div() operator.

My over arching concern is naturally the divergence of the solution, but also, why dividing by a uniform field of 1 causes this divergence.

Thanks in advance, James

PS I'm new to the forums here and have benefited immensely from simply reading related problems! Thanks all

I know I am trying to dig into an old thread, but did you solve the problem? I am working on the same issue.