2.3 Breaking of a dam

Advanced Search

Send us a comment/query

Back to main web site

OpenFOAM guides

- User guide

- Programmer’s guide

Index

Changes from OpenFOAM 2.2 to 2.3

Trademarks in the guides

[next] [prev] [prev-tail] [tail] [up]

2.3 Breaking of a dam

In this tutorial we shall solve a problem of simplified dam break in 2 dimensions using the interFoamThe feature of the problem is a transient flow of two fluids separated by a sharp interface, or free surface. The two-phase algorithm in interFoam is based on the volume of fluid (VOF) method in which a specie transport equation is used to determine the relative volume fraction of the two phases, or phase fraction $g \special {t4ht=$ , in each computational cell. Physical properties are calculated as weighted averages based on this fraction. The nature of the VOF method means that an interface between the species is not explicitly computed, but rather emerges as a property of the phase fraction field. Since the phase fraction can have any value between 0 and 1, the interface is never sharply defined, but occupies a volume around the region where a sharp interface should exist.

The test setup consists of a column of water at rest located behind a membrane on the left side of a tank. At time $t = 0 s \special {t4ht=$ , the membrane is removed and the column of water collapses. During the collapse, the water impacts an obstacle at the bottom of the tank and creates a complicated flow structure, including several captured pockets of air. The geometry and the initial setup is shown in Figure 2.28.

$0.584 m water column 0.584 m 0.4 m 0.048 m 0.024 m 0.2 m 0.092 m \special {t4ht=$

Figure 2.28:

Geometry of the dam break.

2.3.1 Mesh generation

The user should start up FoamX as normal and open the damBreak case in their own $OpenFOAM_RUN/tutorials/interFoam directory. Generate the mesh running blockMesh as described previously. The damBreak mesh consist of 5 blocks; the blockMeshDict entries are given below.

31  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32
33  convertToMeters 0.146;
34
35  vertices
36  (
37      (0 0 0)
38      (2 0 0)
39      (2.16438 0 0)
40      (4 0 0)
41      (0 0.32876 0)
42      (2 0.32876 0)
43      (2.16438 0.32876 0)
44      (4 0.32876 0)
45      (0 4 0)
46      (2 4 0)
47      (2.16438 4 0)
48      (4 4 0)
49      (0 0 0.1)
50      (2 0 0.1)
51      (2.16438 0 0.1)
52      (4 0 0.1)
53      (0 0.32876 0.1)
54      (2 0.32876 0.1)
55      (2.16438 0.32876 0.1)
56      (4 0.32876 0.1)
57      (0 4 0.1)
58      (2 4 0.1)
59      (2.16438 4 0.1)
60      (4 4 0.1)
61  );
62
63  blocks
64  (
65      hex (0 1 5 4 12 13 17 16) (23 8 1) simpleGrading (1 1 1)
66      hex (2 3 7 6 14 15 19 18) (19 8 1) simpleGrading (1 1 1)
67      hex (4 5 9 8 16 17 21 20) (23 42 1) simpleGrading (1 1 1)
68      hex (5 6 10 9 17 18 22 21) (4 42 1) simpleGrading (1 1 1)
69      hex (6 7 11 10 18 19 23 22) (19 42 1) simpleGrading (1 1 1)
70  );
71
72  edges
73  (
74  );
75
76  patches
77  (
78      wall leftWall
79      (
80          (0 12 16 4)
81          (4 16 20 8)
82      )
83      wall rightWall
84      (
85          (7 19 15 3)
86          (11 23 19 7)
87      )
88      wall lowerWall
89      (
90          (0 1 13 12)
91          (1 5 17 13)
92          (5 6 18 17)
93          (2 14 18 6)
94          (2 3 15 14)
95      )
96      patch atmosphere
97      (
98          (8 20 21 9)
99          (9 21 22 10)
100          (10 22 23 11)
101      )
102  );
103
104  mergePatchPairs
105  (
106  );
107
108  // ************************************************************************* //

2.3.2 Boundary conditions

Set the boundary conditions using the same procedure as described before. This should correspond to the description in Figure 2.28 although there is an important issue relating to the walls: the interFoam solver includes modelling of surface tension at the contact point between the interface and wall surface. If the user selects the wall boundary type in FoamX, they will notice that the gamma ( $g \special {t4ht=$ ) field is assigned a gammaContactAngle boundary condition. The user must then specify the following: a static contact angle, theta0 $h0 \special {t4ht=$ ; leading and trailing edge dynamic contact angles, thetaA $hA \special {t4ht=$ and thetaR $hR \special {t4ht=$ respectively; and a velocity scaling function for dynamic contact angle, uTheta.

In this tutorial we would like to ignore surface tension effects between the wall and interface. We can do this by setting the static contact angle, $o h0 = 90 \special {t4ht=$ and the velocity scaling function to 0. However, there is also an option which we shall choose here to select a wallNoSurfaceTension boundary type for the walls. Rather than use the gammaContactAngle boundary condition for gamma, this specifies a zeroGradient type instead.

The top boundary is free to the atmosphere and so is given an atmosphere boundary type; the defaultFaces representing the front and back planes of the 2D problem, is, as usual, an empty type.

2.3.3 Setting initial field

Unlike the previous cases, we shall now specify a non-uniform initial condition for the phase fraction $g \special {t4ht=$ where

{ 1 for the liquid phase g = 0 for the gas phase \special {t4ht=

(2.15)

This will be done by creating a OpenFOAM utility specifically for this purpose. The utility, setGammaDamBreak, has already been created and can be found in the setGammaDamBreak subdirectory of the tutorials/interFoam/damBreak directory. Go into this directory and first of all view the setGammaDamBreak/Make/files file:

1  setGammaDamBreak.C
2
3  EXE = $(OpenFOAM_USER_APPBIN)/setGammaDamBreak

This is simply a statement that the utility is complied using the setGammaDamBreak.C file and the compiled executable is named setGammaDamBreak, stored in the $OpenFOAM_USER_APPBIN directory, which is on the system path.

The code works as follows: the included file createFields.H specifies the gamma field to be read in from file:

1      Info<< "Reading field gamma\n" << endl;
2      volScalarField gamma
3      (
4          IOobject
5          (
6              "gamma",
7              runTime.timeName(),
8              runTime,
9              IOobject::MUST_READ,
10              IOobject::AUTO_WRITE
11          ),
12          mesh
13      );

The IOobject::MUST_READ specifies that the field named gamma will be read from the time directory. The field value of gamma can now be changed. In our case, we shall specify gamma to be equal to 1 in the bottom left part of the mesh and equal to 0 elsewhere. Once this operation is performed, the boundary condition on gamma will be updated and the new field will be written out into the 0 time directory, thereby overwriting the old gamma field. The setGammaDamBreak.C file therefore reads:

1  // The OpenFOAM Project // File: util/setGamma/setGamma.C
2  /*
3  -------------------------------------------------------------------------------
4   =========         | Application
5   \\      /         |
6    \\    /          | Name:   setGamma
7     \\  /           | Family: CFD
8      \\/            |
9      F ield         | OpenFOAM version: 2.3
10      O peration     |
11      A and          | Copyright (C) 2000-2007 Nabla Ltd.
12      M anipulation  |          All Rights Reserved.
13  -------------------------------------------------------------------------------
14  APPLICATION
15      setGamma
16
17  DESCRIPTION
18      Set the initial field of gamma for the dam break problem.
19
20  AUTHOR
21      Hrvoje Jasak.
22
23  -------------------------------------------------------------------------------
24  */
25
26  #include "database.H"
27
28  #include "fvCFD.H"
29  #include "OSspecific.H"
30
31  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32
33  int main(int argc, char *argv[])
34  {
35
36  #   include "setRootCase.H"
37
38  #   include "createDatabase.H"
39  #   include "createMesh.H"
40  #   include "createFields.H"
41
42  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
43
44      const vectorField& centres = mesh.C().internalField();
45
46      gamma.internalField() =
47          neg(centres.component(vector::X) - 0.146)*
48          neg(centres.component(vector::Y) - 0.292);
49
50      gamma.write();
51
52      Info<< "End\n" << endl;
53
54      return(0);
55  }
56
57
58  // ************************************************************************* //

The code is compiled using the wmake compilation script described in detail in subsection 3.2.2. In this case, move to the setGammaDamBreak directory and run wmake by typing

cd $OpenFOAM_RUN/tutorials/interFoam/damBreak/setGammaDamBreak
wmake

It is possible that after compiling their first personalised application, the user may have to update the $PATH by resourcing their $HOME/.bashrc (or $HOME/.cshrc) file, depending on whether they are running bash/ksh or tcsh/csh. More specifically: if the user is running bash or ksh, they would need to type at the command line: . $HOME/.bashrc; if running tcsh or csh, they would need to type source $HOME/.cshrc followed by rehash.

The user should then run setGammaDamBreak on the damBreak case:

setGammaDamBreak $OpenFOAM_RUN/tutorials/interFoam damBreak

At this point the user must now select Read Mesh in FoamX. This will read in the changes to the fields made by executing the setGammaDamBreak externally. Check the gamma field in FoamX to see that it is now nonuniform. Using dxFoam, check that the initial gamma field corresponds to the desired distribution as in Figure 2.29.

$Phase fraction, g 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 \special {t4ht=$

Figure 2.29:

Initial conditions for phase fraction gamma.

2.3.4 Fluid properties

Let us examine the case setup from within FoamX. The transportProperties dictionary contains the material properties for each fluid, separated into two subdictionaries phase1 and phase2. The transport model for each phase is selected by the incompressibleTransportModel keyword. The user should select Newtonian in which case the kinematic viscosity is single valued and specified under the keyword nu. The viscosity parameters for the other models, e.g. CrossPowerLaw , are specified within subdictionaries with the generic name <model>Coeffs, i.e. CrossPowerLawCoeffs in this example. The density is specified under the keyword rho.

The surface tension between the two phases is specified under the keyword sigma. The values used in this tutorial are listed in Table 2.3.

phase1 properties

Kinematic viscosity	$m2 s- 1 \special {t4ht=$	nu	$1.0 × 10 -6 \special {t4ht=$
Density	$-3 kg m \special {t4ht=$	rho	$3 1.0 × 10 \special {t4ht=$

phase2 properties

Kinematic viscosity	$m2 s- 1 \special {t4ht=$	nu	$1.48 × 10-5 \special {t4ht=$
Density	$-3 kg m \special {t4ht=$	rho	1.0

Properties of both phases

Surface tension	$N m - 1 \special {t4ht=$	sigma	0.07

Table 2.3:

Fluid properties for the damBreak tutorial

The environmentalProperties dictionary specifies the gravity acceleration vector which should be set to $-2 (0,9.81,0) m s \special {t4ht=$ for this tutorial.

2.3.5 Time step control

Time step control is an important issue in free surface tracking since the surface-tracking algorithm is considerably more sensitive to the Courant number $Co \special {t4ht=$ than in standard fluid flow calculations. Ideally, we should not exceed an upper limit $Co ~~ 0.2 \special {t4ht=$ in the region of the interface. In some cases, where the propagation velocity is easy to predict, the user should specify a fixed time-step to satisfy the $Co \special {t4ht=$ criterion. For more complex cases, this is considerably more difficult. interFoam therefore offers automatic adjustment of the time step as standard in the controlDict. The user should specify adjustTimeStep to be on and the the maximum $Co \special {t4ht=$ , maxCo to be 0.2. The upper limit on time step maxDeltaT can be set to a value that will not be exceeded in this simulation, e.g. 1.0.

By using automatic time step control, the steps themselves are never rounded to a convenient value. Consequently if we request that OpenFOAM saves results at a fixed number of time step intervals, the times at which results are saved are somewhat arbitrary. However even with automatic time step adjustment, OpenFOAM allows the user to specify that results are written at fixed times; in this case OpenFOAM forces the automatic time stepping procedure to adjust time steps so that it ‘hits’ on the exact times specified for write output. The user selects this with the adjustableRunTime option for writeControl in the controlDict dictionary. The controlDict dictionary entries should be:

31  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32
33  applicationClass interFoam;
34
35  startFrom       startTime;
36
37  startTime       0;
38
39  stopAt          endTime;
40
41  endTime         1;
42
43  deltaT          0.001;
44
45  writeControl    adjustableRunTime;
46
47  writeInterval   0.05;
48
49  cycleWrite      0;
50
51  writeFormat     ascii;
52
53  writePrecision  6;
54
55  writeCompression uncompressed;
56
57  timeFormat      general;
58
59  timePrecision   6;
60
61  runTimeModifiable yes;
62
63  adjustTimeStep  yes;
64
65  maxCo           0.5;
66
67  maxDeltaT       1;
68
69
70  // ************************************************************************* //

2.3.6 Discretisation schemes

The free surface treatment in OpenFOAM does not account for the effects of turbulence. This is a consequence of the fact that the Reynolds averaged approach to turbulence modelling does not match the notion of an infinitesimally thin interface between air and water. As a consequence, all free surface simulations can be viewed as a direct numerical simulation (DNS) of fluid flow. DNS is associated with certain requirements on the mesh size, far beyond the mesh resolution of our test case. We shall therefore stabilise the calculation by using upwind differencing on the convection term in the momentum equation. The numerical diffusion introduced by the upwind scheme will stabilise the calculation. The upwind setting is made in the divSchemes subdictionary of the fvSchemes dictionary; the div(rho*phi,U) keyword, corresponding to the $\~/ •(rfU) \special {t4ht=$ term in the momentum equation, should read Gauss upwind.

The convection terms in the two equations for phase $g \special {t4ht=$ require some attention. We must guarantee boundedness while maintaining reasonable accuracy and therefore use the Gamma scheme for bounded scalars Gamma01 for the interpolation in the convection terms. The Gamma schemes require a coefficient $f \special {t4ht=$ as described in subsection 4.4.1. Here, we opt for best accuracy with $f = 0.2 \special {t4ht=$ on the $\~/ •(fg) \special {t4ht=$ term, represented by the div(phi,gamma) keyword and for the $• \~/ (frbg) \special {t4ht=$ term, represented by the div(phirb,gamma) keyword.

The other discretised terms use commonly employed schemes so that the fvSchemes dictionary entries should therefore be:

31  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32
33  timeScheme      EulerImplicit;
34
35  gradSchemes
36  {
37      default         Gauss linear;
38      grad(U)         Gauss linear;
39      grad(gamma)     Gauss linear;
40  }
41
42  divSchemes
43  {
44      div(rho*phi,U)  Gauss upwind;
45      div(phi,gamma)  Gauss Gamma01 0.2;
46      div(phirb,gamma) Gauss Gamma01 1;
47  }
48
49  laplacianSchemes
50  {
51      default         Gauss linear corrected;
52  }
53
54  interpolationSchemes
55  {
56      default         linear;
57  }
58
59  snGradSchemes
60  {
61      default         corrected;
62  }
63
64  fluxRequired
65  {
66      default         no;
67      pd;
68      pcorr;
69      gamma;
70  }
71
72
73  // ************************************************************************* //

2.3.7 Linear-solver control

In the fvSolution, the PISO subdictionary contains elements that are specific to interFoam. There are the usual correctors to the momentum equation but also correctors to a PISO loop around the $g \special {t4ht=$ phase equation. Of particular interest are the nGammaSubCycles and cGamma keywords. nGammaSubCycles represents the number of sub-cycles within the $g \special {t4ht=$ equation; sub-cycles are additional solutions to an equation within a given time step. It is used to enable the solution to be stable without reducing the time step and vastly increasing the solution time. Here we specify 4 sub-cycles, which means that the $g \special {t4ht=$ equation is solved in $4× \special {t4ht=$ quarter length time steps within each actual time step.

The cGamma keyword is a factor that controls the compression of the interface where: 0 corresponds to no compression; 1 corresponds to conservative compression; and, anything larger than 1, relates to enhanced compression of the interface. From experience, we recommend a value of 1.0 for 2-dimensional problems and 1.5 for 3-dimensions. However, we can set higher values to encourage sharper phase boundaries so long as the solution can maintain stability as in this tutorial where 1.3 is used. The code is now ready to run.

2.3.8 Running the code

Running of the code has been described in detail in previous tutorials. The job can be launched from FoamX or manually. Try the following:

cd $OpenFOAM_RUN/tutorials/interFoam
interFoam . damBreak | tee log

The code will now be run interactively, with a copy of output stored in the log file.

2.3.9 Post-processing

Post-processing of the results can now be done in the usual way. The user can monitor the development of the phase fraction gamma in time; Figure 2.30. If a comparison with experimental data is required, a series of pictures for this case can be found in [?].

$Phase fraction, g 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 \special {t4ht=$

(a) At $t = 0.25 s \special {t4ht=$ .

$Phase fraction, g 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 \special {t4ht=$

(b) At $t = 0.35 s \special {t4ht=$ .

Figure 2.30:

Snapshots of phase $g \special {t4ht=$ .

2.3.10 Running in parallel

The results from the previous example are generated using a fairly coarse mesh. We now wish to increase the mesh resolution and re-run the case. The new case will typically take a few hours to run with a single processor so, should the user have access to multiple processors, we can demonstrate the parallel processing capability of OpenFOAM.

The user should first make a copy of the damBreak case, e.g. using the Clone Case function in the FoamX case browser. The new case should be named damBreakFine. Open the new case and change the blocks description in the blockMeshDict dictionary to

    blocks
    (
        hex (0 1 5 4 12 13 17 16) (46 10 1) simpleGrading (1 1 1)
        hex (2 3 7 6 14 15 19 18) (40 10 1) simpleGrading (1 1 1)
        hex (4 5 9 8 16 17 21 20) (46 76 1) simpleGrading (1 2 1)
        hex (5 6 10 9 17 18 22 21) (4 76 1) simpleGrading (1 2 1)
        hex (6 7 11 10 18 19 23 22) (40 76 1) simpleGrading (1 2 1)
    );

Here, the entry is presented as printed from the blockMeshDict file; in short the user must change the mesh densities, e.g. the 46 10 1 entry, and some of the mesh grading entries to 1 2 1. Once the dictionary is correct, generate the mesh.

As the mesh has now changed from the damBreak example, the user must re-initialise the phase field gamma in the 0 time directory since it contains a number of elements that is inconsistent with the new mesh. Note that there is no need to change the U and p fields since they are specified as uniform which is independent of the number of elements in the field. We wish to initialise the field with a sharp interface, i.e. it elements would have $g = 1 \special {t4ht=$ or $g = 0 \special {t4ht=$ . Updating the field with mapFields may produce interpolated values $0 < g < 1 \special {t4ht=$ at the interface, so it is better to rerun the setGammaDamBreak utility by:

setGammaDamBreak $OpenFOAM_RUN/tutorials/interFoam damBreakFine

As in the damBreak case, if the user is using FoamX, they should select Read Mesh which updates the case server to the changes to gamma caused by running the setGammaDamBreak externally of FoamX.

The method of parallel computing used by OpenFOAM is known as domain decomposition, in which the geometry and associated fields are broken into pieces and allocated to separate processors for solution. The first step required to run a parallel case is therefore to decompose the domain using the decomposePar utility which can be selected from the parallelProcessing menu of Foam Utilities in FoamX. The user must edit the dictionary associated with decomposePar named decomposeParDict which is located in the system directory of the case. The first entry is numberOfSubdomains which specifies the number of subdomains into which the case will be decomposed, usually corresponding to the number of processors available for the case.

In this tutorial, the method of decomposition should be simple and the corresponding simpleCoeffs should be edited according to the following criteria. The domain is split into pieces, or subdomains, in the $x \special {t4ht=$ , $y \special {t4ht=$ and $z \special {t4ht=$ directions, the number of subdomains in each direction being given by the vector $n \special {t4ht=$ . As this geometry is 2 dimensional, the 3rd direction, $z \special {t4ht=$ , cannot be split, hence $n z \special {t4ht=$ must equal 1. The $n x \special {t4ht=$ and $n y \special {t4ht=$ components of $n \special {t4ht=$ split the domain in the $x \special {t4ht=$ and $y \special {t4ht=$ directions and must be specified so that the number of subdomains specified by $nx \special {t4ht=$ and $ny \special {t4ht=$ equals the specified numberOfSubdomains, i.e. $nxny = \special {t4ht=$ numberOfSubdomains. It is beneficial to keep the number of cell faces adjoining the subdomains to a minimum so, for a square geometry, it is best to keep the split between the $x \special {t4ht=$ and $y \special {t4ht=$ directions should be fairly even. The delta keyword should be set to 0.001.

For example, let us assume we wish to run on 4 processors. We would set numberOfSubdomains to 4 and $n = (2,2, 1) \special {t4ht=$ . We close the decomposeParDict and run decomposePar. The screen messages of decomposePar can be monitored and show that the decomposition is distributed fairly even between the processors.

Parallel runs can currently only be executed from the command line. $2 \special {t4ht=$ The user should consult section 3.4 for details of how to run a case in parallel; in this tutorial we merely present an example of running in parallel. We use the LAM implementation of the standard message-passing interface (MPI) which requires the user to create a file listing the host names of the machines on which the case is to be run. Let us name this file machines and place it in the system directory of the damBreakFine case. The machines file contains the names of the machines listed one machine per line. As a test, the user can create a machines file to run on a single node, the local host only, by typing

hostname > $OpenFOAM_TUTORIALS/interFoam/damBreakFine/system/machines

The user should then start LAM by typing

lamboot -v $OpenFOAM_TUTORIALS/interFoam/damBreakFine/system/machines

If LAM starts successfully, the following will be written to screen:

    LAM 7.0.6 - Indiana University

    n-1<PID> ssi:boot:base:linear: booting n0 (<hostName>)
    n-1<PID> ssi:boot:base:linear: finished

The user can terminate LAM at any time by typing

lamhalt -d

If the user has access to more licensed nodes, he/she can add them to the machines file as described in subsubsection 3.4.2.1 and restart LAM. Once LAM is running to the user’s satisfaction, interFoam can be run in parallel by typing:

mpirun -np 4 interFoam $OpenFOAM_RUN/tutorials/interFoam damBreakFine
-parallel < /dev/null >& log &

The case should run in the background and the user can follow its progress by monitoring the log file as usual.

2.3.11 Post-processing a case run in parallel

Once the case has completed running, the decomposed fields and mesh must be reassembled for post-processing using the reconstructPar utility. Simply select the utility from the parallelProcessing menu of Foam Utilities in FoamX or run it from the command line. The results from the fine mesh are shown in Figure 2.31. The user can see that the resolution of interface has improved significantly compared to the coarse mesh.

$Phase fraction, g 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 \special {t4ht=$

(a) At $t = 0.25 s \special {t4ht=$ .

$Phase fraction, g 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0 \special {t4ht=$

(b) At $t = 0.35 s \special {t4ht=$ .

Figure 2.31:

Snapshots of phase $g \special {t4ht=$ with refined mesh.

The user may also post-process a segment of the decomposed domain individually. In the user’s $HOME/.foam2.3/controlDict directory (not the case controlDict directory), add the full path of the damBreakFine case directory to the caseRoots, e.g. the entry:

"$OpenFOAM_RUN/tutorials/interFoam/damBreakFine"

$\special {t4ht=$

Figure 2.32:

Mesh of processor 2 in parallel processed case.

The user should restart dxFoam if it currently running. The damBreakFine directory should be selected from the Root menu of the Read Data window. The user will then be able to select processor0 . . . processor3 from the Case menu. The relevant segments of the domain and fields can be displayed for the given processor. Figure 2.32 shows the mesh from processor 2 following the decomposition of the domain using the simple method.

[next] [prev] [prev-tail] [front] [up]