[rpachaly]

Hello,

I'm trying to simulate a system with several inflows and one outflow that changes the flow rate during the simulation. I want to know the pressurization in three points when the outflow is reduced (see image attached and the U file). There is also some connections to the atmosphere where I am not interested in pressurization. The image attached shows the geometry of the system (don't judge the mesh, it's just an old mesh for demonstration).

I usually get two errors when I'm trying to solve this system: negative temperature and printStack error. What I've tried and not worked:

• Limit Temperature
• Increase the computational domain
• Refine the mesh
• Change fvSchemes and fvSolution

When I consider the atm locations as wall instead of atmosphere BC's, the model runs. But I really don't want this extra pressurization at the points where I have atm (see image).

I want to use the kOmegaSST as turbulence model in a near future. But at this moment I'm turning off the turbulence model.

My BC's for p, p_rgh, T, and U are:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * *	 * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 100000;

boundaryField
{
    walls
    {
        type            calculated;
		value $internalField;
    }

	atm
	{
        type            calculated;
		value $internalField;
	}

    inletTop
    {
        type            calculated;
		value $internalField;
    }

    inletBottom
    {
        type            calculated;
		value $internalField;
    }
	
	inletPipe2
    {
        type            calculated;
		value $internalField;
    }
	
	inletPipe3
    {
        type            calculated;
		value $internalField;
    }
	
	inletPipe1
    {
        type            calculated;
		value $internalField;
    }

    outlet
    {
        type            calculated;
		value $internalField;
    }

}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 1e5;

boundaryField
{

	outlet
    {
        type            zeroGradient;
    }
	
	
    ".*"
    {
        type            fixedFluxPressure;
        value           uniform 1e5;
    }
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 300;

boundaryField
{
    walls
    {
        type            zeroGradient;
    }
	
    atm
    {
        type            inletOutlet;
		inletValue uniform 300;
    }

    inletTop
    {
        type            zeroGradient;
    }

    inletBottom
    {
        type            zeroGradient;
    }
	
	inletPipe1
    {
        type            zeroGradient;
    }

    inletPipe2
    {
        type            zeroGradient;
    }

    inletPipe3
    {
        type            zeroGradient;
    }


    outlet
    {
        type            zeroGradient;
    }

}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    walls
    {
        type            noSlip;
    }

	atm
    {
        type            pressureInletOutletVelocity;
        value           $internalField;
    }

    inletTop
    {
        type            flowRateInletVelocity;
        volumetricFlowRate constant 7.7;
    }
	
	inletBottom
    {
        type            flowRateInletVelocity;
        volumetricFlowRate constant 7.2;
    }
	
	inletPipe1
    {
        type            flowRateInletVelocity;
        volumetricFlowRate constant 3.4;
    }
	
	inletPipe2
    {
        type            flowRateInletVelocity;
        volumetricFlowRate constant 3.4;
    }
	
	inletPipe3
    {
        type            flowRateInletVelocity;
        volumetricFlowRate constant 3.4;
    }
		
	outlet
    {
		type            flowRateOutletVelocity;
		volumetricFlowRate  table
		(
			(0 27.5)
			(59 27.5)
			(60 12)
			(180 12)
			(181 27.5)
			(240 27.5)
		);
	}

}

My fvSolution and fvSchemes are:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      1;
        nAlphaSubCycles 2;
        cAlpha          1;

        MULESCorr       no;
        nLimiterIter    5;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-10;
        relTol          0;
    }

    "pcorr.*"
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-12;
            relTol          0;
            smoother        DICGaussSeidel;
        }
        tolerance       1e-12;
        relTol          0;
        maxIter         200;
    }

    ".*(rho|rhoFinal)"
    {
        solver          diagonal;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-12;
        relTol          0.01;
        smoother        DIC;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-12;
            relTol          0;
            nVcycles        2;
            smoother        DICGaussSeidel;
            nPreSweeps      2;
        }
        tolerance       1e-12;
        relTol          0;
        maxIter         20;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-10;
        relTol          0;
        nSweeps         1;
    }

    "(T|k|B|nuTilda).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-10;
        relTol          0;
    }
	
	"(U|k|omega|s).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        nSweeps         1;
        tolerance       1e-10;
        relTol          0.1;
    };
}



PIMPLE
{
    momentumPredictor no;
    transonic       no;
    nOuterCorrectors 1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
}

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes 
{
   default                  Euler;                           // first order in time
} 

gradSchemes
{
    default             Gauss linear;                
}

divSchemes // Convection Schemes Settings
{
 
   default                                            none;    
   div(phi,alpha)                                  Gauss vanLeer; 
   div(rhoPhi,U)                                   Gauss linearUpwind grad(U);  // UEqn
   div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear; // vanLeer;
   div(rhoPhi,T)                                    Gauss upwind;    // TEqn
   div(rhoPhi,p)                                    Gauss linear;
   div(phi,p)                                        Gauss linear;
   div(rhoPhi,K)                                    Gauss linear;    // TEqn     
   div(phirb,alpha)                                Gauss interfaceCompression 1; // alphaEqn
}

laplacianSchemes
{
   default              Gauss linear corrected;
}

interpolationSchemes
// interpolating cell-centred values to face values
{
   default                      linear;                              
}
snGradSchemes
{
     default                  corrected;
}


wallDist
{
    method meshWave;
}

I believe that the problem is at my BC's. Are they making sense?

Thank you!

[dzordz]

Hi,


I think, based on my own experience, that this negative temperatures crashes are happening because of the algebraic numerical method calculating the interface. Basically the use of MULES and artificial compression. I would suggest using a geometric description of the interface (isoADVECTOR). The compressible solver compressibleInterIsoFoam was just released with the v2006 release. This should greately stabilize your simulations and hopefully avoid any further crashes (while improving your solution accuracy).


Cheers.

[rpachaly]

Thank you so much for your reply. I'm using the interIsoFoam and it is stable right now.

[mikulo]

Originally Posted by rpachaly

Thank you so much for your reply. I'm using the interIsoFoam and it is stable right now.

Hello,

May I ask what is your pMin?