OpenFOAM, residuals, and vortex street

This post appeared due to the question on cfd-online.com OpenFOAM forum. It has been started as yet another post on Kármán vortex street but then (around message #60) original poster revealed an article he was trying to compare with, so I've decided to check influence of different discretisation schemes and convergence criterions on final results.

Geometry and flow conditions

The first version of geometry and boundary conditions can be found in post #11, though they slightly differ from the geometry provided in the article (and as a result, from geometry used in simulations).

So we've got a cylinder with diameter D placed in the channel with a length 16D and a width 8D. At the inlet we have velocity 1 m/s, at the outlet zero pressure and zero gradient in velocity. Initially velocity inside the domain set to 1 m/s. With D = 1 m and fluid kinematic viscosity 100 m2/s, Reynolds number is equal to 100.

First attempt

Well, the first variant of the simulation was just a simulation for its own sake (even Re was 200). Just to produce neat pictures. And these videos:




Flow develops, vortices go.

Second attempt

After reading the original article I've decided to check different meshes, discretisation schemes, and convergence criteria to learn if there's any difference. In the article authors were developing new method in FEM framework, so theirs original mesh was triangular with increased density near a centerline of the channel. Now I think there's a need for another set of simulations: with and without increase of mesh density near the center-line of the channel.

Meshes

I've fallen in love with Gmsh for generating meshes, so these two meshed are produced with it. As Gmsh now has transfinite algorithm, it is more or less like blockMesh with GUI, variables and functions in a mesh file, more convenient way of grading description. Much like in blockMesh one needs to define vertices and then can use GUI to describe lines, surfaces and volumes. Then again switch to text editor for definition of transfinite lines, progressions (in the sense of Gmsh), etc. I will omit iterations during the mesh creation process and just provide GEO files with screenshots of the mesh:

"Coarse" mesh


Coarse cylinder-2D.geo

"Fine" mesh


Fine cylinder-2D.geo

The only difference between coarse and fine mesh is twice increased density in case of "Fine" mesh. So a number of cells in the fine mesh is 4 times higher than in coarse mesh.

fvSchemes

Typical fvSchemes file is provided below. The only line that is changed from simulation to simulation is div(phi,U) Gauss vanLeerV;, I've chosen to compare three different second order discretisation schemes with limiters: van Leer, QUICK, and Gamma.

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
   version     2.0;
   format      ascii;
   class       dictionary;
   location    "system";
   object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
   default         CrankNicolson 1.0;
}

gradSchemes
{
   default         Gauss linear;
}

divSchemes
{
   default         none;
   div(phi,U)      Gauss vanLeerV;
   div((nuEff*dev(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
   default         none;
   laplacian(nuEff,U) Gauss linear corrected;
   laplacian((1|A(U)),p) Gauss linear corrected;
}

interpolationSchemes
{
   default         linear;
}

snGradSchemes
{
   default         corrected;
}

fluxRequired
{
   default         no;
   p               ;
}


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

fvSolution

In case of fvSolution file I changed PIMPLE dictionary:

PIMPLE
{
    momentumPredictor yes;
    nCorrectors       2;
    nOuterCorrectors  20;
    nNonOrthogonalCorrectors 1;

    residualControl
    {
        "(U|p)"
        {
            tolerance 1e-2;
            relTol 0;
        }
    }
}

Keeping everything else the same, I changed tolerance in residualControl dictionary to see if it influence results. I've checked values 1e-2, 1e-4, and 1e-6. This value defines the moment when PIMPLE loop will exit from outer corrector cycle. And usually in the beginning of a simulation number of outer correctors is around 4 (for the case of residual 1e-6 - 8), during the time this number falls to 2 (or 4 for the case of 1e-6). I don't know why all tutorial cases use fixed number of outer corrector iterations, sometimes it can lead to completely incorrect results. Maybe OpenFOAM decided that residualControl is "advanced" technique.

Also there're lots of scary posts on relaxation, so I've decided also to check the influence of relaxation factors on results of a simulation. So instead of standard

relaxationFactors
{
    fields
    {
    }
    equations
    {
       "U.*" 1.0;
       "p.*" 1.0;
    }
}

I used

relaxationFactors
{
    fields
    {
    }
    equations
    {
       "U.*" 0.7;
       "p.*" 0.3;
    }
}

Results

Color contours usually provide only schematic view of the flow, so further (surely after these pictures) I will be comparing oscillations of force coefficients. Sometimes it is possible to plot isolines of pressure and velocity on the same picture (when the phases of the oscillating flow do not differ much). But now, the color contours (time is 150 s):

Pressure
van Leer


QUICK


Gamma


Though values seem to be more or less the same, vortices are shifted in time. The same thing is for velocity components:

Ux
van Leer


QUICK


Gamma


Uy
van Leer


QUICK


Gamma

Discretisation schemes comparison

I've fixed final residual in fvSolution at 1e-6 and just changed discretisation scheme. On the graphs X-axis is time, Y-axis is value. All simulations were done on the 'fine' mesh.

There is phase shift between graphs, I've moved them to compare oscillation amplitude. There is difference but it's not so big.

And as there is phase shift between different scheme I will omit comparison of the pressure and velocity contours. While there is difference between Cdrag oscillation amplitudes, Clift oscillations go rather synchronously:

Final residuals comparison

After that I've decided to check sensitivity of a discretisation scheme on final residual. I've fixed the scheme and run simulations for different final residuals: 1e-2, 1e-4, and 1e-6.
While van Leer and Gamma showed almost no dependence (though there is slight shift when going from 1e-2 to 1e-4 for Gamma):

For QUICK switch from 1e-2 to 1e-4 leads to quite large time shift between oscillations:

But even if the time shift in oscillations seems to be not so big, positions of the vortices are strongly depend on this shift. Below is overlays of Ux isolines for van Leer and Gamma schemes at 150 s:

If for van Leer scheme contours are almost the same, for Gamma results with residual 1e-2 are sufficiently shifted.

Mesh density comparison

Also compared results acquired using meshes with different densities. Number of cells in the meshes differ 4 times as I've reduced densities twice in X and Y directions.
Gamma showed just decrease of Cdrag and Clift values on fine mesh, while with QUICK and van Leer schemes there's also phase shift between results on different meshes.

Cdrag

Clift

To relax or not

Finally I've decided to check if a relaxation is really crucial to results. I've fixed the discretisation scheme (van Leer, though now I think maybe it was better to choose Gamma), the final residual in PIMPLE dictionary and just modified relaxation part as shown in the beginning of the post. Surely due to the relaxation solution converges slowly, i.e. instead of 4-5 outer corrector iterations 9-11 iterations was necessary. Also there is a time shift between results.

Conclusion

So it goes.

As usual case files (if you'd like to reproduce results) are on bitbucket.org.

Let's put something in a tube

Introduction

After I've done simulations of vortex street in 2D, I kept thinking that it would be nice to do more or less the same thing in 3D. And after watching the movie about vortex flow meters by Endress+Hauser, I've decided it would be really funny to put something in a tube and watch vortices going by.

First of all we need to decide on operational parameters of the system. Let's assume that we have 10 m long tube with diameter of 1 m and let's install something at 2 m from inlet plane. Further calculations are based upon the article from Wikipedia about Kármán vortex street.

Re = uD/ν

where u is flow velocity, D is characteristic length, and ν is fluid viscosity. We'll get vortex street for Re values larger than 90. If we take, for example, water with nu around 1000000 it's not hard to get large enough Reynolds number. But then a question of the distance between vortices arises. We also want our vortices inside the tube. Vortex shedding frequency is

f = 0.198*u/d*(1 - 19.7/Re)

where d is cylinder diameter (correct for 250 < Re < 10^5). After playing with tube and cylinder sizes, fluid viscosity and inflow velocities (here is a script) I've decided to go with the following:

D = 1 m
d = 0.2 m
u = 1 m/s
ν = 0.0001 m^2/s

Meshing

For meshing I've used Gmsh and snappyHexMesh. I.e. I've created hexagonal mesh of the tube, then I've created triangular surface mesh for the cylinder and sphere, a nd finally I've used snappyHexMesh to cut small cylinder (or sphere) tube.

Mesh for tube is based upon my previous post. The addition is a parametrization of the GEO file, i.e. one need to set tube leng th and its diameter; gmsh will recalculate positions of the points in the mesh. GEO files for cutting shapes (cylinder, sphere) are even simpler cause we need only triangular surface mesh so there's no need to split a shape into volumes with 6 sides (this operation necessary for transfinite algorithm to work).

One of the problems I've encountered during "cutting" procedure is the need to use more or less cubical mesh for sHM to work. In case of sphere it wan't actually a problem, cause it lies inside "inner" part of the mesh which is rather regular. So one can have graded mesh near the walls of the tube. That was not so in case on cylinder. After several unsuccessful attempts to cut a cylinder hole in graded mesh (though sHM tried to do it but wasn't able to generate acceptable quality mesh, usually non-orthogonality on the new mesh was too high), I've removed grading.
snappyHexMeshDict is quite simple, first I define the surface:

geometry
{
    sphere.stl
    {
        type triSurfaceMesh;
        name sphere;
    }
};

then I refine mesh around the surface

refinementRegions
{
    sphere
    {
        mode distance;
        levels ((0.05 3) (0.1 2) (0.2 1));
    }
}

refinement is based upon distance from the mesh. Other parameters in the dictionary is just a copy from tutorial.

So after running mesh conversion (from gmsh to OpenFOAM) and then running snappyHexMesh, I've got the following meshes for the simulation:

Initial tube mesh

with grading towards the walls

without grading

Mesh with internal cylinder

Mesh with internal sphere

Unlike previous time I've decided not to compare turbulence models and just selected realizable k-epsilon.

Flow movies

After running the simulations I've got the following flow movies:

The cylinder in the tube, vorticity



The cylinder in the tube, velocity



The sphere in the tube, vorticity



The sphere in the tube, velocity



As expected after initial flow development I've got quite nice vortices going along the tube.

Probes

Since the movies are not quite enough I've added probes along the axis of the tube to record values of velocity, pressure, and turbulent viscosity during the simulations (this is a fragment of controlDict):

functions
{
    probes
    {
        type            probes;
        functionObjectLibs ("libsampling.so");
        outputControl   timeStep;
        outputInterval  1;
        probeLocations
        (
            (2.47 0 0)
            (2.90 0 0)
            (3.5 0 0)
            (5 0 0)
            (6 0 0)
        );
        fields
        (
            U
            p
            nut
        );
    }
}

Below is a graphs of the pressure evolution at the location of the first probe (click on the figure for a larger version):

For internal cylinder


For internal sphere


Though evolution has periodic structure but the frequency differs from estimations.

How to run the cases on your own

I've put all necessary files into repository on BitBucket. So you can clone the repo or just download archive of the default branch. There is Cases folder where case files reside. Every case directory has Allprepare script which will do all the necessary utility invocations:

#!/bin/sh
cd ${0%/*} || exit 1

. $WM_PROJECT_DIR/bin/tools/RunFunctions

gmsh - tube.geo
gmsh - cylinder.geo
sed -i -e "s/Created by Gmsh/cylinder/" cylinder.stl
mv cylinder.stl constant/triSurface
gmshToFoam tube.msh
cp -r 0.org 0
changeDictionary
surfaceFeatureExtract
snappyHexMesh -overwrite

The script assumes that you have gmsh in your $PATH and first it creates mesh for tube and cylinder, then corrects STL file created by Gmsh, then places STL files to the directory where sHM will look for it, converts tube mesh to OpenFOAM format, corrects boundary type for walls, and finally runs sHM for final mesh creation. sHM procedure is rather lengthy and needs RAM as the limit on the number of cells is 1500000.