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

Adding porosity to UEqn.H in porousInterFoam

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

Like Tree1Likes
  • 1 Post By jameswilson620

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 21, 2014, 13:04
Default Adding porosity to UEqn.H in porousInterFoam
  #1
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12
jameswilson620 is on a distinguished road
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
jameswilson620 is offline   Reply With Quote

Old   August 21, 2014, 22:33
Default
  #2
Senior Member
 
kmooney's Avatar
 
Kyle Mooney
Join Date: Jul 2009
Location: San Francisco, CA USA
Posts: 323
Rep Power: 18
kmooney is on a distinguished road
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
kmooney is offline   Reply With Quote

Old   August 22, 2014, 22:17
Default
  #3
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12
jameswilson620 is on a distinguished road
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
Attached Images
File Type: jpg Screenshot from 2014-08-22 21:11:48.jpg (8.3 KB, 104 views)
jameswilson620 is offline   Reply With Quote

Old   August 23, 2014, 04:35
Default
  #4
ngj
Senior Member
 
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Copenhagen, Denmark
Posts: 1,903
Rep Power: 37
ngj will become famous soon enoughngj will become famous soon enough
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
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request.
ngj is offline   Reply With Quote

Old   August 23, 2014, 10:07
Default
  #5
Member
 
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12
jameswilson620 is on a distinguished road
Thanks Niels, I'll take a look at that asap.

James
ngj likes this.
jameswilson620 is offline   Reply With Quote

Old   November 25, 2016, 23:19
Default
  #6
Member
 
Anirudh Kulkarni
Join Date: May 2016
Posts: 62
Rep Power: 10
Tempest is on a distinguished road
Quote:
Originally Posted by jameswilson620 View Post
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
&nbsp;

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

surfaceScalarField muEff_porosityf
(
"muEff_porosityf",
muEff/porosityf
&nbsp;


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.
Tempest is offline   Reply With Quote

Reply

Tags
divergence, fvvectormatrix, porosity, porousinterfoam, solver compilation error

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


Similar Threads
Thread Thread Starter Forum Replies Last Post
chtMultiRegionSimpleFoam samiam1000 OpenFOAM Running, Solving & CFD 39 March 31, 2016 09:43
Porosity and permeability in PorousInterFoam jasouza1974 OpenFOAM Pre-Processing 14 August 12, 2015 05:16
help adding the energy equation to porousinterfoam using the enthalpy formulation nadine OpenFOAM Programming & Development 18 June 17, 2014 09:39
Help with chtMultiRegionFoam jbvw96 OpenFOAM Running, Solving & CFD 2 December 26, 2010 18:16
Porosity field in Fluent wojciech FLUENT 1 September 20, 2010 12:19


All times are GMT -4. The time now is 01:38.