DES case - instability with Hexa Mesh

greg.cfd · July 11, 2014, 06:26

Hello!

I am simulating a flow around a wing with DES in OpenFoam. At first, I run a RANS simulation (with k-omega SST turbulence model) to obtain the steady state solution as a starting point for the DES case.

Then I use the DES hybrid method (LES with Spalart-Almaras as turbulence model for near wall regions).

I use a Tetra mesh with layers consisting of Penta elements and everything works fine. The average Courant Number is 0.05 for a time step 1e-04, while the maximum Courant Number is approximately 200. With these settings I get a converged solution (monitoring drag coefficient).

The problem appears when I use a Hexa Mesh for the domain (boundary layers are still Penta elements). The drag coefficient after some time steps fluctuates a lot, and when I visualize the velocity flow field, I get very high velocity values and generally uneven velocity distribution (instabilities) in front of the wing. The solution of course does not converge. I changed the time step to 1e-05 and the situation only got a bit better but instabilities were still there.

The following picture shows what I mean:

https://www.dropbox.com/s/2sutdexajm...Screenshot.png

What I'd like to ask is whether someone knows why a DES case runs without problems at all whith a Tetra mesh, while the same case with a Hexa dominant mesh shows this behaviour.

Here are the case files that I use for the DES simulation (controlDict, fvSchemes, fvSolution, checkMesh log, in case someone notices something wrong)

Thanks a lot in advance!
Greg

controlDict:

Code:

FoamFile
{
	version 2.0;
	format ascii;
	class dictionary;
	location "";
	object controlDict;
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

application pisoFoam;

startFrom startTime;

startTime	0.;

stopAt endTime;

endTime	7.;

deltaT 	0.00001;

writeControl timeStep;

writeInterval	10.;

purgeWrite	0;

writeFormat	ascii;

writePrecision	6;

writeCompression	compressed;

timeFormat	general;

timePrecision	6;

graphFormat	raw;

runTimeModifiable	yes;

functions {
cuttingPlane
	{
    		type            surfaces;
   		functionObjectLibs ("libsampling.so");
		outputControl   timeStep; //outputTime;
		outputInterval 1;

		surfaceFormat   vtk;
		fields          ( p U );

		interpolationScheme cellPoint;

		surfaces
		(
		    yNormal
		    {
		        type            cuttingPlane;
		        planeType       pointAndNormal;
		        pointAndNormalDict
		        {
		            basePoint       (0 0 0);
		            normalVector    (0 1 0);
		        }
		        interpolate     true;
		    }
		);
	}
}

fvSchemes:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default backward;
}

d2dt2Schemes
{
}

gradSchemes
{
    default         Gauss linear;

    grad(nuTilda)   cellLimited Gauss linear 1;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss LUST unlimitedGrad(U);
    //div(phi,U)      Gauss linearUpwind unlimitedGrad(U);
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,nuTilda) Gauss limitedLinear 1;

    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited corrected 0.33;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.33;
}

fluxRequired
{
    default         no;
    p;
}

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

fvSolution:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          0.1;

        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 50;//10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    };

    pFinal
    {
        $p;
        tolerance       1e-6;
        relTol          0;
    };

    "(U|k|B|nuTilda)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-07;
        relTol          0;
    };
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
    "U.*"               1;
    "nuTilda.*"         1;
}


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

checkMesh:

Code:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Enabling all geometry checks.

Time = 0

Mesh stats
    points:           16543491
    faces:            88503364
    internal faces:   86100313
    cells:            36717249
    faces per cell:   4.75536
    boundary patches: 52
    point zones:      0
    face zones:       2
    cell zones:       3

Overall number of cells of each type:
    hexahedra:     4900859
    prisms:        16039478
    wedges:        0
    pyramids:      1893485
    tet wedges:    0
    tetrahedra:    13883427
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking geometry...
    Overall domain bounding box (-11.1939 -3 -1.96509) (10.721 3 3.13891)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (4.75942e-15 1.09539e-16 -1.31463e-15) OK.
    Max cell openness = 9.08051e-15 OK.
    Max aspect ratio = 155.878 OK.
    Minimum face area = 2.0481e-09. Maximum face area = 0.0283969.  Face area magnitudes OK.
    Min volume = 7.41383e-14. Max volume = 0.004096.  Total volume = 339.357.  Cell volumes OK.
    Mesh non-orthogonality Max: 59.9987 average: 17.9918
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.99369 OK.
    Coupled point location match (average 0) OK.
 ***Error in face tets: 2 faces with low quality or negative volume decomposition tets.
  <<Writing 2 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
    Min/max edge length = 3.20285e-05 0.304629 OK.
   *There are 1 faces with concave angles between consecutive edges. Max concave angle = 19.7394 degrees.
  <<Writing 1 faces with concave angles to set concaveFaces
    Face flatness (1 = flat, 0 = butterfly) : min = 0.881177  average = 0.999822
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 0 average: 2.71935
 ***Cells with small determinant (< 0.001) found, number of cells: 139051
  <<Writing 139051 under-determined cells to set underdeterminedCells
    Concave cell check OK.
    Face interpolation weight : minimum: 0.0545556 average: 0.443561
    Face interpolation weight check OK.
    Face volume ratio : minimum: 0.0577036 average: 0.804768
    Face volume ratio check OK.

Failed 2 mesh checks.

End

kingjewel1 · March 14, 2016, 06:45

Hi Greg,

This is interesting and I'm finding similar behaviour. Did you have get to the bottom of it?

Cheers,
M

Quote:

Originally Posted by greg.cfd

Hello!

I am simulating a flow around a wing with DES in OpenFoam. At first, I run a RANS simulation (with k-omega SST turbulence model) to obtain the steady state solution as a starting point for the DES case.

Then I use the DES hybrid method (LES with Spalart-Almaras as turbulence model for near wall regions).

I use a Tetra mesh with layers consisting of Penta elements and everything works fine. The average Courant Number is 0.05 for a time step 1e-04, while the maximum Courant Number is approximately 200. With these settings I get a converged solution (monitoring drag coefficient).

The problem appears when I use a Hexa Mesh for the domain (boundary layers are still Penta elements). The drag coefficient after some time steps fluctuates a lot, and when I visualize the velocity flow field, I get very high velocity values and generally uneven velocity distribution (instabilities) in front of the wing. The solution of course does not converge. I changed the time step to 1e-05 and the situation only got a bit better but instabilities were still there.

The following picture shows what I mean:

https://www.dropbox.com/s/2sutdexajm...Screenshot.png

What I'd like to ask is whether someone knows why a DES case runs without problems at all whith a Tetra mesh, while the same case with a Hexa dominant mesh shows this behaviour.

Here are the case files that I use for the DES simulation (controlDict, fvSchemes, fvSolution, checkMesh log, in case someone notices something wrong)

Thanks a lot in advance!
Greg

controlDict:

Code:

FoamFile
{
	version 2.0;
	format ascii;
	class dictionary;
	location "";
	object controlDict;
}
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

application pisoFoam;

startFrom startTime;

startTime	0.;

stopAt endTime;

endTime	7.;

deltaT 	0.00001;

writeControl timeStep;

writeInterval	10.;

purgeWrite	0;

writeFormat	ascii;

writePrecision	6;

writeCompression	compressed;

timeFormat	general;

timePrecision	6;

graphFormat	raw;

runTimeModifiable	yes;

functions {
cuttingPlane
	{
    		type            surfaces;
   		functionObjectLibs ("libsampling.so");
		outputControl   timeStep; //outputTime;
		outputInterval 1;

		surfaceFormat   vtk;
		fields          ( p U );

		interpolationScheme cellPoint;

		surfaces
		(
		    yNormal
		    {
		        type            cuttingPlane;
		        planeType       pointAndNormal;
		        pointAndNormalDict
		        {
		            basePoint       (0 0 0);
		            normalVector    (0 1 0);
		        }
		        interpolate     true;
		    }
		);
	}
}

fvSchemes:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default backward;
}

d2dt2Schemes
{
}

gradSchemes
{
    default         Gauss linear;

    grad(nuTilda)   cellLimited Gauss linear 1;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss LUST unlimitedGrad(U);
    //div(phi,U)      Gauss linearUpwind unlimitedGrad(U);
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,nuTilda) Gauss limitedLinear 1;

    div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited corrected 0.33;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.33;
}

fluxRequired
{
    default         no;
    p;
}

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

fvSolution:

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          0.1;

        smoother        GaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 50;//10;
        agglomerator    faceAreaPair;
        mergeLevels     1;
    };

    pFinal
    {
        $p;
        tolerance       1e-6;
        relTol          0;
    };

    "(U|k|B|nuTilda)"
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-07;
        relTol          0;
    };
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
    "U.*"               1;
    "nuTilda.*"         1;
}


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

checkMesh:

Code:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Enabling all geometry checks.

Time = 0

Mesh stats
    points:           16543491
    faces:            88503364
    internal faces:   86100313
    cells:            36717249
    faces per cell:   4.75536
    boundary patches: 52
    point zones:      0
    face zones:       2
    cell zones:       3

Overall number of cells of each type:
    hexahedra:     4900859
    prisms:        16039478
    wedges:        0
    pyramids:      1893485
    tet wedges:    0
    tetrahedra:    13883427
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Cell to face addressing OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking geometry...
    Overall domain bounding box (-11.1939 -3 -1.96509) (10.721 3 3.13891)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (4.75942e-15 1.09539e-16 -1.31463e-15) OK.
    Max cell openness = 9.08051e-15 OK.
    Max aspect ratio = 155.878 OK.
    Minimum face area = 2.0481e-09. Maximum face area = 0.0283969.  Face area magnitudes OK.
    Min volume = 7.41383e-14. Max volume = 0.004096.  Total volume = 339.357.  Cell volumes OK.
    Mesh non-orthogonality Max: 59.9987 average: 17.9918
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.99369 OK.
    Coupled point location match (average 0) OK.
 ***Error in face tets: 2 faces with low quality or negative volume decomposition tets.
  <<Writing 2 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
    Min/max edge length = 3.20285e-05 0.304629 OK.
   *There are 1 faces with concave angles between consecutive edges. Max concave angle = 19.7394 degrees.
  <<Writing 1 faces with concave angles to set concaveFaces
    Face flatness (1 = flat, 0 = butterfly) : min = 0.881177  average = 0.999822
    All face flatness OK.
    Cell determinant (wellposedness) : minimum: 0 average: 2.71935
 ***Cells with small determinant (< 0.001) found, number of cells: 139051
  <<Writing 139051 under-determined cells to set underdeterminedCells
    Concave cell check OK.
    Face interpolation weight : minimum: 0.0545556 average: 0.443561
    Face interpolation weight check OK.
    Face volume ratio : minimum: 0.0577036 average: 0.804768
    Face volume ratio check OK.

Failed 2 mesh checks.

End

dickcruz · April 29, 2018, 02:50

Gentlemen, I have a question that's probably ill-informed. I'm not familiar with how one would

Quote:

...use the DES hybrid method (LES with Spalart-Almaras as turbulence model for near wall regions).

Could you help me understand this, please?

Quote:

Originally Posted by kingjewel1

Hi Greg,

This is interesting and I'm finding similar behaviour. Did you have get to the bottom of it?

Cheers,
M

Santiago · May 9, 2018, 11:18

Does your hybrid grid have pyramids in the interface between the boundary layer and the hexahedral region? Have you done checkmesh?

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
sliding mesh problem in CFX	Saima	CFX	46	September 11, 2021 08:38
[snappyHexMesh] SnappyHexMesh for internal Flow	vishwa	OpenFOAM Meshing & Mesh Conversion	24	June 27, 2016 09:54
[ICEM] Hexa mesh, curve mesh setup, bunching law	Anorky	ANSYS Meshing & Geometry	4	November 12, 2014 01:27
[ICEM] Hybrid mesh hexa and tetra on ICEM	PeterDaniel	ANSYS Meshing & Geometry	0	October 22, 2013 13:22
[ICEM] Problem making structured mesh on a surface	froztbear	ANSYS Meshing & Geometry	4	November 10, 2011 09:52

May 9, 2018, 11:18		#4
Santiago Senior Member Santiago Lopez Castano Join Date: Nov 2012 Posts: 354 Rep Power: 16	Does your hybrid grid have pyramids in the interface between the boundary layer and the hexahedral region? Have you done checkmesh?