[JackW]

Hi all,

I have a puzzling issue implementing a solver that incorporates Daly & Harlow's GGDH model in the energy Equation. I have updated the EEqn.H of rhoSimpleFoam to read:

volScalarField fluxTimeScale (
        0.3 * turbulence->mut() / (0.09 * rho * turbulence->k())
    );

    volSymmTensorField alphaEff(
        "alphaEff",
        rho*fluxTimeScale*turbulence->R()
    );

    fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
        //   - fvm::laplacian(turbulence->alphaEff(), he) // linear
      - fvm::laplacian(alphaEff, he) // non-linear
      - fvm::laplacian(turbulence->alpha(), he)

     ==
        fvOptions(rho, he)
    );

I am using the kOmegaSST model, so k(), R() and mut() all return well defined values.

In serial everything runs, but in parallel I get:

[1] 
[1] 
[1] --> FOAM FATAL ERROR: 
[1] Maximum number of iterations exceeded
[1] 
[1]     From function thermo<Thermo, Type>::T(scalar f, scalar T0, scalar (thermo<Thermo, Type>::*F)(const scalar) const, scalar (thermo<Thermo, Type>::*dFdT)(const scalar) const, scalar (thermo<Thermo, Type>::*limit)(const scalar) const) const
[1]     in file /home/jack/OpenFOAM/OpenFOAM-2.4.0/src/thermophysicalModels/specie/lnInclude/thermoI.H at line 76.
[1] 
FOAM parallel run aborting

I am using sensibleEnthalpy, perfectGas and const for transport. As I run the simulation I watch the h residuals rise and I get a crash after ~5-15 iterations. This is obviously the cause as h is proportional to T. Note I have also tried solving for e.

If I look at the flow in paraview, I see small regions (1 or 2 cells) of |T|>3e+8 (should be 297-330).

Note, this behaviour is reasonably fvScheme/fvSolution dependent and if I change the number of procs and decomposition. Also, I have tried from a converged solution from using the standard heat flux model and an initial condition.

Note that I have two timeVaryingMappedFixedValue boundary conditions.

if anyone can help - thanks! I am guessing it is something to do with the processor boundaries and using the turbulence member functions, but I do not know.

Cheers,

Jack

[mechkween]

Hi Jack,


Did you manage to find a solution to this problem? I've also implemented a GGDH in the EEqn.H and my simulation crashes stating a floating point exception.


Cheers