[hxaxtma]

Hi guys,
I am struggeling with following problem. I set up a case with kOmegaSST and rhoSimpleFoam for calculating external flow at Ma=0.6 and rhoSimpleFoam crashes immediatly after Iteration 2. For lower Ma number the solver runs. I do not find any problem in my setup. Maybe you can help me here

1) My Mesh is fine, checkMesh gives only a few nonOrthoCells and skewCells are found.

Checking geometry...
    Overall domain bounding box (-2.5 -1.9997 -1.99995) (4.5 1.99943 1.99999)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.13991e-16 3.30607e-16 6.08633e-16) OK.
    Max cell openness = 1.26282e-15 OK.
    Max aspect ratio = 41.6873 OK.
    Minimum face area = 2.56578e-11. Maximum face area = 0.0112171.  Face area magnitudes OK.
    Min volume = 8.23851e-15. Max volume = 0.0010981.  Total volume = 87.8213.  Cell volumes OK.
    Mesh non-orthogonality Max: 89.4792 average: 10.4901
   *Number of severely non-orthogonal (> 70 degrees) faces: 272.
    Non-orthogonality check OK.
  <<Writing 272 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 6.92056, 30 highly skew faces detected which may impair the quality of the results
  <<Writing 30 skew faces to set skewFaces
    Coupled point location match (average 0) OK.
Failed 1 mesh checks.
End

2) My BCs are the following

Valid fields:
    volScalarField    nut
    volVectorField    U
    volScalarField    k
    volScalarField    alphat
    volScalarField    p
    volScalarField    T
    volScalarField    omega
    volScalarField    epsilon

wall    : geom_3
wall    : geom_2
wall    : geom_1
    scalar        nut        generic
    scalar        k        generic
    scalar        alphat        generic
    scalar        p        zeroGradient
    scalar        T        zeroGradient
    scalar        omega        generic
    scalar        epsilon        generic
    vector        U        fixedValue

patch    : in
    scalar        nut        calculated
    scalar        k        turbulentIntensityKineticEnergyInlet
    scalar        alphat        calculated
    scalar        p        zeroGradient
    scalar        T        fixedValue
    scalar        omega        fixedValue
    scalar        epsilon        generic
    vector        U        fixedValue

wall    : geom_3
    scalar        nut        generic
    scalar        k        generic
    scalar        alphat        generic
    scalar        p        zeroGradient
    scalar        T        zeroGradient
    scalar        omega        generic
    scalar        epsilon        generic
    vector        U        fixedValue

patch    : out
    scalar        nut        calculated
    scalar        k        zeroGradient
    scalar        alphat        calculated
    scalar        p        fixedValue
    scalar        T        zeroGradient
    scalar        omega        zeroGradient
    scalar        epsilon        zeroGradient
    vector        U        zeroGradient

wall    : side
    scalar        nut        calculated
    scalar        k        slip  -> Is lip condition for k, epsilon and omega OK here?
    scalar        alphat        calculated
    scalar        p        slip
    scalar        T        slip
    scalar        omega        slip
    scalar        epsilon        slip
    vector        U        slip

3)Initializing Turbulence Parameters by CFD-Online Tool at u_0=200m/s, nut/nu=10 and intensity of 0.5% results in

k= 1.5
epsilon = 1350
omega = 10000 -> Is This value OK?

4) I started to underrelax the pressure field and the simulation runs further, but crashes just on later iterations. Here are my fvSchemes

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         cellLimited Gauss linear 1;
}

divSchemes
{
    default          bounded Gauss upwind;

    div(phi,U)       bounded Gauss upwind;
    div(phi,e)       bounded Gauss upwind;
    div(phi,epsilon) bounded Gauss upwind;
    div(phi,omega)   bounded Gauss upwind;
    div(phi,k)       bounded Gauss upwind;
    div(phi,Ekp)     bounded Gauss upwind;

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

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;//corrected;
}

wallDist
{
    method meshWave;
}

5) and my fvSolution

    p
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.05;
        smoother        GaussSeidel;
        cacheAgglomeration on;
        nCellsInCoarsestLevel 20;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-08;
        relTol          0.1;
    }

    e
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0.1;
    }

    "(k|epsilon|omega)"
    {
        $U;
        tolerance       1e-07;
        relTol          0.1;
    }

    Phi
    {
        solver          GAMG;
        smoother        DIC;
        cacheAgglomeration on;
        agglomerator    faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels     1;

        tolerance       1e-06;
        relTol          0.01;
    }


}

potentialFlow
{
    nNonOrthogonalCorrectors 10;
}


SIMPLE
{
    nNonOrthogonalCorrectors 10;
    nCorrectors                2;
    rhoMin          0.5;
    rhoMax          2.5;
    pRefCell        0;
    pRefValue        0;

    residualControl
    {
        p               1e-2;
        U               1e-4;
        e               1e-3;

        // possibly check turbulence fields
        "(k|epsilon|omega)" 1e-3;
    }
}

relaxationFactors
{
    fields
    {
        p               0.01; //was 0.3
        rho             0.01;
    }
    equations
    {
        U               0.7;
        "(k|epsilon|omega)"   0.7;
        e               0.5;
    }
}

and this is the error, it always happens at Solving for p:

smoothSolver:  Solving for Ux, Initial residual = 0.515113, Final residual = 0.0248749, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.572182, Final residual = 0.0216566, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.572538, Final residual = 0.0219375, No Iterations 6
smoothSolver:  Solving for e, Initial residual = 0.950916, Final residual = 0.0732206, No Iterations 1
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::calculate() at ??:?
#4  Foam::hePsiThermo<Foam::psiThermo, Foam::pureMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::hConstThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleInternalEnergy> > > >::correct() at ??:?
#5  ? at ??:?
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? at ??:?
Floating point exception

I am struggeling with this case for two weeks now and do not know what I am setting up wrong. Maybe you have an idea.

And last question, If my boundary layer resolution is y+<1, are the upper known BCs still valid?

Thanks in advance

[Tobi]

Hey,

• a) If you use slip, your BC Type should also be slip

• b) BC type generic (never used, so what is it for?)

• c) Use other inlet BC for k, omega & epsilon; use turbulentIntensity... and some mixingLength BC

• d) I am never worked with flows having compressibility phenomena but I know that there is a additional pressure calculation in rhoSimpleFoam denoted by the keyword "transonic" -> fvSolution

  if (simple.transonic())

• e) is your residual Control working? For me the set-up is wrong (check out some tutorials)

Hope this will help,
At last: You wrote me that you do DNS, with rhoSIMPLEFoam you only get a start solution (its not DNS - but I think you know this).

[hxaxtma]

Hi Tobias,

first of all thanks for your hints. I adapted the BCs a bit, using mixingLength etc..! But it does not change at all

For further clarification, I have two inlets, one with Ma=0.6 and a further one with Ma=8!
So the reason of divergence is the big pressure jump in the domain. Therefore I already used the transient simple option (not included in the thread above).
and I am underrelaxing p until 1e-04 as convergence is reached and setting the parameters back to default. At the moment I do not see another proper way for solving this problem.

[Tobi]

What Do you mean "transient simple options"