[petr.f.]

Nice afternoon to all!

what is your experience with cl and cd behaviour during steady (simpleFoam) simulation of external vehicle aerodynamics?

The problem I face is that all the residuals are converged (almost don't change in time for many iterations) but cl and cd coefficients oscillate. And quite a lot - cl in <0.11. 0.195> and cd in <0.305, 0.355> (details in the attached pictures). I would expect some oscillations with k-omega SST model I use, but in 3rd or 4th decimal place, not in 2nd... I'm using snappyHexMesh with 10 prism layers on the car body and wheels, 8 on the chassis and 5 on mudguards with the following checkMesh output

    Overall domain bounding box (-20 -18 -0.2616) (60 1.03886e-32 17.7384)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-3.26662e-16 -4.62421e-16 -8.6276e-14) OK.
    Max cell openness = 4.55546e-15 OK.
    Max aspect ratio = 257.513 OK.
    Minimum face area = 5.435e-10. Maximum face area = 0.453756.  Face area magnitudes OK.
    Min volume = 2.4974e-11. Max volume = 0.302519.  Total volume = 25916.9.  Cell volumes OK.
    Mesh non-orthogonality Max: 83.8692 average: 7.62046
   *Number of severely non-orthogonal (> 70 degrees) faces: 695.
    Non-orthogonality check OK.
  <<Writing 695 non-orthogonal faces to set nonOrthoFaces
    Face pyramids OK.
 ***Max skewness = 7.92292, 99 highly skew faces detected which may impair the quality of the results
  <<Writing 99 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

The whole mesh has ~ 35e6 cells and covers only one half of the car. Average y+ is 5 on the body, 25 on chassis, 25 - 30 on the wheels (so OK for wall functions...) and I'm using pretty standard numerical schemes

ddtSchemes
{
    default             steadyState;
}

gradSchemes
{
    default                   cellLimited leastSquares 1.0;
    grad(U)                  cellLimited leastSquares 1.0;
}

divSchemes
{
    default                 none;
    
    div(phi,U)              bounded Gauss linearUpwindV grad(U);
    
     div(phi,k)              bounded Gauss upwind;
    div(phi,omega)          bounded Gauss upwind;
    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default                 Gauss linear limited 0.5;
}

interpolationSchemes
{
    default                 linear;
}

snGradSchemes
{
    default                 limited 0.5;
}

fluxRequired
{
    default                 no;
    p;
}

and solver settings

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           1e-3; // def. 1e-2
        smoother         DICGaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels      1;
        minIter          2;
        maxIter          50;
    }

    "(U|k|omega).*"
    {
        solver           smoothSolver;
        smoother         DILUGaussSeidel;
        tolerance        1e-8;
        relTol           0;//0.01;// def. 0.1;
        nSweeps          2;
    }

}

SIMPLE
{
    nNonOrthogonalCorrectors 0;

    residualControl
    {
        p               1e-3;
        U               1e-3;
        "(k|omega)"     1e-4;
    }
}

potentialFlow
{
    nNonOrthogonalCorrectors 20;
}

relaxationFactors
{
    fields
    {
        p               0.15;
    }

    equations
    {
        U              0.5;
        k               0.5;
        omega      0.5;
    }
}

cache
{
    grad(U);
}

The leastSquares are used due to mesh non-orthogonality, although in such a small number of non-orthogonal faces the effect is almost negligible. If I change to cellLimited Gauss linear 1; the cl and cd oscillations remain the same...