[yanxiang]

I tested porousInterFoam and interFoam in an empty tube with exactly the same settings, but I got different results. From my understanding, they should behave the same when simulating an empty duct.

What is also interesting is interFoam takes "div((muEff*dev(T(grad(U)))))" while porousInterFoam uses "div((nuEff*dev(T(grad(U)))))". I am running OF222

[wyldckat]

Greetings yanxiang,
The difference in using "mu" and "nu" usually implies that one is using an incompressible model and the other uses compressible. Compare pimpleFoam with rhoPimpleFoam for a clearer idea on the differences and namely regarding which refers to compressible and which refers to incompressible .

As for the results, the initial conditions might be the reason why you're getting such a big difference in results, since one handles fluid as incompressible and the other as compressible. If you provide the two cases, it will make it a lot easier to diagnose the specific reason why you're getting those results. Have a look at this blog post of mine, for more ideas: OpenFOAM: Interesting cases of bad meshes and bad initial conditions

Best regards,
Bruno

[yanxiang]

Dear Bruno,

Thanks for your reply. Indeed, "mu" and "nu" are usually used in the context of compressible and incompressible flows respectively. However, in my case, I use porousInterFoam, which should not be different from interFoam in that regard. Since it's simply an extension of the interFoam with the porous media model, I expect it reduces to the original interFoam when there is no resistance specified in the porosityProperties dictionary.

I am, however, using a dangerous inlet boundary condition for the alpha1 field, which varies periodically with time. This could be the reason why it causes some numerical issues. I will test with a safer boundary condition and see whether they behave the same.

Best,
yanxiang

[wyldckat]

Hi yanxiang,

I took a better look at the cases and solvers in OpenFOAM 2.2.2 and I'm not finding any references to "nuEff" in the tutorials related to porousInterFoam and not even for most of the multiphase tutorials!
And I took a look at the source code for this solver and the line that does use "div(muEff)" is commented out, which means that not even that is being used.
Nonetheless, "muEff" is indeed used for this solver: https://github.com/OpenFOAM/OpenFOAM...terFoam/UEqn.H

Perhaps you're mixing your tests with tutorial cases from past OpenFOAM versions?

Best regards,
Bruno

[yanxiang]

Hi Bruno,

I double checked the two cases, and porousInterFoam doesn't use div(muEff) nor div(nuEff). So no numerical scheme needs to be specified for that.

I also looked into the source code a bit more carefully and noticed that the implementation of UEqn in porousInterFoam is different from that in interFoam not only in the addition of the resistance, but also in the following lines:

interFoam:

   
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(rho, U)
    );

porousInterFoam:

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

Could this subtle difference the reason why these two solvers behave differently when no porous section is specified?

Thanks,
yanxiang

[alexeym]

Hi,

if you take a look at the implementation of turbulence->divDevRhoReff(rho, U) (for example in kEpsilon.H):

tmp<fvVectorMatrix> kEpsilon::divDevRhoReff
(
    const volScalarField& rho,
    volVectorField& U
) const
{
    volScalarField muEff("muEff", rho*nuEff());

    return
    (
      - fvm::laplacian(muEff, U)
      - fvc::div(muEff*dev(T(fvc::grad(U))))
    );
}

so it's more-or-less the same as in porousInterFoam.

Main difference is in two lines:

    fvVectorMatrix UEqn
    (
        ...
     ==
        fvOptions(rho, U)
    );

and

    pZones.addResistance(UEqn);

With fvOptions you can include additional sources (even porosity), and with pZones you define porosity also.

Concerning mu and nu: mu usually used when UEqn has term:

fvm::ddt(rho, U)

nu usually used when UEqn has term:

fvm::ddt(U)

i.e. then UEqn is divided by density.

Concerning your initial question: can you, please, post two archives with your case files? (So people can reproduce the behavior of the solvers, take a closer look at the settings etc.)

[yanxiang]

Hi Alexey,

Thanks a lot for your reply.

I did track down the definition of divDevRhoReff, but could you explain how those last two terms are the same?

fvc::div(muEff*dev(T(fvc::grad(U))))

vs

(fvc::grad(U) & fvc::grad(muEff))

I have attached the two cases: tube_iF is with interFoam and the other with porousInterFoam without the porous zone specified. Please note again that I have used a dangerous boundary condition for the alpha field. It will be interesting to see where the difference is from.

Thanks,
yanxiang

[alexeym]

Hi,

concerning results: till 0.10 s they are the same, after due to BC on alpha solution diverges:

...
MULES: Solving for alpha1
Phase-1 volume fraction = 0.025671  Min(alpha1) = 0  Max(alpha1) = 2.11241
MULES: Solving for alpha1
Phase-1 volume fraction = 0.025671  Min(alpha1) = 0  Max(alpha1) = 2.11232
MULES: Solving for alpha1
Phase-1 volume fraction = 0.025671  Min(alpha1) = 0  Max(alpha1) = 2.11223
MULES: Solving for alpha1
Phase-1 volume fraction = 0.025671  Min(alpha1) = 0  Max(alpha1) = 2.11214
...

Not quite sure codedFixedValue worth using here as BC looks more like uniformFixedValue (it appears in 2.1.0 - http://www.openfoam.org/version2.1.0...conditions.php)

[yanxiang]

Hi Alexey,

Thanks for the tests. Indeed, they diverge after a certain point, and it is likely due to the BC.

Also, thanks for pointing me to the BC's.

Cheers,
yanxiang