# 3.1.1. Harlow and Welch (1965) “Numercial Calculation of Fluid with Free Surface”¶

Harlow, F.H. and Welch, J.E. “Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface” The Physics of Fluids, Volume 8, 1965

## 3.1.1.1. Introduction¶

New technique for numerical solution of free surface problems, designated the “Marker and Cell Method”

Physical assumptions:

Time dependent motion

Viscous

Incompressible

2D with

Cartesian CoordinatesUniform body force (gravity) may act on the fluid

Possible boundary conditions:

Wall boundaries for an

irregular boxLines of reflective symmetry

Prescribed pressure may be applied to the surface -

variable in position and time

Method uses:

Finite difference method

Full Navier Stokes Equations

Dependent variables are pressure and two component velocity

Fluid divided into cells - variables in each cell are single average values

Time divided into increments

Small spatial and temporal steps are needed

Comparison with other methods:

Most of work directed at compressible flows in 1965

Incompressible flow was being directed at weather forecasting

Non-steady viscous flows were being studied - von Karman vortex street behind a cylinder - i.e. Fromm and Harlow (1965) - primary dependent variables were vorticity and the stream function - marker and cell may be easier because BCs are easier to apply and physical significance of modification is easy to visualise

## 3.1.1.2. Description of the Method¶

### 3.1.1.2.1. Outline¶

The derivation of the finite difference method is based on:

The initial velocity field is known and

**divergence free**Coordinates of a set of marker particles is known - showing which regions are fluid and which are not fluid

Pressure field is calculated so that

**the rate of change of the divergence of velocity is zero**, which requires the solution of Poisson’s Equation - which can be done via an SOR methodTwo components of acceleration are calculated

Marker particles move according to the velocity field

Adjustments are made for passage of marker particles across cell boundaries - i.e. the velocity changes when a cell is emptied and filled with fluid

The marker particles are not like the particle-in-cell method, i.e. the marker particles **do not** participate in the calculation, where as in the particle-in-cell method, they do.

### 3.1.1.2.2. Finite Difference Equations¶

Use the Conservative Form of the Navier-Stokes Equations:

Harlow and Welch use \(\psi\) to denote a pressure at a constant density i.e. \(p \over \rho\)

### 3.1.1.2.3. Finite Difference Discretisation¶

Indices i and j count the

**centre**of each cellCell boundaries are labelled with half-integer values

The superscript n+1 is used for time and where it is not used n is implied

The Finite Difference approximations are:

Some of the velocity values are not located at the centre of the mesh diagram, so an average of adjacent values is taken - i.e. **the average of the velocity at the closest two points - which could be the vertical or horizontal direction, depending on which are the adjacent points**

Harlow and Welch now take the divergence of the momentum equation:

Now set:

Applying the principle of continuity:

The divergence of velocity at the new timestep \(\mathbf{u}^{n+1} = 0\)

The divergence at the current timestep \(\mathbf{u}^n \ne 0\)

Re-arrange for a form of the Poisson Equation **at time n**:

Harlow and Welch say that you can neglect viscosity in the Poisson Equation **with strict convergence critera**

The loop in time is like this:

Compute the RHS of the Poisson Equation at n

Compute pressure using Poisson Equation at n

Compute velocities at n+1 using Momentum Equations (from the pressure previously calculated at n)

**The divergence of velocity must also vanish at the boundaries**

### 3.1.1.2.4. Marker Particles¶

Harlow and Welch developed (a now obsolete) particles method for tracking the motion of the free surface

They actually track the motion of particles throughout the fluid (not just at the free surface)

Linear interpolation is performed to calculate the velocity with which a particle should move

The interpolation weighting depends on the distance of the particle from the nearest velocity point

### 3.1.1.2.5. Boundary Conditions at Rigid Walls - Velocity¶

Rigid walls may be of two types:

No-slip

Free-slip (could be a plane of symmetry - or a greased surface)

At the wall \(v = 0\)

On the fluid side \(v \ne 0\)

On the outside \(v' = -v\) (no slip)

On the outside \(v' = +v\) (free slip)

(Similarly for horizontal walls)

In general:

the tangential velocity reverses or remains unchanged for no-slip or free slip

the normal velocity reverses for a free slip wall

the normal velocity is calculated to ensure the divergence vanishes for an external cell for a no-slip wall

### 3.1.1.2.6. Boundary Conditions at Rigid Walls - Pressure¶

The boundary condition for \(\psi\) must be consistent with the velocity

#### 3.1.1.2.6.1. Free slip wall¶

**Vertical wall** From the momentum equation for \(u\) (with the reversal of all normal velocities and no change in the tangential velocity):

\(+ \Rightarrow\) fluid is on the left of the wall

\(- \Rightarrow\) fluid is on the right of the wall

**Horizontal wall** From the momentum equation for \(v\):

\(+ \Rightarrow\) fluid is below the wall

\(- \Rightarrow\) fluid is above the wall

#### 3.1.1.2.6.2. No slip wall¶

Condition is that the divergence of velocity on the outer cell must be zero:

**Vertical wall**

\(+ \Rightarrow\) fluid is on the left of the wall

\(- \Rightarrow\) fluid is on the right of the wall

**Horizontal wall**

\(+ \Rightarrow\) fluid is on the left of the wall

\(- \Rightarrow\) fluid is on the right of the wall

#### 3.1.1.2.6.3. Boundary Conditions at Free Surface - Velocity¶

Condition is that the free surface has zero divergence of velocity

One side facing vacuum: three sides computed normally, fourth side follows from vanished divergence

Two sides facing vacuum: \(\partial u \over \partial x\) and \(\partial v \over \partial y\) must vanish separately - each vacuum dies velocity is set equal to the velocity of the cell across from it

Three sides facing vacuum: rare, but vacuum side opposite fluid side carries the same velocity of that side. The other two vacuum sides which oppose each other are calculated to follow freely the effects of the body force and do not otherwise change

Four sides facing vacuum: similar to three so that this isolated drop follows a free fall trajectory.

#### 3.1.1.2.6.4. Boundary Conditions at Free Surface - Pressure¶

Condition is that the normal stress is zero or of equating it to the applied external pressure

**Neglected the effect of viscous stress in the surface boundary and equated surface pressure directly to applied pressure**(to reduce complexity of simulation)May be valid for small viscosities

## 3.1.1.3. Numerical Stability and Accuracy¶

The method is stable - viscosity is not needed to ensure stability

“Artificial viscosity” term is not needed (which it is in the treatment of shocks in compressible flow and for the elimination of stagnation fluctuations)

Condition 1 for diffusional stability:

Condition 2 for incompressible flow analogy to Courant condition (g is gravity, d is maximum depth):

Harlow and Welsh also did internal consistency checks

## 3.1.1.4. Examples¶

### 3.1.1.4.1. Broken Dam Problem¶

Variant 1: Entire Dam Removed Instantly

Variant 2: Dam Partly Broken (like a sluice gate)

Comparison is good with experiment, but poor with hydrodynamic theory

Free surface is sensitive to grid density, but the effect is negligible

For the partly broken dam problem, the zero-pressure boundary condition applied at the centre of the cells leads to an inward dip of the zero pressure isobar at the earliest time.

## 3.1.1.5. Conclusions¶

Harlow and Walsh conclude by saying:

The technique can be extended to 3D

Other coordinate systems can be used

Other applications tested:

Rayleigh-Taylor instability of a falling free surface

Splash of a falling column of water

Breaking of a dam

Flow under a sluice gate

Flow over an underwater obstacle

The slow creep of a highly viscous blob

Splash of a linear drop

Calculations planned:

Formation of waves by a linear explosion over the surface

Breaking waves on a sloping beach

Modification to cylindrical coordinates allows studies of spherical explosion and the splash of a drop