[ignacio]

Hello Foamers,

I am running a T-Junction case (both water and pipes), using the chtMultiRegionFoam solver.
First I run a steady state and now I am doing a URANS with the steady state output as initial value.
I use k-OmegaSST Low Reynolds model.

I don't know why but the simulation always crashes after 3 iterations... I changed parameters in fvSolution and fvSchemes, but still, it's always after 3 iterations.
My time step is so low that the initial Courant is only 0.1

The problem is at the outflow pipe. Inflow is in another place with much higher velocity but it's working fine. I attach some pictures
No problems at all in checkMesh.

Any idea of how to solve this?
I would really appreciate any kind of comments.


Mesh.jpg

Courant.jpg

[ignacio]

fvSolution

solvers
{
    rho
    {
        solver          PCG
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0;
    }
    rhoFinal
    {
        $rho;
        tolerance       1e-7;
        relTol          0;
    }

    p_rgh
    {        
    solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;

        smoother         DICGaussSeidel;
    nPreSweeps    0;
    nPostSweeps    2;
    nFinestSweeps    2;

        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;
    }
    p_rghFinal
    {
        $p_rgh;
        tolerance        1e-7;
        relTol           0;
    }
    h
    {
        solver           GAMG;
        tolerance        1e-10;
        relTol           0.01;

    smoother    GaussSeidel;
    nPreSweeps    0;
    nPostSweeps    2;
    nFinestSweeps    2;

        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;

    }
    hFinal
    {
        $h;
        tolerance        1e-10;
        relTol           0;
    }
    "(U|k|epsilon|omega)"
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }
    "(U|k|omega|R)Final"
    {
        $U;
        tolerance        1e-7;
        relTol           0;
    }
}

PIMPLE
{
    momentumPredictor        on;    //Accounts for variations by reconstructing the flux
    nCorrectors            3;    //How many times the pressure equation is solved in outer loop
    nOuterCorrectors        2;      //Outer Piso loop
    nNonOrthogonalCorrectors     3;

    residualControl
    {
    h        1e-3;
    }
}
relaxationFactors
{
    fields
    {
        rho             1;
        p_rgh           0.5;
    }
    equations
    {
        U               0.3;
        h               0.5;
        nuTilda         0.5;
        k               0.5;
        epsilon         0.5;
        omega           0.5;
        "ILambda.*"     0.5;
    }
}

fvSchemes

ddtSchemes
{
    default     backward;
}

gradSchemes
{
    default              Gauss linear;
    grad(p)        cellLimited Gauss linear 1;
    grad(U)        cellLimited Gauss linear 1;    
}

divSchemes
{
    default              Gauss limitedLinear 1;
    div(phi,U)           Gauss filteredLinear2V 0.2 0;

    div((muEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default          Gauss linear limited 1;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 0.333;
}

fluxRequired
{
    default         no;
    p_rgh;
}

Boundary Conditions

U
    {
        internalField   uniform (0 0 0);

        boundaryField
        {
            ".*"
            {
                type        fixedValue;
                value        uniform (0 0 0);
            }
            COLDINLETMAXX
            {
                type            timeVaryingMappedFixedValue;
            setAverage      off;
        offset        (0 0 0);
            }
            COLDINLETMINX
            {
                 type           timeVaryingMappedFixedValue;
             setAverage       off;
         offset        (0 0 0);
              }
            HOTINLETMAXZ
            {
                type            timeVaryingMappedFixedValue;
            setAverage      off;
        offset        (0 0 0);
            }
            HOTINLETMINZ
             {
                type            timeVaryingMappedFixedValue;
            setAverage      off;
        offset        (0 0 0);
             }
            OUTLETMAXZ
             {
            type            inletOutlet;
        inletValue    uniform (0 0 0);
        value        uniform (0 0 0);
            }
            OUTLETMINZ
            {
            type            inletOutlet;
        inletValue    uniform (0 0 0);
        value        uniform (0 0 0);
            }
        }
    }

    T
    {
        internalField   uniform 333.15;

        boundaryField
        {
        "WATER_to_.*"
            {
                type                compressible::turbulentTemperatureCoupledBaffleMixed;
                neighbourFieldName  T;
                kappa               fluidThermo;
                kappaName           none;
                value               uniform 333.15;
            }
        HIHGHREFLUX
            {
                type                compressible::turbulentTemperatureCoupledBaffleMixed;
                neighbourFieldName  T;
                kappa               fluidThermo;
                kappaName           none;
                value               uniform 333.15;
            }
            COLDINLETMAXX
            {
            type            fixedValue;
        value        uniform 333.15;
            }
            COLDINLETMINX
            {
            type            fixedValue;
        value        uniform 333.15;
              }
            HOTINLETMAXZ
            {
            type            fixedValue;
        value        uniform 549.15;
            }
            HOTINLETMINZ
             {
            type            fixedValue;
        value        uniform 549.15;
             }
            OUTLETMAXZ
             {
            type            inletOutlet;
        inletValue    uniform 333.15;
        value        uniform 333.15;
            }
            OUTLETMINZ
            {
            type            inletOutlet;
        inletValue    uniform 333.15;
        value        uniform 333.15;
            }
        TOP
            {
                type            zeroGradient;
            }
        BOTTOM
            {
                type            zeroGradient;
            }
        OUTERRAD
            {
                type            zeroGradient;
            }
        HIGHREADIAB
            {
                type            zeroGradient;
            }
        }
    }

    p_rgh
    {
        internalField   uniform 72e5;

        boundaryField
        {
            ".*"
            {
                type            buoyantPressure;
                value           uniform 72e5;
            }
            COLDINLETMAXX
            {
            type            zeroGradient;
            }
            COLDINLETMINX
            {
            type            zeroGradient;
              }
            HOTINLETMAXZ
            {
            type            zeroGradient;
            }
            HOTINLETMINZ
             {
            type            zeroGradient;
             }
            OUTLETMAXZ
             {
            type            fixedValue;
        value        uniform 72e5;
            }
            OUTLETMINZ
            {
            type            fixedValue;
        value        uniform 72e5;
            }      
        }
    }

    p
    {
        internalField   uniform 72e5;

        boundaryField
        {
            ".*"
            {
                type            calculated;
        value        uniform 72e5;
            }
        }
    }
    //------------------------------------------------
    //   Turbulent Low Reynolds k-OmegaSST model
    //------------------------------------------------

    omega
    {
        internalField   uniform 40;

        boundaryField
        {
            ".*"
            {
                type            compressible::omegaWallFunction; 
                value           uniform 1000;
            }
            COLDINLETMAXX
            {
            type            fixedValue;
        value        uniform 40;
            }
            COLDINLETMINX
            {
            type            fixedValue;
        value        uniform 40;
              }
            HOTINLETMAXZ
            {
            type            fixedValue;
        value        uniform 40;
            }
            HOTINLETMINZ
             {
            type            fixedValue;
        value        uniform 40;
             }
            OUTLETMAXZ
             {
            type            inletOutlet;
                inletValue      uniform 40;
            value           uniform 40;
            }
            OUTLETMINZ
            {
            type            inletOutlet;
                inletValue      uniform 40;
            value           uniform 40;
            }
        }
    }

    k
    {
        internalField   uniform 0.0005;

        boundaryField
        {
            ".*"
            {
                type            fixedValue;
        value        uniform 1e-10;
            }
        HIGHREADIAB
            {
                type            compressible::kqRWallFunction;
        value        uniform 1e-10;
            }
        HIGHREFLUX
            {
                type            compressible::kqRWallFunction;
        value        uniform 1e-10;
            }
            COLDINLETMAXX
            {
            type            fixedValue;
        value        uniform 0.0005;
            }
            COLDINLETMINX
            {
            type            fixedValue;
        value        uniform 0.0005;
              }
            HOTINLETMAXZ
            {
            type            fixedValue;
        value        uniform 0.0005;
            }
            HOTINLETMINZ
             {
            type            fixedValue;
        value        uniform 0.0005;
             }
            OUTLETMAXZ
             {
            type            inletOutlet;
                inletValue      uniform 0.0005;
            value           uniform 0.0005;
            }
            OUTLETMINZ
            {
            type            inletOutlet;
                inletValue      uniform 0.0005;
            value           uniform 0.0005;
            }
        }
    }

    mut 
    {
        internalfield   uniform 0;

        boundaryField
        {
            ".*"
            {
                type         mutLowReWallFunction;
        value        uniform 0;
            }
        HIGHREADIAB
            {
                type            mutkWallFunction;
        value        uniform 1e-10;
            }
        HIGHREFLUX
            {
                type            mutkWallFunction;
        value        uniform 1e-10;
            }
            COLDINLETMAXX
            {
            type            calculated;
            }
            COLDINLETMINX
            {
            type            calculated;
              }
            HOTINLETMAXZ
            {
            type            calculated;
            }
            HOTINLETMINZ
             {
            type            calculated;
             }
            OUTLETMAXZ
             {
            type            calculated;;
            }
            OUTLETMINZ
            {
            type            calculated;
            }
    }
    alphat 
    {
        internalfield   uniform 0;

        boundaryField
        {
            ".*"
            {
                type         alphatWallFunction;
        value        uniform 0;
            }
            COLDINLETMAXX
            {
            type            calculated;
        value        uniform 0;
            }
            COLDINLETMINX
            {
            type            calculated;
        value        uniform 0;
              }
            HOTINLETMAXZ
            {
            type            calculated;
        value        uniform 0;
            }
            HOTINLETMINZ
             {
            type            calculated;
        value        uniform 0;
             }
            OUTLETMAXZ
             {
            type            calculated;
        value        uniform 0;
            }
            OUTLETMINZ
            {
            type            calculated;
        value        uniform 0;
            }    
    }

[Lieven]

A few things you can easily try:
1. Set the time step 10x smaller such that courant number < 1.0
2. Take upwind-schemes to increase stability. It might not be interestiging from numerical point of view, but at this moment the priority is not obtaining 'a' result and not the 'best' result. Later-on you can switch to 2nd order schemes with the upwind-solution as initial condition.
3. First run a few timesteps with potentialFoam to improve the initial conditions

Cheers,

L

[ignacio]

Thanks for the response Lieven!

1 and 2 I have already tried them. But I don't see any improvement in the simulation.
I will see what happens with the potentialFoam

I found out a strange thing playing around with the case.
If I change the density at thermphysicalProperties file it's stable! (the commented part is what I changed) Although the min/max temperature start decreasing/increasing with the iteration

thermoType 
{
    type         heRhoThermo;      // Thermophysical model based on sensible enthalpy
    mixture        pureMixture;      // Only water
    transport    polynomial;       // (mu, k) T dependant transport properties. T in Kelvin
    thermo        hPolynomial;      // (Cp)    T dependant thermodynamic properties. T in Kelvin
    equationOfState    icoPolynomial;    // (rho)   Incompressible polynomial eq of state
    specie        specie;
    energy        sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       18;   // Water
    }
    equationOfState
    {
        rhoCoeffs<8>
    (863.76 0 0 0 0 0 0 0);
    //(856.76 1.2436 -0.0026 0 0 0 0 0);     //Density [kg/m3] = 856.76 - 1.2436*T - 0.0026*T^2
    }

Do you know what effect can this change have in stablility?

Thanks a lot

[Lieven]

Ok, best thing to do then is first let the simulation run with
(863.76 0 0 0 0 0 0 0)
next with
(863.76 1.2436 0 0 0 0 0 0)
and finally with
(863.76 1.2436 -0.0026 0 0 0 0 0)

So always start from the previous simulation as initial condition. If the simulation gets unstable, try to figure out why.
.
By the way, you write 1.2436 but according to your correlation it should be -1.2436

[ignacio]

Oh, mistake. The code is correct, the commented equation isn't.

We are now doing the process that you said. I will tell you how it goes.

Thanks again Lieven!