[leke102]

Hello Foamers,

I am trying to make a simulation of atmospheric wind flow over urban areas. I am using handmade blockmesh and my checkMesh evalutation is totally okay. Wind tunnel dimensions are 18,5*4*2 so model also has same dimensions. I am using atmospheric boundary conditions for inlet outlet and top which is named as atmos. For ground I am using atmospheric rough wall function. Ground is divided by 2. Incoming flow is face with higher roughness and housing part has nearly 0 roughness. I am attaching my all files to the post.
My Question is,
simulation is crashing. I am having so huge epsilon,k and time step continuity values and error says bad math happened. I tried to change my whole bc and ic condition and I refined mesh so many times but I couldn't find solution. I am waiting your precious ideas and helps

U File

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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


internalField   uniform (3.28 0 0);

boundaryField
{
    inlet
    {
        type            atmBoundaryLayerInletVelocity;
        #include        "include/ABLConditions" 
    }
    outlet
    {
        type            zeroGradient;
        // type            inletOutlet;
        // inletValue      uniform (0 0 0);
        // value           uniform (6 0 0);
    }

    ground
    {
        type            fixedValue;
        // type            uniformFixedValue;
        // uniformValue    constant (0 0 0);
        value           uniform (0 0 0);
    }
    houses
    {
        type            fixedValue;
        // type            uniformFixedValue;
        // uniformValue    constant (0 0 0);
        value           uniform (0 0 0);
    }


    atmos
    {
        type            atmBoundaryLayerInletVelocity;
        #include        "include/ABLConditions" 
    }
    #include        "include/sideAtmos"


}

P file

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 101000;

boundaryField
{

    inlet
    {
        type            zeroGradient;
        value           uniform 101000;
    }

    outlet
    {
        type            fixedValue;
        meanValue       uniform 101000;
        inletValue      uniform 101000;
        value           uniform 101000;
    }

    ground
    {
        type            zeroGradient;
    }
    houses
    {
        type zeroGradient;
    }

    atmos
    {
        type            zeroGradient;
    }
    #include        "include/sideAtmos"

}

K file

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0.3;

boundaryField
{

    inlet
    {
        type            atmBoundaryLayerInletK;
        #include        "include/ABLConditions"
    }

    outlet
    {
        type            zeroGradient;
        // type            inletOutlet;
        // inletValue      $internalField;
        // value           $internalField;
    }

    ground
    {
        type            kqRWallFunction;
        value           uniform 0.0;
    }
    houses
    {
        type            kqRWallFunction;
        value           uniform 0.0;
    }

    atmos
    {
        type            atmBoundaryLayerInletK;
        #include        "include/ABLConditions"
    }

    #include        "include/sideAtmos"

}

Nut File

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField       uniform 0;

boundaryField
{

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            zeroGradient;
        // type           inletOutlet;
        // inletValue     $internalField;
        // value          $internalField;
    }

    ground
    {
        type            nutkAtmRoughWallFunction;
        z0              uniform 0.0068;
        value          $internalField;


    }
    houses 
    {
        type            nutkAtmRoughWallFunction;
        z0              uniform 0.00001;
        value          $internalField;

    }


    atmos
    {
        type            zeroGradient;
    }
    #include        "include/sideAtmos"

}

Epsilon File

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

//#include        "include/initialConditions"

internalField   uniform 0.016;

boundaryField
{
    #include        "include/sideAtmos"
    outlet
    {
        type            inletOutlet;
        inletValue      uniform 0.001;
        value           uniform 0.001;
    }

    inlet
    {
        type            atmBoundaryLayerInletEpsilon;
        #include        "include/ABLConditions"
    }

    ground
    {
        type            epsilonWallFunction;
        value           uniform 0.0;
    }
    houses
    {
        type            epsilonWallFunction;
        value           uniform 0.0;
    }

    atmos
    {
        type            atmBoundaryLayerInletEpsilon;
        #include        "include/ABLConditions"
    }

    
}

Abl Conditions

Uref                 3.28;
Zref                 0.6;
zDir                 (0 0 1);
flowDir              (1 0 0);
z0                   uniform 0.0068;
zGround              uniform 0.000000001;
value                $internalField;

FvSolution

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{

    p
    {
        solver              GAMG;
        relTol              0.1;
//        tolerance           1e-7;
        smoother            GaussSeidel;

    }


    "(U|k|epsilon)"
    {
        solver smoothSolver;
        smoother GaussSeidel;
        relTol 0.1;
        tolerance       1e-12;
        nSweeps 1;
    }

}

PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors  1;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefPoint   (0.00001 0.000001 300.00001);
    pRefValue   101000.0;

    hydrostaticInitialization yes;
    nHydrostaticCorrectors 5;
}

PISO {
 nCorrectors 2;
 nNonOrthogonalCorrectors 0;
 pRefCell 0;
 pRefValue 0;
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent      no;

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

relaxationFactors
{

    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        "(k|epsilon)"   0.7;
    }

}

FvSchemes

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         steadyState;   //Euler;  
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
    default         none;


    //div(phi,U)          Gauss linear;
    // div(phi,k)          Gauss linear;
    // div(phi,epsilon)    Gauss linear;
 /*  
    div(phi,U)          Gauss fixedBlended 0.9 linear linearUpwind grad(U);
    div(phi,k)          Gauss fixedBlended 0.9 linear linearUpwind grad(k);
    div(phi,epsilon)    Gauss fixedBlended 0.9 linear linearUpwind grad(epsilon);
*/


// first oder upwind
    div(phi,U)      Gauss upwind;   //
    div(phi,k)      Gauss upwind;//
    div(phi,epsilon)      Gauss upwind;//

/*
// mixed
    div(phi,U)          Gauss limitedLinear 0.01;
    div(phi,k)          Gauss limitedLinear 0.01;
*/

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

    
}

laplacianSchemes
{
   // default         Gauss linear corrected;
   default         Gauss linear limited corrected 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
   // default         corrected;
   default         limited corrected 0.5; // 
}
fluxRequired
{
    default         no;
    p_rgh;
}

wallDist
{
    method meshWave;    //meshWave; Poisson; advectionDiffusion;

}

Error is in log file but to make it easier I am adding also here

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in /lib/x86_64-linux-gnu/libc.so.6
#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::lduMa trix::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::dict ionary const&) at ??:?
#7 Foam::fvMatrix<double>::solveSegregatedOrCoupled(F oam::dictionary const&) at ??:?
#8 Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
#9 ? in /opt/OpenFOAM/OpenFOAM-v1806/platforms/linux64GccDPInt32Opt/bin/simpleFoam
#10 __libc_start_main in /lib/x86_64-linux-gnu/libc.so.6
#11 ? in /opt/OpenFOAM/OpenFOAM-v1806/platforms/linux64GccDPInt32Opt/bin/simpleFoam


Best Regards,
Ibrahim

[Svetlana]

Does it crash immediately or after some number of successful iterations?
What happens if you set the inlet velocity to zero?

[leke102]

It crashes after 10 not succesfull iterations. I changed gradschemes from leastsquare to cell limited it is working but no converging.

If I set inlet velocity zero I dont think I will have huge gradients but I wont have correct result also.

Best Regards,
Ibrahim

[Svetlana]

Perhaps simplify it to "flow around a rectangular block in a rectangular wind tunnel" while applying the same BCs, and if the issue still occurs, post the complete case here as an attachment.