[Alasir]

Greetings. I'm a student attempting to simulate a few basic geometries and calculate their drag coefficients.

The specific case is a 3D grid of 1/4 of a sphere (0.1m diameter), floating in air (300K) thats going past with Re=10^4. I've used the 2D airfoil tutorial as a baseline. simpleFoam, Spalart Allmaras equations.

My problem is that after a relatively quick convergence at start, the residuals stops reducing. They fluctuate on the following scales:

U: 10^-5
P: 10^-3
NuTilde: 10^-5

And this is ultimately not enought to give me the expected Cd values.

Things I've tried in order to solve this, mostly from scanning different forum posts:
y+ verified to max 9.
Setting U, P, Nutilda tolerance to 1e-12, reltol to 0.05, 0.01, 0.05.
Setting simple's nNonOrthogonalCorrectors to 4.
Changed nut and nutilda freestream bds and internalstreams to be 1e-4, under the idea that 5*nu (1.57e^-5) would be better than the standard 0.14.
Changed grad(U) and grad(nuTilda) shemes to cellLimited leastSquares 1 from gauss linear.

Image of case:
http://imgur.com/a/oPcWJ

Boundary Conditions: (Note that Geo1 & Geo2 are the sphere's wall patches. Left and Bot are symmetryPlane patches.
nut:

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

internalField   uniform 0.0001;

boundaryField
{
    Inlet
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Outlet
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Top
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Right
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Left
    {
        type            symmetryPlane;
    }

    Bot
    {
        type            symmetryPlane;
    }
    Geo1
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }

    Geo2
    {
        type            nutUSpaldingWallFunction;
        value           uniform 0;
    }
}

Nutilda:

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

internalField   uniform 0.0001;

boundaryField
{
    Inlet
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Outlet
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Top
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Right
    {
        type            freestream;
        freestreamValue uniform 0.0001;
    }

    Left
    {
        type            symmetryPlane;
    }

    Bot
    {
        type            symmetryPlane;
    }
    Geo1
    {
        type            fixedValue;
        value           uniform 0;
    }

    Geo2
    {
        type            fixedValue;
        value           uniform 0;
    }
}

U:

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

internalField   uniform (1.57 0 0);

boundaryField
{
    Inlet
    {
        type            freestream;
        freestreamValue uniform (1.57 0 0);
    }

    Outlet
    {
        type            freestream;
        freestreamValue uniform (1.57 0 0);

    }

    Top
    {
        type            freestream;
        freestreamValue uniform (1.57 0 0);
    }

    Right
    {
        type            freestream;
        freestreamValue uniform (1.57 0 0);
    }

    Left
    {
        type            symmetryPlane;
    }

    Bot
    {
        type            symmetryPlane;
    }
    Geo1
    {
        type            noSlip;
    }

    Geo2
    {
        type            noSlip;
    }

}

P:

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

internalField   uniform 0;

boundaryField
{
    Inlet
    {
        type            freestreamPressure;
    }

    Outlet
    {
        type            freestreamPressure;
    }

    Top
    {
        type            freestreamPressure;
    }

    Right
    {
        type            freestreamPressure;
    }

    Left
    {
        type            symmetryPlane;
    }

    Bot
    {
        type            symmetryPlane;
    }
    Geo1
    {
        type            zeroGradient;
    }

    Geo2
    {
        type            zeroGradient;
    }

}

Transport properties:

transportModel  Newtonian;

rho             [1 -3 0 0 0 0 0] 1.18;

nu              [0 2 -1 0 0 0 0] 1.57e-05;

Turbulence Properties:

simulationType RAS;

RAS
{
    RASModel        SpalartAllmaras;

    turbulence      on;

    printCoeffs     on;
}

ControDict:

application     simpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         5000000;

deltaT          1;

writeControl    timeStep;

writeInterval   2000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

functions
{
    forceCoeffs
    {
        type        forceCoeffs;
        functionObjectLibs ( "libforces.so" );
        outputControl   timeStep;
        timeInterval    1;
        lRef        0.1;        //referanselengde
        log         yes;
        rho             rhoInf;      // Indicates incompressible
        rhoInf      1;       // Redundant for incompressible
        liftDir     (0 0 1);
        dragDir     (1 0 0);
        CofR        (0 0 0);  // Axle midpoint on ground
        pitchAxis   (0 1 0);

        magUInf     1.57;
        patches     ( Geo1 Geo2 );
        Aref        0.001964;        // referanseareal 0.00785
    }
}

fvSheme:

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    grad(U)         cellLimited leastSquares 1;
    grad(nuTilda)   cellLimited leastSquares 1;
    default         Gauss linear corrected;

}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div(phi,nuTilda) bounded Gauss linearUpwind grad(nuTilda);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

fvSolution

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-12;
        relTol          0.01;
        smoother        GaussSeidel;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         2;
        tolerance       1e-12;
        relTol          0.05;
    }

    nuTilda
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        nSweeps         3;
        tolerance       1e-12;
        relTol          0.05;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 4;
    pRefCell        0;
    pRefValue       0;

    residualControl
    {
        p               1e-5;
        U               1e-5;
        nuTilda         1e-5;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        nuTilda         0.7;
    }
}

Anyone got any ideas?

[Alasir]

UPDATE:

I have found a small anamoly in the solution:

http://imgur.com/a/oUYJz

It appears that there is a problem with the mesh here.

I've tried to refine the mesh, but that just gave me 3 such anamolies other places. Anyone got any ideas?

I've used gmsh to generate the mesh:

RefinementLevel=200;


ge=0.05;
gw=ge/RefinementLevel;
l=1.5;
h=0.3;
d=0.05;

//inlet
Point(1)={0, 0, 0, 0.5*ge};
Point(2)={0, h, 0, 1.5*ge};
Point(3)={0, 0, h, 1.5*ge};
Point(4)={0, h, h, 2*ge};

//Outlet
Point(5)={l, 0, 0, ge};
Point(6)={l, h, 0, 1.5*ge};
Point(7)={l, 0, h, 1.5*ge};
Point(8)={l, h, h, 2*ge};

//Geometry
Point(10)={2*l/10, 0, 0, gw};
Point(11)={2*l/10+d, 0, 0, gw};
Point(12)={2*l/10+d, d, 0, gw};
Point(13)={2*l/10+d, 0, d, gw};
Point(14)={2*l/10+2*d, 0, 0, 2*gw};

//Main Lines

//Length
Line(20)={1, 10};
Line(21)={14, 5};
Line(22)={2, 6};
Line(23)={3, 7};
Line(24)={4, 8};

//inlet
Line(25)={1, 3};
Line(26)={3, 4};
Line(27)={4, 2};
Line(28)={2, 1};

//Outlet
Line(29)={5, 7};
Line(30)={7, 8};
Line(31)={8, 6};
Line(32)={6, 5};

//Geometry
Circle(40)={10, 11, 12};
Circle(41)={10, 11, 13};
Circle(42)={12, 11, 13};
Circle(43)={12, 11, 14};
Circle(44)={13, 11, 14};

//Loops
Line Loop(50)={25, 26, 27, 28};//inlet
Line Loop(51)={29, 30, 31, 32};//outlet
Line Loop(52)={25, 23, -29, -21, -44, -41, -20};//bot
Line Loop(53)={-27, 24, 31, -22};//top
Line Loop(54)={28, 20, 40, 43, 21, -32, -22};//left
Line Loop(55)={-26, 23, 30, -24};//right

Line Loop(56)={40, 42, -41};//Geometry 1
Line Loop(57)={42, 44, -43};//Geometry 2

//Planes
Plane Surface(60)={50};
Plane Surface(61)={51};
Plane Surface(62)={52};
Plane Surface(63)={53};
Plane Surface(64)={54};
Plane Surface(65)={55};

Ruled Surface(66)={56};
Ruled Surface(67)={57};

//Naming
Physical Surface("Inlet")={60};
Physical Surface("Outlet")={61};
Physical Surface("Bot")={62};
Physical Surface("Top")={63};
Physical Surface("Left")={64};
Physical Surface("Right")={65};

Physical Surface("Geo1")={66};
Physical Surface("Geo2")={67};

//Volume
Surface Loop(70)={60, 61, 62, 63, 64, 65, 66, 67};
Volume(71)={70};
Physical Volume("Internal")={71};

[piu58]

Without analysing your case en détail:

simpleFoam is a steady state solver. You have a Reynolds number of 10,000. In this Re range we get a Kármán vortex street .

The iteration round in simple may be seen as pseudo time steps. In your case which is in reality transient it may be that the solver doesn't come to an end in your sense. You may select one or more monitoring point and look when the frequency of changing gets stable.