Anders Rynell

CFD with OpenSource Software

OpenFOAM 1.5.dev
Tutorial of interTrackFoam-solver
Peer Reviewed by Anton Berce and Jelena Andric

Department of Applied Mechanics

Division of Fluid Dynamics
Chalmers University of Technology
Göteborg Sweden, 2010

1
2
Contents

1 introduction 1

2 hydrofoil 3

3 domainSettings 5
3.1 blockMesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 7

4 modDamBreak 11
4.1 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2 constant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.2.1 dynamicMeshDict . . . . . . . . . . . . . . . . . . . . 12
4.2.2 faMesh . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.2.3 freeSurfaceProperties . . . . . . . . . . . . . . . . . . 15
4.2.4 polyMesh . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2.5 transportProperties . . . . . . . . . . . . . . . . . . . 20
4.3 system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5 discretisation 23
5.1 blockMesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2 makeFaMesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.3 Automatic mesh motion . . . . . . . . . . . . . . . . . . . . . 24

6 interTrackFoam 27

3
4
Chapter 1

introduction

This tutorial will describe the solver interT rackF oam which is a part of
the OpenFOAM 1.5.dev version. The solver is an incompressible, transient
solver which uses a free surface tracking algorithm. What this really means
will be explained as well as the different choices considering the settings for
a modified case modDamBreak.
First as always when it comes to computational fluid dynamics the domain
has to be discretized into small cells where the properties solved for can be
computed. The equation is solved for each timestep and the computational
mesh is adjusted to the shape of the free surface. This means that the mesh
can not be pre-defined. All the internal points depend strongly on the motion
of the free surface and will move accordingly to the chosen mesh solver. In
order to get reliable results the mesh quality has to be maintained during
the calculations.

1
2
Chapter 2

hydrofoil

The solver interT rackF oam is used when solving the case hydrof oil which
consists of an inclined airfoil situated below a free surface in water. The
pressure distribution around the profile causes waves to propagate on the
water surface. The mesh in turn is adapted to the free surface motion and
will adjust to it. The inlet height is constrained to represent constant inlet
mass flow. The whole domain is shown in Figure 2.1.

Figure 2.1: The whole hydrof oil domain

A closer look at the hydrof oil situated in the refinement mesh area is shown
in Figure 2.2. As can be seen in the figures, hydrof oil consists of thousands
of cells. In order to show the different settings, a different geometry is used,
otherwise the computational time will be unacceptable long. The geometry
from damBreak was chosen since it also includes an obstacle, which will

3
cause a pressure field to propagate, see Figure 2.3. This in turn will cause
the free surface to move.

Figure 2.2: A close-up of the hydrof oil

Figure 2.3: The modDamBreak domain

4
Chapter 3

domainSettings

In order to be able to use the solver interT rackF oam on damBreak some
modifications in the files are necessary. Since the only thing that will change
is the domain, the polyM esh folder from damBreak will be copied into
the hydrof oil folder. Then hydrof oil is simply renamed modDamBreak.
The boundary file is also deleted, since it can cause problem when using
blockM esh. All this is done by typing;

run
cp -r $FOAM_TUTORIALS/interTrackFoam/hydrofoil .
rm -rf hydrofoil/constant/polyMesh
cp -r $FOAM_TUTORIALS/interFoam/damBreak/constant/\
polyMesh hydrofoil/constant/polyMesh
mv hydrofoil modDamBreak
rm -rf modDamBreak/constant/polyMesh/boundary

Now it remains to change the settings for the modified case. Open the
blockMeshDict by typing;

gedit constant/polyMesh/blockMeshDict

3.1 blockMesh
In blockM eshDict file, the patches is changed according to;

patches
(
patch left
(
(0 12 16 4)

5
 (4 16 20 8)
)
patch right
(
(7 19 15 3)
(11 23 19 7)
)
wall bottom
(
(0 1 13 12)
(1 5 17 13)
(5 6 18 17)
(2 14 18 6)
(2 3 15 14)
)
patch freeSurface
(
(8 20 21 9)
(9 21 22 10)
(10 22 23 11)
)
empty frontAndBackPlanes
(
(0 1 5 4)
(4 5 9 8)
(5 6 10 9)
(2 3 7 6)
(6 7 11 10)
(12 13 17 16)
(16 17 21 20)
(17 18 22 21)
(14 15 19 18)
(18 19 23 22)
)
);

mergePatchPairs
(
);

6
The mesh is constructed and the different boundaries are renamed. It is
important that the free surface is named f reeSurf ace otherwise the solver
later on will not be able to track it. No geometrical changes are made but the
left and right sides of the domain are no longer of type ”wall” but ”patch”.
Also the f rontAndBackP lanes is added and is essential for the finite-area
disretization, see f aM eshDef inition.

3.2 boundary Conditions

Also the boundary conditions have to be modified for the velocity, pressure
and displacements. It is important that the names used in blockM eshDict,
i.e. the names for each patch are specified also in the motionU , p and U .
Since the files have been copied from hydrof oil it is important that the
boundary conditions for the hydrof oil is removed and the ”inlet” and ”out-
let” are replaced by ”left” and ”right”. Open the file motionU to change its
corresponding boundary conditions;

gedit 0/motionU

Change the boundaryF ield settings in motionU to

boundaryField
{
bottom
{
type fixedValue;
value uniform (0 0 0);
}

freeSurface
{
type fixedValue;
value uniform (0 0 0);
}

left
{
type slip;
}

right

7
 {
type slip;
}

frontAndBackPlanes
{
type empty;
}
}

Open the file p to change its corresponding boundary conditions

gedit 0/p

Change the boundaryF ield settings in p to

boundaryField
{
bottom
{
type zeroGradient;
}

freeSurface
{
type fixedValue;
value uniform 0.0;
}

right
{
type zeroGradient;
}

left
{
type zeroGradient;
}

frontAndBackPlanes
{
type empty;

8
 }
}

Open the file U to change its corresponding boundary conditions

gedit 0/U

Change the boundaryF ield settings in U to

boundaryField
{
bottom
{
type slip;
}

freeSurface
{
type fixedGradient;
gradient uniform (0 0 0);
}

right
{
type zeroGradient;
}

left
{
type fixedValue;
value uniform (0.8 0 0);
}

frontAndBackPlanes
{
type empty;
}
}

The finite area discretization of the f reeSurf ace also needs some settings,
therefore change the boundary in f aM eshDef inition. f aM eshDef inition
can be found in the f aM esh folder located in the constant directory.

9
boundary
{
left
{
type patch;
ownerPolyPatch freeSurface;
neighbourPolyPatch left;
}

right
{
type patch;
ownerPolyPatch freeSurface;
neighbourPolyPatch right;
}

frontAndBackPlanes
{
type empty;
ownerPolyPatch freeSurface;
neighbourPolyPatch frontAndBackPlanes;
}
}
Open f reeSurf aceP roperties and change the f ixedF reeSurf aceP atches1(inlet)
to f ixedF reeSurf aceP atches1(lef t);. The modified damBreak can now be
used. Type blockM esh, makeF aM esh and interT rackF oam. In order to
save time, modify the Allrun file so it looks like
#!/bin/sh

. $WM_PROJECT_DIR/bin/tools/RunFunctions

application="interTrackFoam"

runApplication blockMesh
runApplication makeFaMesh
runApplication $application
Despite the fact that the domain has been changed, the computational is
still time consuming. The constrain on the f reeSurf ace curvature (no mass
flux) makes the computation sensitive to large displacements.

10
Chapter 4

modDamBreak

All OpenFOAM cases contains at its minimum three directories in order

to run the application. This directories are 0, constant and system. This
is where the settings are made. A primary explanation of these are now
presented.

4.1 0
In the time directory, 0;, the initial conditions are specified for the properties
solved for, namely

motionU p U

where motionU ,p and U are the displacements, pressure and velocity, respec-
tively. The settings for each has already been setup and can be visualized in
Section 3.2. The settings for the hydrofoil case differs from the ones presented
here but the topology looks the same. Information of the various settings
can be found at the [1].
Dependent on the settings which are made in controlDict, more files with
results will be written and placed here. The names for each is based on the
simulated time at which the data is written.

4.2 constant
Constant directory contains the following;

dynamicMeshDict faMesh freeSurfaceProperties polyMesh

transportProperties

11
In short, in the constant directory the mesh settings are made, both the
geometry and the settings which correspond to the mesh and its motion
during the solving process. A further look into these settings is presented
below.

4.2.1 dynamicMeshDict
twoDMotion yes;

solver laplaceFaceDecomposition;

diffusivity patchEnhanced;

distancePatches 1 (freeSurface);

frozenDiffusion yes;

pseudoSolid
{
poissonsRatio 0.3;
nCorrectors 3;
convergenceTolerance 1e-9;
};

The various choices of settings in this files affect the characteristics of the
mesh. If twoDM otion is chosen only two dimensional movements are allowed
and twoDP ointCorrectors[2] is used during the calculations. twoDP oint−
Correctors prevents the mesh from twisting. It remains to determine the mo-
tion of the mesh points. The solver-types available are laplaceF aceDecomposition,
pseudoSolidF aceDecomposition and RBF M otionSolver. The laplaceF aceDecomposition
uses a Laplacian equation to solve the motion of the mesh points emanating
from the motion of the freeSurface. Since the magnitude of the point motion
is not uniformly distributed through the whole domain, the diffusivity type
chosen will help to conserve the mesh quality. One advantageous way is to
confine the largest deformation to the internal part of the mesh. Table 4.1
shows the various diffusivity options.
In the quality-based methods, the diffusion field is a function of the cell
quality measurements while the distance-based methods uses the distance to
the nearest boundary, specified by distanceP atches [3]. The other option
f ile uses a diffusivity field specified in a file. Since the free surface mo-

12
 diffusivity type
distancebased linear
quadratic
exponental
patchEnhanced
qualitybased unif orm
distortionEnergy
def ormationEnergy
other f ile

Table 4.1: The diffusivity options

tion in this case is very smooth and the fact that the mesh is coarse, the
choice of diffusivity type will not affect the solution. An example where
the diffusivity options makes a huge difference can be illustrated in [4].
f rozenDif f usion means that the initial diffusion rate is ”frozen” and does
not change during the calculation process. The main difference between the
pseudoSolidF aceDecomposition and
laplaceF aceDecomposition is that the former also deals with rotation. Since
the pseudo-solid solver requires a greater storage because of the coupled mo-
tion vector components, it is a matter of computational effort. If the mesh
quality is not substantially improved its choice as a solver, can not be jus-
tified. Figure 4.2 illustrates the mesh rotation obtained from solving the
pseudo-solid equation while Figure 4.1 is the laplacian.
The RBF M otionSolver[5] uses a different approach. It uses radial basis
function (therof the abbreviation RBF) to interpolate the motion of the mov-
ing boundary. ”A radial basis function (RBF) is a real-valued function whose
value depends only on the distance from the origin” [6]. A subset of control-
points along the f reeSurf ace are taken, which are used in the interpolation
process. This is an algebraic formulation compared to the partial equations
of laplacian and pseudo-solid type. For more information about its settings,
a deeper look into the solver icoDyM F oam[7] is an alternative.

13
 Figure 4.1: Laplace mesh motion solver

Figure 4.2: Pseudo-solid mesh motion solver

4.2.2 faMesh
The folder f aM esh consists of

boundary.gz faceLabels.gz faMeshDefinition

14
The two former are compressed files which simply describes the f reeSurf ace
and are constructed when executing makeF aM esh. The conditions for the
edges of f reeSurf ace are listed in f aM eshDef inition.
polyMeshPatches 1( freeSurface );

boundary
{
left
{
type patch;
ownerPolyPatch freeSurface;
neighbourPolyPatch left;
}

right
{
type patch;
ownerPolyPatch freeSurface;
neighbourPolyPatch right;
}

frontAndBackPlanes
{
type empty;
ownerPolyPatch freeSurface;
neighbourPolyPatch frontAndBackPlanes;
}
}
Only one polyMeshPatch is chosen, namely the f reeSurf ace. As an example
the first definition means that the condition for the inlet edge is of type patch,
belongs to the f reeSurf ace and is neighbour with the inlet. In the same
way the three other edges are specified.

4.2.3 freeSurfaceProperties
The f reeSurf aceP roperties defines the conditions at the f reeSurf ace.
Some of the settings involve the nature of surfactants, therefore a brief ex-
planation of what surfactants are and its given equations are given.

15
twoFluids no;

normalMotionDir no;

freeSurfaceSmoothing no;

cleanInterface yes;

muFluidA muFluidA [ 1 -1 -1 0 0 0 0 ] 0;

muFluidB muFluidB [ 1 -1 -1 0 0 0 0 ] 1.5e-5;

rhoFluidA rhoFluidA [ 1 -3 0 0 0 0 0 ] 1000.0;

rhoFluidB rhoFluidB [ 1 -3 0 0 0 0 0 ] 1.0;

surfaceTension surfaceTension [ 1 -2 0 0 0 0 0 ] 0.0;

g g [ 0 1 -2 0 0 0 0 ] (0 -9.81 0);

fixedFreeSurfacePatches 1 ( inlet );

surfactantProperties
{
bulkConc bulkConc [ 0 -3 0 0 1 0 0 ] 1.0e-2;

saturatedConc saturatedSurfConc [ 0 -2 0 0 1 0 0 ] 5.0e-6;

adsorptionCoeff adsorptionCoeff [ 0 3 -1 0 -1 0 0 ] 40.0;

desorptionCoeff desorptionCoeff [ 0 -3 0 0 1 0 0 ] 8.93e-2;

bulkDiffusion bulkDiffusion [ 0 2 -1 0 0 0 0 ] 1.0e-9;

diffusion diffusion [ 0 2 -1 0 0 0 0 ] 1.0e-9;

temperature temperature [ 0 0 0 1 0 0 0 ] 293.0;

}

From the top to the bottom, twoF luids means the possibility to simulate
two fluids and their interaction along the f reeSurf ace. In order to do so, a

16
separate mesh has to be created for each fluid and proper boundary, which
bounds the fluids[8] is needed. Also the properties corresponding for each
fluid has to be specified, as in Table 4.2.

muFluidA muFluidA [ 1 -1 -1 0 0 0 0 ] 0
muFluidB muFluidB [ 1 -1 -1 0 0 0 0 ] 1.5e-5
rhoFluidA rhoFluidA [ 1 -3 0 0 0 0 0 ] 1000.0
rhoFluidB rhoFluidB [ 1 -3 0 0 0 0 0 ] 1.0

Table 4.2: Viscosity and density for each fluid

normalM otionDir deals with the mesh motion, and if used, the discretized
cells are free to move, both horizontally and vertically. Figure 4.3 is without
normal motion direction, i.e. motion is only allowed in the vertical direction
while Figure 4.4 shows the case where motion in both direction is allowed.

Figure 4.3: Original mesh motion

f reeSurf aceSmoothing, smooth the f reeSurf ace curvature, taking surface

tension into account. One of the most important settings is the cleanInterf ace
option. If the free surface is not chemically clean, surfactant (surface active
agents) chemicals concentrate on the free surface. Surfactants are compounds
that alter the interfacial tension and are convected and diffused in the bulk
fluid, as well as at the interface. This, in turn, will influence the local surface

17
 Figure 4.4: normalM otionDir is used

tension and significantly influence the behaviour of the system. The con-
centration on the free surface acts as a boundary condition for the volume
transport. The surfactant properties are shown in Table 4.3.

bulkConc Φ/C bulk surfactant concentration

saturatedConc Φ∞ saturated surfactant concentration
absorptionCoeff ka parameter of absorption kinetics
desorptionCoeff βa parameter of absorption kinetics
bulkDiffusion D bulk surfactant diffusion coefficient
diffusion Ds surfactant diffusion coefficient along the surface
Temperature T Temperature

Table 4.3: Surfactant properties

The dynamic condition, which states that the forces acting on the interface
are in equilibrium depends on the surface tension coefficient. This coefficient
is calculated using
Φ
σ = σ0 + RT Φ∞ ln(1 − ) (4.1)
Φ∞
where σ0 is the surface tension of a clean surface (surf aceT ension given
above), R is the universal gas constant and the other parameters are pre-
sented in Table 4.3. The transport of surfactant in the bulk fluid is given

18
by
Z I I
d
CdV + n · (v − vs )CdS = n · (D∇C)dS (4.2)
dt V S S

Transport of surfactant along an arbitrary surface

Z I I Z
d
ΦdS + m · (vt − bt )ΦdL = m · (Ds ∇s Φ)dL + sΦ dS (4.3)
dt S ∂S ∂S S

where sΦ is the source/sink of surfactant per unit area due to absorption and
desorption and is given by

sΦ = ka [Cs (Φ∞ − Φ) − βΦ] (4.4)

More information can be found[8]. surf aceT ension is the tension of the
f reeSurf ace when no surfactants are present and g is the gravitational ac-
celeration. To ensure constant mass flow f ixedF reeSurf aceP atches1(lef t)
option is used. It locks the left boundary so it can not move. If the boundary
also should be part of the solution the domain will look like the domain in
the background of Figure 4.5 and 4.6, respectively, instead of the original one
shown in foreground.

Figure 4.5: Inlet height change

19
 Figure 4.6: Inlet height change

4.2.4 polyMesh
As in all cases in OpenFOAM, the polyM esh directory contains the poly-
mesh. The domain for hydrof oil has already been setup and does not include
a blockM eshDict file but rather

boundary cells.gz faces.gz neighbour.gz owner.gz points.gz

However, the modDamBreak includes this file as shown in Section 3. The

polymesh is based around faces with internal faces connecting two cells and
boundary faces addressing a cell and a boundary patch. The owner and
neighbour cell labels therefore describe the connection between cells. Open-
Foam needs all these files in order to build up a domain/mesh[1].

4.2.5 transportProperties
transportP roperties specifies the viscosity for the application and contains
only it, namely

nu nu [0 2 -1 0 0 0 0] 1.5e-3;

20
4.3 system
In the system dictionary, the settings necessary in order to control the solu-
tion procedure are set. For the modDamBreak case these include;

controlDict faSchemes faSolution fvSchemes fvSolution

tetFemSolution

In controlDict, solver (interT rackF oam for this case) and time (start/end
time as well as timestep and writing controls of the output) parameters are
set, for example startT ime, deltaT and writeControl. The f aSchemes and
f vSchemes are taking care of the different types of discretisation schemes
used in the solution and f aSolution as well as f vSolution includes the
equation solvers, tolerances and other algorithm controls. It is important
to know that solvers in this case mean the ”linear-solver”, which is used
when solving the sets of discretised equation and not the application solver
(interT rackF oam), which is used to solve a specific case (modDamBreak,
hydrof oil etc). Algorithm controls refer to the settings for SIMPLE algo-
rithm used when implicity solving for pressure. Another file tetF emSolution
is also part of system and it is necessary in order to solve the equations cor-
responding to the given mesh motion.

21
22
Chapter 5

discretisation

5.1 blockMesh
The blockM esh utility provided in OpenFOAM is used in order to discre-
tise the computational space into a finite number of convex polyhedral cells
bounded by convex polygons. By ”convex” it is meant that the normals
point in the same direction. Also the computational domain is split into a
finite number of time-steps, which means that the equations also are solved
in time[9]. The discretised domain modDamBreak is shown in Figure 5.1.

Figure 5.1: Finite-volume discretization

23
5.2 makeFaMesh
As has already been mentioned, the settings concerning the finite-area dis-
cretisation is done in f aM eshDef initon (Section 4.2.2). The reason why the
makeF aM esh utility is used, is the fact that the f reeSurf ace class requires
it. The f reeSurf ace is shown in Figure 5.2.

Figure 5.2: Finite-area discretization

5.3 Automatic mesh motion

The polyhedral cells are split into tetrahedrals[10]. The mesh motion equa-
tion (Section 4.2.1) is discretised on the tetrahedral decomposition using
standard second-order finite element method. The finite element matrix ob-
tained is then solved using an iterative linear equation solver (Section 4.3,
tetF emSoltution). The motion of the boundaries are specified as bound-
ary conditions on the motionU field, and then used in the solution to up-
date point position. It is important to note that, the interface (i.e. the
f reeSurf ace) is moved according to the calculated displacements after the
iteration process for each time steps is carried out. This means that the
f reeSurf ace is fixed during the calculation for each timestep and moved
afterwards. If this was not the case the mesh-motion swept would be wrong

24
at the interface. The mass flux at the interface should be close to zero in
order to conserve the mass.

25
26
Chapter 6

interTrackFoam

An explanation of the solver interT rackF oam with all its features will now
be presented. First of all, as always in OpenFoam code, header files are
included in order to include classes which will be used in the computational
process.

#include "fvCFD.H"
#include "motionSolver.H"
#include "freeSurface.H"
#include "OFstream.H"

int main(int argc, char *argv[])

{
# include "setRootCase.H"
# include "createTime.H"
# include "createMeshNoClear.H"
# include "createFields.H"
# include "initContinuityErrs.H"

The f vCF D.H file includes the finite volume class definitons, motionSolver.H
includes for instance the twoDP ointCorrector mentioned in Section 4.2.1,
f reeSurf ace.H is used in the implementation of the freeSurface tracking
method considering a moving mesh and OF stream.H among others han-
dle compressed files. Header files are mainly used to make the code struc-
ture easier to use. The initial conditions are then specified (setRootCase.H,
createT ime.H, createM esh−
N oClear.H,createF ields.H and initContinuityErrs.H).
Info << "\nStarting time loop\n" << endl;

27
 for (runTime++; !runTime.end(); runTime++)
{
Info << "Time = " << runTime.value() << endl;

# include "readInterfaceSIMPLEControls.H"
# include "CourantNo.H"

interface.updateDisplacementDirections();

interface.moveMeshPointsForOldFreeSurfDisplacement();

interface.smooth();

The iteration process starts and will continue until the time set in controlDict
is reached, Section 4.3. The non-orthogonal controls are read
(readInterf aceSIM P LEControls.H) from the f vSolution file in order to
control the correction of the pressure field later on. The maximum as well
as the medium Courant number(CourantN o.H) are also calculated. When
the fields was created in the header file createF ields.H an interface of class
f reeSurf ace was created.
The subsequent tables (Table 6.1,6.2 and 6.3, respectively) will now explain
the commands that are used in the solver.

interf ace.updateDisplacementDirections()
The interf ace defined in createF ields.H is updated according to the
displacement direction. This is only done if normalM otionDir (Section
4.2.3) is allowed.

Table 6.1: Update the displacement directions

interf ace.moveM eshP ointsF orOldF reeSurf Displacement()

All internal mesh points are moved according to the displacement
of the interface (motionU ) in the previous timestep.

Table 6.2: Move internal mesh points according to motionU

28
 interf ace.smooth()
Smooth the interface if this option is specified, i.e. f reeSurf aceSmoothing
(Section 4.2.3). Also the twoDM otion has to be set to ”no”.

Table 6.3: Surfacesmoothing

The fluid flow equations are solved using a segregated SIMPLE procedure[11],
taking into account the kinematic and dynamic condition of the interface, as
well as the surface tension. The no-flux condition is satisfied in an iterative
sequence, providing the boundary condition for mesh motion on the free
surface (motionU ). This is shown in the following code lines

// --- SIMPLE loop

for (label timeCorr=0; timeCorr<=nTimeCorr; timeCorr++)
{
p.storePrevIter();

interface.correctBoundaryConditions();

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phiNet, U)
- fvm::laplacian(mu, U)
);

UEqn().relax();

solve(UEqn() == - fvc::grad(p));

volScalarField AU = UEqn().A();

U = UEqn().H()/AU;
U.correctBoundaryConditions();

UEqn.clear();

phi = (fvc::interpolate(U) & mesh.Sf());

The pressure from previous iteration is stored using p.storP revIter() and the
boundary conditions at the interface is corrected using interf ace.correct−

29
BoundaryConditions. It is important to note that the .correct−
BoundaryConditions does a couple of things, namely

void freeSurface::correctBoundaryConditions()
{
correctVelocity();
correctVelocityGradient();
correctSurfactantConcentration();
correctPressure();
}

A temporary matrix of of type f vV ectorM atrix is used to solve for the

velocity U and it looks like
→
−
∂ρ U →
− →
−
+ ∇ · φ U − ∇ · µ∇ U (6.1)
∂t
By using a relaxation factor (UEqn.relax()) defined in f vSolution (section
4.3) the iteration process is limited but less sensitive. The relaxed U Eqn is
then solved by coupling the pressure obtained from previous iteration to it,
i.e. →
−
∂ρ U →
− →
−
+ ∇ · φ U − ∇ · µ∇ U = −∇p (6.2)
∂t
The central coefficients from the U Eqn matrix (UEqn().A()) are calculated
and the velocity U . The U.correctBoundaryConditions command is used to
update the boundary conditions for the new velocity field while U Eqn.clear()
clears the vector to reduce peak memory. It is important to note that the
velocity distribution is stored without being affected by the pressure. The
velocities for each cell faces (phi) are obtained by interpolation. This is done
and necessary in order to later on calulate the fluxes at the faces.
Next step is to calculate the pressure and then correct it. Since an explicit
equation for the pressure is not available, one common way is to derive it sim-
ply by using the divergence of the momentum equation and by substituting
it in the continuity equation.

// Non-orthogonal pressure corrector loop

for (label nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
tmp<fvScalarMatrix> pEqn
(
fvm::laplacian(1.0/AU, p) == fvc::div(phi)
);

30
 pEqn().setReference(0, 0.0);
pEqn().solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn().flux();
}
}

# include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector

p.relax();

// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();

// Move mesh
interface.movePoints();

// Update motion fluxes

phiNet = fvc::interpolate(rho)*(phi - fvc::meshPhi(rho, U));

# include "freeSurfaceContinuityErrs.H"

Solve the pressure equaton iteratively and correct the fluxes at the cell faces
using control points[8]. Number of non-orthogonal corrector steps is set in
f vSolution, nN onOrthCorr and the iteration procedure is continued un-
til it is reached. As before for the U Eqn, a temporary matrix of type
f vScalarM atrix is used to define the pressure equation pEqn. By using
continuityErrs.H the continuity errors are calculated. Under-relax the pres-
sure for the momentum corrector and then correct the boundary conditions
for U . Move the mesh considering the new boundary conditions and update
the motion fluxes.

runTime.write();

Info << "ExecutionTime = "

<< scalar(runTime.elapsedCpuTime())

31
 << " s\n" << endl << endl;
}

Info << "End\n" << endl;

return(0);
}

The last thing done for each time step is to print the ”ExecutionTime”.

32
Bibliography

[1] OpenFOAM documentation, “User guide.” http://openfoam.com/

docs/user/, October 2010.

[2] OpenFOAM, “twodpointcorrector.” ~/OpenFOAM-1.5-dev/src/

meshTools/twoDPointCorrector, October 2010.

[3] A. O. Gonzlez, “Mesh motion alternatives in openfoam,” tech. rep.,

Institution of Applied Mechanics at Chalmers Technical University, Gte-
borg, Sweden, 2009. http://www.tfd.chalmers.se/~hani/kurser/
OS_CFD_2009/AndreuOliverGonzalez/ProjectReport_Corrected.
pdf.

[4] H. Jasak, “Dynamic mesh handling in openfoam,” tech. rep., Faculty

of Mechanical Engineering and Naval Architecture, University of Za-
greb, Croatia, London, England, 2009. http://powerlab.fsb.hr/ped/
kturbo/openfoam/papers/dynamicMesh_AIAA2009.pdf.

[5] H. Jasak and Z. Tukovic, “Dynamic handling in openfoam applied to

fluid-structure interaction simulations,” tech. rep., Faculty of Mechan-
ical Engineering and Naval Architecture, University of Zagreb, Croa-
tia, London, England, 2010. http://web.univ-ubs.fr/limatb/EG2M/
Disc_Seminaire/ECCOMAS-CFD2010/papers/01178.pdf.

[6] Wikipedia, “Wikipedia.” http://en.wikipedia.org/wiki/Radial_

basis_function, October 2010.

[7] OpenFOAM, “icodymfoam.” ~/OpenFOAM-1.5-dev/applications/

solvers/incompressible/icoDyMFoam, October 2010.

[8] Z. Tukovic and H. Jasak, “Simulation of free-rising bubble with soluble

surfactant using moving mesh finite volume/area method,” tech. rep.,
Faculty of Mechanical Engineering and Naval Architecture, University
of Zagreb, Croatia, London, England, 2008. http://powerlab.fsb.hr/
ped/kturbo/openfoam/papers/CFD2008.pdf.

33
 [9] Z. Tukovic and H. Jasak, “Updated lagrangian finite volume
solver for large deformation dynamic response of elastic body,”
tech. rep., Faculty of Mechanical Engineering and Naval Ar-
chitecture, University of Zagreb, Croatia, London, England,
2007. http://powerlab.fsb.hr/ped/kturbo/openfoam/papers/
TukovicJasak_NonLinElastodynamics_FAMENA_18-06-2007.pdf.

[10] Z. Tukovic and H. Jasak, “Automatic mesh motion for the unstruc-
tured finite volume method,” tech. rep., Faculty of Mechanical Engi-
neering and Naval Architecture, University of Zagreb, Croatia, London,
England, 2004. http://powerlab.fsb.hr/ped/kturbo/openfoam/
papers/MeshMotionJCPManuscript.pdf.

[11] OpenFOAMWiki, “Openfoamwiki.” http://openfoamwiki.net, Octo-

ber 2010.

34