[lcbuijs]

Dear Foamers,

I'm simulating a small Rayleigh-Benard cell filled with water under Boussinesq conditions. I've already been running successful simulations using buoyantBoussinesqPimpleFoam and now I want to generate results for the same case using buoyantPimpleFoam, in order to implement temperature dependent properties after validation.

The problem I'm currently facing is that solving using buoyantPimpleFoam appears to generate much higher Courant numbers, resulting in a simulation time that is ridiculously long. I have no idea why this is happening, because all my simulation parameters seem to be the same as when using buoyantBoussinesqPimpleFoam. In short, the only changes that I've made are using a different solver, and using polynomial values for the transport coefficients (which I've set to the same constants as before for the time being).

Can you tell me what I'm doing wrong? As you can probably tell, I'm not an expert at finite volume methods. My previous simulations yielded perfectly fine results however.

Many thanks in advance and best regards,
Luuk

fvSolution:

solvers
{
    "rho.*"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       0;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-05;
        relTol          0;
        smoother        DIC;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-5;
    }

    "(U|h|e|k|epsilon|R)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-5;
    }

    "(U|h|e|k|epsilon|R)Final"
    {
        $U;
        tolerance       1e-5;
    }
}

PIMPLE
{
    momentumPredictor yes;
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

relaxationFactors
{
    "(U|T|k|epsilon|R)" 1;
    "(U|T|k|epsilon|R)Final" 1;
}

fvSchemes:

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
    div(phi,h)      Gauss linear;
    div(phi,e)      Gauss linear;
    div(phi,k)      Gauss linear;
    div(phi,epsilon) Gauss linear;
    div(phi,R)      Gauss linear;
    div(phi,K)      Gauss linear;
    div(phi,Ekp)    Gauss linear;
    div(R)          Gauss linear;
    div((muEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         uncorrected;
}

fluxRequired
{
    default         no;
    p_rgh;
}

thermophysicalProperties:

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnthalpy;
}


mixture
{
    // coefficients for water

    specie
    {
        nMoles          1;
        molWeight       18.01528;
    }
    equationOfState
    {
        rhoCoeffs<8>    ( 998 0 0 0 0 0 0 0 );
    }
    thermodynamics
    {
        Hf              0;
        Sf              0;
        CpCoeffs<8>     ( 4180 0 0 0 0 0 0 0 );
    }
    transport
    {
        muCoeffs<8>     ( 8.94e-4 0 0 0 0 0 0 0 );
        kappaCoeffs<8>  ( 0.58 0 0 0 0 0 0 0 );
    }
}

[Svensen]

Try to decrease the relaxation factor for
"(U|T|k|epsilon|R)" 1;
to about 0.3;

[lcbuijs]

Originally Posted by Svensen

Try to decrease the relaxation factor for
"(U|T|k|epsilon|R)" 1;
to about 0.3;

Thank you for your reply! I have given this a try, but it didn't fix the problem.

I have found something strange that could be the cause of the differences though. It seems that buoyantPimpleFoam and buoyantBoussinesqPimpleFoam use different dimensions for the p and p_rgh fields: bBPF uses [0 2 -2 ... ] and bPF uses [1 -1 -2 ... ]. I don't really understand this difference in definition, which makes it hard to compare. Changing the dimensions in bPF to those used in bBPF causes the solver to fail, as can be expected.

I've also noticed that the values of the internal fields for U, p and p_rgh differ greatly between bPF and bBPF after one timestep. bPF yields much higher U values (resulting in a higher Courant number), and the difference in definition of the pressure fields is of course persistent throughout the simulation. Also, the values of p_rgh in bPF vary slightly more than in bBPF.

I'm suspecting that all of these differences are correlated but I'm having a hard time keeping a clear view of what is what. Can anybody tell me what is going on or where I should look to find out?

Again many thanks in advance and best regards!
Luuk

[Svensen]

"bBPF uses [0 2 -2 ... ] and bPF uses [1 -1 -2 ... ]"

In the first case, the [0 2 -2 ...] you need to divide your Pressure values in Pascals by the value of density and you will obtain the correct value.
In the seconds case, you need no calculations and you can use here you Pressure values in Pascals

[lcbuijs]

Originally Posted by Svensen

"bBPF uses [0 2 -2 ... ] and bPF uses [1 -1 -2 ... ]"

In the first case, the [0 2 -2 ...] you need to divide your Pressure values in Pascals by the value of density and you will obtain the correct value.
In the seconds case, you need no calculations and you can use here you Pressure values in Pascals

You are completely right. Luckily I figured this out shortly after my last post. In the meantime, I seem to have figured out where (part of) the problem is. bPF seems to be having trouble with the rhoConst entry in the thermophysicalProperties dict. Setting this to perfectGas fixes the Courant number problem. My explanation is that setting constant density confuses the solver because there is no buoyancy, but that would not directly explain the high velocity values. Can you tell me whether or not this explanation is correct?

Also, if this is true, I would still like to compare with my previous simulations in bBPF. Is perfectGas a good way to compare with the Boussinesq approximation or should I use another entry for equationOfState?

Many thanks!
Luuk

EDIT: It seems I was wrong. Setting equationOfState to perfectGas only fixes the high Courant value for a number of timesteps. After a while the values still increase and the timesteps converge towards a floating point exception... I'm still thinking the problem lies somewhere within the density settings. Any clues?

[lcbuijs]

Hi again guys,

I've printed out the density values as a scalar field and everything seems to be fine. I was expecting that uniform density would imply no buoyancy and therefore a uniform zero velocity field. The solver solves for the gradient in p_rgh though, which is of course nonzero. The velocity values I end up with are quite high as I said before. Can anyone tell me if my reasoning is incorrect and why?

Many thanks!