Session 3. Numerical problem C. Fernandes, L.L. Ferr as, J.M. N - - PowerPoint PPT Presentation

Session 3. Numerical problem

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Institute for Polymers and Composites (i3N), University of Minho

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

CloudPYME Project

El proyeto CloudPYME (id: 0682 CLOUDPYME2 1 E) est´ a cofi- nanciado por la Comisi´

• n Europea a trav´

es de el Fondo Europeo de Desarrollo Regional (FEDER), dentro de la tercera convocatoria de proyectos del Programa Operativo de Cooperaci´

• n Transfronteriza

Espa˜ na-Portugal 2007-2013 (POCTEP).

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Problem description

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 Fig. 1 in which all the boundaries of the square are walls. The top wall moves

• n the x-direction at a speed of 1 m/s while the other 3 are station-

ary.

Figure 1: Geometry of the lid driven cavity

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Pre-processing

Cases are setup in OpenFOAM by editing case files. A case being simulated involves data for mesh, fields, properties, control pa- rameters, etc. In OpenFOAM this data is stored in a set of files within a case directory. To access the cavity case, the user should change to the case directory: >> cd $FOAM RUN/cavity

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

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

Figure 2: Block structure of the mesh for the cavity

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

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

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

Figure 3: BlockMeshDict file for the cavity case

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

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

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

The file log blockMesh, Fig. 4, contains the information related to the generation of the mesh:

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

Figure 4: blockMesh output

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

To check if any error is obtained during the mesh generation we can type: >> checkMesh > log checkMesh

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

The file log checkMesh, Fig. 5, contains the information related to the mesh statistics and topology:

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh generation

Figure 5: checkMesh output

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Boundary and initial conditions

Once the mesh generation is complete, the user can look at initial fields set up for this case. The case is set up to start at time t = 0 s, so the initial field data is stored in a 0 sub-directory of the cavity

• directory. The 0 sub-directory contains 2 files, p and U, one for

each of the pressure (p) and velocity (U) fields whose initial values and boundary conditions must be set.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Boundary and initial conditions

Let us examine file p, Fig. 6:

Figure 6: Pressure initial conditions

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Boundary and initial conditions

For those using the incompressible solvers (icoFoam, simpleFoam, pisoFoam, pimpleFoam), we need to remind that the pressure used by these solvers is the modified pressure or P = p/ρ. So do not forget to multiply the pressure by the density in order to get the physical pressure.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Boundary and initial conditions

Let us examine file U, Fig. 7:

Figure 7: Velocity initial conditions

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Physical properties

The physical properties for the case are stored in dictionaries whose names are given the suffix . . . Properties, located in the constant directory tree. For this case, the only property that must be speci- fied is the kinematic viscosity which is stored from the transport- Properties dictionary.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Physical properties

The keyword for kinematic viscosity is nu (ν). This case will be run with a Reynolds number of 10, where the Reynolds number is defined as: Re = d|U| ν Here d = 0.1 m, |U| = 1m/s, so that for Re = 10, ν = 0.01m2/s. The correct file entry, Fig. 8, for kinematic viscosity is thus specified below:

Figure 8: Transport properties file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Control

Input data relating to the control of time and reading and writing

• f the solution data are read in from the controlDict dictionary.
• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Control

To achieve temporal accuracy and numerical stability when running icoFoam, a Courant number of less than 1 is required. The Courant number is defined for one cell as: Co = δt|U| δx 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: δt = Coδx |U| = 1 × 0.005 1 = 0.005s

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Control

For this case, the entries in the controlDict, Fig. 9, are shown below:

Figure 9: Control dictionary file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

The user specifies the choice of finite volume discretisation schemes in the fvSchemes dictionary in the system directory. The specifica- tion of the linear equation solvers and tolerances and other algorithm controls is made in the fvSolution dictionary, similarly in the system directory.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

For this case, the entries in the fvSchemes, Fig. 10, are shown below:

Figure 10: fvSchemes dictionary file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

Other discretization schemes are available, Fig. 11:

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

Figure 11: fvSchemes options

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

For this case, the entries in the fvSolution, Fig. 12, are shown below:

Figure 12: fvSolution dictionary file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

• Fig. 13 shows other options for the fvSolution file:

Figure 13: fvSolution options

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Discretisation and linear-solver settings

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.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Viewing the mesh

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

• r

>> touch cavity.foam >> paraview cavity.foam

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Viewing the mesh

• Fig. 14 shows the mesh visualized in Paraview.

Figure 14: Mesh

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Running an application

OpenFOAM applications can be run in two ways: as a foreground process, i.e. one in which the shell waits until the command has fin- ished before giving a command prompt; as a background process,

• ne which does not have to be completed before the shell accepts

additional commands.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Running an application

On this occasion, we will run icoFoam in the foreground. The icoFoam solver is executed by entering the case directory and typing >> icoFoam Later on this tutorial we will run icoFoam in background process.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Running an application

Sometimes it is difficult to determine all the available run time selection classes. To make it easier we can use the banana trick. For example, in the pressure matrix solver change the pressure solver to banana, like in Fig. 15.

Figure 15: Banana trick

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Running an application

Running the icoFoam application again will produce the following error, Fig. 16.

Figure 16: Banana trick output

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Paraview

Try to export pressure and velocity profiles. For example, velocity field with vectors, Fig. 17, can be visualized on Paraview.

Figure 17: Velocities in the cavity case

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Instead, if it is necessary to compare the numerical results with analytical ones we should extract some field values from Open- FOAM results files. In that case we can use sample utility to do post-processing.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

A good set of data for comparison is the data of Ghia et al. (1982), since it includes tabular results for various of Reynolds numbers beginning in Re = 100. The case solved so far has had a Reynolds number of 10. We will now increase the Reynolds number to 100 so we can compare against data from Ghia et al. (1982). The user should make a copy of the cavity case and name it cavityHighRe by typing: >> cd $FOAM RUN/ >> cp -r cavity cavityHighRe

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Since the Reynolds number is required to be increased by a factor

• f 10, decrease the kinematic viscosity by a factor of 10, i.e.

to 0.001 m2/s. We can now run this case by restarting from the solution at the end

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

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Run icoFoam for this case from the case directory and view the run time information. To run the job in the background type >> icoFoam > log simulation & Now the information of the simulation in stored in the log file. We can assess this information typing >> cat log simulation

• r

>> tail -n 500 log simulation

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

In the previous run 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 10−6), the run has effectively converged and can be stopped once the field data has been written out to a time directory.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

For example, from the log file, Fig. 18, of the cavityHighRe case appears that the velocity has already converged after 1.43s and ini- tial pressure residuals are small; No Iterations 0 indicates that the solution of U does not have to be improved.

Figure 18: Log file for high Reynolds number cavity case

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

We can check the convergence of the simulation using the utility gnuplot to draw the residuals of pressure and velocity profiles. Fig. 19 shows the content of the Residuals file

Figure 19: Residuals file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

The command to see the graphic is >> gnuplot Residuals

• Fig. 20 shows the residuals behavior for the cavityHighRe case.

Figure 20: Monitoring simulation residuals

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Now we will sample the velocity values along the y-direction in the middle location of x-direction. To do that we need to use a sampleDict file to give information of which fields we want to extract and the respective position. First we should find an sam- pleDict file to do our changes. >> cd $FOAM TUTORIALS >> grep -r “sampleDict” *

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

With the last commands we can see that inside the /heatTrans- fer/buoyantSimpleFoam/buoyantCavity/ tutorial there is a sample- Dict file that we can use. Next, we will copy that file to our system cavityHighRe folder >> cd $FOAM RUN >> cp $FOAM TUTORIALS/heatTransfer/ buoyantSimpleFoam/buoyantCavity/system/sampleDict cavity- HighRe/system

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Now we should change the contents of the sampleDict file, Fig. 21, to match the following

Figure 21: SampleDict file

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Finally, inside the cavityHighRe case we type >> sample –latestTime And it will be created a folder named postProcessing with the subfolder sets and inside this the subfolder 2 with the file containing velocity components on our specific location. Copy the file U OpenFOAM U.xy to the cavityHighRe case >> cp postProcessing/sets/2/U OpenFOAM U.xy .

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

We can now plot this OpenFOAM numerical values against Ghia et al. values. For that we need to create one file with the Ghia et

• al. values for Re = 100, Fig. 22, and one gnuplot file, Fig. 23, to

plot numerical and Ghia et al. values simultaneously. The two files should look

Figure 22: Ghia et al. values

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Figure 23: Numerical comparison cavity file

Finally, to obtain the graphic comparing the results from Ghia et al. and OpenFOAM type >> gnuplot cavity

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Sample

Here we have plotted u/U on the vertical line x/L = 0.5, Fig. 24. We have reasonably good agreement with the previous results from Ghia et al.

Figure 24: Ghia et al. vs. OpenFOAM cavity results

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

The mesh resolution will now be increased by a factor of two in each direction. The results from the coarser mesh will be mapped

• nto the finer mesh to use as initial conditions for the problem. The

results will be compared again to the ones of Ghia et al.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

First, we need to create a new case named cavityHighReFine from the cavityHighRe. The commands are the following: >> cd $FOAM RUN/ >> mkdir cavityHighReFine >> cp -r cavityHighRe/constant cavityHighReFine >> cp -r cavityHighRe/system cavityHighReFine >> cd cavityHighReFine

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

Now we will increase the number of cells in the mesh by editing the blockMeshDict file. The block specification inside this file should be changed. We should change (20 20 1) to (40 40 1). Then, the new refined mesh should be created by running >> blockMesh

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

We will use now the mapFields utility to map the pressure and velocity fields of the cavityHighRe onto the cavityHighReFine. Since geometry and boundary types of both source and target fields are identical we can use this function, otherwise we should run the case from the beginning. Since the final results of the cavityHighRe case are stored in the 2.0 directory, we should set the startTime variable in the controlDict file to 2.0 and startFrom to startTime. Finally, we can map the fields between the two cases running >> mapFields ../cavityHighRe -consistent

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

Before running the case, we should make some modifications in the controlDict file. To maintain a Courant number less than 1, the step size must be halved since the size of all cells has halved. Therefore, we should set deltaT to 0.0025 s. In addition, we can specify data output at fixed intervals of time. For that, we should change writeinterval to 0.1 and writeControl to runTime. Finally, as we use the final solution of the coarser mesh we can set endTime to 2.5 s because its expected a short period to achieve convergence.

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

From the cavityHighReFine directory we can run the case as >> icoFoam > log simulation & When the simulation finish repeat the same steps to plot Ghia et

• al. results versus the ones obtained with the finer mesh. We should
• btain a graphic, Fig. 25, similar to
• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem

Mesh refinement study

Figure 25: Ghia et al. vs. OpenFOAM cavity results with refined mesh

• C. Fernandes, L.L. Ferr´

as, J.M. N´

• brega

Session 3. Numerical problem