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 turbFoam solver will be used for turbulent, isothermal, incompressible flow.

2.1.1 Pre-processing

The pre-processing in OpenFOAM is done using FoamX, a JAVA GUI tool for managing OpenFOAM cases. The use of FoamX is described in more detail in chapter 5, but the command descriptions given in the text below should be sufficient to guide the user through the tutorial.

First the user must start up FoamX. Here, we are assuming that the FoamX host browser is being run on the local machine; otherwise, to run on a remote machine the user should consult chapter 5. Simply type

The user should now open the $OpenFOAM_RUN/tutorials/icoFoam directory by a double click on the root directory icon. This opens the directory revealing all the case directories located at the $OpenFOAM_RUN/tutorials/icoFoam path. The user should select the cavity case, again by a double click on the case icon.

The case opens in the same panel of the FoamX window and presents the user with a directory tree containing the description of the case: Dictionaries of input data and control parameters; a list of the Fields required in the problem e.g. velocity; and, the Mesh.

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

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 6.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 FoamX. This is done by clicking the right mouse button with the cursor over the case name cavity at the top of the directory tree. A menu opens to allow the user to select blockMesh from the meshUtilities -> generation sub-menu of the Foam Utilities menu as shown in Figure 2.4.

2.1.1.2 Boundary and initial conditions

Once the mesh generation is complete, click the right mouse button on Mesh in the case directory tree and select the Read Mesh & Fields function to load the mesh into FoamX, as shown in Figure 5.16. The names of the patches will appear in a folder named Patches. The user must click on each patch in turn to specify the physical boundary conditions as shown in Figure 5.17. Clicking on a patch brings up a window requesting the physical boundary types. For the cavity, the wall type should be selected for the fixedWalls and movingWall patches. The frontAndBack patch represents the front and back planes of the 2D case and therefore must be set as empty. Notice that as each boundary type is selected, the window displays the patch conditions for the primitive solution variables for each physical type, e.g. fixedValue for U and zeroGradient for p for the wall type.

Next select the Fields to set the internal and boundary fields. Clicking on U brings up a window requesting internal field, fields on the wall boundaries and the reference level as shown in Figure 2.5.

To enter an entry with multiple values, e.g. a vector, simply click on the entry and a button will appear on the right hand side of the that window; clicking on that button brings up a table containing the individual values of the entry that can be edited as normal as shown in Figure 2.5.

Clicking on p brings up a window requesting internal field and reference values; both should be set to 0. Clicking on U brings up a window requesting internal field and reference values and values on the wall boundaries. The internal and reference fields should both be set to . The velocity should be set to on the fixedWalls and on the movingWall.

All values of the field are stored relative to the reference level; in most circumstances, including here, it is set to zero but can be set to another value to improve accuracy in certain cases.

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 must ensure that the kinematic viscosity is set correctly by double clicking on 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 ,  1 , so that for  10,  0.01 . The correct dimensions for kinematic viscosity are specified with its value, in SI units. The dimension is described in terms of powers of SI base units of measurement , in this case , or

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. Firstly, the user must set the start/stop times and the time step for the run. 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 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 must 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 ; in fact, with hindsight, we discover that 0.5  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 . With a time step of 0.005 , 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 , 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 remaining entries in the controlDict are shown in Figure 2.6.

2.1.1.5 Discretisation and linear-solver settings

The user specifies the choice of finite volume discretisation schemes in the fvSchemes dictionary. The specification of the linear equation solvers and tolerances and other algorithm controls is made in the fvSolution dictionary. The user is free to view these dictionaries but we do not need to discuss them at this stage as they are pre-configured to enable the tutorials to run satisfactorily.

2.1.1.6 Saving data to file

The mesh has now been generated, the boundary conditions and fields have been initialised and the control parameters and material properties have been set. This data must be saved to the case files by clicking the menu button with the floppy disk icon at the top of the FoamX window.

For the remainder of the manual:
There are menu buttons at the top of the FoamX window. To check which function a button performs, hold the cursor over the button for 1  and a text description will appear.

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 dxFoam, the post-processing tool supplied with OpenFOAM. The dxFoam post-processing is started by typing at a command prompt

The bounding box for the geometry appears in the Image window. Reset the view by selecting Reset from the Options menu in the Image window; this will display the front view in the - plane. Now select the Mesh window from the menu in the Window Control window. In the Mesh window, increase the grid Opacity to, say, 0.5. The mesh appears in the Image window as shown in Figure 2.8.

The mesh view can be manipulated from the View Control window selected from the Options menu in the Image window. The user can experiment with the controls, most of which are fairly self-explanatory but which are explained in full in Table 7.2. The user should select Off Front from the Set View menu in order to verify that the mesh is, in fact, 3-dimensional with 1 cell through the thickness.

To exit dxFoam, go the the Display window and select Quit from the File menu, clicking No in the Save Confirmation window.

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. There are two ways that this can be done: either by clicking the Start Calculation Now button () in FoamX; or by typing at a command prompt:

2.1.4 Post-processing

As soon as results are written to time directories, they can be viewed using dxFoam. Open dxFoam as before and select the cavity case as described in subsection 2.1.2. In order to visualise field data more easily, it is advisable to switch off the mesh Grid by setting grid opacity to 0.0 in the Mesh window, accessible from the Window Control window as shown in Figure 2.9.

2.1.4.1 Contour plots

To view a contour plot of pressure, open the Field 1 visualisation window from the Window Control window.

Figure 2.10 displays the operations to display contour plots described in the following text. Select volScalarField from the Type menu and p from the Name menu. The plane in which the contour plot is drawn must be coincident with the 2D mesh plane. Therefore open the Map Plane, by clicking the button at the top of the Field 1 window, and ensure the plane 1 normal is set to .

Now return to the Field 1 visualisation window and check the GO box. The initial pressure field will appear inside the square domain, all blue since the initial field is 0 everywhere. Now go into the Read Data window, select the Time from stepper option and increment the time steps by clicking the Right Arrow button of Stepper. The pressure field solution at 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.11.

By default, none of the buttons are selected in the Options panel which means that the contours are plotted on the plane within the domain volume, i.e. across the complete plane in this case. The contour plot grades the colours between data points to give a smooth distribution in color across the domain.

Alternatively the user may fill complete volumes with colour rather than colouring a map plane by switching on the Volume button in the Options panel. If the Cell button is also selected, a single value will be attributed to each cell so that each cell will be denoted by a single colour with no grading. The user should experiment with these buttons in the Options panel.

The user may find that the image colours are improved by selecting the Lighting window from the Window Control window and switching the spot light Intensity to 0.0.

In the Field 1 visualisation window, set the opacity to a low value, e.g. 0.4 and ensure the contours are switched on in the Contour Attributes panel. The contours appear in the image; they can be switched off in the Contour Attributes panel, or exclusively switched on by selecting only. The contours can be further adjusted in the Map Plane window. In the bottom panel, the user can either choose that a contour line is drawn for a single value of pressure or that several contours -- specified in the number of C box -- are drawn at equal intervals across the range of pressure. Try increasing the number of contours and see the effect on the image.

A colour bar appears on the screen that can be adjusted from the Color range window also selected by clicking the button at the top right of the Field 1 window. The range of the colour bars is adjusted to scale automatically to the maximum and minimum values of the field. The user can disable this by checking the Limit colors box for a particular field and setting appropriate maximum and minimum values. The orientation of the legend can be switched between horizontal and vertical.

2.1.4.2 Vector plots

Before we start to plot the vectors of the flow velocity, first switch off the isoline plot of pressure field by checking the disable box in the Options panel of the Field 1 visualisation window. It may also help to reduce the Opacity of the Bounding Box in the Mesh window to 0.2 so that the vectors can be seen clearly on the cavity boundaries.

Now open the Field 2 visualisation window and choose volVectorField from the Type menu, U from the Name menu and check the GO box. The vectors may either be plotted on a specified plane or within the cell volumes. For complex geometries, it can be difficult to generate a clear image if a vector is plotted for every cell, in which case plotting on a plane is often preferred. Here, we can try both and, in preparation, the user should set the plane for this field to be coincident with the 2D mesh plane by entering for the plane 2 normal in the Map Plane window.

The velocity field should now load and provided the time step is greater than 0.0, some vectors should appear in the Image. The user will need to modify the appearance of the vectors in the Glyph Attributes panel of the the Field 2 visualisation window. Firstly, the glyph scale should be reduced to around 0.004. The glyph type can also be changed between: arrow2D, a 2-dimensional arrow; speedy, a 2-dimensional line; and standard and spiffy, both 3-dimensional arrows. The user is encouraged to experiment with these options with the view set to Off Front in the View Control window in order to see clearly the difference between the types of vector. It may also be useful to select normalise glyphs when the glyphs are all displayed with the same length with only the colour to distinguish the magnitude. This is particularly useful where there is a significant variation in the magnitude of the vector field across the domain, as in this example. A vector plot with normalised glyphs is shown in Figure 2.12.

By default, none of the buttons are selected in the Options panel which means that the vectors are plotted on cell edges/vertices that intersect the plane. In order to plot the vectors at the cell centres both the the Volume and Cell buttons should be checked. The user should experiment with the buttons in Options panel.

2.1.4.3 Streamline plots

dxFoam has the function for drawing streamlines of velocity fields selected in the Field 1 visualisation window. Therefore, in preparation of drawing streamlines of the flow, switch off the velocity vector plot for Field 2 by disabling the GO button if the Field 2 visualisation window.

Now open the Field 1 visualisation window and select volVectorField from the Type menu and U from the Name menu. Check the GO box and check the disable box in the Options panel. It may also help to increase the Opacity of the Bounding Box in the Mesh window to 1.0 again.

Now open the Streamlines window from the Window Control window. The user must set the parameters as shown in Figure 2.13. In particular, the user should check the Show Grid button in order to view the streamline rake; the meaning of streamline rake and its setting is discussed further in subsection 7.2.6.

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

Close the cavity case in FoamX by clicking the Close Case button (). 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. In the FoamX case browser window, simply place the cursor over the cavity case and click the right mouse button to bring up a menu from which Clone Case can be selected. The user should enter the root path, the name of the new case, cavityFine, and the icoFoam application class. The times option specifies which time directories are copied into the cloned case. In this example, we require no time directory from cavity, otherwise after refining the mesh, the size of the copied fields will be inconsistent with the size of the mesh. Therefore we select noTime. By clicking OK, the new case is created and presented in the case browser window. The user should open the cavityFine case.

2.1.5.2 Creating the finer mesh

We now wish to increase the number of cells in the mesh by using blockMesh. As before, select blockMesh from the meshUtilities -> generation sub-menu of the Foam Utilities menu by clicking the right mouse button with the cursor over the case name at the top of the directory tree. The blockMesh window appears in which the user should press the Edit Dictionary button. The user should then select the value cell of the the blocks item and in the blocks window, select (the only) block 0. This produces a new window with a table of entries for the block. The first entry, hex, describes the type of block, a hexahedron -- in fact the only option in blockMesh-- which is specified by an ordered list of vertex labels. The second entry, cellDensity, is the mesh density: 20 cells in the -direction; 20 in the -direction; and 1 in the -direction, since it is a 2 dimensional problem. The third entry, expansionRatio specifies the mesh grading which is discussed in subsection 2.1.6.

In this case we wish to increase the mesh density to 41 cells in the and directions. Select the cellDensity entry and edit the elements, replacing 20 20 1 by 41 41 1. Then close the dictionary editing window and the dictionary will save automatically. Press the Start Run button and the mesh is generated with the running status of blockMesh reported in the terminal window as before. To view the mesh, use dxFoam as described in subsection 2.1.2. The user may be prompted to read the new mesh into FoamX. The user must do so if he/she wishes, because the number of faces in the the boundary patch definitions have changed.

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  in the controlDict dictionary and startFrom should be set to startTime. The user should save these changes.

The case is ready to run mapFields. Click the right mouse button with the cursor over the case name cavityFine at the top of the directory tree. A menu opens to allow the user to select mapFields from the preProcessing sub-menu of the Foam Utilities menu. A window containing arguments to the mapFields utility opens as shown in Figure 2.15. The name of the target case, with full root path, into which the results are being mapped are set correctly by default in <rootAndCase>. The user must enter the root path and name of the source case in <sourceRootAndPath>, e.g. as the example shows in Figure 2.15.

The user must select the -consistent option from the arguments by selecting on in the [consistent] value box. The mapFields utility will run when the user clicks Start Run. The progress messages in the terminal window should tell the user that the data has been interpolated from the cavity case to the cavityFine case.

2.1.5.4 Control adjustments

To maintain a Courant number of less that 1, as discussed in subsubsection 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 to 0.0025  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 the 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 . Make sure these settings are correct and then save the case.

2.1.5.5 Running the code as a background process

The user should experience running icoFoam as a background process. Press the Start Calculation button () and a window appears. With the background button checked, click Start Run. The case runs in the background and is complete when a line of text appears in the terminal window beginning with Finished doing.... The user can then view the log file by clicking View Log.

2.1.5.6 Vector plots with a user grid

The results on the finer mesh can now be compared with those of the coarser mesh. Select dxFoam and open the cavityFine case. Go to the final results from the run, i.e. those in the 0.7 time directory. Display the velocity vectors in the Field 1 window, making sure that they are not normalised -- the normalise glyphs box should be unchecked in the Glyph Attributes panel of the window. By default dxFoam plots one velocity vector per cell; since we want to compare results from two cases with different numbers of cells, it is useful to fix the number of vectors. Therefore the user grid box should be checked in the Glyph Attributes panel of the window and the density set to a value that displays a reasonable number of vectors, e.g. 500. The user can now switch from the cavityFine case to the cavity case to compare the velocity fields.

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 and is released with some standard utilities for this purpose, namely Ucomponents and magU. When Ucomponents is run on a case, say 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 the Ucomponents utility on both cavity and cavityFine cases. Like all utilities, Ucomponents can be executed from the menus opened by clicking the right mouse button with the cursor over the case name in FoamX. Ucomponents is located in the postProcessing -> velocityField sub-menu; execute it on cavityFine. For cavity the user may wish to execute Ucomponents from the command line by the following command:

Let us plot for the cavityFine case along a vertical line through the centre of the cavity. First, load the case and select Ux from volScalarField in the Field 1 window. Next return to the Graph window; the Grid Point of the graph line should coincide with the centre of the geometry, i.e. . The Vector has magnitude of half the cavity height and is in the vertical direction, i.e. . Set these values and check the Reveal Line box as shown in Figure 2.16; in this figure we have selected a White Background from the Write Control window, opened from the Window Control. A graph line should appear, dividing the cavity along its vertical centreline. Check the Do Plot box and a plot of against should appear in the Display window similar to that shown in Figure 2.16. Note that the Display window may be hidden behind other windows and the graph can be switched off by deselecting the Do Plot box.

Set Plot Labels accordingly. The user can now switch between the cavityFine and cavity case to compare plots.

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 subsubsection 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 $OpenFOAM_RUN/tutorials/icoFoam directory.

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

Once familiar with the blockMeshDict file for this case, the user can execute blockMesh from the command line using a command of the form:

The graded mesh can be viewed as before using dxFoam as described in subsection 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 subsubsection 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 first and last 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 .

On this occasion we shall demonstrate the fact that any changes can be made to the case dictionaries by simply editing the relevant file. Here we wish to edit the time and control information which is stored in the cavityGrade/system/controlDict file. Open this file in an editor of your choice. As stipulated in the previous paragraph ensure that the time step deltaT is set to 2.5e-3 and the writeInterval is set to 40.

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 , i.e. the endTime should be 0.8.

2.1.6.3 Mapping fields, running the code and post-processing

As in subsubsection 2.1.5.3, use mapFields to map the final results from case cavityFine onto the mesh for case cavityGrade. The mapFields utility is executed by a line of the form

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 no secondary vortices in the flow. The Reynolds number will now be raised to 50, at which point secondary vortices should start to appear in the corners of the domain if the mesh is fine enough. The coarsest mesh in case cavity will be used initially. The user should make a copy of the cavity case and name it cavityHighRe. The user can use the Clone Case function in FoamX as described in subsubsection 2.1.5.1 or simply copy the cavity case directory by typing:

2.1.7.1 Pre-processing

Let us return to use FoamX for managing the cavityHighRe case. If the user has created the cavityHighRe case by making a copy of cavity as described above, the case will probably not appear in the case directory tree panel of the case browser window. To make it appear, the user must select Refresh Case Browser by either: pressing the right mouse button with the cursor over the Licenced Hosts icon at the top of the case directory tree; or, by selecting the Refresh Case Browser button from the menu buttons.

From the case directory tree, open the cavityHighRe case and edit the transportProperties dictionary. Since the Reynolds number is required to be increased by a factor of 5, decrease the kinematic viscosity by a factor of 5, i.e. to  . We can 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 1.5 . Save the case.

2.1.7.2 Running the code

Run icoFoam for this case from the case directory and view the run time information.

The velocity has already converged after 0.925  and initial pressure residuals are small. The code will next write out results at 1  so execution must continue until then. The user may stop a run either from the command line or from FoamX. To stop the run from the command line, first type ps to find out the process identification (PID) number of icoFoam (or ps -e if icoFoam is not displayed in this list, for example if it is typed in a different terminal window). To stop icoFoam, use the kill command followed by the PID number of icoFoam using the -9 option if necessary:

The alternative of monitoring jobs is by using the process editor in FoamX. The process editor is opened using the Start Process Editor function () in the Case Browser window. The user can view the jobs that are currently running in the runningJobs table; the table can be refreshed by clicking the read button. If the user wishes to stop a job that is running, they should highlight the job in the table and click the kill button. They will be prompted to confirm the killing of the job.

2.1.7.3 Post-processing

View the results in dxFoam. Display the velocity vectors using the Normalise Glyphs option in the Glyph Attributes panel of the Field 1 window and set the Glyph Scale to so that the vectors are reasonable lengths. Look for secondary vortices in the corners of the cavity; for , there should be no secondary vortices. Go back to FoamX and increase the Reynolds number further by decreasing viscosity. Save, re-run the case and look for secondary vortices in dxFoam, e.g. Figure 2.18 shows vortices appearing for a Reynolds Number of 100. Note that as the Reynolds number increases the time to convergence also increases. The user must extend the endTime accordingly.

While displaying the vectors, increment the time sequence from the original low Reynolds number solution at 0.5 . Watch the higher Reynolds number flow develop to a converged solution.

2.1.7.4 Exercise

Continue increasing the Reynolds number and watch more vortices appear. As the number of vortices increases the mesh resolution around them will need to increase in order to resolve the more complicated flow patterns. At what Reynolds number does the flow become unsteady?

2.1.8 High Reynolds number flow

subsection 2.1.7 demonstrates that solving higher Reynolds number flows becomes increasingly expensive if an accurate solution is required. 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 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 turbFoam.

2.1.8.1 Pre-processing

In FoamX, open the cavity case in the $HOME_RUN/tutorials/turbFoam directory (N.B: the $OpenFOAM_RUN/tutorials/turbFoam directory). Generate the mesh by running blockMesh from within FoamX 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.

Set the boundary conditions using FoamX as described in subsubsection 2.1.1.2. The selection of the wall type boundary condition assigns a zeroGradient boundary condition to both and . To set the initial conditions, select the fields as described previously. The initial conditions for and are and 0 respectively as before. However, positive values for and must be given to avoid division by 0 in the solution algorithm. We can specify reasonable initial conditions for and in terms of 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 20% of the box width, 0.1 , then and are given by:

Set these initial conditions for and .

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.

To select the turbulence model open the turbulenceProperties dictionary. The turbulence model is selected by the turbulenceModel entry. It gives a long list of available models that are listed in Table 3.9. The user should select the kEpsilon which is is the standard model; the user should also ensure that turbulence calculation is switched on. The coefficients relating to the model are stored in a standard dictionary under kEpsilonCoeffs; the model also uses the wallFunctionCoeffs.

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

2.1.8.2 Running the code

Execute turbFoam using any of the methods described earlier in this tutorial. 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  and, on this occasion, set startFrom to latestTime. This instructs turbFoam to read the start data from the latest time directory, i.e. 10.0. The endTime should be set to 20  since the run converges a lot slower than the laminar case. Restart the run as before and monitor the convergence of the solution.

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 open in FoamX the case cavityClipped in the icoFoam directory which 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 and read the mesh into the case browser as described previously. Ensure that the patch types are set correctly: the fixedWalls and lid patches should be set to wall, the lid patch being the movingWall patch of the cavity case renamed for sake of clarity. Now save the case, which saves the field data into time directory 0.5 since the startTime is set to 0.5 in the controlDict. View the geometry and fields at 0.5 .

Now we wish to map the velocity and pressure fields from cavity onto the new fields of cavityClipped. Simply open the mapFields utility as before in FoamX and edit the arguments: set the root and case of the source and target but make sure the -consistent option is switched off.

Return to the mapFields utility window and select Edit Dictionary. A window opens containing the arguments followed by 2 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 we will include the fixedWalls to demonstrate the interpolation process.

Now the user should run mapFields, either from FoamX or from the command line:

The user should now select the Read Mesh function to load the new fields into FoamX . The user can view the mapped field as shown in Figure 2.19. The boundary patches have inherited values from the source case as we expected. Having demonstrated this, however, we actually wish to reset the velocity on the fixedWalls patch to . Open the U field in FoamX, select the fixedWalls patch and change the field from nonuniform to uniform . Now run the case with icoFoam. The solution at 0.6 is shown in Figure 2.20. Vortices are seen forming in the bottom corners, which would be expected for the modified geometry.