[bentkj]

Hi everyone, I am attempting to run a simulation of a commerical jet (NASA common research model) at 37,000ft where the flow is mildly compressible with a freestream Mach = 0.2, but local regions are observed accelerate up to M~0.5 based on simpleFoam results of the same setup. Because of this local acceleration in the flow field exceeding M= 0.3, I am trying to get a solution to account for compressibility with the rhoSimpleFoam solver. However, after setting it up and attempting to run it, rhoSimpleFoam returns an error:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 4.1
Exec   : rhoSimpleFoam
Date   : Jan 22 2018
Time   : 16:10:40
Host   : "ben-ubuntu1604"
PID    : 42862
Case   : /home/ben/OpenFOAM/ben-4.1/run/CRM_WingWakeImpingement/CRM-Unstruct-CLmax-12deg-RANS-SA-rhoSimpleFoam
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

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

Create mesh for time = 0


SIMPLE: convergence criteria
    field p	 tolerance 1e-07
    field U	 tolerance 1e-07
    field nuTilda	 tolerance 1e-07
    field e	 tolerance 1e-07

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 SpalartAllmaras
Selecting patchDistMethod meshWave
SpalartAllmarasCoeffs
{
    sigmaNut        0.66666;
    kappa           0.41;
    Cb1             0.1355;
    Cb2             0.622;
    Cw2             0.3;
    Cw3             2;
    Cv1             7.1;
    Cs              0.3;
}

No MRF models present

No finite volume options present


Starting time loop

wallShearStress wallShearStress:
    processing all wall patches

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 0.989842, Final residual = 0.0636575, No Iterations 12
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 0.0615403, No Iterations 12
smoothSolver:  Solving for Uz, Initial residual = 0.996507, Final residual = 0.0711557, No Iterations 12
smoothSolver:  Solving for e, Initial residual = 0.757366, Final residual = 0.054747, No Iterations 3
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0389688, No Iterations 9
time step continuity errors : sum local = 1.55654e-08, global = -2.36362e-18, cumulative = -2.36362e-18
rho max/min : 0.543587 0.121896
smoothSolver:  Solving for nuTilda, Initial residual = 0.999812, Final residual = 0.0697661, No Iterations 4
ExecutionTime = 217.24 s  ClockTime = 218 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.00316191, Final residual = 0.000227228, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.0822842, Final residual = 0.0060974, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.0018152, Final residual = 0.000137586, No Iterations 4
smoothSolver:  Solving for e, Initial residual = 0.00229383, Final residual = 0.0001295, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.952563, Final residual = 0.0461374, No Iterations 3
time step continuity errors : sum local = 3.33308e-07, global = -1.90393e-18, cumulative = -4.26755e-18
rho max/min : 17.3624 0
smoothSolver:  Solving for nuTilda, Initial residual = 0.00043298, Final residual = 3.11907e-05, No Iterations 4
ExecutionTime = 322.94 s  ClockTime = 323 s

Time = 3

smoothSolver:  Solving for Ux, Initial residual = 0.00067038, Final residual = 5.0315e-05, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.00196303, Final residual = 0.000150027, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.000252499, Final residual = 1.9195e-05, No Iterations 4

#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::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#4  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#5  ? at ??:?
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? at ??:?

I understand that there seems to be an attempt to divide by zero somewhere hinting an incorrect BC setup but I cannot uncover what's wrong with my set up. I've also attached my 0/, fvSchemes, fvSolutions, and polyMesh/boundary file below. Because of the 5file upload limit I've consolidate the 0/* files as one 0folder.txt.

To avoid confusion, I've named by boundaries to utilize the wildcard (".*") for the boundaries and walls of the geometry. Unfortunately, my mesh file is large (>1gb) and sharing it is difficult. The mesh was created with ICEM (unstructured) and imported into openfoam with the fluent3Dmeshtofoam utility. checkMesh returns a 'Mesh OK' result and simpleFoam runs it well.

If it matters, I am also running OFv4.1 on Ubuntu1604.

Any input is really much appreciated as this stage is vital for my research! Thank you very much in advance!!

[rcastilla]

Hello,
In my experience, rhoSimpleFoam is quite unstable even for moderate Mach numbers. I would suggest to try with rhoCentralFoam although it is a transient solver. It is density based and it is much more stable for high compressibility.

Best regards

Robert

[jurich]

Hello Ben,

If you look at Time = 2, density max/min go to 17.3 and 0. That's your divide by zero. You can add rhoMax and rhoMin parameters to your fvSolution file to limit this and let the solution get started, but it indicates an error in your problem set up.

In your U file you have the internalField velocity set to the same value as your free stream boundaries. In the first couple of solver steps, that will produce high gradients as it interacts with your geometry. Especially since you have second order schemes. You may want to cut that to ~1/10 of the value, or even 0.

Take a look at the rhoSimpleFoam squareBend tutorial for solver settings; your settings seem more like simpleFoam settings.

Best of luck.

[bentkj]

Originally Posted by jurich

Hello Ben,

If you look at Time = 2, density max/min go to 17.3 and 0. That's your divide by zero. You can add rhoMax and rhoMin parameters to your fvSolution file to limit this and let the solution get started, but it indicates an error in your problem set up.

In your U file you have the internalField velocity set to the same value as your free stream boundaries. In the first couple of solver steps, that will produce high gradients as it interacts with your geometry. Especially since you have second order schemes. You may want to cut that to ~1/10 of the value, or even 0.

Take a look at the rhoSimpleFoam squareBend tutorial for solver settings; your settings seem more like simpleFoam settings.

Best of luck.

Thanks for your suggestion Mr Joseph. Like you mentioned, limiting the density (rho) does help progress the solution by a few more steps, but the problem lies with my setup. I'm still unable to discover what's wrong but will continue looking (squarebend tutorial runs fine).

Cutting the internalfield velocity to (0 0 0) also helps the simulation progress and changes the error message to:

Time = 3

smoothSolver:  Solving for Ux, Initial residual = 0.546285, Final residual = 0.00166209, No Iterations 8
smoothSolver:  Solving for Uy, Initial residual = 0.350791, Final residual = 0.0015506, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.356911, Final residual = 0.00232054, No Iterations 6
smoothSolver:  Solving for e, Initial residual = 0.923695, Final residual = 0.000999545, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.720858, Final residual = 0.00235791, No Iterations 5
time step continuity errors : sum local = 6.81735e-05, global = 1.32415e-18, cumulative = 7.3745e-19
rho max/min : 1 0.1
smoothSolver:  Solving for nuTilda, Initial residual = 0.0552017, Final residual = 0.000316008, No Iterations 4
ExecutionTime = 444.17 s  ClockTime = 445 s

Time = 4

smoothSolver:  Solving for Ux, Initial residual = 0.404296, Final residual = 0.00108116, No Iterations 8
smoothSolver:  Solving for Uy, Initial residual = 0.242708, Final residual = 0.00102003, No Iterations 6
smoothSolver:  Solving for Uz, Initial residual = 0.207327, Final residual = 0.00175429, No Iterations 6
smoothSolver:  Solving for e, Initial residual = 0.553096, Final residual = 0.000721111, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.697383, Final residual = 0.00574708, No Iterations 5
time step continuity errors : sum local = 0.000101322, global = 1.23099e-19, cumulative = 8.6055e-19
rho max/min : 1 0.1
smoothSolver:  Solving for nuTilda, Initial residual = 0.421538, Final residual = 0.00308192, No Iterations 4
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  log in "/lib/x86_64-linux-gnu/libm.so.6"
#4  Foam::nutUWallFunctionFvPatchScalarField::calcYPlus(Foam::Field<double> const&) const at ??:?
#5  Foam::nutUWallFunctionFvPatchScalarField::calcNut() const at ??:?
#6  Foam::nutWallFunctionFvPatchScalarField::updateCoeffs() at ??:?
#7  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::Boundary::evaluate() at ??:?
#8  Foam::RASModels::SpalartAllmaras<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correctNut(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#9  Foam::RASModels::SpalartAllmaras<Foam::EddyDiffusivity<Foam::ThermalDiffusivity<Foam::CompressibleTurbulenceModel<Foam::fluidThermo> > > >::correct() at ??:?
#10  ? at ??:?
#11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#12  ? at ??:?

From what it looks like seems there's some problems with the turbulence model setup (and possibly wallfunction?). For reference, I've initialised the turbulence parameters nut and nuTilda according to: https://turbmodels.larc.nasa.gov/spalart.html. I will need to double check.

[bentkj]

Originally Posted by rcastilla

Hello,
In my experience, rhoSimpleFoam is quite unstable even for moderate Mach numbers. I would suggest to try with rhoCentralFoam although it is a transient solver. It is density based and it is much more stable for high compressibility.

Best regards

Robert

Thank you for your suggestion Robert. is rhoCentralFoam valid even for mildly compressible (M=0.3 to 0.6) and turbulent flows? I will try it out. Thank you.

[rcastilla]

Originally Posted by bentkj

Thank you for your suggestion Robert. is rhoCentralFoam valid even for mildly compressible (M=0.3 to 0.6) and turbulent flows? I will try it out. Thank you.

Yes, at least I have used it for a wide range of Ma and turbulent flow.
Good luck
Robert

[maddy11]

Hi Robert,


I'm an Openfoam beginner. I'm modelling transonic flow around RAE2822 aerofoil. Rhocentralfoam seems to be running for almost 44 hours now. The deltaT seems to show e-09 when I actually gave e-04. Can you please tell me why this is happening?