OpenFOAM v9 User Guide

[version 11][version 10][version 9][version 8][version 7][version 6]

[chapter contents][previous topic][next topic][index]

2.3 Breaking of a dam

In this tutorial we shall solve a problem of simpliﬁed dam break in 2 dimensions using the interFoam.The feature of the problem is a transient ﬂow of two ﬂuids separated by a sharp interface, or free surface. The two-phase algorithm in interFoam is based on the volume of ﬂuid (VOF) method in which a specie transport equation is used to determine the relative volume fraction of the two phases, or phase fraction $\text{[math]}$ , 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 ﬁeld. Since the phase fraction can have any value between 0 and 1, the interface is never sharply deﬁned, 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 $\text{[math]}$ , 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 ﬂow structure, including several captured pockets of air. The geometry and the initial setup is shown in Figure 2.21 .

$0.584 m water column 0.584 m 0.292m 0.048 m 0.1461 m 0.1459 m 0.024 m \relax \special {t4ht=$

Figure 2.21: Geometry of the dam break.

2.3.1 Mesh generation

The user should go to their run directory and copy the damBreak case from the $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak directory, i.e.

run
cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/damBreak/damBreak .

Go into the damBreak case directory and generate the mesh running blockMesh as described previously. The damBreak mesh consist of 5 blocks; the blockMeshDict entries are given below.

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

2.3.2 Boundary conditions

The user can examine the boundary geometry generated by blockMesh by viewing the boundary ﬁle in the constant/polyMesh directory. The ﬁle 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 speciﬁed. 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 generic patch, the wall type contains no geometric or topological information about the mesh and only diﬀers from the plain patch in that it identiﬁes the patch as a wall, should an application need to know, e.g. to apply special wall surface modelling. For example, the interFoam solver includes modelling of surface tension and can include wall adhesion at the contact point between the interface and wall surface. Wall adhesion models can be applied through a special boundary condition on the alpha ( $\text{[math]}$ ) ﬁeld, e.g. the constantAlphaContactAngle boundary condition, which requires the user to specify a static contact angle, theta0.

In this tutorial we would like to ignore surface tension eﬀects between the wall and interface. We can do this by setting the static contact angle, $\text{[math]}$ . However, rather than using the constantAlphaContactAngle boundary condition, the simpler zeroGradient can be applied to alpha on the walls.

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

totalPressure which is a ﬁxedValue condition calculated from speciﬁed total pressure p0 and local velocity U;
pressureInletOutletVelocity, which applies zeroGradient on all components, except where there is inﬂow, in which case a ﬁxedValue condition is applied to the tangential component;
inletOutlet, which is a zeroGradient condition when ﬂow outwards, ﬁxedValue when ﬂow is inwards.

At all wall boundaries, the ﬁxedFluxPressure boundary condition is applied to the pressure ﬁeld, which adjusts the pressure gradient so that the boundary ﬂux matches the velocity boundary condition for solvers that include body forces such as gravity and surface tension.

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

2.3.3 Setting initial ﬁeld

Unlike the previous cases, we shall now specify a non-uniform initial condition for the phase fraction $\text{[math]}$ where

{ 1 for the water phase αwater = 0 for the air phase \relax \special {t4ht=

(2.15)

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

16
17defaultFieldValues
18(
19    volScalarFieldValue alpha.water 0
20);
21
22regions
23(
24    boxToCell
25    {
26        box (0 0 -1) (0.1461 0.292 1);
27        fieldValues
28        (
29            volScalarFieldValue alpha.water 1
30        );
31    }
32);
33
34
35// ************************************************************************* //

The defaultFieldValues sets the default value of the ﬁelds, i.e. the value the ﬁeld takes unless speciﬁed otherwise in the regions sub-dictionary. That sub-dictionary contains a list of subdictionaries containing fieldValues that override the defaults in a speciﬁed region. The region 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 deﬁne the set of cells of the water region. The phase fraction $\text{[math]}$ is deﬁned as 1 in this region.

The setFields utility reads ﬁelds from ﬁle and, after re-calculating those ﬁelds, will write them back to ﬁle. In the damBreak tutorial, the alpha.water ﬁeld is initially stored as a backup named alpha.water.orig. A ﬁeld ﬁle with the .orig extension is read in when the actual ﬁle does not exist, so setFields will read alpha.water.orig but write the resulting output to alpha.water (or alpha.water.gz if compression is switched on). This way the original ﬁle is not overwritten, so can be reused.

The user should therefore execute setFields like any other utility by:

setFields

Using paraFoam, check that the initial alpha.water ﬁeld corresponds to the desired distribution as in Figure 2.22.

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

Figure 2.22: Initial conditions for phase fraction alpha.water.

2.3.4 Fluid properties

Let us examine the transportProperties ﬁle in the constant directory. The dictionary ﬁrst contains the names of each ﬂuid phase in the phases list, here water and air. The material properties for each ﬂuid are then separated into two dictionaries water and air. 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 speciﬁed under the keyword nu. The viscosity parameters for other models, e.g.CrossPowerLaw , would otherwise be speciﬁed as described in section 7.3 . The density is speciﬁed under the keyword rho.

The surface tension between the two phases is speciﬁed by the keyword sigma . The values used in this tutorial are listed in Table 2.3 .

water properties

Kinematic viscosity	$\text{[math]}$	nu	$\text{[math]}$
Density	$\text{[math]}$	rho	$\text{[math]}$

air properties

Kinematic viscosity	$\text{[math]}$	nu	$\text{[math]}$
Density	$\text{[math]}$	rho	1.0

Properties of both phases

Surface tension	$\text{[math]}$	sigma	0.07

Table 2.3: Fluid properties for the damBreak tutorial

Gravitational acceleration is uniform across the domain and is speciﬁed in a ﬁle named g in the constant directory. Unlike a normal ﬁeld ﬁle, e.g. U and p, g is a uniformDimensionedVectorField and so simply contains a set of dimensions and a value that represents $\text{[math]}$ for this tutorial:

16
17dimensions [0 1 -2 0 0 0 0];
18value (0 -9.81 0);
19
20
21// ************************************************************************* //

2.3.5 Turbulence modelling

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

16
17simulationType laminar;
18
19
20// ************************************************************************* //

2.3.6 Time step control

Time step control is an important issue in transient simulation and the surface-tracking algorithm in interface capturing solvers. The Courant number $\text{[math]}$ needs to be limited depending on the choice of algorithm: with the “explicit” MULES algorithm, an upper limit of $\text{[math]}$ for stability is typical in the region of the interface; but with “semi-implicit” MULES, speciﬁed by the MULESCorr keyword in the fvSolution ﬁle, there is really no upper limit in $\text{[math]}$ for stability, but instead the level is determined by requirements of temporal accuracy.

In general it is diﬃcult to specify a ﬁxed time-step to satisfy the $\text{[math]}$ criterion, so interFoam oﬀers automatic adjustment of the time step as standard in the controlDict. The user should specify adjustTimeStep to be on and the the maximum $\text{[math]}$ for the phase ﬁelds, maxAlphaCo , and other ﬁelds, maxCo , to be 1.0. 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 ﬁxed 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 ﬁxed times; in this case OpenFOAM forces the automatic time stepping procedure to adjust time steps so that it ‘hits’ on the exact times speciﬁed for write output. The user selects this with the adjustableRunTime option for writeControl in the controlDict dictionary. The controlDict dictionary entries should be:

16
17application     interFoam;
18
19startFrom       startTime;
20
21startTime       0;
22
23stopAt          endTime;
24
25endTime         1;
26
27deltaT          0.001;
28
29writeControl    adjustableRunTime;
30
31writeInterval   0.05;
32
33purgeWrite      0;
34
35writeFormat     binary;
36
37writePrecision  6;
38
39writeCompression off;
40
41timeFormat      general;
42
43timePrecision   6;
44
45runTimeModifiable yes;
46
47adjustTimeStep  yes;
48
49maxCo           1;
50maxAlphaCo      1;
51
52maxDeltaT       1;
53
54
55// ************************************************************************* //

2.3.7 Discretisation schemes

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

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 ( $\text{[math]}$ ), denoted by the div(rhoPhi,U) keyword, uses Gauss linearUpwind grad(U) to produce good accuracy. Here, we have opted for best stability with $\text{[math]}$ . The $\text{[math]}$ term, represented by the div(phi,alpha) keyword uses a specialist interfaceCompression scheme where the speciﬁed coeﬃcient 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 adopt a value of 1.0 which is employed in this example.

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

16
17ddtSchemes
18{
19    default         Euler;
20}
21
22gradSchemes
23{
24    default         Gauss linear;
25}
26
27divSchemes
28{
29    div(rhoPhi,U)  Gauss linearUpwind grad(U);
30    div(phi,alpha)  Gauss interfaceCompression vanLeer 1;
31    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
32}
33
34laplacianSchemes
35{
36    default         Gauss linear corrected;
37}
38
39interpolationSchemes
40{
41    default         linear;
42}
43
44snGradSchemes
45{
46    default         corrected;
47}
48
49
50// ************************************************************************* //

2.3.8 Linear-solver control

In the fvSolution ﬁle, the alpha.water sub-dictionary in solvers contains elements that are speciﬁc to interFoam. Of particular interest are the nAlphaSubCycles and cAlpha keywords. nAlphaSubCycles represents the number of sub-cycles within the $\text{[math]}$ 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 $\text{[math]}$ equation is solved in $\text{[math]}$ half length time steps within each actual time step.

2.3.9 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 ﬁles:

cd $FOAM_RUN/damBreak
interFoam | tee log

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

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

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

(a) At $\text{[math]}$ .

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

(b) At $\text{[math]}$ .

Figure 2.23: Snapshots of phase $\text{[math]}$ .

2.3.11 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 ﬁrst clone the damBreak case, e.g. by

run
foamCloneCase damBreak damBreakFine

Enter the new case directory 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 ﬁle; 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 by running blockMesh.

As the mesh has now changed from the damBreak example, the user must re-initialise the phase ﬁeld 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 ﬁelds since they are speciﬁed as uniform which is independent of the number of elements in the ﬁeld. We wish to initialise the ﬁeld with a sharp interface, i.e. it elements would have $\text{[math]}$ or $\text{[math]}$ . Updating the ﬁeld with mapFields may produce interpolated values $\text{[math]}$ at the interface, so it is better to rerun the setFields utility.

The mesh size is now inconsistent with the number of elements in the alpha.water ﬁle in the 0 directory, so the user must delete that ﬁle so that the original alpha.water.orig ﬁle is used instead.

rm 0/alpha.water
setFields

The method of parallel computing used by OpenFOAM is known as domain decomposition, in which the geometry and associated ﬁelds are broken into pieces and allocated to separate processors for solution. The ﬁrst 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 speciﬁc utility, i.e. in $FOAM_UTILITIES/parallelProcessing/decomposePar for this case.

The ﬁrst entry is numberOfSubdomains which speciﬁes 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 $\text{[math]}$ , $\text{[math]}$ and $\text{[math]}$ directions, the number of subdomains in each direction being given by the vector $\text{[math]}$ . As this geometry is 2 dimensional, the 3rd direction, $\text{[math]}$ , cannot be split, hence $\text{[math]}$ must equal 1. The $\text{[math]}$ and $\text{[math]}$ components of $\text{[math]}$ split the domain in the $\text{[math]}$ and $\text{[math]}$ directions and must be speciﬁed so that the number of subdomains speciﬁed by $\text{[math]}$ and $\text{[math]}$ equals the speciﬁed numberOfSubdomains, i.e. $\text{[math]}$ numberOfSubdomains. It is beneﬁcial 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 $\text{[math]}$ and $\text{[math]}$ 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 $\text{[math]}$ . The user should run decomposePar with:

decomposePar

The terminal output shows that the decomposition is distributed fairly even between the processors.

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

mpirun -np 4 interFoam -parallel > log &

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

$\relax \special {t4ht=$

Figure 2.24: Mesh of processor 2 in parallel processed case.

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

(a) At $\text{[math]}$ .

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

(b) At $\text{[math]}$ .

Figure 2.25: Snapshots of phase $\text{[math]}$ with reﬁned mesh.

2.3.12 Post-processing a case run in parallel

Once the case has completed running, the decomposed ﬁelds and mesh can be reassembled for post-processing using the reconstructPar utility. Simply execute it from the command line. The results from the ﬁne mesh are shown in Figure 2.25 . The user can see that the resolution of interface has improved signiﬁcantly compared to the coarse mesh.

The user may also post-process an individual region 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

paraFoam -case processor1

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

OpenFOAM v9 User Guide - 2.3 Breaking of a dam