[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...

[Aurelien Thinat]

Hi petr.f.,

Your mesh presents high value of cell's non-orthogonality. You should try adding some nNonOrthogonalCorrectors (something like 2 or 3 instead of 0).
You can also try to decrease a bit your schemes : limited 0,333 instead of limited 0,5.

Then if it fails, averaging your solution over few hundreds of iteration should do the trick.

If you kept the output log file, post it on the forum. It would be easier to spot any numerical problem (continuity error...).

[petr.f.]

Hi Aurelien,

thanks for the advice. I've tried the averaging and I got a reasonable result with standard deviation 0.006, which is not that bad. However, if I try to simulate the same case with StarCCM+, the oscillations are in the magnitude of 1e-3 (which is one order lower than what happens with simpleFoam).

I'll try the non-ortho correctors and decreasing the schemes as well.

The continuity error is quite OK (around 2e-7 in the last 5000 iterations).

Regards, P.

[RodriguezFatz]

Originally Posted by petr.f.

However, if I try tu simulate the same case with StarCCM+, the oscillations are in the magnitude of 1e-3 (which is one order lower than what happens with simpleFoam).

I guess that shows that StarCCM introduces ridiculous amounts of numerical dissipation.

[ssss]

Maybe your simulation is purely unsteady, so even if you use steady solvers you will have some sort of oscillations.

[Blanco]

Isn't the pressure residual a little bit too high?

Other residuals looks good, but the pressure res seems high in my opinion...can you say it is converged?

I had similar problems with external aerodyn. cases and I ofted had numerical oscillations in Cd and Cl when the nonOrthogonality of cells was high...and in those cases, I've always obtained high pressure residual.