Terminal Falling Velocity of a Sand Grain

The first stop for polluted water entering a water work is normally a large tank, where large particles are left to settle. More generally, gravity settling is an economical method of separating particles. If the fluid in the tank is moving at a controlled low velocity, the particles can be sorted in separate containers according to the time it takes for them to reach the bottom.

This application simulates a spherical sand grain falling in water. The grain accelerates from standstill and rapidly reaches its terminal velocity. The results agree with experimental studies. The model is an axially symmetric fluid-flow simulation in a moving coordinate system, coupled to an ordinary differential equation (ODE) describing the grain’s motion.

Model Definition

The model couples the flow simulation in cylindrical coordinates with an ODE for the force balance of the particle. Due to axial symmetry, you can model the flow in 2D instead of 3D. The figure below shows the modeling domain.

Domain Equations

The fluid flow is described by the Navier–Stokes equations

where ρ denotes the density (kg/m3), u the velocity (m/s), η the viscosity (Ns/m2), and p the pressure (Pa). The fluid is water with a viscosity of 1.51·10−3 Ns/m2 and a density of 1000 kg/m3. The model uses the accelerating reference system of the sand grain. This means that the volume force density F is given by:

where a (m/s2) is the acceleration of the grain and g = 9.81 m/s2 is the acceleration due to gravity. The ODE that describes the force balance is:

where m (kg) denotes the mass of the particle, x (m) the position of the particle, Fg (N) the gravitational force, and Fz the z-component of the force that the water exerts on the sand grain. The gravitational force is given by:

where Vgrain (m3) is the volume of the sand grain and ρgrain (kg/m3) its density. The force that the water exerts on the grain is calculated by integrating the normal component of the stress tensor over the surface of the particle. Because the model is axially symmetric, only the force’s z-component is nonzero:

where r (m) is the radial coordinate and n is the normal vector on the surface of the grain.

The initial values for position and velocities are

Boundary Conditions

At the sphere’s surface, the fluid velocity relative the sphere is zero, that is u = 0 — a situation described by the no slip wall condition. At the inlet of the fluid domain the velocity equals the falling velocity:

. At the outer boundary of the water domain, the normal velocity and the tangential shear stress both vanish, which means that a symmetry condition applies. Furthermore, a neutral condition, n · [−pI + η(∇u + (∇u)T)] = 0, describes the outlet, and an axial symmetry condition models the symmetry axis at r = 0.

Results

Figure 1 shows the velocity field at the final simulation time, t = 1 s, when the particle has reached steady state.

Figure 1: The velocity field at steady state. Note that the velocities are plotted in the reference system of the sand grain.

The following series of snapshots (Figure 2) displays velocity field from a moment just after the sand grain is released until it is approaching steady state. Notice the recirculation forming above the grain.

Figure 2: The velocity field around the sand grain at a series of times (t = 0.025 s, 0.05 s, 0.075 s, 0.1 s, 0.15 s, 0.2 s, 0.3 s, 0.5 s, and 0.75 s). For a color legend, see Figure 1.

Figure 3 shows the falling velocity of the grain as a function of time.

Figure 3: Falling velocity (m/s) of the grain versus time. After the solution time of 1 s, the velocity approaches the terminal velocity.

The terminal velocity equals 0.288 m/s. When this state is reached, the gravity and the forces from the water cancel out. Figure 4 shows the forces on the sand grain.

Figure 4: The forces on the sand grain. The force that the water exerts on the sphere (blue line) increases as the grain gains speed. The gravity force (red line) remains the same, and the total force (black line) tends toward zero as the solution approaches steady state.

Several approximate equations have been proposed for the terminal velocity of a sphere falling in a fluid. Ref. 1 cites the following expression for the total force that the fluid exerts on the sphere, as a function of the velocity:

Here d (m) is the diameter of the sphere, ρ (kg/m3) is the fluid density, v (m/s) is the velocity, and Re = (ρvd)/μ is the Reynolds number, with μ (Ns/m2) being the viscosity of the fluid. The gravity force is given analytically as Fg = πd3(ρ − ρs)g/6, where ρs (kg/m3) is the density of the sphere. Equating the two forces and introducing the values used in the simulation gives an approximate terminal velocity of 0.284 m/s.

The same reference discusses correction factors for nonspherical particles. You can easily adapt the model to hold for a general axially symmetric object (by redrawing the geometry) or even an arbitrarily shaped object (by modeling in 3D).

Reference

1. J.M. Coulson and J.F. Richardson, Chemical Engineering vol. 2, 4th ed., 1993.

COMSOL_Multiphysics/Fluid_Dynamics/falling_sand

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

A set of global parameters is provided in a text file.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file falling_sand_parameters.txt.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 6e-3.

In the text field, type 14e-3.

Locate the section. In the text field, type -6e-3.

Click

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 1e-3.

Click

Difference 1 (dif1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Find the subsection. Click to select the

toggle button.

Select the object only.

Click

This completes the model geometry.

Definitions

Define a nonlocal integration coupling and a variable using this coupling. Later, you will use this variable to calculate the drag force.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundaries 6 and 7 only.

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Description
F_z	intop1(-w_lm[Pa])	Drag force

Here, w_lm is a Lagrange multiplier and represents the reaction flux when imposing the z-component of the fluid velocity at the wall by means of weak constraints.

Laminar Flow (spf)

Fluid Properties 1

In the window, under click .

In the window for , locate the section.

From the ρ list, choose . In the associated text field, type rho_water.

From the μ list, choose . In the associated text field, type mu_water.

Volume Force 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. Specify the F vector as


0	r
-rho_water*(Xdott+g_const)	z

Next, use weak constraints to impose the wall boundary condition on the grain. Add a node for the equation of motion for the falling grain, and constraint the Lagrange multiplier corresponding to the horizontal component of the flux at the extreme points of the wall so that it respects the symmetry of the problem. By default, these feature are not visible.

Click the