[shock77]

Hi,


I am trying to run a simulation with openfoam4 and rhoCentralFoam. I am trying to model the supersonic outflow from a pressure vessel and I am interested in transient effects.


I have attached an image of the residuals. They rise slighty and then basically explode. I have tracked the min and max values of U, p and T during runtime. Just before the crash, the temperature drops close to 0 in some areas and rises to 1e200 in other parts. I have tried to limit the temperature but it didnt work (I dont know why):

temperatureLimit
{
    type            limitTemperature;
    active          true;

        limitTemperatureCoeffs
        {
            selectionMode   all;
            Tmin     10;
            Tmax     10000;
        }
}

It works if I run the simulation first order, but I want it to run second order for more accuracy. Here is my fvSchemes dict:

fluxScheme      Kurganov;

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         cellLimited Gauss linear 1;
}

divSchemes
{
    default         Gauss linear;
    div(tauMC)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
    reconstruct(rho) vanAlbada;
    reconstruct(U)   vanAlbadaV;
    reconstruct(T)   vanAlbada;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

I also tried CrankNicolson with varying coefficient and tried different Courant numbers (0.1 - 0.5). Nothing helped so far.


My Boundary Conditions are:


Inlet: Outlet:



p - uniformTotalPressure; p - waveTransmissive;

T - uniformTotalTemperature; T - inletOutlet;

U - pressureInletOutletVelocity; U - inletOutlet;



Walls are all slip.


Maybe someone has some good advice.

[shock77]

I have noticed, that fvOptions is not included in rhoCentralFoam.


So I have added a few lines into a custom solver and included fvOption.H and createFvOptions.H:

if (!inviscid)
        {
    fvScalarMatrix EEqnRCF
    (
     fvm::ddt(rho, e) - fvc::ddt(rho, e)
        - fvm::laplacian(turbulence->alphaEff(), e) ==
        fvOptions(rho, e)
    );

    fvOptions.constrain(EEqnRCF);
        EEqnRCF.solve();

        fvOptions.correct(e);
            thermo.correct();
            rhoE = rho*(e + 0.5*magSqr(U));
        }

limitTemperature seems to work now.