4.1 Breaking of a dam

Tutorial path:

• $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak

This example introduces the following OpenFOAM features for the first time:

• multiphase flow using VOF interface capturing method
• non-uniform initial conditions setup using setFields utility
• running case in parallel
• post-processing a case in parallel

4.1.1 Problem specification

In this tutorial we shall solve a problem of simplified dam break in 2 dimensions using the interFoam.The 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 , 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 , 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 4.1.

---

---

4.1.2 Mesh generation

The user should go to the damBreak case in their $FOAM_RUN/tutorials/multiphase/interFoam/laminar directory. Generate the mesh running blockMesh as described previously. The damBreak mesh consist of 5 blocks; the blockMeshDict entries are given below.

4.1.3 Boundary conditions

The user can examine the boundary geometry generated by blockMesh by viewing the boundary file in the constant/polyMesh directory. The file contains a list of 5 boundary patches: leftWall, rightWall, lowerWall, atmosphere and defaultFaces. The user should notice the type of the patches. The atmosphere is a standard patch, i.e. has no special attributes, merely an entity on which boundary conditions can be specified. The defaultFaces patch is empty since the patch normal is in the direction we will not solve in this 2D case. The leftWall, rightWall and lowerWall patches are each a wall. Like the plain patch, the wall type contains no geometric or topological information about the mesh and only differs from the plain patch in that it identifies the patch as a wall, should an application need to know, e.g. to apply special wall surface modelling.

A good example is that the interFoam solver includes modelling of surface tension at the contact point between the interface and wall surface. The models are applied by specifying the alphaContactAngle boundary condition on the alpha.water () field. With it, the user must specify the following: a static contact angle, theta0 ; leading and trailing edge dynamic contact angles, thetaA and thetaR 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, and the velocity scaling function to 0. However, the simpler option which we shall choose here is to specify a zeroGradient type on alpha.water, rather than use the alphaContactAngle boundary condition.

The top boundary is free to the atmosphere so needs to permit both outflow and inflow according to the internal flow. We therefore use a combination of boundary conditions for pressure and velocity that does this while maintaining stability. They are:

• totalPressure which is a fixedValue condition calculated from specified total pressure p0 and local velocity U;
• pressureInletOutletVelocity, which applies zeroGradient on all components, except where there is inflow, in which case a fixedValue condition is applied to the tangential component;
• inletOutlet, which is a zeroGradient condition when flow outwards, fixedValue when flow is inwards.

At all wall boundaries, the fixedFluxPressure boundary condition is applied to the pressure field, which calculates the normal gradient from the local density gradient.

The defaultFaces patch representing the front and back planes of the 2D problem, is, as usual, an empty type.

4.1.4 Setting initial field

Unlike the previous cases, we shall now specify a non-uniform initial condition for the water phase fraction, , where

(4.1)

This is achieved by running the setFields utility. It requires a setFieldsDict dictionary, located in the system directory, whose entries for this case are shown below.

The defaultFieldValues sets the default value of the fields, i.e. the value the field takes unless specified otherwise in the regions sub-dictionary. That sub-dictionary contains a list of subdictionaries containing fieldValues that override the defaults in a specified region. The region is expressed in terms of a topoSetSource that creates a set of points, cells or faces based on some topological constraint. Here, boxToCell creates a bounding box within a vector minimum and maximum to define the set of cells of the liquid region. The phase fraction is defined as 1 in this region.

The setFields utility reads fields from file and, after re-calculating those fields, will write them back to file. Because the files are then overridden, it is recommended that a backup is made before setFields is executed. In the damBreak tutorial, the alpha.water field is initially stored as a backup only, named alpha.water.orig. Before running setFields, the user first needs to copy alpha.water.orig to alpha.water, e.g. by typing:

The user should then execute setFields as any other utility is executed. Using paraFoam, check that the initial alpha.water field corresponds to the desired distribution as in Figure 4.2.

---

---

4.1.5 Fluid properties

Let us examine the transportProperties file in the constant directory. Its 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 transportModel 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 4.1.

Gravitational acceleration is uniform across the domain and is specified in a file named g in the constant directory. Unlike a normal field file, e.g. U and p, g is a uniformDimensionedVectorField and so simply contains a set of dimensions and a value that represents for this tutorial:

4.1.6 Turbulence modelling

As in the cavity example, the choice of turbulence modelling method is selectable at run-time through the simulationType keyword in turbulenceProperties dictionary. In this example, we wish to run without turbulence modelling so we set laminar:

4.1.7 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 than in standard fluid flow calculations. Ideally, we should not exceed an upper limit 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 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 yes and the maximum for the phase fields, maxAlphaCo, and other fields, maxCo, to be 0.5. 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:

4.1.8 Discretisation schemes

The interFoam solver uses the multidimensional universal limiter for explicit solution (MULES) method, created by OpenCFD, to maintain boundedness of the phase fraction independent of underlying numerical scheme, mesh structure, etc. The choice of schemes for convection are therefore not restricted to those that are strongly stable or bounded, e.g. upwind differencing.

The convection schemes settings are made in the divSchemes sub-dictionary of the fvSchemes dictionary. In this example, the convection term in the momentum equation (), denoted by the div(rhoPhi,U) keyword, uses Gauss linearUpwind grad(U) to produce good accuracy. The term, represented by the div(phi,alpha) keyword uses the vanLeer scheme. The term, represented by the div(phirb,alpha) keyword, can similarly use the vanLeer scheme, but generally produces smoother interfaces using the linear scheme.

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

4.1.9 Linear-solver control

In the fvSolution, the PISO sub-dictionary 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 phase equation. Of particular interest are the nAlphaSubCycles and cAlpha keywords. nAlphaSubCycles represents the number of sub-cycles within the 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 2 sub-cycles, which means that the equation is solved in half length time steps within each actual time step.

The cAlpha 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. We generally recommend a value of 1.0 which is employed in this example.

4.1.10 Running the code

Running of the code has been described in detail in previous tutorials. Try the following, that uses tee, a command that enables output to be written to both standard output and files:

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

4.1.11 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 alpha.water in time, e.g. see Figure 4.3.

---

---

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

Enter the new case directory and change the blocks description in the blockMeshDict dictionary to 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 alpha.water 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_rgh 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 or . Updating the field with mapFields may produce interpolated values at the interface, so it is better to rerun the setFields utility. There is a backup copy of the initial uniform field named 0/alpha.water.org that the user should copy to 0/alpha.water before running setFields:

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. There is a dictionary associated with decomposePar named decomposeParDict which is located in the system directory of the tutorial case; also, like with many utilities, a default dictionary can be found in the directory of the source code of the specific utility, i.e. in $FOAM_UTILITIES/parallelProcessing/decomposePar for this 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 , and directions, the number of subdomains in each direction being given by the vector . As this geometry is 2 dimensional, the 3rd direction, , cannot be split, hence must equal 1. The and components of split the domain in the and directions and must be specified so that the number of subdomains specified by and equals the specified numberOfSubdomains, i.e.  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 and 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 . When running decomposePar, we can see from the screen messages that the decomposition is distributed fairly even between the processors.

The user should consult User Guide section ?? 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 openMPI implementation of the standard message-passing interface (MPI). As a test here, the user can run in parallel on a single node, the local host only, by typing:

The user may run on more nodes over a network by creating a file that lists the host names of the machines on which the case is to be run as described in User Guide section ??. The case should run in the background and the user can follow its progress by monitoring the log file as usual.

---

---

---

---

4.1.13 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 execute it from the command line. The results from the fine mesh are shown in Figure 4.5. The user can see that the resolution of interface has improved significantly compared to the coarse mesh.

The user may also post-process a segment of the decomposed domain individually by simply treating the individual processor directory as a case in its own right. For example if the user starts paraFoam by

then processor1 will appear as a case module in ParaView. Figure 4.4 shows the mesh from processor 1 following the decomposition of the domain using the simple method.