OpenFOAM v12 User Guide

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

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

2.2 Breaking of a dam

In this example we shall solve a problem of a simpliﬁed dam break in 2 dimensions using the incompressibleVoF modular solver. 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 incompressibleVoF is based on the volume of ﬂuid (VoF) method in which a phase 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 phases 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 precisely 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.15 .

$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.15: Geometry of the dam break.

2.2.1 Mesh generation

The user should go to their run directory and copy the damBreakLaminar case from the $FOAM_TUTORIALS/incompressibleVoF directory, i.e.

run
cp -r $FOAM_TUTORIALS/incompressibleVoF/damBreakLaminar .

Go into the damBreakLaminar case directory and generate the mesh running blockMesh as described previously. The damBreakLaminar mesh consist of ﬁve 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

55defaultPatch

56{

57    type empty;

58}

59

60boundary

61(

62    leftWall

63    {

64        type wall;

65        faces

66        (

67            (0 12 16 4)

68            (4 16 20 8)

69        );

70    }

71    rightWall

72    {

73        type wall;

74        faces

75        (

76            (7 19 15 3)

77            (11 23 19 7)

78        );

79    }

80    lowerWall

81    {

82        type wall;

83        faces

84        (

85            (0 1 13 12)

86            (1 5 17 13)

87            (5 6 18 17)

88            (2 14 18 6)

89            (2 3 15 14)

90        );

91    }

92    atmosphere

93    {

94        type patch;

95        faces

96        (

97            (8 20 21 9)

98            (9 21 22 10)

99            (10 22 23 11)

100        );

101    }

102);

103

104

105// ************************************************************************* //

The mesh is written into a set of ﬁles in a polyMesh directory in the constant directory. The user can list the contents of the directory to reveal the set of ﬁles/

ls constant/polyMesh

The ﬁles include: points, a list of the cell vertices; faces, a list of the cell faces; owner and neighbour, containing the indices of cells connected to a given face; boundary, a description of the boundary patches.

2.2.2 Boundary conditions

The boundary ﬁle can be read and understood by the user. The user should take a look at its contents, either by opening it in a ﬁle editor or printing out in the terminal window using the cat utility.

cat constant/polyMesh/boundary

The ﬁle contains a list of ﬁve boundary patches: leftWall, rightWall, lowerWall, atmosphere and defaultFaces. The user should notice the type of the patches. Firstly, the atmosphere is a standard patch, i.e. has no special attributes, merely an entity on which boundary conditions can be speciﬁed. Then, the defaultFaces patch is formed of block faces that are omitted from the boundary sub-dictionary in the blockMeshDict ﬁle. Those block faces form a patch whose properties are speciﬁed in a defaultPatch sub-dictionary in the blockMeshDict ﬁle. In this case, the default type is set to 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. This is required by some modelling, e.g. turbulent wall functions and some turbulence models which include the distance to the nearest wall in their calculations.

With VoF speciﬁcally, surface tension models 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 alphaContactAngle boundary condition, which requires the user to specify a static contact angle, theta0.

This example ignores surface tension eﬀects between the wall and interface. This can be done can do this by setting the static contact angle, $\text{[math]}$ , but a simpler approach is to apply the zeroGradient condition 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:

prghTotalPressure, applied to the pressure ﬁeld, minus the hydrostatic component, $\text{[math]}$ , given by Equation 6.4 ;
pressureInletOutletVelocity, applied to velocity $\text{[math]}$ , which sets zeroGradient on all components of $\text{[math]}$ , except where there is inﬂow, in which case a ﬁxedValue condition is applied to the tangential component;
inletOutlet applied to other ﬁelds, 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.2.3 Phases

The ﬂuid phases are speciﬁed in the phaseProperties ﬁle in the constant directory as follows:

16

17phases          (water air);

18

19sigma           0.07;

20

21

22// ************************************************************************* //

It lists two phases, water and air. Equations for phase fraction are solved for the phases in the list, except the last phase listed, i.e. air in this case. Since there are only two phases, only one phase fraction equation is solved in this case, for the water phase fraction $\text{[math]}$ , speciﬁed in the ﬁle alpha.water in the 0 directory.

The phaseProperties ﬁle also contains an entry for the surface tension between the two phases, speciﬁed by the keyword sigma in units $\text{[math]}$ .

2.2.4 Setting initial ﬁelds

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.5)

This is 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    volScalarFieldValue tracer.water 0

21    volScalarFieldValue tracer.air 0

22);

23

24regions

25(

26    boxToCell

27    {

28        box (0 0 -1) (0.1461 0.292 1);

29        fieldValues

30        (

31            volScalarFieldValue alpha.water 1

32        );

33    }

34    boxToCell

35    {

36        box (0 0 -1) (0.1461 0.146 1);

37        fieldValues

38        (

39            volScalarFieldValue tracer.water 1

40        );

41    }

42    boxToCell

43    {

44        box (0 0.292 -1) (0.1461 0.438 1);

45        fieldValues

46        (

47            volScalarFieldValue tracer.air 1

48        );

49    }

50);

51

52

53// ************************************************************************* //

The defaultFieldValues sets the default value of the ﬁelds, i.e. the value the ﬁeld takes unless speciﬁed otherwise in the regions list. That list contains sub-dictionaries with fieldValues that override the defaults in a speciﬁed region. Regions are deﬁned by functions that deﬁne a set of points, faces, or cells. Here boxToCell creates a box deﬁned by minimum and maximum bounds that deﬁnes the set of cells of the water region. The phase fraction $\text{[math]}$ is speciﬁed 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 damBreakLaminar case, the alpha.water ﬁeld is initially stored in its original form with the name 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 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.16.

$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.16: Initial conditions for phase fraction alpha.water.

2.2.5 Fluid properties

The physical properties for the air and water phases are speciﬁed in physicalProperties.air and physicalProperties.water ﬁles, respectively, in the constant directory. Physical properties describe characteristics of the ﬂuid in a absence of ﬂow. Each ﬁle speciﬁes the viscosity model through the viscosityModel keyword, which is set to constant to indicate the value is unchanging. The viscosity is then speciﬁed by the nu keyword in units $\text{[math]}$ . The density of each ﬂuid is also speciﬁed by the keyword rho in units $\text{[math]}$ . The physicalProperties.air ﬁle is shown below as an example:

16

17viscosityModel  constant;

18

19nu              1.48e-05;

20

21rho             1;

22

23

24// ************************************************************************* //

If the viscosity does change according to the ﬂow, e.g. as in non-Newtonian or visco-elastic ﬂuids, then those models are speciﬁed through the momentumProperties ﬁle, as described in section 8.3 .

2.2.6 Gravity

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

16

17dimensions      [0 1 -2 0 0 0 0];

18value           (0 -9.81 0);

19

20

21// ************************************************************************* //

2.2.7 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.2.8 Time step control

The simulation of the damBreakLaminar case is fully transient so the time step requires attention. The Courant number $\text{[math]}$ is an important consideration relating to time step. It is dimensionless parameter which can be deﬁned for each cell as:

Co = δt|U-| δx \relax \special {t4ht=

(2.6)

where $\text{[math]}$ is the time step, $\text{[math]}$ is the magnitude of the velocity through that cell and $\text{[math]}$ is the cell size in the direction of the velocity. With explicit solution, stability requires the maximum $\text{[math]}$ at least; stricter limits exist depending on the choice of advection scheme. Implicit solutions do not have the same stability limit of the maximum $\text{[math]}$ , but temporal accuracy becomes more relevant as $\text{[math]}$ increased beyond 1.

Time step control is particularly important with interface-capturing. The incompressibleVoF solver module uses the multidimensional universal limiter for explicit solution (MULES), created by Henry Weller, to maintain boundedness of the phase fraction. $\text{[math]}$ needs to be limited depending on the choice of MULES algorithm. With the original explicit MULES algorithm, an upper limit of $\text{[math]}$ for stability is typically required. However, there is also the semi-implicit version of MULES, speciﬁed by the MULESCorr switch in the fvSolution ﬁle. For semi-implicit MULES, 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 since $\text{[math]}$ is changing from cell to cell during the simulation. Instead, automatic adjustment of the time step is speciﬁed in the controlDict by switching adjustTimeStep to on and specifying the maximum $\text{[math]}$ for the phase ﬁelds, maxAlphaCo , and other ﬁelds, maxCo . In this example, the maxAlphaCo and maxCo are set to 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 with automatic time step adjustment, results can be written at ﬁxed times using the adjustableRunTime option for writeControl in the controlDict dictionary. With this option, the automatic time stepping procedure further adjusts their time steps so that it ‘hits’ on the exact times speciﬁed by the writeInterval, set to 0.05 in this example. The controlDict dictionary entries are shown below.

16

17application     foamRun;

18

19solver          incompressibleVoF;

20

21startFrom       startTime;

22

23startTime       0;

24

25stopAt          endTime;

26

27endTime         1;

28

29deltaT          0.001;

30

31writeControl    adjustableRunTime;

32

33writeInterval   0.05;

34

35purgeWrite      0;

36

37writeFormat     ascii;

38

39writePrecision  6;

40

41writeCompression off;

42

43timeFormat      general;

44

45timePrecision   6;

46

47runTimeModifiable yes;

48

49adjustTimeStep  yes;

50

51maxCo           1;

52maxAlphaCo      1;

53

54maxDeltaT       1;

55

56// ************************************************************************* //

2.2.9 Discretisation schemes

The MULES method, used by the incompressibleVoF modular solver, maintains boundedness of the phase fraction independently of the underlying numerical scheme, mesh structure, etc. The choice of schemes for convection are therefore not restricted to those that are strongly stable or bounded, such as 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.

The $\text{[math]}$ term, represented by the div(phi,alpha) keyword uses a bespoke 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 use a value of 1.0, as in this example.

The other discretised terms use commonly employed schemes so that the fvSchemes dictionary entries is as follows.

16

17ddtSchemes

18{

19    default         Euler;

20}

21

22gradSchemes

23{

24    default         Gauss linear;

25}

26

27divSchemes

28{

29    div(phi,alpha)  Gauss interfaceCompression vanLeer 1;

30

31    div(rhoPhi,U)   Gauss linearUpwind grad(U);

32

33    div(phi,k)      Gauss upwind;

34    div(phi,epsilon) Gauss upwind;

35    div(phi,omega)  Gauss upwind;

36

37    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;

38

39    div(alphaPhi,tracer) Gauss vanLeer;

40}

41

42laplacianSchemes

43{

44    default         Gauss linear uncorrected;

45}

46

47interpolationSchemes

48{

49    default         linear;

50}

51

52snGradSchemes

53{

54    default         uncorrected;

55}

56

57

58// ************************************************************************* //

2.2.10 Linear-solver control

In the fvSolution ﬁle, the sub-dictionary in solvers for alpha.water contains elements that are speciﬁc to the MULES algorithm as shown below.

"alpha.water.*"

{

   nAlphaCorr      2;

   nAlphaSubCycles 1;



   MULESCorr       yes;

   nLimiterIter    5;



   solver          smoothSolver;

   smoother        symGaussSeidel;

   tolerance       1e-8;

   relTol          0;

}

MULES calculates two limiters to keep the phase fraction within the lower and upper bounds of 0 and 1. The limiter calculation is iterative, with the number of iterations speciﬁed by nLimiterIter. As a rule of thumb, nLimiterIter can be set to 3 + $\text{[math]}$ to maintain boundedness.

The semi-implicit version of MULES is activated by the MULESCorr switch. It ﬁrst calculates an implicit, upwind solution before applying MULES as a higher-order correction. The linear solver must be conﬁgured for the implicit, upwind solution, through the solver, smoother, tolerance and relTol parameters.

The nAlphaCorr keyword which controls the number of iterations of the phase fraction equation within a solution step. The iteration is used to overcome nonlinearities in the advection which are present in this case due to the interfaceCompression scheme.

2.2.11 Running the code

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

foamRun | tee log

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

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

$PICT\relax \special {t4ht=$

Figure 2.17: Phase fraction $\text{[math]}$ at $\text{[math]}$ = 0.25 s (left) and 0.50 s (right).

2.2.13 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. Using a ﬁner mesh, we can then demonstrate the parallel processing capability of OpenFOAM.

The user should ﬁrst clone the damBreakLaminar case, e.g. by

run
foamCloneCase damBreakLaminar damBreakLaminarFine

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

The user should then rerun the setFields utility. However, 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 setFields uses the original alpha.water.orig ﬁle.

rm 0/alpha.water
setFields

Parallel computing uses 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. Also, sample dictionaries can be found within the etc directory in the OpenFOAM installation, which can be copied to the case directory by running the foamGet script, e.g. (you do not need to do this):

foamGet decomposeParDict

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.

This example uses the simple method of decomposition. It requires the simpleCoeffs to be conﬁgured 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. For example, let us assume we wish to run on 4 processors. We would set numberOfSubdomains to 4 and n to (2, 2, 1). The user should run decomposePar by the following command.

decomposePar

The terminal output shows that the decomposition is distributed evenly between the sub-domains. The decomposition writes the mesh and ﬁelds of each sub-domain into separate sub-directories named processor<N>, where N is the sub-domain ID, e.g. 0, 1, 2, etc. The user should list the ﬁles in the case directory to conﬁrm that four directories processor0, processor1, processor2 and processor3 exist.

This example presents running in parallel with the openMPI implementation of the standard message-passing interface (MPI). The following command runs on 4 cores of a local multi-processor CPU.

mpirun -np 4 foamRun -parallel

The user can consult section 3.4 for more details of how to run a case in parallel. For example, 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 .

2.2.14 Post-processing a case run in parallel

When the case runs in parallel, the results are written into time directories within the processor<N> sub-directories. The user can conﬁrm this by listing the time directories for the processor0 directory.

ls processor0

It is possible to post-process an individual sub-domain by treating the individual processor directory as a case in its own right. For example, to view the processor1 sub-domain in ParaView, the user can launch paraFoam by running the following command.

paraFoam -case processor1

Figure 2.18 shows the mesh from this sub-domain, following the decomposition of the domain using the simple method.

$\relax \special {t4ht=$

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

Alternatively, the decomposed ﬁelds and mesh can be reassembled for post-processing in serial. The reconstructPar utility takes the ﬁeld ﬁles from time directories for the processor sub-domain and builds equivalent ﬁeld ﬁles for the complete domain. The user should run the utility as follows.

reconstructPar

The ﬁelds are reconstructed are written to solution time directories in the case directory. These ﬁelds can be visualised as normal in ParaView. The results from the ﬁne mesh are shown in Figure 2.19. The user can see that the resolution of interface has improved signiﬁcantly compared to the coarse mesh.

$PICT\relax \special {t4ht=$

Figure 2.19: Phase fraction $\text{[math]}$ at $\text{[math]}$ = 0.25 s (left) and 0.50 s (right).

OpenFOAM v12 User Guide - 2.2 Breaking of a dam