OpenFOAM v9 User Guide - 2.1 Lid-driven cavity flow

2.1 Lid-driven cavity flow

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

---

---

2.1.1 Pre-processing

Cases are setup in OpenFOAM by editing case files. Users should select an editor of choice with which to do this, such as emacs, vi, gedit, 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 the 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. This tutorial consists of a set of cases located in $FOAM_TUTORIALS/incompressible/icoFoam/cavity, the first of which is simply named cavity. As a first step, the user should copy the cavity case directory to their run directory.

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 in the - plane. A uniform mesh of 20 by 20 cells will be used initially. The block structure is shown in Figure 2.2 .

---

---

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

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:
 

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:

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  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 () and velocity () fields whose initial values and boundary conditions must be set. Let us examine file p:

There are 3 principal entries in field data files:

dimensions

specifies the dimensions of the field, here kinematic pressure, i.e.  (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 noSlip condition. The top surface moves at a speed of 1 m/s in the -direction so requires a fixedValue condition with value of uniform (1 0 0).

2.1.1.3 Physical properties

The physical properties for the case are stored in dictionary files in the constant directory. 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 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:

(2.1)

where and are the characteristic length and velocity respectively and is the kinematic viscosity. Here  0.1 m,  1 m/s, so that for  10,  0.01 . The correct file entry for kinematic viscosity is thus specified below:

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.4 . In this tutorial we wish to start the run at time 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:

(2.2)

where is the time step, is the magnitude of the velocity through that cell and is the cell size in the direction of the velocity. The flow velocity varies across the domain and we must ensure everywhere. We therefore choose based on the worst case: the maximum 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 will occur next to the lid where the velocity approaches 1 . The cell size is:

(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:

(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 th time step where the value 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:

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 ParaView, the post-processing tool supplied with OpenFOAM. The ParaView post-processing is conveniently launched on OpenFOAM case data by executing the paraFoam script from within the case directory.

Any UNIX/Linux executable 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, which allows the shell to accept additional commands while it is still running. Since it is convenient to keep ParaView open while running other commands from the terminal, we will launch it in the background using the & operator by typing

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

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 scroll down to 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. in the Coloring section, select Solid Color;
2. click Edit (in Coloring) and select an appropriate colour e.g. black (for a white background);
3. select Wireframe from the Representation menu. The background colour can be set in the View Render panel below the Display panel in the Properties window.

---

---

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 Camera Parallel Projection is selected near the bottom of the View Render panel.

2.1.3 Running an application

Like any UNIX/Linux executable, OpenFOAM applications can be run either in the foreground or background. On this occasion, we will run icoFoam in the foreground. The icoFoam solver is executed either by entering the case directory and typing

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

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.

---

---

---

---

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 uncheck Cache Mesh and click Refresh Times at the top Properties window (scroll up the panel if necessary). When The time data will be loaded into ParaView.

In order to view the solution at  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 at the top of the ParaView window, as shown in Figure 6.4 .

2.1.4.1 Colouring surfaces

To view pressure, the user should go to 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 :

1. select Surface from the Representation menu;
2. select in Coloring
3. click the Rescale button to set the colour scale to the data range, if necessary.

The pressure field should appear as shown in Figure 2.5 , with a region of low pressure at the top left of the cavity and one of high pressure at the top right of the cavity.

With the point icon () the pressure field is interpolated across each cell to give a continuous appearance. Instead if the user selects the cell icon, , from the Coloring 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 legend can be added by either by clicking the Toggle Color Legend Visibility button in the Active Variable Controls toolbar or the Show button in the Coloring section of the Display panel. The legend can be located in the image window by drag and drop with the mouse. The Edit button, either in the Active Variable Controls toolbar or in the Coloring panel of the Display panel, opens the Color Map Editor window, as shown in Figure 2.6 , where the user can set a range of attributes of the colour scale and the color bar.

---

---

In particular, ParaView defaults 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 the Choose Preset button (with the heart icon) in the Color Scale Editor. The conventional color scale for CFD is Blue to Red Rainbow which is only listed if the user types the name in the Search bar or checks the gearwheel to the right of that bar.

After selecting Blue to Red Rainbow and clicking Apply and Close, the user can click the Save as Default button at the absolute bottom of the panel (file save symbol) so that ParaView will always adopt this type of colour bar.

The user can also edit the color legend properties, such as text size, font selection and numbering format for the scale, by clicking the Edit Color Legend Properties to the far right of the search bar, as shown in Figure 2.6 .

2.1.4.2 Cutting plane (slice)

If the user rotates the image, by holding down the left mouse button in the image window and moving the cursor, they can see that they have now coloured the complete geometry surface by the pressure. In order to produce a genuine 2-dimensional contour plot the user should first create a cutting plane, or ‘slice’. With the cavity.OpenFOAM module highlighted in the Pipeline Browser, the user should select the Slice filter from the Filters menu in the top menu of ParaView (accessible at the top of the screen on some systems). The Slice filter can be initially found in the Common sub-menu, but once selected, it moves to the Recent sub-menu, disappearing from the Common sub-menu. The cutting plane should be centred at and its normal should be set to (click the Z Normal button).

2.1.4.3 Contours

Having generated the cutting plane, contours can be created using by applying the Contour filter. With the Slice module highlighted in the Pipeline Browser, the user should select the Contour filter. In the Properties panel, the user should select pressure from the Contour By menu. Under Isosurfaces, the user could delete the default value with the minus button, then add a range of 10 values. The contours can be displayed with a Wireframe representation if the Coloring is solid or by a field, e.g. pressure.

---

---

---

---

2.1.4.4 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->Common menu. The Properties window panel should appear as shown in Figure 2.7 . This outdated image omits Orientation Array which must be set to U.

We suggest popular filters are added to the Filters->Favourites menu. In the resulting Properties panel, the velocity field, U, must be selected from the vectors menu. The user should set the Scale Array for the glyphs to be No scale array, with Set Scale Factor set to 0.005. Under Glyph Mode, the user should select All Points.

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 can also select Show Color Legend in Edit Color Map. The output is shown in Figure 2.8, 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. However, it is clear 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 -direction when the glyph scaling off and the velocity magnitude is 0.

---

---

---

---

2.1.4.5 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.9 . The Seed points should be specified along a High Resolution Line Source running vertically through the centre of the geometry, i.e. from to . For the image in this guide we used: a point Resolution of 21; Maximum Step Length of 0.5; Initial Step Length of 0.2; and, Integration Direction BOTH. The Runge-Kutta 4/5 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.10 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 go to the run directory, by typing

Note that there is also a convenient alias, named run, that reproduces the command above to change directory to $FOAM_RUN, simply by typing run.

The cavityFine case can be created by making a new case directory and copying the relevant directories from the cavity case.

The user can then prepare to run the new case by changing into the case directory.

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 the system directory 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 and 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 within the cavityFine directory by

The utility should run with output to the terminal including:

2.1.5.4 Control adjustments

To maintain a Courant number of less that 1, as discussed in section 2.1.1.4 , the time step must now be halved since the size of all cells has halved. Therefore deltaT should be set to 0.0025 s in the controlDict dictionary. Field data is currently written out at an interval of a fixed number of time steps. Here we demonstrate how to specify data output at fixed intervals of time. Under the writeControl keyword in controlDict, instead of requesting output by a fixed number of time steps with the timeStep entry, a fixed amount of run time can be specified between the writing of results using the runTime entry. In this case the user should specify output every 0.1  and therefore should set writeInterval to 0.1 and writeControl to runTime. Finally, since the case is starting with a solution obtained on the coarse mesh we only need to run it for a short period to achieve reasonable convergence to steady-state. Therefore the endTime should be set to 0.7 s. Make sure these settings are correct and then save the file.

2.1.5.5 Running the code as a background process

The user should experience running icoFoam as a background process, redirecting the terminal output to a log file that can be viewed later. From the cavityFine directory, the user should execute:

2.1.5.6 Vector plot with the refined mesh

The user can open multiple cases simultaneously in ParaView; essentially because each new case is simply another module that appears in the Pipeline Browser. There is an inconvenience when opening a new OpenFOAM case in ParaView because it expects that case data is stored in a single file which has a file extension that enables it to establish the format. However, OpenFOAM stores case data in multiple files without an extension in the name, within a specific directory structure. The ParaView reader module works on the basis that, when opening case data in OpenFOAM format, it is passed a dummy (empty) file with the .OpenFOAM extension that resides in the case directory. The paraFoam script automatically creates this file — hence, the cavity case module is called cavity.OpenFOAM.

If the user wishes to open a second case directly from within ParaView, they need to create such a dummy file. They can do this ‘by hand’ or, more simply, use the paraFoam script with the option -touch. For the cavityFine example, that involves executing from the case directory:

Now the cavityFine case can be loaded into ParaView by selecting Open from the File menu, and having navigated to the cavityFine directory, opening cavityFine.OpenFOAM. The user can now make a vector plot of the results from the refined mesh in ParaView. The plot can be compared with the cavity case by enabling glyph images for both case simultaneously.

---

---

2.1.5.7 Plotting graphs

The user may wish to visualise the results by extracting some scalar measure of velocity and plotting 2-dimensional graphs along lines through the domain. OpenFOAM is well equipped for this kind of data manipulation. There are numerous utilities that do specialised data manipulations, and the postProcess utility that includes a broad range of generic post-processing functionality. The functions within postProcess can be listed by typing:

The components and mag functions provide useful scalar measures of velocity. When the components function is executed on a case, e.g. cavity, it reads in the velocity vector field from each time directory and, in the corresponding time directories, writes scalar fields Ux, Uy and Uz representing the , and components of velocity.

The user can run postProcess with the components function on both cavity and cavityFine cases. For example, for the cavity case the user should go into the cavity directory and execute postProcess as follows:

The individual components can be plotted as a graph in ParaView. It is quick, convenient and has reasonably good control over labelling and formatting, so the printed output is a fairly good standard. However, to produce graphs for publication, users may prefer to write raw data and plot it with a dedicated graphing tool, such as gnuplot or Grace/xmgr. To do this, we recommend using the sampling tools, described in section 6.3.2 and section 2.2.3 .

Before commencing plotting, the user needs to load the newly generated Ux, Uy and Uz fields into ParaView. To do this, the user should click the Refresh Times at the top of the Properties panel for the cavity.OpenFOAM module which will cause the new fields to be loaded into ParaView and appear in the Volume Fields window. Ensure the new fields are selected and the changes are applied, i.e. click Apply again if necessary. Also, data is interpolated incorrectly at boundaries if the boundary regions are selected in the Mesh Parts panel. Therefore the user should deselect the patches in the Mesh Parts panel, i.e. movingWall, fixedWall and frontAndBack, and apply the changes.

Now, in order to display a graph in ParaView the user should select the module of interest, e.g.cavity.OpenFOAM and apply the Plot Over Line filter from the Filter->Data Analysis menu. This opens up a new XY Plot window below or beside the existing 3D View window. A PlotOverLine module is created in which the user can specify the end points of the line in the Properties panel. In this example, the user should position the line vertically up the centre of the domain, i.e. from to , in the Point1 and Point2 text boxes. The Resolution can be set to 100.

On clicking Apply, a graph is generated in the XY Plot window. In the Display panel, the user should set Attribute Mode to Point Data. The Use Data Array option can be selected for the X Axis Data, taking the arc_length option so that the x-axis of the graph represents distance from the base of the cavity.

The user can choose the fields to be displayed in the Line Series panel of the Display window. From the list of scalar fields to be displayed, it can be seen that the magnitude and components of vector fields are available by default, e.g. displayed as U_X, so that it was not necessary to create Ux using the components function. Nevertheless, the user should deselect all series except Ux (or U_x). A square colour box in the adjacent column to the selected series indicates the line colour. The user can edit this most easily by a double click of the mouse over that selection.

In order to format the graph, the user should modify the settings below the Line Series panel, namely Line Color, Line Thickness, Line Style, Marker Style and Chart Axes.

Also the user can click one of the buttons above the top left corner of the XY Plot. The third button, for example, allows the user to control View Settings in which the user can set title and legend for each axis, for example. Also, the user can set font, colour and alignment of the axes titles, and has several options for axis range and labels in linear or logarithmic scales.

Figure 2.12 is a graph produced using ParaView. The user can produce a graph however he/she wishes. For information, the graph in Figure 2.12 was produced with the options for axes of: Standard type of Notation; Specify Axis Range selected; titles in Sans Serif 12 font. The graph is displayed as a set of points rather than a line by activating the Enable Line Series button in the Display window. Note: if this button appears to be inactive by being “greyed out”, it can be made active by selecting and deselecting the sets of variables in the Line Series panel. Once the Enable Line Series button is selected, the Line Style and Marker Style can be adjusted to the user’s preference.

---

---

2.1.6 Introducing mesh grading

The error in any solution will be more pronounced in regions where the form of the true solution differ widely from the form assumed in the chosen numerical schemes. For example a numerical scheme based on linear variations of variables over cells can only generate an exact solution if the true solution is itself linear in form. The error is largest in regions where the true solution deviates greatest from linear form, i.e. where the change in gradient is largest. Error decreases with cell size.

It is useful to have an intuitive appreciation of the form of the solution before setting up any problem. It is then possible to anticipate where the errors will be largest and to grade the mesh so that the smallest cells are in these regions. In the cavity case the large variations in velocity can be expected near a wall and so in this part of the tutorial the mesh will be graded to be smaller in this region. By using the same number of cells, greater accuracy can be achieved without a significant increase in computational cost.

A mesh of cells with grading towards the walls will be created for the lid-driven cavity problem and the results from the finer mesh of section 2.1.5.2 will then be mapped onto the graded mesh to use as an initial condition. The results from the graded mesh will be compared with those from the previous meshes. Since the changes to the blockMeshDict dictionary are fairly substantial, the case used for this part of the tutorial, cavityGrade, is supplied in the $FOAM_TUTORIALS/incompressible/icoFoam/cavity directory. The user should copy the cavityGrade case into the run directory, then follow the steps below.

2.1.6.1 Creating the graded mesh

The mesh now needs 4 blocks as different mesh grading is needed on the left and right and top and bottom of the domain. The block structure for this mesh is shown in Figure 2.13 .

---

---

The user can view the blockMeshDict file in the system subdirectory of cavityGrade; for completeness the key elements of the blockMeshDict file are also reproduced below. Each block now has cells in the and directions and the ratio between largest and smallest cells is .

Once familiar with the blockMeshDict file for this case, the user can execute blockMesh from the command line. The graded mesh can be viewed as before using paraFoam as described in section 2.1.2 .

2.1.6.2 Changing time and time step

The highest velocities and smallest cells are next to the lid, therefore the highest Courant number will be generated next to the lid, for reasons given in section 2.1.1.4 . It is therefore useful to estimate the size of the cells next to the lid to calculate an appropriate time step for this case.

When a nonuniform mesh grading is used, blockMesh calculates the cell sizes using a geometric progression. Along a length , if cells are requested with a ratio of between the last and first cells, the size of the smallest cell, , is given by:

(2.5)

where is the ratio between one cell size and the next which is given by:

(2.6)

and

(2.7)

For the cavityGrade case the number of cells in each direction in a block is 10, the ratio between largest and smallest cells is and the block height and width is 0.05 . Therefore the smallest cell length is 3.45 . From Equation 2.2 , the time step should be less than 3.45  to maintain a Courant of less than 1. To ensure that results are written out at convenient time intervals, the time step deltaT should be reduced to 2.5  and the writeInterval set to 40 so that results are written out every 0.1 s. These settings can be viewed in the cavityGrade/system/controlDict file.

The startTime needs to be set to that of the final conditions of the case cavityFine, i.e. 0.7. Since cavity and cavityFine converged well within the prescribed run time, we can set the run time for case cavityGrade to 0.1 s, i.e. the endTime should be 0.8.

2.1.6.3 Mapping fields

First,

As in section 2.1.5.3 , use mapFields to map the final results from case cavityFine onto the mesh for case cavityGrade. From the cavityGrade directory, execute mapFields by:

Now run icoFoam from the case directory and monitor the run time information. View the converged results for this case and compare with other results using post-processing tools described previously in section 2.1.5.6 and section 2.1.5.7 .

2.1.7 Increasing the Reynolds number

The cases solved so far have had a Reynolds number of 10. This is very low and leads to a stable solution quickly with only small secondary vortices at the bottom corners of the cavity. We will now increase the Reynolds number to 100, at which point the solution takes a noticeably longer time to converge. The coarsest mesh in case cavity will be used initially. The user should clone the cavity case and name it cavityHighRe. Rather than copying individual directories (system, constant, etc.) as described previously, the foamCloneCase can be used, which copies the relevant directories in one step. By default the 0 time directory is copied, but here the user can use the -latestTime option to copy the latest time directory, 0.5, which can be used as the initial field data for our simulation. The example also uses the run alias as a quick way to change to the run directory.

2.1.7.1 Pre-processing

Go into the cavityHighRe case and edit the transportProperties dictionary in the constant directory. Since the Reynolds number is required to be increased by a factor of 10, decrease the kinematic viscosity by a factor of 10, i.e. to  . We now run this case by restarting from the solution at the end of the cavity case run. To do this we can use the option of setting the startFrom keyword to latestTime so that icoFoam takes as its initial data the values stored in the directory corresponding to the most recent time, i.e. 0.5. The endTime should be set to 2 s.

2.1.7.2 Running the code

Run icoFoam for this case from the case directory and view the run time information. When running a job in the background, the following UNIX commands can be useful:

nohup

enables a command to keep running after the user who issues the command has logged out;

nice

changes the priority of the job in the kernel’s scheduler; a niceness of -20 is the highest priority and 19 is the lowest priority.

This is useful, for example, if a user wishes to set a case running on a remote machine and does not wish to monitor it heavily, in which case they may wish to give it low priority on the machine. In that case the nohup command allows the user to log out of a remote machine he/she is running on and the job continues running, while nice can set the priority to 19. For our case of interest, we can execute the command in this manner as follows:

In previous runs you may have noticed that icoFoam stops solving for velocity U quite quickly but continues solving for pressure p for a lot longer or until the end of the run. In practice, once icoFoam stops solving for U and the initial residual of p is less than the tolerance set in the fvSolution dictionary (typically ), the run has effectively converged and can be stopped once the field data has been written out to a time directory. For example, at convergence a sample of the log file from the run on the cavityHighRe case appears as follows in which the velocity has already converged after 1.395 s and initial pressure residuals are small; No Iterations 0 indicates that the solution of U has stopped:

2.1.8 High Reynolds number flow

View the results in paraFoam and display the velocity vectors. The secondary vortices in the corners have increased in size. The user can then increase the Reynolds number further by decreasing the viscosity and then rerun the case. The number of vortices increases so the mesh resolution around them will need to increase in order to resolve the more complicated flow patterns. In addition, as the Reynolds number increases the time to convergence increases. The user should monitor residuals and extend the endTime accordingly to ensure convergence.

The need to increase spatial and temporal resolution then becomes impractical as the flow moves into the turbulent regime, where problems of solution stability may also occur. Of course, many engineering problems have very high Reynolds numbers and it is infeasible to bear the huge cost of solving the turbulent behaviour directly. Instead Reynolds-averaged simulation (RAS) turbulence models are used to solve for the mean flow behaviour and calculate the statistics of the fluctuations. The standard model with wall functions will be used in this tutorial to solve the lid-driven cavity case with a Reynolds number of . Two extra variables are solved for: , the turbulent kinetic energy; and, , the turbulent dissipation rate. The additional equations and models for turbulent flow are implemented into a OpenFOAM solver called pisoFoam.

2.1.8.1 Pre-processing

Go back to the run directory and copy the cavity case in the $FOAM_RUN/tutorials/incompressible/pisoFoam/RAS directory (N.B: the pisoFoam/RAS directory), renaming it cavityRAS to avoid a clash with the existing cavity tutorial. Go into the new case directory.

Generate the mesh by running blockMesh as before. Mesh grading towards the wall is not necessary when using the standard model with wall functions since the flow in the near wall cell is modelled, rather than having to be resolved.

A range of wall function models is available in OpenFOAM that are applied as boundary conditions on individual patches. This enables different wall function models to be applied to different wall regions. The choice of wall function models are specified through the turbulent viscosity field, in the 0/nut file:

This case uses standard wall functions, specified by the nutWallFunction type on the movingWall and fixedWalls patches. Other wall function models include the rough wall functions, specified through the nutRoughWallFunction keyword.

The user should now open the field files for and (0/k and 0/epsilon) and examine their boundary conditions. For a wall boundary condition, is assigned a epsilonWallFunction boundary condition and a kqRwallFunction boundary condition is assigned to . The latter is a generic boundary condition that can be applied to any field that are of a turbulent kinetic energy type, e.g. , or Reynolds Stress . The initial values for and are set using an estimated fluctuating component of velocity and a turbulent length scale, . and are defined in terms of these parameters as follows:

where is a constant of the model equal to 0.09. For a Cartesian coordinate system, is given by:

(2.10)

where , and are the fluctuating components of velocity in the , and directions respectively. Let us assume the initial turbulence is isotropic, i.e. , and equal to 5% of the lid velocity and that , is equal to 5% of the box width, 0.1 , then and are given by:

These form the initial conditions for and . The initial conditions for and are and 0 respectively as before.

Turbulence modelling includes a range of methods, e.g. RAS or large-eddy simulation (LES), that are provided in OpenFOAM. The choice of turbulence modelling method is selectable at run-time through the simulationType keyword in momentumTransport dictionary (known as turbulenceProperties prior to OpenFOAM v8). The user can view this file in the constant directory:

The options for simulationType are laminar, RAS and LES. With RAS selected in this case, the choice of RAS modelling is specified in a RAS subdictionary. The turbulence model is selected by the model entry from a long list of available models that are listed in Section 7.2.1.1 . The kEpsilon model should be selected which is the standard model; the user should also ensure that turbulence calculation is switched on.

The coefficients for each turbulence model are stored within the respective code with a set of default values. Setting the optional switch called printCoeffs to on will make the default values be printed to standard output, i.e. the terminal, when the model is called at run time. The coefficients are printed out as a sub-dictionary whose name is that of the model name with the word Coeffs appended, e.g. kEpsilonCoeffs in the case of the kEpsilon model. The coefficients of the model, e.g. kEpsilon, can be modified by optionally including (copying and pasting) that sub-dictionary within the RAS sub-dictionary and adjusting values accordingly.

The user should next set the laminar kinematic viscosity in the transportProperties dictionary. To achieve a Reynolds number of , a kinematic viscosity of   is required based on the Reynolds number definition given in Equation 2.1 .

Finally the user should set the startTime, stopTime, deltaT and the writeInterval in the controlDict. Set deltaT to 0.005 s to satisfy the Courant number restriction and the endTime to 10 s.

2.1.8.2 Running the code

Execute pisoFoam by entering the case directory and typing “pisoFoam” in a terminal. In this case, where the viscosity is low, the boundary layer next to the moving lid is very thin and the cells next to the lid are comparatively large so the velocity at their centres are much less than the lid velocity. In fact, after  100 time steps it becomes apparent that the velocity in the cells adjacent to the lid reaches an upper limit of around 0.2  hence the maximum Courant number does not rise much above 0.2. It is sensible to increase the solution time by increasing the time step to a level where the Courant number is much closer to 1. Therefore reset deltaT to 0.02 s and, on this occasion, set startFrom to latestTime. This instructs pisoFoam to read the start data from the latest time directory, i.e.10.0. The endTime should be set to 20 s since the run converges a lot slower than the laminar case. Restart the run as before and monitor the convergence of the solution. View the results at consecutive time steps as the solution progresses to see if the solution converges to a steady-state or perhaps reaches some periodically oscillating state. In the latter case, convergence may never occur but this does not mean the results are inaccurate.

2.1.9 Changing the case geometry

A user may wish to make changes to the geometry of a case and perform a new simulation. It may be useful to retain some or all of the original solution as the starting conditions for the new simulation. This is a little complex because the fields of the original solution are not consistent with the fields of the new case. However the mapFields utility can map fields that are inconsistent, either in terms of geometry or boundary types or both.

As an example, let us copy the cavityClipped case from the tutorials directory in the user’s run directory, and change into the cavityClipped directory:

The case consists of the standard cavity geometry but with a square of length removed from the bottom right of the cavity, according to the blockMeshDict below:

Generate the mesh with blockMesh. The patches are set accordingly as in previous cavity cases. For the sake of clarity in describing the field mapping process, the upper wall patch is renamed lid, previously the movingWall patch of the original cavity.

In an inconsistent mapping, there is no guarantee that all the field data can be mapped from the source case. The remaining data must come from field files in the target case itself. Therefore field data must exist in the time directory of the target case before mapping takes place. In the cavityClipped case the mapping is set to occur at time 0.5 s, since the startTime is set to 0.5 s in the controlDict. Therefore the user needs to copy initial field data to that directory, e.g. from time 0:

Before mapping the data, the user should view the geometry and fields at 0.5 s.

Now we wish to map the velocity and pressure fields from cavity onto the new fields of cavityClipped. Since the mapping is inconsistent, we need to edit the mapFieldsDict dictionary, located in the system directory. The dictionary contains 2 keyword entries: patchMap and cuttingPatches. The patchMap list contains a mapping of patches from the source fields to the target fields. It is used if the user wishes a patch in the target field to inherit values from a corresponding patch in the source field. In cavityClipped, we wish to inherit the boundary values on the lid patch from movingWall in cavity so we must set the patchMap as:

---

---

---

---

The cuttingPatches list contains names of target patches whose values are to be mapped from the source internal field through which the target patch cuts. In this case, the fixedWalls patch is a noSlip condition so the internal values cannot be interpolated to the patch. Therefore the cuttingPatches list can simply be empty:

If the user does wish to interpolate internal values from the source case to the fixedWalls patch in the target case, a fixedValue boundary condition needs to be specified on the patch, whose value can then be overridden during the mapping process; the fixedWalls patch then needs to be included in the cuttingPatches list.

The user should run mapFields, from within the cavityClipped directory:

The user can view the mapped field as shown in Figure 2.14 . The fixedWalls patch has not inherited values from the source case as we expected. The user can then run the case with icoFoam.

2.1.10 Post-processing the modified geometry

Velocity glyphs can be generated for the case as normal, first at time 0.5 s and later at time 0.6 s, to compare the initial and final solutions. In addition, we provide an outline of the geometry which requires some care to generate for a 2D case. The user should select Extract Block from the Filter menu and, in the Parameter panel, highlight the patches of interest, namely the lid and fixedWalls. On clicking Apply, these items of geometry can be displayed by selecting Wireframe in the Display panel. Figure 2.15 displays the patches in black and shows vortices forming in the bottom corners of the modified geometry.