Compressible nozzle flow - no problem for fluent but impossible with openFOAM?

flar.t · November 23, 2021, 11:45

Hello everybody,

after struggling with that case for a long while now, I need some help with openFOAM. Perhaps there is someone out there who could give me some guidance.
I'm working on a 2D axisymmetric model of the compressible flow through a converging-diverging nozzle. After the flow leaves the nozzle it passes a backwards facing step:
Bild1.png

Solving for the same geometry and mesh fluent converges smoothly within a few thousand iterations. This holds for the realizable k-epsilon (non-equilibrium wall function) and the k-omega turbulence model both with viscous heating enabled. The solver configuration is as follows:

Pressure-Velocity Coupling:
- Scheme: SIMPLE
- Flux Type: Rhi-Cho: distance based
Spatial Discretization:
- Gradient: Least Squares Cell based
- Pressure: Second Order
- Density, Momentum, Turbulent Kinetic Energy, Turbulent Dissipation Rate, Energy: Second Order Upwind
Under-Relaxation Factors
- Pressure: 0.3
- Density, Body Forces: 1
- Momentum: 0.7
- Turbulent Kinetic Energy, Turbulent Dissipation Rate: 0.8
- Turbulent Viscosity: 1
- Energy: 1
Solution Limits:
- Temperature: 1-5000K
- Pressure: 1-5e10 Pa
Advanced Solution Controls
- Multigrid Solver with Termination 0.1 and no Stabilization Method
- Max Coarse Levels: 40
- Smoother: Gauss-Seidel
- Spatial Discretization Limiter: Standard, Cell to Face Limiting

Boundary conditions are:

Inlet:
- Mass Flow Rate: 0.003 kg/s
- Initial Gauge Pressure: 0
- Turbulence Intensity: 5%
- Turbulence Viscosity Ration: 10
Outlet:
- Gauge Pressure : 101325 Pa
- Backflow Pressure Specification: Total Pressure
- Backflow Turbulence Intensity: 5%
- Backflow Turbulence Viscosity Ration: 10
Wall:
- no Slip
- Standard Roughness Model
Axis

However, when I tried to rebuild the case with openFOAM I'm unable to get a converged solution. With the realizable k-epsilon model the solver at least doesn't crash but with the k-omega it will blow up within a couple of iterations. In both cases the time step continuity errors are huge while the linear solvers for the different fields show small residuals.

What I did so far:

Started from a openFOAM tutorials > incompressible > simpleFoam > pitzDaily since there is no similar case for rhoSimpleFoam available.
Read the corresponding parts of the OpenFOAM: User Guide v2012 and OpenFOAM: API Guide v2012 to understand how to make a similar setup like I had in fluent. Additionally I read a lot of forum posts concerning the solution and schemes settings for similar cases.
Prepared a fine mesh with blockMesh and the makeAxialMesh utility (http://openfoamwiki.net/index.php/Co...akeAxialMesh); checkMesh returns no errors:

Code:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Time = 0

Mesh stats
    points:           432246
    internal points:  0
    faces:            861240
    internal faces:   428995
    cells:            215250
    faces per cell:   5.99412311266
    boundary patches: 7
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     213985
    prisms:        1265
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
    Patch               Faces    Points   Surface topology                  
    inlet               50       101      ok (non-closed singly connected)  
    outlet              240      481      ok (non-closed singly connected)  
    upperWall           1455     2912     ok (non-closed singly connected)  
    lowerWall           0        0        ok (empty)                        
    frontAndBack        0        0        ok (empty)                        
    frontAndBack_pos    215250   216756   ok (non-closed singly connected)  
    frontAndBack_neg    215250   216756   ok (non-closed singly connected)  

Checking faceZone topology for multiply connected surfaces...
    No faceZones found.

Checking basic cellZone addressing...
    No cellZones found.

Checking geometry...
    Overall domain bounding box (0 0 -0.000436193873653) (0.08 0.01 0.000436193873653)
    Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Wedge frontAndBack_pos with angle 2.49762357054 degrees
    Wedge frontAndBack_neg with angle 2.49762357054 degrees
    All edges aligned with or perpendicular to non-empty directions.
    Boundary openness (4.32643646952e-18 2.00893658027e-17 -8.27156673543e-16) OK.
    Max cell openness = 3.03021770314e-16 OK.
    Max aspect ratio = 10.4202379507 OK.
    Minimum face area = 3.28108080426e-11. Maximum face area = 1.34703394904e-07.  Face area magnitudes OK.
    Min volume = 1.00959120606e-15. Max volume = 1.27411305502e-11.  Total volume = 2.53751454593e-07.  Cell volumes OK.
    Mesh non-orthogonality Max: 23.0894947956 average: 3.34709052281
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.568549620689 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

Tried different settings for the gradient schemes (according to rhoSimpleFoam crash help or https://www.openfoam.com/documentati...t-example.html and https://www.openfoam.com/documentati...e-example.html)

Tried different settings for the solvers and under-relaxation factors

GAMG solver (similar to the fluent multigrid setting) with different nCellsInCoarsestLevel (https://wiki.openfoam.com/OneraM6_by_Michael_Alletto or https://www.researchgate.net/post/Is...-fine-hex-mesh) and different tolerances (Questions about fvsolution, Why OF set relTol to control convergence?)

SIMPLE with "consistent no" and relaxation or with "consistent yes" and no relaxation (rhoSimpleFoam crash help)

Changed boundary conditions

Calculated the boundary conditions for turbulence parameters similar to fluent (turbulentIntensity for k and mixingLength for epsilon)

Checked the yplus (<3) and ratio of integral length scale to grid length scale (>5) according to http://www.wolfdynamics.com/training..._tutorials.pdf . Nevertheless, I used y+ insensitive wall functions.

Tried small modifications for the boundary conditions (e.g. U@outlet: inletOutlet vs. zeroGradient) as discussed in SimpleFoam high time-step continuity error on a simple geometry

Tried to "ramp" the inlet flow rate with a table after the flow has converged for a very small flow rate as initialization. Also tried initialization with potentialFoam

Checked for working boundary conditions with an incompressible simpleFoam run but even when the time step continuity errors are quite small (1e-4) compared to the compressible case there is no real convergence for the field variables

There are many similar questions in this forum, e.g.:

SimpleFoam high time-step continuity error on a simple geometry

rhoSimpleFoam time step continuity errors

What happens with my k and epsilon after a few timesteps

time step continuity error, fvScheme

But most of the suggested solutions concern the above already mentioned modifications: bounded schemes vs. unbounded, small under relaxation factors or no relaxation at all, problems with the mesh, or problems with the boundary conditions. Using the bounded schemes stops the negative values for the turbulence variables. But I did not get rid of the high time step continuity errors. Perhaps my boundray conditions are wrong because the incompressible case did not converge, too.

I would really appreciate any help in solving this problem.

Best regards,
Ralf

These are my case files:

0/alphat:

Code:

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

internalField   uniform 1e-3;

boundaryField
{
    inlet
    {
        type            calculated;
        value        $internalField;
    }
    outlet
    {
        type            calculated;
        value        $internalField;
    }
    upperWall
    {
        type            compressible::alphatWallFunction;
        value           $internalField;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

0/epsilon:

Code:

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 200;

boundaryField
{
    inlet
    {
        type            turbulentMixingLengthDissipationRateInlet;
        mixingLength    5.96e-3;
        value           uniform 200;

    }

    outlet
    {
        type            zeroGradient;
        value        uniform 200;
    }

    upperWall
    {
        type            epsilonWallFunction;
        value        uniform 200;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

0/k:

Code:

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

internalField   uniform 1;

boundaryField
{
    inlet
    {
        type            turbulentIntensityKineticEnergyInlet;
        intensity       0.05;
        value           uniform 1;
    
    }
    outlet
    {
        type            zeroGradient;
        value           uniform 1;
    }
    upperWall
    {
        type         kqRWallFunction;
        value         uniform 1;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

0/nut:

Code:

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            calculated;
    value        $internalField;
    }
    outlet
    {
        type            calculated;
    value        $internalField;
    }
    upperWall
    {
        type            nutkWallFunction;
        value           $internalField;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

0/p:

Code:

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


internalField   uniform 101325;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
           type            fixedValue;
        value           uniform 101325;        
    }
    upperWall
    {
        type            zeroGradient;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }   
}

0/T:

Code:

dimensions      [0 0 0 1 0 0 0];


internalField   uniform 298;

boundaryField
{
    inlet
    {
        type         fixedValue;
        value        $internalField;
    }
    outlet
    {
        type            inletOutlet;
        inletValue    $internalField;
        value        $internalField;
    }
    upperWall
    {
        type            zeroGradient;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

0/U:

Code:

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

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type flowRateInletVelocity;
        massFlowRate constant 4e-5;
        //massFlowRate tableFile;
        //file    "system/mdottable";
        //outOfBounds    clamp;
        rhoInlet    1;
    }
    outlet
    {
        type            inletOutlet;
        inletValue    uniform (0 0 0);
        value        uniform (0 0 0);
    }
    upperWall
    {
        type        noSlip;
    }
    lowerWall
    {
        type            symmetryPlane;
    }
    frontAndBack
    {
        type            empty;
    }
    frontAndBack_pos
    {
        type            wedge;
    }
    frontAndBack_neg
    {
        type            wedge;
    }
}

constant/thermophysicalProperties:

Code:

thermoType
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       28.966;
    }
    thermodynamics
    {
        Cp              1006.43;
        Hf              2.544e+6;
    }
    transport
    {
        mu              1.7894e-05;
        Pr              0.737;
    }
}

constant/turbulenceProperties:

Code:

simulationType  RAS;

RAS
{
   RASModel        realizableKE;

   turbulence      on;

   printCoeffs     on;
}

system/controlDict:

Code:

application     rhoSimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         25000;

deltaT          1;

writeControl    timeStep;

writeInterval   500;

purgeWrite      0;

writeFormat     ascii;

writePrecision  12;

writeCompression off;

timeFormat      general;

timePrecision   12;

runTimeModifiable true;

functions
{
    #includeFunc mag(U)
    #includeFunc flowRatePatch(name=inlet)
    #includeFunc flowRatePatch(name=outlet)
    #includeFunc patchAverage(name=inlet,p)
    minmaxdomain
    {
    type fieldMinMax;

    libs        (fieldFunctionObjects);

    enabled true; //true or false

    mode component;

    writeControl outputTime;

    log true;

    fields (p U ); //k epsilon omega nut);
    }

    yplus
    {
        type yPlus;
        libs        (fieldFunctionObjects);
        enabled true;
        writeControl outputTime;
    }
    
    writecellvol
    {
        type writeCellVolumes;
        libs        (fieldFunctionObjects);
        enabled true;
        writeControl outputTime;
    }

    integralLenghtScale
    {
    // Mandatory entries (unmodifiable)
    type        turbulenceFields;
    libs        (fieldFunctionObjects);

    // Mandatory entries (runtime modifiable)
    // Either field or fields entries
    field      L;
    enabled true;
    writeControl outputTime;
    }
}

system/fvOptions:

Code:

temperature_constraints
{
    type            limitTemperature;

    selectionMode   all;
    active        true;
    
    limitTemperatureCoeffs
    {
       selectionMode   all;
       min  1;
       max  5000;
    }
}

system/fvSchemes:

Code:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         cellLimited Gauss linear 1; //faceLimited leastSquares 1.0;

}

divSchemes
{
    default                 none;

    div(phi,U)              bounded Gauss linearUpwind default; //Gauss limitedLinearV 1.0;

    div(phi,k)              bounded Gauss linearUpwind default; //Gauss limitedLimitedLinear 1.0 1e-8 1e10;
    div(phi,epsilon)         bounded Gauss linearUpwind default; //Gauss limitedLimitedLinear 1.0 1e-8 1e10;
    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
    
    div(phi,e)              bounded Gauss linearUpwind default; //bounded Gauss upwind;    
    div(phi,Ekp)            bounded Gauss linearUpwind default; //bounded Gauss upwind;

    div(phid,p)             Gauss limitedLinear 1.0;

    div((phi|interpolate(rho)),p)  Gauss limitedLinear 1.0;
}

laplacianSchemes
{
    default         Gauss linear limited 1;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 1;
}

system/fvSolution:

Code:

solvers
{
   p
   {
       solver          GAMG;
       tolerance       1e-07;
       relTol          0;
       smoother        GaussSeidel;
       nCellsInCoarsestLevel 400;
   }
   "(U|e|k|epsilon|nuTilda)"
   {
       solver          GAMG;
       tolerance       1e-07;
       relTol          0;
       smoother        GaussSeidel;
       nCellsInCoarsestLevel 400;
   }
}

SIMPLE
{
   nNonOrthogonalCorrectors 0;
   transonic       yes;
   consistent      no;

   residualControl
   {
       p               1e-5;
       U               1e-6;
       e               1e-5;
       "(k|epsilon|nuTilda)" 1e-5;
   }
}

relaxationFactors
{
   fields
   {
       p               0.3;
   }
   equations
   {
       U               0.9;
       e               0.8;
       k               0.9;
       epsilon         0.9;
       nuTilda         0.8;
   }
}

dlahaye · November 24, 2021, 09:38

Greetings,

rhoSimpleFoam solves conservation of mass, momentum and energy. The boundary conditions (inlet especially) should be compatible with the coupling of pressure, velocity and temperature. This makes a compressible rhoSimpleFoam case harder to set up than an incompressible simpleFoam case.

The ideal gas law defines a relation between temperate, pressure and density (and thus velocity through mass flow rate). Fixing two quantities, the third one can be determined. This means that temperature, pressure and density can *not* be set (or user defined) independently from each other.

In my case, e.g., it proved to be valuable to fix inlet pressure (e.g p = 101325 Pascal) and inlet temperature (e.g. T = 773.15 K) to arrive at density rho = p / (R_{specific} * T) = 101325 / (287.058 * 773.15) = 0.45654 kg/m^3 . Given density and area of the inlet, the inlet velocity can be determined.

In your case set-up, I fail to understand the inlet Gauge pressure set to 0.

My understanding is that you are very close to a working solution. Good luck.

flar.t · November 24, 2021, 12:19

Thank you very much for your reply.
The "Initial Gauge Pressure" is only specified in the fluent case. As far as I understand it is used only for initialization. The openFOAM case uses a zeroGradient boundary condition at the inlet. I thought with that BC the pressure itself is not specified and the ideal gas law is used as you described to calculate it from the fixed velocity and temperature.

In the meanwhile I rebuilt an incompressible case from scratch and the simpleFoam solver converges really smoothly within 2000 iterations for the same BCs on pressure, velocity and turbulence as mentioned above. The moment I change from incompressible to compressible(*) the solver crashes.

(*) by

- adding thermophysicalProperties,
- changing pressure from m^2/s^2 to kg/(m*s^2) in 0/p,
- adding rho to massFlowRate in 0/U,
- adding 0/T and 0/alphat,
- modifying divScheme for the shear-rate tensor (rho*...),
- adding the schemes and solver for the energy eq.

I noticed that the convergence strongly depends on the gradient schemes for the incompressible case. While kEpsilon is fine with "Gauss linear" the kOmega works only when gradSchemes is set to "cellLimited leastSquares 1". Perhaps the crashes result from the wrong schemes or solver settings.. Are there any recommendations for compressible internal flows?

dlahaye · November 24, 2021, 12:50

My bad: pressure indeed has zeroGradient on the inlet patch as you suggest.

What do you mean by "solver crashing"? What happens exactly? Can you share solver log file?

What happens is case that rho on inlet is increased from 1 to 100, say (increasing density should decrease velocity at same mass flow rate)? Does this lowering of the inlet velocity make the case easier to converge?

flar.t · November 24, 2021, 17:06

By crashing I mean the solver stops more or less rapidly.
I think these are floating point exceptions visible in this example but I have no idea where they come from.

Code:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _7bdb509494-20201222 OPENFOAM=2012
Arch   : "LSB;label=32;scalar=64"
Exec   : rhoSimpleFoam
Date   : Nov 24 2021
Time   : 21:40:22
Host   : xxxx
PID    : 6509
I/O    : uncollated
Case   : run/darm22_RhoSimpleFoam_mdot_kOmegaSST3
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: convergence criteria
    field p     tolerance 0.01
    field U     tolerance 0.001
    field "(k|e|epsilon|omega|f|v2)"     tolerance 0.001

Reading thermophysical properties

Selecting thermodynamics package 
{
    type            hePsiThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleInternalEnergy;
}

Reading field U

Reading/calculating face flux field phi

Creating turbulence model

Selecting turbulence model type RAS
Selecting RAS turbulence model kOmegaSST
Selecting patchDistMethod meshWave
RAS
{
    RASModel        kOmegaSST;
    turbulence      on;
    printCoeffs     on;
    alphaK1         0.85;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.856;
    gamma1          0.555555555556;
    gamma2          0.44;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    a1              0.31;
    b1              1;
    c1              10;
    F3              false;
    decayControl    false;
    kInf            0;
    omegaInf        0;
}

No MRF models present

Creating finite volume options from "system/fvOptions"

Selecting finite volume options type limitTemperature
    Source: temperature_constraints
    - selecting all cells
    - selected 215250 cell(s) with volume 2.53751454593e-07

Starting time loop

surfaceFieldValue flowRatePatch(name=inlet):
    operation     = sum


surfaceFieldValue flowRatePatch(name=outlet):
    operation     = sum


surfaceFieldValue patchAverage(name=inlet,p):
    operation     = areaAverage


turbulenceFields integralLenghtScale: storing fields:
    turbulenceProperties:L

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 0.0728821995915, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Uz, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for e, Initial residual = 0.999955154113, Final residual = 0.0982901196004, No Iterations 3
limitTemperature temperature_constraints Lower limited 0 (0%) of cells
limitTemperature temperature_constraints Upper limited 0 (0%) of cells
limitTemperature temperature_constraints Unlimited Tmax 298
Unlimited Tmin 298
GAMG:  Solving for p, Initial residual = 0.999999999373, Final residual = 0.0610521226431, No Iterations 21
time step continuity errors : sum local = 5.14879974762, global = -0.0250572123018, cumulative = -0.0250572123018
smoothSolver:  Solving for omega, Initial residual = 0.163249893053, Final residual = 0.0138597003377, No Iterations 3
bounding omega, min: -31743.4509378 max: 10700467.3544 average: 79932.7997844
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.0721647330784, No Iterations 5
ExecutionTime = 4.68 s  ClockTime = 5 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.307889634693, Final residual = 0.0206100973894, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.429556159428, Final residual = 0.0429469947432, No Iterations 3
smoothSolver:  Solving for Uz, Initial residual = 0.502502639153, Final residual = 0.0112927535803, No Iterations 1
smoothSolver:  Solving for e, Initial residual = 0.999956608923, Final residual = 0.0255575610827, No Iterations 2
limitTemperature temperature_constraints Lower limited 0 (0%) of cells
limitTemperature temperature_constraints Upper limited 0 (0%) of cells
limitTemperature temperature_constraints Unlimited Tmax 298.02092684
Unlimited Tmin 298
GAMG:  Solving for p, Initial residual = 0.00282266766787, Final residual = 0.000277825629673, No Iterations 12
time step continuity errors : sum local = 3102.24138499, global = -1592.12762623, cumulative = -1592.15268344
smoothSolver:  Solving for omega, Initial residual = 0.051531840972, Final residual = 0.00407122098254, No Iterations 1
bounding omega, min: -1743823.00544 max: 4006080691.61 average: 89426.2874374
smoothSolver:  Solving for k, Initial residual = 0.327474312565, Final residual = 0.0122084139505, No Iterations 2
bounding k, min: -1.09712825705 max: 14.4920020548 average: 0.389117493712
ExecutionTime = 7.67 s  ClockTime = 8 s

Time = 3

smoothSolver:  Solving for Ux, Initial residual = 0.117821454106, Final residual = 0.00901320173028, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0927055396049, Final residual = 0.0089339016751, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.375478299784, Final residual = 0.0141813272293, No Iterations 1
smoothSolver:  Solving for e, Initial residual = 0.279167325572, Final residual = 0.0118560063607, No Iterations 2
limitTemperature temperature_constraints Lower limited 8 (0.00371660859466%) of cells
limitTemperature temperature_constraints Upper limited 0 (0%) of cells
limitTemperature temperature_constraints Unlimited Tmax 415.915872131
Unlimited Tmin 176.302097449
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in /lib64/libpthread.so.0
#3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
#4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
#7  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
#8  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#9  Foam::fvMatrix<double>::solve() at ??:?
#10  ? at ??:?
#11  __libc_start_main in /lib64/libc.so.6
#12  ? at ??:?

I managed to get a solution for very very low massFlowRates and by setting the relTol near to zero. The latter seems to reduce the problem of high time step continuity errors I saw before. I will try to ramp the flowrate slowly by using a table.

Furthermore I'm thinking about the outlet boundary conditions. I guess there is backflow in some parts of the outlet. So I set the type to inletOutlet in 0/U. But I have no idea which inletValue I can use. For k and omega the zeroGradient at the outlet might also be problematic.

dlahaye · November 24, 2021, 17:15

Do you require a waveTransmissive boundary condition for the pressure at the outlet? See https://openfoamwiki.net/index.php/H...dary_condition

DevilX · November 25, 2021, 09:41

Hey,

you tried a lot, thats true. RhoSimpleFoam is kind of special and very sensitive. Give a transient calculation with rhoPimple a try, it's much more robust and versatile.

Your log file indicates, that your timeStepErrors are really high ( around 10^7 bigger then they should), also you have high bounding values in k and omega. What velocities (maximium) are you expecting in your domain? Simplefoam should give you a value there. Is it really near the transonic region of your fluid (fvSolution -transonic - yes should always be off most times).

What are your timeStepErrors and Residuals in SimpleFoam? Btw: this solver is so robust, it runs with most settings.

flar.t · November 25, 2021, 12:08

I'm really glad to get your feedback - thank you. This helps a lot.

In the meanwhile I was able to check two things:
1. I was unsure about the right usage of the waveTransmissive boundary condition. So I added a really long duct at the outlet to move the place for the zeroGradient boundary condition far away from the nozzle throat - expecting the flow to be more homogeneous there. And indeed this results in a very smooth convergence for lower mass flow rates (about 1/10 of the actual flow rate).

2. When I then ramp the mass flow rate to higher values this seems to be okay up to 1/4 of the actual flow rate. Then the residuals rise again, turbulence and temperature become bounded and at last I get a floating point exception. Here are the last steps from the log:

Code:

Time = 10007

smoothSolver:  Solving for Ux, Initial residual = 0.682938315105, Final residual = 0.000287364778599, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.682938327772, Final residual = 0.000335502672822, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.956672521947, Final residual = 0.0287743274877, No Iterations 1000
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 0.00171898979327, No Iterations 1000
limitTemperature temperature_constraints Lower limited 4332 (2.01254355401%) of cells
limitTemperature temperature_constraints Upper limited 10694 (4.96817653891%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.248165319778, Final residual = 1.69988878639e-05, No Iterations 18
time step continuity errors : sum local = 9.08666371988e+32, global = -7.09410752815e+18, cumulative = -7.09410752815e+18
smoothSolver:  Solving for omega, Initial residual = 0.555683310243, Final residual = 7.65079881608e-06, No Iterations 1000
bounding omega, min: -2.22543731148e+24 max: 3.73801862533e+24 average: 6.25428322118e+19
smoothSolver:  Solving for k, Initial residual = 0.190530340693, Final residual = 8.86466389855e-09, No Iterations 99
bounding k, min: -45650036586.1 max: 1.09206304413e+12 average: 19063766.1197
ExecutionTime = 20.13 s  ClockTime = 20 s

Time = 10008

smoothSolver:  Solving for Ux, Initial residual = 0.337580144242, Final residual = 0.000100646177355, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.337548080513, Final residual = 9.32237768003e-05, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.307343200808, Final residual = 0.00266242950238, No Iterations 1000
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 5.65715128871e-08, No Iterations 1000
limitTemperature temperature_constraints Lower limited 115602 (53.7059233449%) of cells
limitTemperature temperature_constraints Upper limited 3259 (1.51405342625%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.952214969982, Final residual = 6.07324943602e-05, No Iterations 9
time step continuity errors : sum local = 1.24414784332e+47, global = 2.24296712762e+33, cumulative = 2.24296712762e+33
smoothSolver:  Solving for omega, Initial residual = 0.000916815161875, Final residual = 7.5135316827e-09, No Iterations 33
bounding omega, min: -7.36369332318e+32 max: 7.45087842067e+34 average: 8.45192123834e+29
smoothSolver:  Solving for k, Initial residual = 0.0652221006548, Final residual = 9.86389746847e-09, No Iterations 27
bounding k, min: -7.95004417902e+19 max: 2.54187707242e+21 average: 2.49264192669e+16
ExecutionTime = 27.47 s  ClockTime = 28 s

Time = 10009

smoothSolver:  Solving for Ux, Initial residual = 0.707100864244, Final residual = 8.17380806564e-09, No Iterations 34
smoothSolver:  Solving for Uy, Initial residual = 0.707102511597, Final residual = 8.09615585829e-09, No Iterations 34
smoothSolver:  Solving for Uz, Initial residual = 0.677345826642, Final residual = 3.47383376751e-07, No Iterations 1000
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 4.84816984568e-09, No Iterations 18
limitTemperature temperature_constraints Lower limited 99354 (46.1574912892%) of cells
limitTemperature temperature_constraints Upper limited 10946 (5.08524970964%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.652473647695, Final residual = 5.7463368848e-05, No Iterations 16
time step continuity errors : sum local = 1.15109831238e+56, global = -9.67393479366e+41, cumulative = -9.67393477123e+41
smoothSolver:  Solving for omega, Initial residual = 0.71639185423, Final residual = 5.03738174817e-09, No Iterations 14
bounding omega, min: -9.40492911923e+29 max: 9.37100865931e+33 average: 4.10924824413e+29
smoothSolver:  Solving for k, Initial residual = 0.742849536333, Final residual = 6.31077894166e-09, No Iterations 13
bounding k, min: -9.91457445321e+18 max: 2.28501407243e+21 average: 3.36552987043e+16
ExecutionTime = 29.96 s  ClockTime = 30 s

Time = 10010

smoothSolver:  Solving for Ux, Initial residual = 0.629389455891, Final residual = 3.12801623638e-06, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.629388913105, Final residual = 1.98095241143e-06, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.925952525366, Final residual = 2.71850304045e-08, No Iterations 1000
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 8.7328752983e-09, No Iterations 31
limitTemperature temperature_constraints Lower limited 5227 (2.42833914053%) of cells
limitTemperature temperature_constraints Upper limited 149022 (69.2320557491%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.539866259316, Final residual = 4.9318115972e-05, No Iterations 9
time step continuity errors : sum local = 6.43187006542e+68, global = 1.16146455737e+54, cumulative = 1.16146455737e+54
smoothSolver:  Solving for omega, Initial residual = 0.000345869080579, Final residual = 9.58730961687e-09, No Iterations 6
bounding omega, min: -1.51533821371e+32 max: 1.49080821324e+34 average: 8.06873195342e+29
smoothSolver:  Solving for k, Initial residual = 0.000147693344574, Final residual = 7.64781028605e-09, No Iterations 5
bounding k, min: -3.03356951287e+19 max: 1.75336888299e+21 average: 6.53704323548e+16
ExecutionTime = 35.49 s  ClockTime = 36 s

Time = 10011

smoothSolver:  Solving for Ux, Initial residual = 0.66610414697, Final residual = 0.000621669998517, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.666103047247, Final residual = 0.000410915153699, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.632811666534, Final residual = 72.2468595028, No Iterations 1000
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 6.21515996553e-09, No Iterations 18
limitTemperature temperature_constraints Lower limited 57763 (26.8353077816%) of cells
limitTemperature temperature_constraints Upper limited 66936 (31.0968641115%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.39167792921, Final residual = 2.5286286111e-05, No Iterations 12
time step continuity errors : sum local = 5.46160355531e+83, global = 1.16422002872e+68, cumulative = 1.16422002872e+68
smoothSolver:  Solving for omega, Initial residual = 3.70809948994e-05, Final residual = 4.5090165572e-09, No Iterations 10
bounding omega, min: -3.26955250853e+39 max: 1.38304174851e+41 average: 1.06914348437e+36
smoothSolver:  Solving for k, Initial residual = 1.48407726826e-05, Final residual = 3.2271897245e-09, No Iterations 11
bounding k, min: -7.45554271837e+23 max: 8.11191480185e+28 average: 5.8911348066e+23
ExecutionTime = 41.1 s  ClockTime = 41 s

Time = 10012

smoothSolver:  Solving for Ux, Initial residual = 0.793459788924, Final residual = 7.71857459656e-09, No Iterations 37
smoothSolver:  Solving for Uy, Initial residual = 0.793459473542, Final residual = 7.71856327155e-09, No Iterations 37
smoothSolver:  Solving for Uz, Initial residual = 0.560225739981, Final residual = 7.75970306657e-09, No Iterations 39
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 2.5198167101e-07, No Iterations 1000
limitTemperature temperature_constraints Lower limited 118319 (54.9681765389%) of cells
limitTemperature temperature_constraints Upper limited 42348 (19.6738675958%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.712697558163, Final residual = 7.10567897157e-05, No Iterations 19
time step continuity errors : sum local = 1.18719938884e+91, global = -6.57209032417e+76, cumulative = -6.57209031253e+76
smoothSolver:  Solving for omega, Initial residual = 0.0454299588808, Final residual = 6.46035024839e-09, No Iterations 23
bounding omega, min: -1.84987790127e+39 max: 1.00064926973e+43 average: 9.17667453186e+37
smoothSolver:  Solving for k, Initial residual = 0.0385242789961, Final residual = 8.17381323816e-09, No Iterations 9
bounding k, min: -1.23972360477e+32 max: 1.08563089711e+33 average: 8.83847936331e+27
ExecutionTime = 43.64 s  ClockTime = 44 s

Time = 10013

smoothSolver:  Solving for Ux, Initial residual = 0.497624650435, Final residual = 4.87058969705e-09, No Iterations 29
smoothSolver:  Solving for Uy, Initial residual = 0.497903225512, Final residual = 4.86876696485e-09, No Iterations 29
smoothSolver:  Solving for Uz, Initial residual = 0.570478907314, Final residual = 9.52400499797e-09, No Iterations 30
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 8.24709699828e-09, No Iterations 26
limitTemperature temperature_constraints Lower limited 119153 (55.3556329849%) of cells
limitTemperature temperature_constraints Upper limited 65867 (30.600232288%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.743150135483, Final residual = 4.20822840176e-05, No Iterations 5
time step continuity errors : sum local = 9.2987226444e+95, global = 6.03048038874e+81, cumulative = 6.03041466784e+81
smoothSolver:  Solving for omega, Initial residual = 1.4758490331e-06, Final residual = 6.68083715593e-09, No Iterations 7
bounding omega, min: -3.6782010718e+42 max: 9.50602386695e+43 average: 1.17393202904e+39
smoothSolver:  Solving for k, Initial residual = 1.04960751006e-08, Final residual = 4.77378562102e-09, No Iterations 4
bounding k, min: -2.68770735334e+32 max: 4.51128912912e+33 average: 6.12559229196e+28
ExecutionTime = 44.5 s  ClockTime = 45 s

Time = 10014

smoothSolver:  Solving for Ux, Initial residual = 0.544102085423, Final residual = 8.9839774487e-09, No Iterations 23
smoothSolver:  Solving for Uy, Initial residual = 0.544070041353, Final residual = 8.98402588636e-09, No Iterations 23
smoothSolver:  Solving for Uz, Initial residual = 0.549530349327, Final residual = 5.34064788471e-09, No Iterations 23
smoothSolver:  Solving for e, Initial residual = 1, Final residual = 0.00211127002229, No Iterations 1000
limitTemperature temperature_constraints Lower limited 96275 (44.7270615563%) of cells
limitTemperature temperature_constraints Upper limited 118290 (54.9547038328%) of cells
limitTemperature temperature_constraints Unlimited Tmax 5000
Unlimited Tmin 1
GAMG:  Solving for p, Initial residual = 0.244574783005, Final residual = 1.80656793475e-05, No Iterations 19
time step continuity errors : sum local = 1.50876019421e+102, global = -1.32173121908e+87, cumulative = -1.32172518866e+87
smoothSolver:  Solving for omega, Initial residual = 1.43058998943e-05, Final residual = 6.0894017599e-09, No Iterations 4
bounding omega, min: -7.92551170334e+42 max: 3.92247524649e+43 average: 6.2522736871e+38
smoothSolver:  Solving for k, Initial residual = 3.55770947532e-07, Final residual = 4.62807131238e-09, No Iterations 2
bounding k, min: -9.74285631755e+32 max: 4.12723898701e+34 average: 2.58919841857e+29
ExecutionTime = 46.95 s  ClockTime = 47 s

Time = 10015

smoothSolver:  Solving for Ux, Initial residual = 0.51964422401, Final residual = 0.00660365318865, No Iterations 1000
smoothSolver:  Solving for Uy, Initial residual = 0.519658858041, Final residual = 0.00783449220059, No Iterations 1000
smoothSolver:  Solving for Uz, Initial residual = 0.319326144599, Final residual = 0.00035545345465, No Iterations 1000
[2] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[2] #1  Foam::sigFpe::sigHandler(int) at ??:?
[2] #2  ? in /lib64/libpthread.so.0
[2] #3  Foam::multiply(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
[2] #4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator*<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
[2] #5  Foam::fv::boundedConvectionScheme<double>::fvcDiv(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
[2] #6  ? at ??:?
[2] #7  ? at ??:?
[2] #8  __libc_start_main in /lib64/libc.so.6
[2] #9  ? at ??:?

@DevilX:
I tried the simpleFoam incompressible solver. This one gives nice convergence and time step continuity errors around 1e-4 as you can see in this log:

Code:

Time = 1497

smoothSolver:  Solving for Ux, Initial residual = 0.000100500888787, Final residual = 8.90568075304e-06, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.00100435845354, Final residual = 8.95753961279e-05, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.000578782900718, Final residual = 5.17059123177e-05, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.000307317720265, Final residual = 2.73690798359e-05, No Iterations 3
time step continuity errors : sum local = 0.0314680286054, global = 0.000215929883069, cumulative = 17.0819715067
smoothSolver:  Solving for omega, Initial residual = 0.000189409035749, Final residual = 1.53940593495e-05, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.000239689487637, Final residual = 2.09473499003e-05, No Iterations 6
ExecutionTime = 888.39 s  ClockTime = 889 s

Time = 1498

smoothSolver:  Solving for Ux, Initial residual = 0.000100341349694, Final residual = 8.89145140224e-06, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.00100288646037, Final residual = 8.94435032349e-05, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.000577551294907, Final residual = 5.16052069827e-05, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.000307355450312, Final residual = 2.80029736407e-05, No Iterations 3
time step continuity errors : sum local = 0.0321770726825, global = -0.000101536466035, cumulative = 17.0818699703
smoothSolver:  Solving for omega, Initial residual = 0.000186158849077, Final residual = 1.57109669246e-05, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.000239410721694, Final residual = 2.09196646412e-05, No Iterations 6
ExecutionTime = 888.97 s  ClockTime = 890 s

Time = 1499

smoothSolver:  Solving for Ux, Initial residual = 0.000100190439595, Final residual = 8.8780677202e-06, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.00100123785007, Final residual = 8.92963811531e-05, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.000576308406416, Final residual = 5.14990578231e-05, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.000306761427373, Final residual = 2.71315079717e-05, No Iterations 3
time step continuity errors : sum local = 0.0311601766532, global = 0.000192217775192, cumulative = 17.082062188
smoothSolver:  Solving for omega, Initial residual = 0.000192195856679, Final residual = 1.54182053449e-05, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.000239151556127, Final residual = 2.08927250343e-05, No Iterations 6
ExecutionTime = 889.55 s  ClockTime = 890 s

Time = 1500

smoothSolver:  Solving for Ux, Initial residual = 0.000100031217945, Final residual = 8.86382202319e-06, No Iterations 6
smoothSolver:  Solving for Uy, Initial residual = 0.000999777848919, Final residual = 8.9165314016e-05, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.000575137221727, Final residual = 5.13952295765e-05, No Iterations 6
GAMG:  Solving for p, Initial residual = 0.000307531956787, Final residual = 2.76136481807e-05, No Iterations 3
time step continuity errors : sum local = 0.0316947674859, global = -5.90868822675e-05, cumulative = 17.0820031012
smoothSolver:  Solving for omega, Initial residual = 0.00019415877079, Final residual = 1.5487326225e-05, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 0.00023881617953, Final residual = 2.08642842142e-05, No Iterations 6
ExecutionTime = 890.13 s  ClockTime = 891 s


SIMPLE solution converged in 1500 iterations

Decreasing the relTol will make them even smaller.

From the simpleFoam solution I would expect max. velocities around 600 m/s. Fluent calculates a max. velocity of 542 m/s for the compressible case with the same flow rate. So the transonic option might be justified.

I will try the rhoPimpleFoam next and the waveTransmissive bc if I find some useful values for the parameters.

dlahaye · November 25, 2021, 14:51

We have found below to work well.
Here fieldInf 101325 is the initial value for the pressure on the domain.
Good luck!

Outlet
{
type waveTransmissive;
psi thermo

si;
lInf 1;
field p;
gamma 1.4;
fieldInf 101325;
value $internalField;
}

DevilX · November 26, 2021, 06:47

Okay, so SimpleFoam works good and accurate as necessary. The timeStepErrors are related to the pressure tolerance, so squeezing this will result in lower values.

For your problem with rhoSimple:
Your log FIle indicates, that your residuals are really high after 10000 Iterations, the turbulence values are extraordinary, so BC´s maybe the real issue here.
Take into consideration, that your BC´s are not perfect for these very high velocities, but i am not sure, which to choose. Give these fromDomenico a try.

My guess is, that at these velocities, the steady solver just comes to a problem. This case is so dependend on time, that neglectimg these terms actually causes divergence.

youknowme · March 30, 2022, 02:21

HI,
Could you please post your BCs and system folder of simpleFoam case?
i am too trying to simulate a pressure driven flow and i am kinda stuck because of bounding of k/epsilon
TIA

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Compressible Flow Through Rocket Nozzle	rootbeer_cigarettes	FLUENT	2	March 12, 2019 10:42
OpenFOAM v3.0+ ??	SBusch	OpenFOAM	22	December 26, 2016 15:24
Compressible Nozzle Flow	sebastian	OpenFOAM Running, Solving & CFD	14	September 21, 2016 11:47
Having Problem solving 2D supersonic flow around a plug nozzle	chrislloyd	FLUENT	7	July 22, 2015 14:09
Finding someone to help with some OpenFOAM code for compressible, nonisothermal flow?	bzz77	CFD Freelancers	2	December 15, 2014 15:36

November 24, 2021, 09:38		#2
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 802 Blog Entries: 1 Rep Power: 18	Greetings, rhoSimpleFoam solves conservation of mass, momentum and energy. The boundary conditions (inlet especially) should be compatible with the coupling of pressure, velocity and temperature. This makes a compressible rhoSimpleFoam case harder to set up than an incompressible simpleFoam case. The ideal gas law defines a relation between temperate, pressure and density (and thus velocity through mass flow rate). Fixing two quantities, the third one can be determined. This means that temperature, pressure and density can not be set (or user defined) independently from each other. In my case, e.g., it proved to be valuable to fix inlet pressure (e.g p = 101325 Pascal) and inlet temperature (e.g. T = 773.15 K) to arrive at density rho = p / (R_{specific} * T) = 101325 / (287.058 * 773.15) = 0.45654 kg/m^3 . Given density and area of the inlet, the inlet velocity can be determined. In your case set-up, I fail to understand the inlet Gauge pressure set to 0. My understanding is that you are very close to a working solution. Good luck.

November 24, 2021, 12:19		#3
flar.t New Member Join Date: Feb 2020 Posts: 6 Rep Power: 6	Thank you very much for your reply. The "Initial Gauge Pressure" is only specified in the fluent case. As far as I understand it is used only for initialization. The openFOAM case uses a zeroGradient boundary condition at the inlet. I thought with that BC the pressure itself is not specified and the ideal gas law is used as you described to calculate it from the fixed velocity and temperature. In the meanwhile I rebuilt an incompressible case from scratch and the simpleFoam solver converges really smoothly within 2000 iterations for the same BCs on pressure, velocity and turbulence as mentioned above. The moment I change from incompressible to compressible() the solver crashes. () by - adding thermophysicalProperties, - changing pressure from m^2/s^2 to kg/(ms^2) in 0/p, - adding rho to massFlowRate in 0/U, - adding 0/T and 0/alphat, - modifying divScheme for the shear-rate tensor (rho...), - adding the schemes and solver for the energy eq. I noticed that the convergence strongly depends on the gradient schemes for the incompressible case. While kEpsilon is fine with "Gauss linear" the kOmega works only when gradSchemes is set to "cellLimited leastSquares 1". Perhaps the crashes result from the wrong schemes or solver settings.. Are there any recommendations for compressible internal flows?

November 24, 2021, 12:50		#4
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 802 Blog Entries: 1 Rep Power: 18	My bad: pressure indeed has zeroGradient on the inlet patch as you suggest. What do you mean by "solver crashing"? What happens exactly? Can you share solver log file? What happens is case that rho on inlet is increased from 1 to 100, say (increasing density should decrease velocity at same mass flow rate)? Does this lowering of the inlet velocity make the case easier to converge?

November 24, 2021, 17:15		#6
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 802 Blog Entries: 1 Rep Power: 18	Do you require a waveTransmissive boundary condition for the pressure at the outlet? See https://openfoamwiki.net/index.php/H...dary_condition

November 25, 2021, 09:41		#7
DevilX Member Daniel Join Date: Jan 2021 Posts: 39 Rep Power: 5	Hey, you tried a lot, thats true. RhoSimpleFoam is kind of special and very sensitive. Give a transient calculation with rhoPimple a try, it's much more robust and versatile. Your log file indicates, that your timeStepErrors are really high ( around 10^7 bigger then they should), also you have high bounding values in k and omega. What velocities (maximium) are you expecting in your domain? Simplefoam should give you a value there. Is it really near the transonic region of your fluid (fvSolution -transonic - yes should always be off most times). What are your timeStepErrors and Residuals in SimpleFoam? Btw: this solver is so robust, it runs with most settings.

November 25, 2021, 14:51		#9
dlahaye Senior Member Domenico Lahaye Join Date: Dec 2013 Posts: 802 Blog Entries: 1 Rep Power: 18	We have found below to work well. Here fieldInf 101325 is the initial value for the pressure on the domain. Good luck! Outlet { type waveTransmissive; psi thermosi; lInf 1; field p; gamma 1.4; fieldInf 101325; value $internalField; }

November 26, 2021, 06:47		#10
DevilX Member Daniel Join Date: Jan 2021 Posts: 39 Rep Power: 5	Okay, so SimpleFoam works good and accurate as necessary. The timeStepErrors are related to the pressure tolerance, so squeezing this will result in lower values. For your problem with rhoSimple: Your log FIle indicates, that your residuals are really high after 10000 Iterations, the turbulence values are extraordinary, so BC´s maybe the real issue here. Take into consideration, that your BC´s are not perfect for these very high velocities, but i am not sure, which to choose. Give these fromDomenico a try. My guess is, that at these velocities, the steady solver just comes to a problem. This case is so dependend on time, that neglectimg these terms actually causes divergence. dlahaye and hogsonik like this.

March 30, 2022, 02:21		#11
youknowme New Member anonymous Join Date: Mar 2022 Posts: 2 Rep Power: 0	HI, Could you please post your BCs and system folder of simpleFoam case? i am too trying to simulate a pressure driven flow and i am kinda stuck because of bounding of k/epsilon TIA