[gian93]

Hi Guys,
I'm simulating the unsteady flow in an exhaust pipe with rhoPimpleFoam.
The simulation is a pressure driven simulation. I impose at the inlet a variable bc for the total pressure and a variable static temperature.
The problem is that, unless checkMesh gives positive results with no errors,
simulations blows up in few steps.

I've imposed this BC

inlet (there are multiple inlets since the engine has more than one cylinder)

P >> uniformTotalPressure (tableFile over 2 entire engine cycles)
T >> uniformFixedValue (tableFile over 2 entire engine cycles)
U >> pressureInletOutletVelocity

Outlet (one single outlet)

p >> fixedValue (downstream static pressure before turbine)
T >> inletOutlet
U >>pressureInletOutletVelocity

k-epsilon model with wall function activated

I meshed only the exhaust duct. No additional tube or volumes to stabilize the pressure was connected downstream the x in 1 junction.

I've tried different set up with constant temperature and constant inlet total pressure... no way..

Are those BC suitable for this case? Is rhoPimpleFoam suitable for this type of problems? I would expect pressure waves in the channel.

[gian93]

these are fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1906 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default Euler;
}

gradSchemes
{

default cellMDLimited Gauss linear 1.0;
grad(p) Gauss linear;
grad(U) cellLimited Gauss linear 1;
grad(k) cellLimited Gauss linear 1;
grad(epsilon) cellLimited Gauss linear 1;


}

divSchemes
{


default none;
div(phiU,p) Gauss upwind;

div(phi,U) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,K) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,h) Gauss upwind;
div(U) Gauss upwind;



}

laplacianSchemes
{
default Gauss linear limited 0.5;//0.33;
}

interpolationSchemes
{
default linear;
}

snGradSchemes // surface normal gradient schemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}

wallDist
{
method meshWave;
}



// ************************************************** *********************** //




and fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1906 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{

rho
{
solver PCG;
preconditioner DIC;
minIter 0;
maxIter 2000;
tolerance 1e-08;
relTol 0;
}

rhoFinal
{
solver PCG;
preconditioner DIC;
minIter 0;
maxIter 2000;
tolerance 1e-08;
relTol 0;
}


p


{

solver GAMG;
tolerance 1e-8;
relTol 0.1;
smoother GaussSeidel;
nPreSweeps 1;
nPostSweeps 2;
nFinestSweeps 3;
scaleCorrection true;
directSolveCoarsestLevel false;
cacheAgglomeration true;
nCellsInCoarsestLevel 500;
agglomerator faceAreaPair;
mergeLevels 1;


}

pFinal
{
$p;
tolerance 1e-8;
relTol 0;

maxIter 2000;
}


"(U|h|k|epsilon|omega)"
{
solver PBiCGStab;
preconditioner DILU;
tolerance 1e-08; //1e-06;
relTol 0.1;
}

"(U|h|k|epsilon|omega)Final"
{
$U;
tolerance 1e-08;
relTol 0;
}


}

SIMPLE
{
nNonOrthogonalCorrectors 4;
}

PIMPLE
{

nOuterCorrectors 6;
nCorrectors 4;



rhoMin rhoMin [ 1 -3 0 0 0 ] 1e-5;
rhoMax rhoMax [ 1 -3 0 0 0 ] 1e5;


momentumPredictor on;
nNonOrthogonalCorrectors 8;

transonic on;
turbOnFinalIterOnly off;
}

fluxRequired
{
p;
}


relaxationFactors
{

equations
{
U 0.8;

k 0.2;
epsilon 0.2;


}
fields
{
p 0.2;
}

}


// ************************************************** *********************** //

[gian93]

nobody can give me a hint?
my doubts are related to the fact that rhoPimpleFoam is effectively suitable for this type of simulation in which pressure waves and high Mach number are expected and on the correctness of the BC that I have imposed.

Any help are appreciated

[HPE]

If I am not wrong, you have stated this: 'unless checkMesh fails', you do not observe this issue. Bad mesh quality almost always suggests problems of any sort.

[gian93]

Hi HPE ,
I apologize about that. I've make a mistake! I meant that checkMesh do not gives errors. The only thing is that max non-orthogonality is a little bit high..

Checking geometry...
Overall domain bounding box (-0.191584 -0.653172 0.0211068) (0.189914 -0.241768 0.272275)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-3.00775e-16 2.09832e-16 6.18816e-16) OK.
Max cell openness = 4.04177e-16 OK.
Max aspect ratio = 15.7391 OK.
Minimum face area = 8.36371e-09. Maximum face area = 2.91039e-06. Face area magnitudes OK.
Min volume = 8.35547e-13. Max volume = 4.19577e-09. Total volume = 0.00170685. Cell volumes OK.
Mesh non-orthogonality Max: 68.5039 average: 5.93194
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 1.65157 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End

I've probably missed something in the BC . What do you think about them?
this is the last part of the log of the simulation … (here reported some slices)
Cumulative error is increasing after few iterations.

time step continuity errors : sum local = 9.2315e-05, global = 1.14419e-11, cumulative = 0.0279472
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.0258494, Final residual = 6.79523e-09, No Iterations 4
DILUPBiCGStab: Solving for k, Initial residual = 0.0571322, Final residual = 3.3752e-09, No Iterations 4
ExecutionTime = 1656.01 s ClockTime = 1669 s

yPlus yPlus write:
writing field yPlus
patch Cyl5Pipe y+ : min = 1.28177, max = 9.87123, average = 2.03242
patch Cyl6Pipe y+ : min = 1.29158, max = 11.7239, average = 1.9949
patch Cyl7Pipe y+ : min = 1.20546, max = 35.8124, average = 2.57857
patch Cyl8Pipe y+ : min = 1.29889, max = 5.18343, average = 1.86761
patch fourInOne y+ : min = 0.718431, max = 3.11383, average = 1.87058
patch restVolume y+ : min = 1.1881, max = 3.04107, average = 1.80754
Courant Number mean: 0.0135969 max: 15.6855
deltaT = 4.4049e-06
Time = 2.5391e-05

PIMPLE: iteration 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCGStab: Solving for Ux, Initial residual = 0.106514, Final residual = 0.00156979, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.0532714, Final residual = 0.000641332, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.0618325, Final residual = 0.000705202, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.0670228, Final residual = 0.00569434, No Iterations 1
GAMG: Solving for p, Initial residual = 0.101657, Final residual = 0.00369779, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00604541, Final residual = 0.000449623, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00075055, Final residual = 4.70863e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 9.16001e-05, Final residual = 6.7737e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 1.25512e-05, Final residual = 1.10542e-06, No Iterations 3
…….
……
…..
time step continuity errors : sum local = 0.000834054, global = 0.000767087, cumulative = 0.0371216
GAMG: Solving for p, Initial residual = 0.0524838, Final residual = 0.00242037, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00382485, Final residual = 0.000336143, No Iterations 2
GAMG: Solving for p, Initial residual = 0.000514833, Final residual = 4.01419e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 7.08663e-05, Final residual = 6.27203e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 1.07244e-05, Final residual = 1.06986e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 1.87011e-06, Final residual = 1.13159e-07, No Iterations 4
GAMG: Solving for p, Initial residual = 3.24669e-07, Final residual = 2.38007e-08, No Iterations 3
GAMG: Solving for p, Initial residual = 9.80524e-08, Final residual = 7.67282e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 6.5788e-08, Final residual = 4.72473e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.92555e-08, Final residual = 4.13328e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.74656e-08, Final residual = 4.01796e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.68362e-08, Final residual = 4.00461e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.65253e-08, Final residual = 4.02891e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.63408e-08, Final residual = 4.07174e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.61609e-08, Final residual = 4.12716e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.60003e-08, Final residual = 4.198e-09, No Iterations 3
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.000832999, global = 0.000766636, cumulative = 0.0378882
GAMG: Solving for p, Initial residual = 0.0524932, Final residual = 0.00242304, No Iterations 2
GAMG: Solving for p, Initial residual = 0.00382904, Final residual = 0.000336586, No Iterations 2
GAMG: Solving for p, Initial residual = 0.0005156, Final residual = 4.01908e-05, No Iterations 3
GAMG: Solving for p, Initial residual = 7.09604e-05, Final residual = 6.27803e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 1.07348e-05, Final residual = 1.07055e-06, No Iterations 3
GAMG: Solving for p, Initial residual = 1.86949e-06, Final residual = 1.1331e-07, No Iterations 4
GAMG: Solving for p, Initial residual = 3.23173e-07, Final residual = 2.38251e-08, No Iterations 3
GAMG: Solving for p, Initial residual = 9.76716e-08, Final residual = 7.64007e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 6.52732e-08, Final residual = 4.67573e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.87994e-08, Final residual = 4.07805e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.69443e-08, Final residual = 3.96139e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.62847e-08, Final residual = 3.94669e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.59683e-08, Final residual = 3.97036e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.57785e-08, Final residual = 4.01183e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.55885e-08, Final residual = 4.06501e-09, No Iterations 3
GAMG: Solving for p, Initial residual = 5.54168e-08, Final residual = 5.33879e-10, No Iterations 5
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 0.000832747, global = 0.000766515, cumulative = 0.0386548
DILUPBiCGStab: Solving for epsilon, Initial residual = 0.00598927, Final residual = 9.42791e-06, No Iterations 1
DILUPBiCGStab: Solving for k, Initial residual = 0.0124009, Final residual = 1.91494e-05, No Iterations 1
PIMPLE: iteration 4
DILUPBiCGStab: Solving for Ux, Initial residual = 0.0398893, Final residual = 0.000846991, No Iterations 1
DILUPBiCGStab: Solving for Uy, Initial residual = 0.0191814, Final residual = 0.000304033, No Iterations 1
DILUPBiCGStab: Solving for Uz, Initial residual = 0.0217133, Final residual = 0.000332122, No Iterations 1
DILUPBiCGStab: Solving for h, Initial residual = 0.0167912, Final residual = 0.00025334, No Iterations 1
-------------------------------------------------------
Primary job terminated normally, but 1 process returned
a non-zero exit code.. Per user-direction, the job has been aborted.
-------------------------------------------------------
--------------------------------------------------------------------------
mpirun detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

Process name: [[56472,1],12]
Exit code: 144

[cryabroad]

I normally don't use that many nNonOrthogonalCorrectors (you have 8 of them), I either have 1 or 0, occsinally I use 2. If you think you need that many nNonOrthogonalCorrectors, maybe you need to create a better mesh. Also, your nCorrectors is set to 4, which is also not common. 2 is probably good enough. I would use much more nOuterCorrectors, set it to like 50 and see if the simulation converges. By the way, relaxing the pressure field could be helping. And, why is the transonic option set to on? The flow is trasonic?

Something that can really help you figure out the problem is, since you are doing RANS, run a steady-state simulation and see if that converges or not. If yes, you can play with the settings specifically for PIMPLE(nOuterCorrectors, etc.). If not, that means your boundary conditions are likely to be not good.

[gian93]

Originally Posted by cryabroad

I normally don't use that many nNonOrthogonalCorrectors (you have 8 of them), I either have 1 or 0, occsinally I use 2. If you think you need that many nNonOrthogonalCorrectors, maybe you need to create a better mesh. Also, your nCorrectors is set to 4, which is also not common. 2 is probably good enough. I would use much more nOuterCorrectors, set it to like 50 and see if the simulation converges. By the way, relaxing the pressure field could be helping. And, why is the transonic option set to on? The flow is trasonic?

Something that can really help you figure out the problem is, since you are doing RANS, run a steady-state simulation and see if that converges or not. If yes, you can play with the settings specifically for PIMPLE(nOuterCorrectors, etc.). If not, that means your boundary conditions are likely to be not good.

hi ! thanks for your reply, I will try to follow your hints for the set up of the Pimple. I think the Mach number of the flow involved is really close the unit in some point of the domain (especially near the engine exhaust port, my inlet patches).
I've seen that without turbulence activated the simulation runs, in addition I can say that the problem should be related to the strong pressure gradient that arouse cyclically inside each exhaust collector. I've tried to smooth the pressure tableFile initial values but when pressure increases the simulation blows up soon.

[cryabroad]

Sorry but I'm not that familiar with trasonic flows and how to solve them properly, I'm sure someone else on this forum can help.

You mentioned that the simulation runs without activating turbulence, that sounds interesting. I would think that without turbulence activated your case should blow up pretty quickly, but it's the opposite. Also, you have multiple inlets (probably with different pressures and temperatures) and a single outlet, and you are close to Mach 1. It' sounds to be not an easy case even for a steady-state simulation, not to mention transient processes. Still, I would set up a quick test, which is a steady-state simulation with different pressure and temperature values at the different inlets. See if it converges or not.

[gian93]

Originally Posted by cryabroad

You mentioned that the simulation runs without activating turbulence, that sounds interesting. I would think that without turbulence activated your case should blow up pretty quickly

yes your' re right. I've checked again and it is obvious.
I cannot get rid of the problem of continuity error!

I've tried a simple channel flow (maintaining rhoPimpleFoam) with reduced inlet total pressure and kOmegaSST TURB MODEL but nothing change: sum local, global and cumulative error increase at any iteration.. here are my bc for the simplified test:

U

PHP Code:

walls    {   type            fixedValue;        value           uniform (0 0 0);                 }Inlet    {       type            pressureInletOutletVelocity;        phi        phi;        rho        rho;        value           uniform (0 0 0);    }outlet    {        type            pressureInletOutletVelocity;        phi        phi;        rho        rho;        value           uniform (0 0 0);      }} 


T

PHP Code:

internalField   uniform 1000; //mean valueboundaryField{        Inlet    {       type            totalTemperature;       gamma           1.4;       T0              uniform 1160;       U               U;       phi             phi;     }    outlet    {        type       inletOutlet;       inletValue      uniform 1220;        value           uniform 1220;    }   wall    {        type            fixedValue ;    value         uniform 1000;//zeroGradient;    }} 


p

PHP Code:

internalField   uniform 361022;boundaryField{            wall    {        type            zeroGradient;    }    Inlet    {                type           totalPressure;        U               U;        phi             phi;        rho             rho;        psi             none;        gamma           1.4;    p0              uniform 380000;    value           uniform 380000;        }outlet    {        type            fixedValue;        value           uniform 361022;    }} 


nut

PHP Code:

boundaryField{   wall    {        type            nutUSpaldingWallFunction;    value           uniform 0;    }    Inlet    {        type            calculated;        value           uniform 0;    }  outlet    {        type            calculated;        value           uniform 0;    }} 


k

PHP Code:

boundaryField{ wall    {        type            kqRWallFunction;        value           uniform 4.2;    }  Inlet    {        type            turbulentIntensityKineticEnergyInlet;        intensity       0.05;        value           uniform 4.2;         /* type            inletOutlet;        inletValue      uniform 4.2;        value           uniform 4.2;*/    } outlet    {         type            inletOutlet;        inletValue      uniform 4.2;        value           uniform 4.2;    }} 


Omega

PHP Code:

internalField   uniform 388;boundaryField{  wall            {                 type            omegaWallFunction;                                 value           uniform 388;            }  Inlet            {               type            turbulentMixingLengthFrequencyInlet;                mixingLength    0.023;                value           uniform 388;            /*type            inletOutlet;                inletValue      uniform 388;                value           uniform 388;*/            }     outlet            {                type            inletOutlet;                inletValue      uniform 388;                value           uniform 388;           }} 


alphat

PHP Code:

internalField   uniform 0;boundaryField{walls    {        type            compressible::alphatWallFunction;        Prt             0.85;        value           uniform 0;    }Inlet    {        type            calculated;         value           uniform 0;    }outlet    {        type            calculated;        value           uniform 0;    }}