[arashgmn]

I'm trying to simulate a rotating 2D airfoil in the wind. For that, I have defined a rotating region that encompasses my airfoil and a static region for the rest of the domain. To couple the two regions, I used cyclicAMI.

I managed to set up a working case using k-omega SST for turbulent flow. The result of p,U and k are reasonable. Yet, I see a large discontinuity in the coupled interface for omega and nut regardless of mesh refinement level. It seems that for these two parameters the solution in the static region is not really influenced by rotating and these fields (more or less) preserve their initial values Look at the pictures below.

p field (it's fine)
U field (it's fine)
k field (it's fine)
nut (it's NOT fine)
omega (it's NOT fine)


This happens when I use my own mesh generated by gmsh and not with snappyHexMesh utility. To put in more contrast, I made a mesh using snappyHexMesh and observed that this discontinuity doesn't occur. This showed me that whatever is causing this issue involves the create Baffle/Patchs.


UPDATE:
In the snappy case, the mesh had no boundary layer (for speedup). After I made a lighter case of the same problem using gmsh but with no boundary layer, the problem resolved somehow! it's just an observation. I'm not really suggesting that the boundary layer is causing the problem here.



How I make my mesh in gmsh?
To be able to define the sliding interface, I define a physical surface in gmsh (called AMI). This internal face has to be later baffled. The full list of physical surfaces defined in gmsh are

- AMI
- Inlet (type --> inlet)
- Outlet (type --> outlet)
- Airfoil (type --> wall)
- TopBottom
- FrontEnd (type --> empty)

I use gmshToFoam to make polyMesh folder. I later, edit the boundary file such that the patch types coincide with the expressions in the parenthesis above.

How I set up my BafflesDict?
Here is my baffle dict to transform the AMI into two patches.

internalFacesOnly true;

baffles
{	
    cyclicAMI
    {
        //- Use predefined faceZone to select faces and orientation.
        type        faceZone;
        zoneName    AMI;

        patches
        {
            master //fix
            {
                //- Master side patch
                name            AMI1;
                type            cyclicAMI;
                matchTolerance  0.0001;
                neighbourPatch  AMI2;
                transform       noOrdering;
                // Switch to zero-gradient for low weights
                lowWeightCorrection 0.2;
		inGroups(Rotating);
	     }

            slave //rotating
            {
                //- Slave side patch
                name            AMI2;
                type            cyclicAMI;
                matchTolerance  0.0001;
                neighbourPatch  AMI1;
                transform       noOrdering;
                lowWeightCorrection 0.2;
		inGroups(Stationary);
            }
        }
    }
}

I then run the following commands to apply this baffling:

createBaffles -overwrite; mergeOrSplitBaffles -split -overwrite

How the initial conditions are set up?
I first run a steady-state case with simpleFoam and later use its converged values as the initial condition for the transient-state solver (pimpleFoam). Initial conditions for the steady-state case (before baffling) have the following structure (this is actually for k, the rest are similar)

internalField   uniform 24;

boundaryField
{
    TopBottom
    {
        type            inletOutlet;
        inletValue      uniform 24;
        value           uniform 24;
    }
    Inlet
    {
        type            freestream;
        freestreamValue uniform 24;
        value           uniform 24;
    }
    Outlet
    {
        type            inletOutlet;
        inletValue      uniform 24;
        value           uniform 24;
    }
    FrontBack
    {
        type            empty;
    }
    Airfoil
    {
        type            kqRWallFunction;
        value           uniform 24;
    }
    AMI
    {
        type            slip;
    }
}

which after baffling becomes

internalField   uniform 24;

boundaryField
{
    FrontBack
    {
        type            empty;
    }
    TopBottom
    {
        type            inletOutlet;
        inletValue      uniform 24;
        value           uniform 24;
    }
    Inlet
    {
        type            freestream;
        freestreamValue uniform 24;
        value           uniform 24;
    }
    Outlet
    {
        type            inletOutlet;
        inletValue      uniform 24;
        value           uniform 24;
    }
    Airfoil
    {
        type            kqRWallFunction;
        value           uniform 24;
    }
    AMI1
    {
        type            cyclicAMI;
        value           uniform 0;
    }
    AMI2
    {
        type            cyclicAMI;
        value           uniform 0;
    }
}

Other relevant info

openFoam version: v1812

fvScheme

ddtSchemes
{
    default Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
    grad(U)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;//limitedLinear 1;
    div(phi,omega)  Gauss upwind;//limitedLinear 1;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited corrected 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

wallDist
{
    method meshWave;
}

fvSolution

solvers
{
    "pcorr.*"
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother	 GaussSeidel;
    }

    p
    {	
        $pcorr;
	tolerance        1e-7;
        relTol           0.01;
    }

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

    "(U|k|omega)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0.1;
    }

    "(U|k|omega)Final"
    {
        $U;
        tolerance       1e-08;
        relTol          0;
    }
/*
    cellDisplacement
    {
        solver          GAMG;
        tolerance       1e-5;
        relTol          0;
        smoother        GaussSeidel;
    }
*/
}


PIMPLE
{
    correctPhi               yes;
    nOuterCorrectors    100; 
    nCorrectors            2;  
    nNonOrthogonalCorrectors 2;

    residualControl
    {
    p
 	{
        relTol          0.001;
        tolerance       1e-7;
	}
    "(U|k|omega)"
	{
	relTol		0.01;
	tolerance	1e-7;
	}
    }
}

relaxationFactors
{
    fields
    {
        p               0.9;
    }
    equations
    {	
	"(U|k|omega)"   		0.7;
        "(U|k|omega)Final"		1.0;
    }
}

cache
{
    grad(U);
}

controlDict

application     pimpleFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         1;

deltaT          0.001;

writeControl    adjustableRunTime;

writeInterval   0.01;

purgeWrite      0;

writeFormat     ascii;

writePrecision  7;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

adjustTimeStep  yes;

maxCo           5;

dynamicMeshDict

dynamicFvMesh   dynamicMotionSolverFvMesh;

motionSolver    solidBody;

cellZone        Rotating;

solidBodyMotionFunction  oscillatingRotatingMotion;

origin        (0.5 0 0);
amplitude     (0 0 28); // degree [it's weird though that this is not in rad]
omega         1.884955592; // 2*Pi*f_mean  rad/s

I appreciate your help to fix that explains the reason and the remedy for the discontinuity.

Arash

[HPE]

>> This happens when I use my own mesh generated by gmsh and not with snappyHexMesh utility.

Might be related to the fact that the type of cells is tetrahedron in GMSH? I think your observation suggests that. If you can create a hexahedron-based mesh in GMSH and test it - you can isolate the matter further.

[arashgmn]

Originally Posted by HPE

>> This happens when I use my own mesh generated by gmsh and not with snappyHexMesh utility.

Might be related to the fact that the type of cells is tetrahedron in GMSH? I think your observation suggests that. If you can create a hexahedron-based mesh in GMSH and test it - you can isolate the matter further.

Meanwhile, I tried using a coarse (again generated with gmsh) but without a boundary layer for faster simulation. Turned out the problem disappears! I checked it in both v1912 and v1812. Also, now I remember I turned off the boundary layer in snappyHexMesh as well (again for speedup) and saw no problem.

My observations hint that the boundary layer somehow is involved but I don't really see any link! I'm very confused!

[arashgmn]

So I found the reason and the solution.

Only omega and nut seem wrong. This means they affect each other somehow. As nut is not computed indirectly in k-omega SST, we can conclude that the computation of omega is wrong.

I realized the omega field outside of the inner circle, is constant by time. It didn't evolve. This led me to the fact that maybe the tolerance for omega is too high which was the case. After setting a tolerance of 1e-9 for omega, everything seemed to work normally.

Take home message
In turbulent fluids, the omega will change mostly around the wall. So if the wall boundaries are only a small fraction of your total cells, the residuals might be dominated by the almost-constant areas of the domain. Hence, it's probably a good practice to choose a small tolerance for omega.