[mmkr825]

Hi Foamers,
I have been stuck with the Floating point exception error. After hacking old discussions about this error, i fixed "unset FOAM_SIGFPE", but still I am running into another Problem.

Error:

Time = 3nl [0 0 -1 0 0 0 0] 0.005626
DILUPBiCG: Solving for Ux: solution singularity
DILUPBiCG: Solving for Uy: solution singularity
DILUPBiCG: Solving for T: solution singularity
DICPCG: Solving for p: solution singularity
time step continuity errors : sum local = nan, global = nan, cumulative = nan
ExecutionTime = 0.33 s ClockTime = 0 snl [0 0 -1 0 0 0 0] 0.005626

--> FOAM FATAL IO ERROR:
wrong token type - expected Scalar, found on line 3 the word 'nan'

file: /home/malli_reddy/OpenFOAM/malli_reddy-2.1.1/run/caseFile/system/data::solverPerformance::U at line 3.

From function operator>>(Istream&, Scalar&)
in file lnInclude/Scalar.C at line 91.

FOAM exiting

Can any one help to solve this problem. Thanks in advance.

Regards
M Reddy

[Bernhard]

Floating point errors are often a sign for problems with your case or set-up. Ignoring the error does not make the prolem go away. You should try to find out what the problem is with your case.

[mmkr825]

Hi Bernhard,
Thanks very much for quick reply. I am dealing with suspension flow (particles present in liquid media) and i prepared my solver by editing simpleFoam solver. I think my case and set up is OK, because when i give good initial guess for velocity (by running my case with another solver which is having only one liquid phase), the case is running very well and the results are accurate.
But when i run directly for suspension fluid with out giving good initial guess, i got the above mentioned error. Could you please tell me how to fix this problem.

Thanks
M Reddy.

[Bernhard]

My first guess would be, that somewhere in your own model, you divide by derivatives of velocities (strain rate or something). Without a good initialization, these derivatives are zero, causing the FPE. Of course when initializing with a real solution for the velocity field, these derivatives will be non-zero.

[mmkr825]

Hi Bernhard,
Your guess is correct. In my solver i am using strain rate term in my solver. Can you please help how to fix this without giving good initial guess. Should i change any fvschemes.

Thanks
M Reddy.

[Bernhard]

You cannot solve this via fvSchemes.

You have three options.
1. Start from a solved velocity field as initial guess. This is a perfectly fine approach.
2. Adopt your solve, such the flow field gets solved once, before the first iteration.
3. Where you divide by the strain rate, at a very small number. 1/1e-20 is fine, while 1/0 is not.

[mmkr825]

Hi Bernhard,
In my solver i have already added a small non-local term(a small constant) to local shear rate. Actually there are two models which gives the flow Physics of suspension fluid. For the model-1 your approach is working fine, but for model-2 it is not working. Even though the good initial guess in velocity has no effect.
Should i add anything to .bashrc to ignore such operation (divide by zero).

Thanks
m Reddy

[Bernhard]

Adjusting your .bashrc will certainly not be a solution to everything.

You should give more details on the models to see how the problem can be solved. You can easily find out which line is giving you the error (if not by the error itself, then by the infamous 'Info << "Foo 1" << endl;')

[mmkr825]

Hi Bernhard,
Thank you Very much for your kindness. I am attaching the links for the solver and case file.

Solver : https://www.dropbox.com/s/tm8i4y2u0a...M25Foam.tar.gz


Casefile: https://www.dropbox.com/s/rxro15ga6v...aseFile.tar.gz


My case is simple 2D rectangular pipe. Can you suggest me any tutorial to fix this problem.

Thanks
M Reddy

[Bernhard]

Maybe you can paste relevant pieces of the code here as well? Just a few line, that would make things easier.

[mmkr825]

Hi Bernhard,
very much thanks for your valuable time. In my solver, totally there are three fields (U, p and T).

Momentum Equation:

tmp<fvVectorMatrix> UEqn
    (
      - fvm::laplacian(2*nu*ns, U)
      + fvc::div(Epnn1)
    );

    UEqn().relax();

    sources.constrain(UEqn());

    solve(UEqn() == -fvc::grad(p));

Continuity Equation :

{
    p.boundaryField().updateCoeffs();
    #include "TEqn.H"
    volScalarField rAU(1.0/UEqn().A());
    U = rAU*UEqn().H();
    UEqn.clear();
    
    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
    adjustPhi(phi, U, p);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAU, p) == fvc::div(phi)
        );

        pEqn.setReference(pRefCell, pRefValue);

        pEqn.solve();

        if (simple.finalNonOrthogonalIter())
        {
            phi -= pEqn.flux();
        }
    }

    #include "continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U -= rAU*fvc::grad(p);
    U.correctBoundaryConditions();
    sources.correct(U);
}

T prediction:

solve
    (
        fvm::div(phi, T)
     + fvc::div(J)
    );

The parameters are:

label leftWall = mesh.boundaryMesh().findPatchID("leftWall");
// leftWall
forAll(mesh.boundaryMesh()[leftWall],i)
{
//std::cout<<i<<"\n";
Q.boundaryField()[leftWall][i].component(tensor::XX) = 1; //or any other values
Q.boundaryField()[leftWall][i].component(tensor::XY) = 0;
Q.boundaryField()[leftWall][i].component(tensor::XZ) = 0;
Q.boundaryField()[leftWall][i].component(tensor::YX) = 0;
Q.boundaryField()[leftWall][i].component(tensor::YY) = 0.8;
Q.boundaryField()[leftWall][i].component(tensor::YZ) = 0;
Q.boundaryField()[leftWall][i].component(tensor::ZX) = 0;
Q.boundaryField()[leftWall][i].component(tensor::ZY) = 0;
Q.boundaryField()[leftWall][i].component(tensor::ZZ) = 0.5;
}


label fixedWalls = mesh.boundaryMesh().findPatchID("fixedWalls");
// fixedWalls
forAll(mesh.boundaryMesh()[fixedWalls],i)
{
//std::cout<<i<<"\n";
Q.boundaryField()[fixedWalls][i].component(tensor::XX) = 1; //or any other values
Q.boundaryField()[fixedWalls][i].component(tensor::XY) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::XZ) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::YX) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::YY) = 0.8;
Q.boundaryField()[fixedWalls][i].component(tensor::YZ) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::ZX) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::ZY) = 0;
Q.boundaryField()[fixedWalls][i].component(tensor::ZZ) = 0.5;
}


label rightWall = mesh.boundaryMesh().findPatchID("rightWall");
// rightWallvvv
forAll(mesh.boundaryMesh()[rightWall],i)
{

Q.boundaryField()[rightWall][i].component(tensor::XX) = 1; //or any other values
Q.boundaryField()[rightWall][i].component(tensor::XY) = 0;
Q.boundaryField()[rightWall][i].component(tensor::XZ) = 0;
Q.boundaryField()[rightWall][i].component(tensor::YX) = 0;
Q.boundaryField()[rightWall][i].component(tensor::YY) = 0.8;
Q.boundaryField()[rightWall][i].component(tensor::YZ) = 0;
Q.boundaryField()[rightWall][i].component(tensor::ZX) = 0;
Q.boundaryField()[rightWall][i].component(tensor::ZY) = 0;
Q.boundaryField()[rightWall][i].component(tensor::ZZ) = 0.5;
}

// assign Q on internal Field
/////////////////////////////////////////////////////////
forAll(mesh.cells(),j)
{
Q[j].component(tensor::XX) = 1;
Q[j].component(tensor::XY) = 0;
Q[j].component(tensor::XZ) = 0;
Q[j].component(tensor::YX) = 0;
Q[j].component(tensor::YY) = 0.8;
Q[j].component(tensor::YZ) = 0;
Q[j].component(tensor::ZX) = 0;
Q[j].component(tensor::ZY) = 0;
Q[j].component(tensor::ZZ) = 0.5;
}
//////////////////////////////////////////////////////////*/

//Viscosity of the suspending fluid dimensionlessed by suspending fluid viscosity (nu)

volScalarField ns
    (
        IOobject
        (
            "ns",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        pow( (1.0- (T/Tm)), -1.82)
    );

//////////////////////////////////////////////////////////

//Particle shear stress np=ns-1;

volScalarField np
    (
        IOobject
        (
            "np",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        ns-1
    );

//////////////////////////////////////////////////////////

//calculating normal stress viscosity nN=(Kn*pow(T/Tm,2)*pow(1-T/Tm,-2));

volScalarField nN
    (
        IOobject
        (
            "nN",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        Kn*pow(T/Tm,2)*pow(1-T/Tm,-2)
    );
//////////////////////////////////////////////////////////

//calculating hindrance mobility function f(T)=(1-T/Tm)*pow(1-T, c-1), where c=[2,4]

volScalarField f
    (
        IOobject
        (
            "f",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        (1-T/Tm)*pow(1-T, c-1)
    );

/////////////////////////////////////////////////////////


//calculating local shear rate g(gamma dot) given by sqrt(2E:E)

//volScalarField g(1.41421*sqrt(magSqr(symm(fvc::grad(U))))+nl);

// where nl is the non local shear rate

volScalarField g
    (
        IOobject
        (
            "g",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        sqrt(magSqr(symm(fvc::grad(U))))+nl
    );

/////////////////////////////////////////////////////////

//calculating particle normal stress diagonal tensor Epnn1=nu*nN*g*Q;

volTensorField Epnn1
    (
        IOobject
        (
            "Epnn1",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        nu*nN*g*Q
    );
/////////////////////////////////////////////////////////


//volVectorField J(2*a*a/(9*nu)*pow(1-T,c)*(fvc::div(Ep)));

volVectorField J
    (
        IOobject
        (
            "J",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        -2*a*a/(9*nu)*f*fvc::div(Epnn1)+fvc::laplacian(2*a*a/(9*nu)*f*2*nu*np, U)
    );

I prepared my solver by modifying simpleFoam solver. Thanks again.

Regards
M Reddy