2.1 Lid-driven cavity flow

From en.wiki.laduga.ru
Jump to: navigation, search


This tutorial will describe how to pre-process, run and post-process a case involving isothermal, incompressible flow in a two-dimensional square domain. The geometry is shown in Figure 2.1 in which all the boundaries of the square are walls. The top wall moves in the x-direction at a speed of 1 m/s while the other 3 are stationary. Initially, the flow will be assumed laminar and will be solved on a uniform mesh using the icoFoam solver for laminar, isothermal, incompressible flow. During the course of the tutorial, the effect of increased mesh resolution and mesh grading towards the walls will be investigated. Finally, the flow Reynolds number will be increased and the pisoFoam solver will be used for turbulent, isothermal, incompressible flow.


Fig2.1.png

Figure 2.1: Geometry of the lid driven cavity.


2.1.1 Pre-processing

Cases are setup in OpenFOAM by editing case files. Users should select an xeditor of choice with which to do this, such as emacs, vi, gedit, kate, nedit, etc. Editing files is possible in OpenFOAM because the I/O uses a dictionary format with keywords that convey sufficient meaning to be understood by even the least experienced users.

A case being simulated involves data for mesh, fields, properties, control parameters, etc. As described in section 4.1, in OpenFOAM this data is stored in a set of files within a case directory rather than in a single case file, as in many other CFD packages. The case directory is given a suitably descriptive name, e.g. the first example case for this tutorial is simply named cavity. In preparation of editing case files and running the first cavity case, the user should change to the case directory


   cd $FOAM_RUN/tutorials/incompressible/icoFoam/cavity

2.1.1.1 Mesh generation

OpenFOAM always operates in a 3 dimensional Cartesian coordinate system and all geometries are generated in 3 dimensions. OpenFOAM solves the case in 3 dimensions by default but can be instructed to solve in 2 dimensions by specifying a ‘special’ empty boundary condition on boundaries normal to the (3rd) dimension for which no solution is required.

The cavity domain consists of a square of side length d = 0.1 m in the x-y plane. A uniform mesh of 20 by 20 cells will be used initially. The block structure is shown in Figure 2.2.


Fig2.2.png

Figure 2.2: Block structure of the mesh for the cavity.


The mesh generator supplied with OpenFOAM, blockMesh, generates meshes from a description specified in an input dictionary, blockMeshDict located in the constant/polyMesh directory for a given case. The blockMeshDict entries for this case are as follows:

   1  /*--------------------------------*- C++ -*----------------------------------*\
   2  | =========                 |                                                 |
   3  | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
   4  |  \\    /   O peration     | Version:  2.3.0                                 |
   5  |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
   6  |    \\/     M anipulation  |                                                 |
   7  \*---------------------------------------------------------------------------*/
   8  FoamFile
   9  {
   10      version     2.0;
   11      format      ascii;
   12      class       dictionary;
   13      object      blockMeshDict;
   14  }
   15  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
   16  
   17  convertToMeters 0.1;
   18  
   19  vertices
   20  (
   21      (0 0 0)
   22      (1 0 0)
   23      (1 1 0)
   24      (0 1 0)
   25      (0 0 0.1)
   26      (1 0 0.1)
   27      (1 1 0.1)
   28      (0 1 0.1)
   29  );
   30  
   31  blocks
   32  (
   33      hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
   34  );
   35  
   36  edges
   37  (
   38  );
   39  
   40  boundary
   41  (
   42      movingWall
   43      {
   44          type wall;
   45          faces
   46          (
   47              (3 7 6 2)
   48          );
   49      }
   50      fixedWalls
   51      {
   52          type wall;
   53          faces
   54          (
   55              (0 4 7 3)
   56              (2 6 5 1)
   57              (1 5 4 0)
   58          );
   59      }
   60      frontAndBack
   61      {
   62          type empty;
   63          faces
   64          (
   65              (0 3 2 1)
   66              (4 5 6 7)
   67          );
   68      }
   69  );
   70  
   71  mergePatchPairs
   72  (
   73  );
   74  
   75  // ************************************************************************* //

The file first contains header information in the form of a banner (lines 1-7), then file information contained in a FoamFile sub-dictionary, delimited by curly braces ({…}).

For the remainder of the manual: For the sake of clarity and to save space, file headers, including the banner and FoamFile sub-dictionary, will be removed from verbatim quoting of case files

The file first specifies coordinates of the block vertices; it then defines the blocks (here, only 1) from the vertex labels and the number of cells within it; and finally, it defines the boundary patches. The user is encouraged to consult section 5.3 to understand the meaning of the entries in the blockMeshDict file.

The mesh is generated by running blockMesh on this blockMeshDict file. From within the case directory, this is done, simply by typing in the terminal:

   blockMesh

The running status of blockMesh is reported in the terminal window. Any mistakes in the blockMeshDict file are picked up by blockMesh and the resulting error message directs the user to the line in the file where the problem occurred. There should be no error messages at this stage.

2.1.1.2 Boundary and initial conditions

Once the mesh generation is complete, the user can look at this initial fields set up for this case. The case is set up to start at time t = 0 s, so the initial field data is stored in a 0 sub-directory of the cavity directory. The 0 sub-directory contains 2 files, p and U, one for each of the pressure (p) and velocity (U) fields whose initial values and boundary conditions must be set. Let us examine file p:

   17  dimensions      [0 2 -2 0 0 0 0];
   18  
   19  internalField   uniform 0;
   20  
   21  boundaryField
   22  {
   23      movingWall
   24      {
   25          type            zeroGradient;
   26      }
   27  
   28      fixedWalls
   29      {
   30          type            zeroGradient;
   31      }
   32  
   33      frontAndBack
   34      {
   35          type            empty;
   36      }
   37  }
   38  
   39  // ************************************************************************* //

There are 3 principal entries in field data files:

dimensions
specifies the dimensions of the field, here kinematic pressure, i.e. $ m^2 s^{-2} $ (see section 4.2.6 for more information);
internalField
the internal field data which can be uniform, described by a single value; or nonuniform, where all the values of the field must be specified (see section 4.2.8 for more information);
boundaryField
the boundary field data that includes boundary conditions and data for all the boundary patches (see section 4.2.8 for more information).

For this case cavity, the boundary consists of walls only, split into 2 patches named: (1) fixedWalls for the fixed sides and base of the cavity; (2) movingWall for the moving top of the cavity. As walls, both are given a zeroGradient boundary condition for p, meaning "the normal gradient of pressure is zero". The frontAndBack patch represents the front and back planes of the 2D case and therefore must be set as empty.

In this case, as in most we encounter, the initial fields are set to be uniform. Here the pressure is kinematic, and as an incompressible case, its absolute value is not relevant, so is set to uniform 0 for convenience.

The user can similarly examine the velocity field in the 0/U file. The dimensions are those expected for velocity, the internal field is initialised as uniform zero, which in the case of velocity must be expressed by 3 vector components, i.e. uniform (0 0 0) (see section 4.2.5 for more information).

The boundary field for velocity requires the same boundary condition for the frontAndBack patch. The other patches are walls: a no-slip condition is assumed on the fixedWalls, hence a fixedValue condition with a value of uniform (0 0 0). The top surface moves at a speed of 1 m/s in the $ x $-direction so requires a fixedValue condition also but with uniform (1 0 0).

2.1.1.3 Physical properties

The physical properties for the case are stored in dictionaries whose names are given the suffix …Properties, located in the Dictionaries directory tree. For an icoFoam case, the only property that must be specified is the kinematic viscosity which is stored from the transportProperties dictionary. The user can check that the kinematic viscosity is set correctly by opening the transportProperties dictionary to view/edit its entries. The keyword for kinematic viscosity is nu, the phonetic label for the Greek symbol $ \nu $ by which it is represented in equations. Initially this case will be run with a Reynolds number of 10, where the Reynolds number is defined as:

$ Re= \frac{d |U|}{\nu} $ (2.1)

where $ d $ and $ |U | $ are the characteristic length and velocity respectively and $ \nu $ is the kinematic viscosity. Here $ d = 0.1 m , |U | = 1 ms^{ -1} $, so that for $ Re = 10, \nu = 0.01 m^2 s^{ -1} $. The correct file entry for kinematic viscosity is thus specified below:

17   
18  nu              nu [ 0 2 -1 0 0 0 0 ] 0.01; 
19   
20   
21  // ************************************************************************* //

2.1.1.4 Control

Input data relating to the control of time and reading and writing of the solution data are read in from the controlDict dictionary. The user should view this file; as a case control file, it is located in the system directory.

The start/stop times and the time step for the run must be set. OpenFOAM offers great flexibility with time control which is described in full in section 4.3. In this tutorial we wish to start the run at time $ t = 0 $ which means that OpenFOAM needs to read field data from a directory named 0 — see section 4.1 for more information of the case file structure. Therefore we set the startFrom keyword to startTime and then specify the startTime keyword to be 0.

For the end time, we wish to reach the steady state solution where the flow is circulating around the cavity. As a general rule, the fluid should pass through the domain 10 times to reach steady state in laminar flow. In this case the flow does not pass through this domain as there is no inlet or outlet, so instead the end time can be set to the time taken for the lid to travel ten times across the cavity, i.e. 1 s; in fact, with hindsight, we discover that 0.5 s is sufficient so we shall adopt this value. To specify this end time, we must specify the stopAt keyword as endTime and then set the endTime keyword to 0.5.

Now we need to set the time step, represented by the keyword deltaT. To achieve temporal accuracy and numerical stability when running icoFoam, a Courant number of less than 1 is required. The Courant number is defined for one cell as:

$ Co = \frac {\delta t |U|} {\delta x} $ (2.2)

where $ \delta t $ is the time step, $ |U | $ is the magnitude of the velocity through that cell and $ \delta x $ is the cell size in the direction of the velocity. The flow velocity varies across the domain and we must ensure $ Co < 1 $ everywhere. We therefore choose $ \delta t $ based on the worst case: the maximum $ Co $ corresponding to the combined effect of a large flow velocity and small cell size. Here, the cell size is fixed across the domain so the maximum $ Co $ occur next to the lid where the velocity approaches $ 1 m s^{-1} $.

The cell size is:

$ \delta x = \frac {d}{n} = \frac {0.1}{20} = 0.005 m $ (2.3)

Therefore to achieve a Courant number less than or equal to 1 throughout the domain the time step deltaT must be set to less than or equal to:

$ \delta t = \frac {Co \delta x}{|U|} = \frac {1 \times 0.0005}{1} = 0.005 s $ (2.4)

As the simulation progresses we wish to write results at certain intervals of time that we can later view with a post-processing package. The writeControl keyword presents several options for setting the time at which the results are written; here we select the timeStep option which specifies that results are written every $ n $th time step where the value $ n $ is specified under the writeInterval keyword. Let us decide that we wish to write our results at times 0.1, 0.2,…, 0.5 s. With a time step of 0.005 s, we therefore need to output results at every 20th time time step and so we set writeInterval to 20.

OpenFOAM creates a new directory named after the current time, e.g. 0.1 s, on each occasion that it writes a set of data, as discussed in full in section 4.1. In the icoFoam solver, it writes out the results for each field, U and p, into the time directories. For this case, the entries in the controlDict are shown below:

17   
18  application     icoFoam; 
19   
20  startFrom       startTime; 
21   
22  startTime       0; 
23   
24  stopAt          endTime; 
25   
26  endTime         0.5; 
27   
28  deltaT          0.005; 
29   
30  writeControl    timeStep; 
31   
32  writeInterval   20; 
33   
34  purgeWrite      0; 
35   
36  writeFormat     ascii; 
37   
38  writePrecision  6; 
39   
40  writeCompression off; 
41   
42  timeFormat      general; 
43   
44  timePrecision   6; 
45   
46  runTimeModifiable true; 
47   
48   
49  // ************************************************************************* //

2.1.1.5 Discretisation and linear-solver settings

The user specifies the choice of finite volume discretisation schemes in the fvSchemes dictionary in the system directory. The specification of the linear equation solvers and tolerances and other algorithm controls is made in the fvSolution dictionary, similarly in the system directory. The user is free to view these dictionaries but we do not need to discuss all their entries at this stage except for pRefCell and pRefValue in the PISO sub-dictionary of the fvSolution dictionary. In a closed incompressible system such as the cavity, pressure is relative: it is the pressure range that matters not the absolute values. In cases such as this, the solver sets a reference level by pRefValue in cell pRefCell. In this example both are set to 0. Changing either of these values will change the absolute pressure field, but not, of course, the relative pressures or velocity field.

2.1.2 Viewing the mesh

Before the case is run it is a good idea to view the mesh to check for any errors. The mesh is viewed in paraFoam, the post-processing tool supplied with OpenFOAM. The paraFoam post-processing is started by typing in the terminal from within the case directory

   paraFoam

Alternatively, it can be launched from another directory location with an optional -case argument giving the case directory, e.g.

   paraFoam -case $FOAM_RUN/tutorials/incompressible/icoFoam/cavity

This launches the ParaView window as shown in Figure 6.1. In the Pipeline Browser, the user can see that ParaView has opened cavity.OpenFOAM, the module for the cavity case. Before clicking the Apply button, the user needs to select some geometry from the Mesh Parts panel. Because the case is small, it is easiest to select all the data by checking the box adjacent to the Mesh Parts panel title, which automatically checks all individual components within the respective panel. The user should then click the Apply button to load the geometry into ParaView.

The user should then open the Display panel that controls the visual representation of the selected module. Within the Display panel the user should do the following as shown in Figure 2.3: (1) set Color By Solid Color; (2) click Set Ambient Color and select an appropriate colour e.g. black (for a white background); (3) in the Style panel, select Wireframe from the Representation menu. The background colour can be set by selecting View Settings… from Edit in the top menu panel.


Fig2.3.png

Figure 2.3: Viewing the mesh in paraFoam.


Especially the first time the user starts ParaView, it is recommended that they manipulate the view as described in section 6.1.5. In particular, since this is a 2D case, it is recommended that Use Parallel Projection is selected in the General panel of View Settings window selected from the Edit menu. The Orientation Axes can be toggled on and off in the Annotation window or moved by drag and drop with the mouse.

2.1.3 Running an application

Like any UNIX/Linux executable, OpenFOAM applications can be run in two ways: as a foreground process, i.e. one in which the shell waits until the command has finished before giving a command prompt; as a background process, one which does not have to be completed before the shell accepts additional commands.

On this occasion, we will run icoFoam in the foreground. The icoFoam solver is executed either by entering the case directory and typing

   icoFoam

at the command prompt, or with the optional -case argument giving the case directory, e.g.

   icoFoam -case $FOAM_RUN/tutorials/incompressible/icoFoam/cavity

The progress of the job is written to the terminal window. It tells the user the current time, maximum Courant number, initial and final residuals for all fields.


Fig2.4.png

Figure 2.4: Displaying pressure contours for the cavity case.



Fig2.5.png

Figure 2.5: Pressures in the cavity case.


2.1.4 Post-processing

As soon as results are written to time directories, they can be viewed using paraFoam. Return to the paraFoam window and select the Properties panel for the cavity.OpenFOAM case module. If the correct window panels for the case module do not seem to be present at any time, please ensure that: cavity.OpenFOAM is highlighted in blue; eye button alongside it is switched on to show the graphics are enabled;

To prepare paraFoam to display the data of interest, we must first load the data at the required run time of 0.5 s. If the case was run while ParaView was open, the output data in time directories will not be automatically loaded within ParaView. To load the data the user should click Refresh Times in the Properties window. The time data will be loaded into ParaView.

2.1.4.1 Isosurface and contour plots

To view pressure, the user should open the Display panel since it controls the visual representation of the selected module. To make a simple plot of pressure, the user should select the following, as described in detail in Figure 2.4: in the Style panel, select Surface from the Representation menu; in the Color panel, select Color by $ \cdot p $ and Rescale to Data Range. Now in order to view the solution at $ t = 0.5 s $, the user can use the VCR Controls or Current Time Controls to change the current time to 0.5. These are located in the toolbars below the menus at the top of the ParaView window, as shown in Figure 6.4. The pressure field solution has, as expected, a region of low pressure at the top left of the cavity and one of high pressure at the top right of the cavity as shown in Figure 2.5.

With the point icon ($ \cdot p $) the pressure field is interpolated across each cell to give a continuous appearance. Instead if the user selects the cell icon, $ \Box p $, from the Color by menu, a single value for pressure will be attributed to each cell so that each cell will be denoted by a single colour with no grading.

A colour bar can be included by either by clicking the Toggle Color Legend Visibility button in the Active Variable Controls toolbar, or by selecting Show Color Legend from the View menu. Clicking the Edit Color Map button, either in the Active Variable Controls toolbar or in the Color panel of the Display window, the user can set a range of attributes of the colour bar, such as text size, font selection and numbering format for the scale. The colour bar can be located in the image window by drag and drop with the mouse.

New versions of ParaView default to using a colour scale of blue to white to red rather than the more common blue to green to red (rainbow). Therefore the first time that the user executes ParaView, they may wish to change the colour scale. This can be done by selecting Choose Preset in the Color Scale Editor and selecting Blue to Red Rainbow. After clicking the OK confirmation button, the user can click the Make Default button so that ParaView will always adopt this type of colour bar.

If the user rotates the image, they can see that they have now coloured the complete geometry surface by the pressure. In order to produce a genuine contour plot the user should first create a cutting plane, or ‘slice’, through the geometry using the Slice filter as described in section 6.1.6.1. The cutting plane should be centred at $ (0.05,0.05,0.005) $ and its normal should be set to $ (0,0,1) $ (click the Z Normal button). Having generated the cutting plane, the contours can be created using by the Contour filter described in section 6.1.6.

2.1.4.2 Vector plots

Before we start to plot the vectors of the flow velocity, it may be useful to remove other modules that have been created, e.g. using the Slice and Contour filters described above. These can: either be deleted entirely, by highlighting the relevant module in the Pipeline Browser and clicking Delete in their respective Properties panel; or, be disabled by toggling the eye button for the relevant module in the Pipeline Browser.

We now wish to generate a vector glyph for velocity at the centre of each cell. We first need to filter the data to cell centres as described in section 6.1.7.1. With the cavity.OpenFOAM module highlighted in the Pipeline Browser, the user should select Cell Centers from the Filter->Alphabetical menu and then click Apply.

With these Centers highlighted in the Pipeline Browser, the user should then select Glyph from the Filter->Alphabetical menu. The Properties window panel should appear as shown in Figure 2.6. In the resulting Properties panel, the velocity field, U, is automatically selected in the vectors menu, since it is the only vector field present. By default the Scale Mode for the glyphs will be Vector Magnitude of velocity but, since the we may wish to view the velocities throughout the domain, the user should instead select off and Set Scale Factor to 0.005. On clicking Apply, the glyphs appear but, probably as a single colour, e.g. white. The user should colour the glyphs by velocity magnitude which, as usual, is controlled by setting Color by U in the Display panel. The user should also select Show Color Legend in Edit Color Map. The output is shown in Figure 2.7, in which uppercase Times Roman fonts are selected for the Color Legend headings and the labels are specified to 2 fixed significant figures by deselecting Automatic Label Format and entering %-#6.2f in the Label Format text box. The background colour is set to white in the General panel of View Settings as described in section 6.1.5.1.

Note that at the left and right walls, glyphs appear to indicate flow through the walls. On closer examination, however, the user can see that while the flow direction is normal to the wall, its magnitude is 0. This slightly confusing situation is caused by ParaView choosing to orientate the glyphs in the x-direction when the glyph scaling off and the velocity magnitude is 0.

2.1.4.3 Streamline plots

Again, before the user continues to post-process in ParaView, they should disable modules such as those for the vector plot described above. We now wish to plot streamlines of velocity as described in section 6.1.8.

With the cavity.OpenFOAM module highlighted in the Pipeline Browser, the user should then select Stream Tracer from the Filter menu and then click Apply. The Properties window panel should appear as shown in Figure 2.8. The Seed points should be specified along a Line Source running vertically through the centre of the geometry, i.e. from (0.05, 0, 0.005) to (0.05, 0.1, 0.005). For the image in this guide we used: a point Resolution of 21; Max Propagation by Length 0.5; Initial Step Length by Cell Length 0.01; and, Integration Direction BOTH. The Runge-Kutta 2 IntegratorType was used with default parameters.

On clicking Apply the tracer is generated. The user should then select Tube from the Filter menu to produce high quality streamline images. For the image in this report, we used: Num. sides 6; Radius 0.0003; and, Radius factor 10. The streamtubes are coloured by velocity magnitude. On clicking Apply the image in Figure 2.9 should be produced.

2.1.5 Increasing the mesh resolution

The mesh resolution will now be increased by a factor of two in each direction. The results from the coarser mesh will be mapped onto the finer mesh to use as initial conditions for the problem. The solution from the finer mesh will then be compared with those from the coarser mesh.

2.1.5.1 Creating a new case using an existing case

We now wish to create a new case named cavityFine that is created from cavity. The user should therefore clone the cavity case and edit the necessary files. First the user should create a new case directory at the same directory level as the cavity case, e.g.

cd $FOAM RUN/tutorials/incompressible/icoFoam
mkdir cavityFine

The user should then copy the base directories from the cavity case into cavityFine, and then enter the cavityFine case.

cp -r cavity/constant cavityFine
cp -r cavity/system cavityFine
cd cavityFine

2.1.5.2 Creating the finer mesh

We now wish to increase the number of cells in the mesh by using blockMesh. The user should open the blockMeshDict file in an editor and edit the block specification. The blocks are specified in a list under the blocks keyword. The syntax of the block definitions is described fully in section 5.3.1.3; at this stage it is sufficient to know that following hex is first the list of vertices in the block, then a list (or vector) of numbers of cells in each direction. This was originally set to (20 20 1) for the cavity case. The user should now change this to (40 40 1) and save the file. The new refined mesh should then be created by running blockMesh as before.

2.1.5.3 Mapping the coarse mesh results onto the fine mesh

The mapFields utility maps one or more fields relating to a given geometry onto the corresponding fields for another geometry. In our example, the fields are deemed ‘consistent’ because the geometry and the boundary types, or conditions, of both source and target fields are identical. We use the -consistent command line option when executing mapFields in this example.

The field data that mapFields maps is read from the time directory specified by startFrom/startTime in the controlDict of the target case, i.e. those into which the results are being mapped. In this example, we wish to map the final results of the coarser mesh from case cavity onto the finer mesh of case cavityFine. Therefore, since these results are stored in the 0.5 directory of cavity, the startTime should be set to 0.5 s in the controlDict dictionary and startFrom should be set to startTime.

The case is ready to run mapFields. Typing mapFields -help quickly shows that mapFields requires the source case directory as an argument. We are using the -consistent option, so the utility is executed from withing the cavityFine directory by

mapFields ../cavity -consistent

The utility should run with output to the terminal including:

Source: ".." "cavity"
Target: "." "cavityFine"
Create databases as time
Case
: ../cavity
nProcs : 1
Source time: 0.5
Target time: 0.5
Create meshes
Source mesh size: 400 Target mesh size: 1600
Consistently creating and mapping fields for time 0.5
Creating mesh-to-mesh addressing ...
Overlap volume: 0.0001
Creating AMI between source patch movingWall and target patch movingWall ...
interpolating p
interpolating U
End

2.1.5.4 Control adjustments