第 24 回数値流体力学シンポジウム
E3-4

有限要素法を用いた自由境界を持つ流れの二次元シミュレーション
Simulation of Flow Problems with Free Boundaries in 2D using the Finite Element Method.
°Shi Han NG, 電通大, 東京都調布市調布ヶ丘 1-5-1, E-mail：shihanng@matuttis.mce.uec.ac.jp
Hans-Georg MATUTTIS，電通大, 東京都調布市調布ヶ丘 1-5-1, E-mail：hg@mce.uec.ac.jp
Shi Han NG, The University of Electro-Communications,
Graduate School of Informatics and Engineering, Chofu, Chofugaoka 1-5-1, Tokyo 182-8585.
Hans-Georg MATUTTIS, The University of Electro-Communications,
Graduate School of Informatics and Engineering, Chofu, Chofugaoka 1-5-1, Tokyo 182-8585.

We simulate flow problems with the Finite Element method for free boundaries whose motion is obtained by time
integration of the surface’s flow velocities. The volume conservation is reasonable, and for the low-viscosity-range,
we obtain good aggreement with the lubrication approximation. For the decay of a square fluid column, we findthat for high viscosity, the corner of the square travels towards the wave front, while for low viscosity, it stays
where it is, consistent with experimental results.

1.

Introduction
Many technical problem in the dynamics of complex
ﬂuids, like multi-phase-ﬂow and realistic ﬂow in porous
media and granular materials deal with the interaction
between ﬂuids and interfaces, rather than with the dynamics of ﬂuid alone. Our long-term aim is to simulate
multi-phase ﬂow problems in granular materials ”microscopically”, i.e. by simulating the particles individuallyas polyhedra via the discrete-element simulation, andthe ﬂuid ﬂow inside the ﬂuid space via the ﬁnite element
method. The ”ﬁrst principles” approach of simulating
the microscopic interactions allows to test the validityof approximations common in the ﬁeld.
The implementation of the following algorithms is in
MATLAB1 . We give no CPU-usage data and timings, as we are more interested to ”get it right” than”get it fast”. Once we have a satisfying algorithm, we
can focus on the performance optimization. Though weare aware that this approach may be computationallycostly, the CPU-performance is still increasing exponentially, and with the advent of GPGPU-cards and Multicore-PC-processors, it seems that substantial sustained
computing power becomes available in the hands of researchers far away from supercomputing centers with
their cumbersome application procedures. Our intention is, rather then focusing on massive-parallel implementations with time-consuming parallelizations, to develop a methodology where we can reduce the degrees of
freedom without lowering the accuracy of the geometrical description, and to use time-integration methods forwhich the time-step is only limited by the time-scale ofthe physical phenomena, rather than by the lamentablestability requirements of easy-to implement algorithms.

2.

Simulation method for the fluid part
Because we would like to run the simulation on arbitrary geometries, especially on pore-spaces which area result of the dynamics of granular particles, as wellas the surfaces, Finite Elements (FE) are the methodof choice for the discretization, as they can be usedalso with triangular meshes. Meshes in the pore-space
may become arbitrarily small or at least arbitrarily intractable for quadrilateral grids (see Fig. 1), without actually contributing much to the ﬂow ﬁeld. Therefore,
we would like not to be limited with our step-size by
the mesh-size due to the Neumann-stability conditions
for explicit integrators, so an implicit integrator seems
preferable. Most of the ﬂow problems involving suspended material will be incompressible, which from thepoint of theoretical approaches means that the pressures turn out as Lagrange-multipliers2 which inhibit
the velocities from compressing the ﬂuid. In numericalanalysis, these types of problems are called Diﬀerential
Algebraic Equations (DAE) and can be written in thefollowing form:
(
)(
)
(
)
M
GT (q)u0
f (q, u)
=
(. 1)
G(q)
0
λ
−gqq (u, u)This is an equation of motion for the masses M, theJacobian of the constraints G, the time derivatives ofthe velocities u0 and the external forces f, as well as the
Hessian of the constraint (rewritten with some vector
analytical identities) gqq (u, u). The DAE-form for the
Navier-Stokes equation by the FEM-formulation with
Backward-Diﬀerence Formula of second order (BDF2)
by Gresho and Sani4 is
[
](
)
1
M
+
K
+
N
(u
)
C
u
n+1
n+1
∆tn
=
CT
0
Pn+1
(
)
1
∆tn M un + fn+1
,
(2)
gn+1
where M is the mass-matrix of the problem, K is thestiﬀness matrix, N is matrix of the non-linear terms,
fn+1 are the external forces , the un+1 are the ﬂow
velocities and the Pn+1 are the pressures. We do not use
any up-winding. Comparison with eq. (1) shows thatthe pressures in eq. (2) indeed turn out as the Lagrange

Fig. 1: Snapshot of a Discrete-Element method forhexagonal particles with triangular grid in which a
pore-space which could hardly be discretized via aquadrilateral meshes.

1c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウム
E3-4”Economy of thought”: ”Ockham’s razor” is the principle that unnecessary principles should not be introduced to solve a problem: When our grid deﬁnes the boundary already, why introduce another”grid” for the boundary. In particular, we arerather averse to introduce an entity which does notexist in nature: The ﬂuid surface should also be
its physical boundary, so the surface information
should also be suﬃcient to model the boundary.

velocities andpressures
only velocities

Fig.
2: Friedrichs-Keller grid (”FEM-grid”) for
P2P1-elements, where the velocities are given at thesites with the crosses and the circles, and the pressures
are deﬁned on the sites with the crosses only.

Mechanic Impedance: When we introduce, additionally to the boundary, a geometric entity with which
the ﬂuid interacts, there is a risk that this will alter
the mechanic impedance (the way that mechanic
signal propagate over the boundary, or how they
are reﬂected there) without the possibility to control or evaluate this eﬀect.

parameters. Therefore, instantaneous changes of thepressure are possible, there is no necessity, as in the casefor the Marker-and-Cell method, for a smooth iteration,
which behaves eﬀectively as an ”equation of motion”
for the pressures and leads to restrictions with respect
to the timestep and the accuracy. Though pressures
which ﬂuctuate strongly are an indication of a too large
timestep, it is reassuring to have method which indicates
too large timesteps without blowing up.
Eq. (2) is implicit, so that the values of the newtimestep n + 1 must be compute via Newton-Raphson
iteration from the previous values, where the ”Matrix”
is the Jacobian of the velocities. Usually we needabout three iteration for ﬁxed boundary conditions. As
the BDF2-method is not self-starting, i.e. at the ﬁrst
timestep, additionally to the initial conditions data from”before” the initial conditions are needed, and additionally, consistent initial conditions are needed for the pressure, so that the incompressibility condition is not violated. Therefore, we start our solution from a stationary state, computed via Newton-Raphson iteration forthe stationary Navier-Stokes equation. The approachof integrating out the time dimension via a solver for
ordinary diﬀerential equations (ODE), instead of discretizing the time direction also via Finite Elements issometimes referred to as ”semi-discretization.” For the
spatial discretization, we use P1 P2 -elements, i.e. thepressures are approximated with aﬃne elements, whilethe velocities are approximated with quadratic elements(see Fig. 2).

Practicability: In the future, we want to simulate
porous media with particles and many thousand
ﬂuid interstices etc. Introducing such an immense
number of additional surfaces, including the corresponding overlaps or intersections may lead to
non-unique algorithmic choices.
Analogy: Finite Elements methods in structural mechanics don’t need additional data structures to describe the surfaces. So shouldn’t we be able to dowithout additional constructs als in ﬂuid modeling?
Another aspect is that our semi-discrete implementation is already a formulation of the FEM-equations
with ordinary diﬀerential equations. Because the FEMpart of the simulation yields the velocities of the ﬂuid atevery lattice point, we should therefore be able to calculate new positions of the lattice points at the surface,
and use the new position of the lattice points without
any additional data structure. So if the position of thesurface at the discrete timestep t = n would be given(i,j)
for the lattice site (i, j) as xn , and we had obtained
the velocity by our FEM-procedure on the free surface(i,j)
as vn (for a sketch, see Fig. 4) we could compute thenew position of the lattice point via an Euler-integration
for a time-step τ as3.

Philosophy of the surface modeling
Conventional Eulerian, i.e. grid based approaches
leave the grid ”as it is” and introduce additional data
structures to trace the surface. In the Marker-and-Cell(MAC) method5, 6 , these ”markers” are used to interpolate the surface between the grid points for which theﬂow is computed. More recent methods solve partial
diﬀerential equations to determine the position of thesurface with various approaches (level set method7 , volume of ﬂuid method8 , advection equation9 ). As we havea Finite Element simulation at hand which allows to restructure and re-mesh the grid in every timestep10 , wewould rather stick with our FEM-grid without introducing additional data structures for a variety of reasons:v(i,j)(3)

Of course, due to the bad stability- and accuracyproperties, it will be advisable to look for ODE-solvers
with better properties. Because we do not want to deform our P2 P1 -elements, we move only the corners of
the triangles and interpolate the center points of theedges accordingly.
With respect to the boundary conditions one modiﬁcation in comparison to conventional ﬂuid simulations
becomes necessary: Imagine a gravity wave (a wave

Expected Boundary
movement
Expected
movementv(i+1,j)surfacesurfaceboundaryv(i−1,j)(i,j)+ vn(i,j) τ.
xn+1 = x(i,j)
n

Fluid
fluid

fluid

Fig. 4: Gravity wave (left) and schematic discretization near the surface in the inset, to the right. Grid
points with normal velocities and tangential velocity set
to zero (◦), all velocity components velocities computed(×) and the grid points with normal velocity set to zero
while the tangential velocity is computed (¤).

Fig. 3: Gravity wave and discretization near the surface(inset) with normal velocities set to zero (circles) andvelocities computed (crosses).
2

c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウム
E3-4

which moves under the inﬂuence of gravity, not a oscillatory solution of the ﬁeld equations of general relativity, of course) in a container of length l and wave
length 2l, with a phase so that the extrema are at theboundary, see Fig. 4. Obviously, the surface must move
down at the maximum, and must move up at the minimum. While the tangential velocities for the elementson the boundary away from the surface must be set to0 ”inside” the computational domain, for the elementon the surface the tangential component of the velocity
must be computed, else no movement of the surface ofthe boundary is possible. At the same time, the normal velocity for this element must be zero. In an earlier
stage of the program development, we forgot and obtained curious, elongated triangles. This formulation
with moving elements has the advantage that, as eq. (2)allows the inclusion of arbitrary external forces fn+1 ,additional models for the surface tension can be implemented in Finite-Element formulation, making use ofwell-established techniques from structural mechanics
like the beam models and equations.
4

4

3

3

2

2

1

1

0

0

2

4

6

8

0

4

4

3

3

2

2

1

1

0

0

2

4

6

8

0

0

2

4

6

8

0

2

4

6

8

available only at discrete times, so one-step methods(Runge-Kutta-type methods which use several function
evaluations per timestep) cannot be implemented. For
the problem in Fig. 5, we tentatively implemented an
Euler-Step, as well as a second order Adams-Bashforthmethod (AB2), which computes the positions at thenext timestep n + 1from the current positions xn , thecurrent velocities vn and the previous velocities vn−1 as(
xn+1 = xn + τ

4.

Test implementation for a gravity wave
To test our implementation, we start with the simulation of an actual gravity wave. First of all, wewant to know how well the volume is conserved for different integrators. We compute the initial stationary
state by setting the pressure at the upper free surface
to zero. The Newton-Iteration for the stationary stateyields the hydrostatic pressure, higher pressures below
the maximum of the wave, and lower pressure of theminimum. Because we have to make use of the velocities at the intervals which are determined by our BDF2integration of the internal ﬂow, we have the velocity
Euler
2nd Order Adams-BashforthWithout interpolation
With interpolation

31.99931.99631.99831.99431.992

Area

Area

(4)

As a tentative parameter of the verisimilitude of thesimulation, we take the conservation of the volume. For
our gravity wave, the viscosity was chosen as unity (inSI-units), the system size width 8 × (average) height
4 (in units of [m]). We prefer to avoid dimensionlessunits, because with free surfaces we have not only inertia eﬀects from the ﬂow, but also from the gravity, andin the future, another scale may enter due to the surface
tension. For easier comparison with the experiment, we
would like to stay with SI-units. The initial pressure
proﬁle we obtained with zero pressure on the surface,which lead to a hydrostatic pressure distribution, where
the pressure was higher below ”hills”, and lower under”valleys” of the wave. For the gravity wave in this section, the re-meshing is done is such a way that we wecompute the new position of the new mesh points on thesurface. The new position for the the new mesh points
below the surface is chosen in the horizontal direction
vertically below the surface points. The new position in
the vertical direction is computed equidistantly along
the vertical. As can be seen in Fig. 5, the elements at
the left and right boundary are shortened in this process, while the elements which are what was initially the
maximum of the wave, are elongated.
As an initial conﬁrmation, we made a test-run witha ﬂat surface and found that in the absence of the motion of a surface, the volume is constant. The result for
the volume for τ = 0.008 and a spatial discretizationof Nelem = 4760 elements can be seen in Fig. 6: Obviously, the volume conservation is better for the higher
integration order integrator. When the motion of thesurfaces slows down, the increase of the error also slows
down. Actually, we have up to now discussed only how
in eqs. (3) and (4) the old velocities to obtain the newpositions, but we have to discuss at which position weevaluate the velocities. Because we have a Eulerian formulation, we should not just use the old positions at
diﬀerent timesteps: As the philosophy of our method
is, rather than trace the movement of the lattice points
in a Lagrangian way, just to predict the new position ofFig. 5: Snapshots from the time evolution of the gravity wave at the initial condition (t = 0.0, upper left),(t = 1.2, upper right), t = 3.6 (lower left) and t = 5.6(lower right).

31.998)13
vn − vn−1 .
2
2

31.99
31.98831.99731.99631.99531.98631.99431.984
0

1

2

3
Time

4

5

0

Fig. 6: Time evolution of the volume for the integration
with Euler (×) and Adams Bashforth of second order(◦) for τ = 0.008 and 4760 grid-points. Symbols are notdrawn for every timestep to keep the ﬁgure readable.

1

2

3
Time

4

5

Fig. 7: Time evolution of the volume for the second order Adams Bashforth method with (-) and without (×)
interpolation for stepsize of τ = 0.008 and Nelem = 4760
elements.
3

c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウムE3-43231.99831.99531.996
Area

Area

31.99
31.99431.992
31.99

31.98

31.98831.97531.986
0

1

2

3
Time

4

31.97

5

0

Fig. 8: Time evolution of the volume with AB2 forthe stepsize of τ = 0.008 and diﬀerent numbers ofmesh points Nelem = 1088 (×), Nelem = 2600 (◦) and
Nelem = 4760 (4).
the surface, we choose for AB2()
3
1
xn+1 = xn + τ
vn (xn ) − vn−1 (xn ) ,
2
2

(5)

4

3

Area

2

1

1

0

1

0

2

3
Time

4

5

4

2 5.

2

Simulation of the collapse of a water column
Now we can turn to more interesting problems forthe sake of verifying the algorithm. A popular test1 case, especially for Lagrangian1 (particle-based) methods
is the time evolution of a water-”step”. While boundary conditions are not a great concern for particle methfor our FEM-method we have to choose the ”right”
0 ods,
0 node on the surface, theboundary conditions: At every
0pressures
2 6 are set2to4zero,
8 at every
0 4 node
6 which does
2 6 not
8
touch a wall, the velocities are computed, and at every node which touches a wall, the normal velocities

31.99
31.9884from the time evolution of the gravity wave with one(left) and two maxima (right) at t = 1.6.

31.9942

3Time

3 Fig. 11: Zoom into the corners
3 of the the snapshots

31.99631.9922

accuracy of the volume conservation becomes (as is to
be expected) better with decreasing mesh size, i.e. increasing number of grids. In the next step, we want toverify the eﬀect of the timestep for a system with 1088
elements. As can be seen in Fig. 9, the improvementfrom τ = 0.06 to τ = 0.008 is monotonous. We have toconclude that the remaining error (the deviation from32) is rather due to the lattice discretization than due
to the time discretization, especially if we compare withthe result for 4760 grid points in Fig. 8.
According to Panton12 (p. 536), the ﬂow of a gravity
wave must be irrotational. By inspection of the ﬂow
ﬁeld, we could verify that for various viscosities thiswas the case. It would it is now interesting how ouralgorithm behaves if we have more wave maxima. Fora system with the same amplitude, but half the wave
length with our highest resolution of Nelem = 4760 anda stepsize of τ = 0.008, we can see (Fig. 10) that theerror is signiﬁcantly larger than for a single wave. In
improvement in the accuracy will be only possible witha smaller mesh-size. Nevertheless, for the time being,
let us stick with the error in the conservation without
any additional tricks like rescaling of elements etc, because the error in the volume may have its uses, as alarge error may indicate deﬁciencies in the spatial discretization. The number of Newton-Raphson iteration
varies between four and two iterations, it depends on
the ﬂow velocities, but not on the resolution of the grid.
Actually, the source of the error in the volume conservation seems to be largely due to the behavior at the solidboundary: Fig. 11 shows that the out-most elements on
the surface have a rather weird angle, compared to therest of the surface, which is rather smooth.

That our argument about the interpolation is justiﬁed
can be seen from Fig. 7: The error in the volume conservation is larger without interpolation, though not signiﬁcantly. So in the following, we will use interpolated
velocities if the surface moves downward. If the surface moves upward, we would need extrapolation, but
we have to decided to do without, as neither is extrapolation an exact science, nor, with respect to the volumeconservation for interpolation, can be expect signiﬁcant
improvements.
Let us next study the eﬀect of the mesh size for theconservation of the volume. For the same stepsize ofτ = 0.008, we have investigated meshes with 1088, 2600
and 4760 mesh-points.
As can be seen in Fig.
4
4 8, the3

1

Fig. 10: Time evolution of the volume for Nelem = 4760
and stepsize τ = 0.008 for integration with AB2 for asingle (×) and double (°) wave.

i.e. the new position of a grid-point of the surface iscomputed from the current position of the grid-point
and the previous ﬂow velocities at the current position.
If at step n − 1 the grid-point was not at xn yet, wehave to obtain it by interpolation. This is easy in ourFEM-formulation if xn is lower than the previous surface height, i.e. if the surface is falling. Unfortunately,
if the surface is rising, we would have to extrapolate topoints where there was no ﬂuid at all, so we have toabandon the use of eq. (5) and just use the velocities of
the grid point at the previous timestep in AB2,()
31
xn+1 = xn + τ
vn (xn ) − vn−1 (xn−1 ). .
(6)22

31.998Single wave
Double wave

31.985

50Fig. 9: Time evolution of the volume for 4760 elements02
04
and various timesteps for integration with AB2 for (frombelow) τ = 0.06 (×), τ = 0.03 (◦), τ = 0.015 (4) andτ = 0.008 (−).
4

c 2010 by JSFM
Copyright °4

第 24 回数値流体力学シンポジウム
E3-4
2.5

0.05
0.04

2
Z /L

0.03
1.5

0.02
0.01

1
0.5
0

0

Huppert
FEM (Small Meshsize)
FEM (Large Meshsize)
0

50

100
p
t 2g /L

150

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.05
0.04

200

0.03

Fig. 12: Comparison of time evolution of the advancement of a wave front under the lubrication approximation with our simulation with large (circles, 400 triangles) and ﬁne (crosses, 3000 triangles) grid for 1/ν = 10
and stepsize τ = 0.01.

0.02
0.01
0

Fig. 14: Snapshot at t = 8.5[s] for the ﬁne (above) andrough mesh (below) for the collapse of the squre watercolumn at 1/ν = 10 and stepsize τ = 0.01.

are set to zero, and the tangential velocities are computed. (Choosing the pressure hydrostatic on the vertical face of the step would inhibit any evolution, as
the pressure ﬁeld would be the same as in the presenceof a wall). The next issue concerns the choice of thegrid: Because the deformations are not small, and cannot only be taken into account by rescaling the gridin the vertical direction, we have to re-mesh the grid
globally. For this, we use the algorithm which we originally developed for grains in ﬂuid 13 . As it was initiallywritten for convex domains, some modiﬁcations havebecome necessary.
Huppert14 derived the time dependence the advancement of the front Z with time t as
( 3 )1/5
gq
Z(t) = 1.411
t1/5
(7)
3ν

wavefront looks rather ragged both for ﬁne and roughmeshes. For this, our meshing algorithm responsible
which adds new nodes if the distance between nodes becomes larger than 1.5 the choosen characteristic length
for the nodes. It seems that in the next version of theremeshing algorithm, we also have to take care of thestructure of the ﬂow ﬁeld. Near the left wall, one can seethat when we supply a good resolution, actually a thinﬁlm develops (Fig. 14, above), while for too rough grids,(Fig. 14, below), the ﬁlm is actually cut oﬀ, and only ameniscus remains near the wall. This shape of the ﬂuid
surface is purely an eﬀect of the viscosity, as we have notyet modeled cohesive or adhesive forces in the simulation. On the one hand, it is conforting to see that we canrun our algorithm both with physically small resolution(to capture even thin ﬁlms), or with rough resolution if
we are only interested in a more large scale dynamics.
Certainly is not possible for many algorithms to be run
beyond multiples of the actual resolution, so we consider our implementation as a success (due, probably,to the stiﬀ BDF2-integrator for the ﬂow ﬁeld). On theother hand, we have no error indicator if we are actually
above the resolution scale of physical phenomena, so no
automatic error detection for the spatial resolution ispossible (for the time step, the Newton-iteration does
not converge if the timestep is choosen too large.)
For higher Reynolds numbers, we want to compare
our simulation with experimental data rather than withother simulations or theories. Especially with the experiment by Martin and Moyce16 is of interest, as theseauthors argue that for their dimensions of the vessels,
air resistance and surface tension can be neglected, exactly the conditions which we have in our simulation.

for an initial area q, gravity g and viscosity ν under
the lubrication approximation, i.e. under the assumption that eﬀects due to surface tension etc. could be
neglected and the height h of the water column is muchsmaller than its width l15 . Thought the latter assumption is rather problematic in our case, our data in Fig 12
compare rather well with eq. (7) at 1/ν = 10. Changesof ν over a reasonable range of parameters did not leadto any changes in the curve. We leave it to the reader
to interpret the good correspondence as either a conﬁrmation of our simulation method or the lubrication
approximation. As the agreement with the lubrication
approximation becomes better for larger mesh sizes (circles in Fig 12) , but the error in the volume in Fig 13
increases up to about 1% for the ﬁnd grid, and 1.5% for
the rough mesh, it is worth having a look at the actual
wave fronts. We see in Fig. 14 that at t = 8.5[s] the
-3

x 10

InitialFilling

3.3

0.074s

Area

3.29

0.10s

0.24s

3.28
3.27

Exact Area
FEM (Small Meshsize)
FEM (Large Meshsize)

3.26

Fig. 15: Advancing of the ﬂuid front in the experimentof Martin and Moyce16 for the collapse of a ﬂuid columnof 57 × 57 mm. The line-width is approximately thewidth of the shadows in the original frames, graphics
is mirrored compare to Martin and Moyce, because forthe simulation it is more convenient to have the origin
on the left.

3.25
0

2

4

6
Time

8

10

Fig. 13: Error in the volume for the collapse of thesqure water column at 1/ν = 10 and stepsize τ = 0.01.
5c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウム
E3-4

upper surface of the column is convex, and their experimental data are not clear enough to determine theshape of the column. In our simulation, where we startwith a singular surface (the gradient has a jump at theright corner), it turns out that the singularity is preserved. Though there is a certain uncertainty, as wehave to investigate the eﬀects of the grid discretizationin more detail, if we assume that our simulation gives
the correct solution, we have to conclude that such a
singularity in the surface does not vanish on hydrodynamical grounds, but due to surface tension, which is
not yet implemented in our code.

t=0
t=0.5
t=3.5
t=7.0
t=10.5

0.05
0.04
0.03
0.02
0.01
0
0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

t=0
t=0.5
t=3.5
t=7.0
t=10.5

0.05
0.04
0.03

6.

Summary and Conclusions
We have implemented a Finite-Element method withfree surfaces where we integrate out the motion of thesurface elements according to the velocity data obtained
from the FEM-method on the surface. Compared toconventional eﬀorts, which try to solve partial diﬀerential equations for the motion of the surface, the additional eﬀort in our method with respect to new data
structures etc. is negligible. The method shows structures (development of a non-convex horizontal surface)which is lacking in particle-based simulations for thesame problem. This gives us hope that the method canbe extended (with better re-meshing algorithms and theinclusion of surface tension) so that it may yield a higher
degree of realism (not to say precision) than other methods which rely on unphysical surface- or particle dynamics. Certainly, the results are aﬀected by remeshingand mesh reﬁnement, and certainly we could identify aproblem in the remeshing of rapidly advancing fronts.
Beyond that, we want to implement surface tension,modeled by beam equations for which the curvature isevaluated at the ﬂuid surface, as well as better gridding
algorithms with a ﬁner mesh towards the surface thanin the bulk of the ﬂuid.

0.02
0.01
0
0

0.02

0.04

0.06

0.08

0.1

0.12

0.14

Fig. 16: Time evolution for the fronts in for for theﬁne (above) and rough mesh (below) from Fig. 14 forthe collapse of the squre water column at 1/ν = 10 andstepsize τ = 0.01. The corner travels towards the front.

Number of Iteration

We have extracted the outlines of the graphics and associated them with the corresponding with the time data
data in Fig. 15. In the following, we focus on the shapefor short time scales, which we consider a more meaningful for the verisimilitude of the simulation than thelong-term time evolution, in which geometric details are
smeared out and momentum-conservation and energydecay dominate.
What is striking in the high-speed pictures (300frames/s) of Martin and Moyce16 is, that the initial
sharp step-shape between initial water level and verticalboundary (dotted oval in Fig. 15) is, within the limits
set by the shadows in the ﬁgures, rather well preserved
during the collapse. Moreover, the upper surfaces looses
convexity both on the left and on the right side. The resilience of the ”corner” on the right is rather surprising,
not to say counter-intuitive, but as the pressures in upper right corner, are negligible, there are no forces which
could cause the decay of the angle. For low inverse viscosity 1/ν = 10 in Fig. 16, the cusp from the initialsquare proﬁls traveled fast towards the front. Within
the deadline for this preprint, we were not able to run
the simulation at the necessary inverse viscosity (about
1/ν = 10000), but for the largest inverse viscosities we
could reach ( 1/ν = 500), the cusp indeed stayed where
it was, see Fig. 17 As we have no surface tension implemented, the cusp becomes sharper and sharper, as in
the case of triangular ocean waves above a shear ﬂow,
which also steepen and become sharper until they break.
In the particle-simulations by Koshizuka et al.17 , the

0.04

0.08

3
Time

4

5

Using at the nth timestep the velocity vn−2 from twotimesteps earlier in eq. (8) may not contribute to theimprovement of the accuracy in the case of strong movement of the surface. We therefore replaced vn−2 by using the centered-diﬀerence formula of the acceleration

0
0.06

2

Appendix: Other Integrators
Beyond AB2, we have implemented a variant of thethird order Adams-Bashforth method (AB3)
(
)
23
4
5
xn+1 = xn + τ
vn − vn−1 + vn−2 .
(8)
12
3
12

0.01

0.04

1

Fig. 18: Time evolution of the volume for the integration with Euler (×), Adams Bashforth of second (◦) andthird (full line) order for τ = 0.008. Symbols are notdrawn for every timestep to keep the ﬁgure readable.

0.02

0.02

31.985

0

0.03

0

31.99

31.98

t=0
t=0.0675
t=0.1350
t=0.2025
t=0.2700

0.05

31.9950.1

Fig. 17: Time evolution for the fronts in for for theﬁne (above) and rough mesh (below) from Fig. 14 forthe collapse of the squre water column at 1/ν = 10 andstepsize τ = 0.01. The corner travels towards the front.

an+1 =
6vn − vn−2
,
2τ

(9)

c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウム
E3-4

which is available as the ”right hand side” from eq. (2)so thatvn−2 = vn − 2τ an+1 .
(10)

third order, so will be the results for the xi . As our velocities and acceleration are only computed up to second
order with our BDF2-integration in eq. (2), we have toconclude that the third order formula is meaningless anddoes not lead to improved accuracy. We could conﬁrm
using manipulated output from the BDF2-integration:
When we intentionally introduced an error there, theAB3-integration was much more strongly aﬀected thanthe AB2-integration.
When we were already about it, we decided to try out
other methods: The midpoint-rule

In that case, our variant of the AB3-method looks like
(
)
7
4
5
xn+1 = xn + τ
vn − vn−1 − an−1 .
(11)
3
3
6
In analogy to eq. (6), we used an−1 (xn ) if interpolation was possible, and an−1 (xn−1 ) where extrapolation would have become necessary. Initially, the result
looked promising, as can be seen in Fig. 18: The volumeconservation was better than for the ﬁrst and second order method.

xn+1 = xn + vn τ + an τ 2 ,named after the fact that it can be derived from centereddiﬀerences which can be derived from eq. (10) . Further,
we wanted to try out the Adams-Moulton-formula inthird order

32.00832.006Area32.004xn+1 = xn + τ

32.00232

31.99631.994
1

2

3
Time

4

5

Fig. 19: Time evolution of the volume for 4760 elements and various timesteps for integration with AB3for τ = 0.03 (◦) τ = 0.015 (4) τ = 0.008 (−) τ = 0.004(¤) and τ = 0.002 (¦).
Unfortunately, when we decreased the stepsize, itturned out that the algorithm did not converge in thirdorder, but only linearly, as can be seen in Fig. 19. In
principle, numerical algorithms may not converge whendecreasing the timestep due to the fact that a decreasein timestep, while decreases the discretization error,
but at the same time increases the rounding error (see.
Fig. 20). Nevertheless, in our case, that cannot be thereason, as AB2 converges without problems with thesame timestep.
Error
Total error
Worse
accuracyxn+1 = xn + τ(13)1(8vn + 10τ an + 4vn−1 ) .
12

(14)

In Fig. 21, we have contrasted the results for diﬀerent
methods: The only method which gives more accurate
results than AB2 for the timestep used is AB3, andthat only because it is not yet fully converged: With
smaller timestep (necessary for larger velocities and/or
smaller viscosities), its error will become larger than in
Fig. 21, so for our input data (computed with BDF2),
Adams-Bashforth in second order gives the most satisfying results.
31.99831.996Area

31.99431.99231.99



optimal31.98831.986

Rounding errorN

31.9840

Betteraccuracy o1(5vn+1 + 8vn − vn−1 ) .
12

Because using implicit integrators with the ”new” velocities vn+1 makes the re-computation of the whole
BDF2-step necessary, we tried to avoid this by modifying the Adams-Moulton formula in eq. (13) to an ”explicit formula” using the centered diﬀerences in eq. (10)for timestep n so that

31.9980(12)Truncation error
1/N m +1

1

2

3
Time

4

5

Fig. 21: Numerical error in the conservation of the volume for several methods 4760 mesh points and timestep
τ = 0.008: Euler (×), Midpoint (¤), 2nd Order
Adams-Bashforth (◦), 3rd Order Adams-Bashforth (4)and 3rd Order Adams-Moulton (−).

total number of timesteps N
LargerSmallerFig. 20: Schematic representation of the numerical error
as a sum of the truncation error and the discretization
error: The smallest timestep does not give the mostaccurate result, but there is an optimal timestep.
The introduction of the acceleration an−1 in eq. (10)instead of vn−2 is not the reason for the loss of accuracy
either: For a time integration of x = −ẍ with ”exact”
input-data, eq. (11) lead to a smaller error than witheq. (11). In In fact, we had overlooked the fact that forn-th order integration with Adams-Bashforth, all variables which enter in the formula must be obtained to nth order. If the vi , ai used in eqs. (8,11) are of lower thanReferences
1. The MathWorks, Inc. (2009), MATLAB Version
7.9.0.529 (R2009b).
2. Lanczos, C., The Variational Principles of Mechanics, 4th Edition (Dover Publications, 1986).
3. Hairer, E. and Wanner, G., Solving Ordinary Diﬀerential Equations II: Stiﬀ and Diﬀerential-Algebraic
Problems, 2nd Edition (Springer-Verlag, 1996).
7

c 2010 by JSFM
Copyright °第 24 回数値流体力学シンポジウム
E3-4

4. Gresho, P. M. and Sani, R. L., Incompressible Flowand the Finite Element Method Volume 2: Isothermal Laminar Flow (John Wiley and Sons, Ltd.,
2005), pp. 805–806.
5. Browne, L. W. B., The Marker and Cell Technique,in Numerical Simulation of Fluid Motion, ed. Noye,J. (North-Holland Publishing Company, 1978), pp.
223–247.
6. Harlow, F. H. and Welch, J. E., “Numerical Calculation of Time-Dependent Viscous Incompressible
Flow of Fluid with Free Surface”, The Physics of
Fluids, 8(12), (1965), pp. 2182–2189.
7. Osher, S. and Fedkiw, R., Level Set Methods andDynamic Implicit Surfaces (Springer, 2003).
8. Hirt, C. W. and Nichols, B. D., “Volume of Fluid(VOF) Method for the Dynamics of Free Boundaries”, Journal of Computational Physics, 39(1),(1981), pp. 201–225.
9. Caboussat, A., Maronnier, V., Picasso, M., andRappaz, J., Numerical Simulation of Three Dimensional Free Surface Flows with Bubbles, in Lecture
Notes in Computational Science and Engineering,Vol. 35 (Springer-Verlag, 2003), pp. 69–86.
10. Sakai, Y. (2007), “Dynamic Remeshing for Finite
Element Method (in Japanese)”, Undergraduate’s
thesis, The University of Electro-Communications,
Department of Mechanical Engineering and Intelligent Systems.
11. Engeln-Müllges, G. and Uhlig, F., Numerical Algorithms with Fortran (Springer, 1996).
12. Panton, R. L., Incompressible Flow Second Edition(John Wiley & Sons, Inc., 1996).
13. Ng, S. H. (2009), “Generation of triangular disordered Grids for Fluid Dynamics Simulation using
the Finite Element Method”, Undergraduate’s thesis, The University of Electro-Communications, Department of Mechanical Engineering and Intelligent
Systems.
14. Huppert, H. E., “The propagation of twodimensional and axisymmetric viscous gravity currents over a rigid horizontal surface”, Journal of
Fluid Mechanics, 121, (1982), pp. 43–58.
15. Gratton, J., Mahajan, S. M., and Minotti, F., Non
Newtonian Gravity Creeping Flow (Trieste: International Centre for Theoretical Physics, 1988), pp.
1–18.
16. Martin, J. C. and Moyce, W. J., “Part IV. An Experimental Study of the Collapse of Liquid Columns
on a Rigid Horizontal Plane”, Philosophical Transactions of the Royal Society of London. Series
A, Mathematical and Physical Sciences, 244(882),(1952), pp. 312–324.
17. Koshizuka, S. and Oka, Y., “Moving-Particle SemiImplicit Method for Fragmentation of Incompressible Fluid”, Nuclear Science and Engineering, 123,
(1996), pp. 421–434.

8

c 2010 by JSFM
Copyright °