Created in COMSOL Multiphysics 5.

Two-Phase Flow Modeling of a Dense Suspension

This model is licensed under the COMSOL Software License Agreement 5.4.
All trademarks are the property of their respective owners. See www.comsol.com/trademarks.
Introduction
Liquid-solid mixtures (suspensions) are important in a variety of industrial fields, such as
oil and gas refinement, paper manufacturing, food processing, slurry transport, and
wastewater treatment. Several different modeling approaches have been developed,
ranging from discrete, particle-based methods to macroscopic, semi-empirical two-phase
descriptions. Particle-based methods are suitable when there is a limited number of solid
particles. When, on the other hand, there are many particles, it is better to use a
macroscopic, or averaged, model that tracks the volume fractions of the phases.

The following example illustrates how you can set up a macroscopic two-phase flow model
in COMSOL Multiphysics using the Mixture Model, Laminar Flow interface. The model
is based on the “diffusive flux” model described in Ref. 1, Ref. 2, and Ref. 3, suitable for
liquid-solid mixtures with high concentrations of solid particles. It accounts for not only
buoyancy effects but also shear-induced migration; that is, the tendency of particles to
migrate toward regions of lower shear rates.

The model simulates the flow of a dense suspension consisting of light, solid particles in a
liquid placed between two concentric cylinders. The inner cylinder rotates while the outer
one is fixed.

Model Definition
A suspension is a mixture of solid particles and a liquid. The dynamics of a suspension can
be modeled by a momentum transport equation for the mixture, a continuity equation,
and a transport equation for the solid phase volume fraction. The Mixture Model, Laminar
Flow interface automatically sets up these equations. It uses the following equation to
model the momentum transport:

∂j T
ρ + ρ ( j ⋅ ∇ )j + ρ c ε ( j slip ⋅ ∇ )j = – ∇ p – ∇ ⋅ μ [ ∇j + ∇j ] + ρg
∂t
T
– ∇ ⋅ [ ρ c ( 1 – φ c ε )u slip j slip ] – ρ c ε ( j ⋅ ∇ )j slip

where j is the volume-averaged mixture velocity (m/s), p denotes the pressure (Pa), g
refers to the acceleration of gravity (m/s2), ε is the reduced density difference, uslip gives
the relative velocity between the solid and the liquid phases (m/s), and jslip denotes the
slip flux (m/s). Further, ρ = ( 1 – φ s )ρ f + φ s ρ s is the mixture density, where ρf and ρs are
the pure-phase densities (kg/m3) of liquid and solids, respectively, and φ s is the solid-
phase volume fraction (m3/m3). Finally, μ represents the mixture viscosity (Ns/m2)
according to the Krieger-type expression

2 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

φ s – 2,5φ max
μ = μ f  1 – ------------ (1)
 φ max

where μf is the dynamic viscosity of the pure fluid and φ max is the maximum packing
concentration.

The mixture model uses the following form of the continuity equation

∇⋅j = 0 (2)

The transport equation for the solid-phase volume fraction is

∂φ s
+ ∇ ⋅ ( φs us ) = 0 (3)
∂t

The solid-phase velocity, us, is given by us = u + ( 1 − cs ) uslip, where cs = φsρs/ρ is the

dimensionless particle mass fraction and u is the mass-averaged mixture velocity (m/s).
Consequently, Equation 3 is equivalent to

∂φ s
+ ∇ ⋅ ( φ s u + φ s ( 1 – c s )u slip ) = 0 (4)
∂t

Rao and others (Ref. 2) formulate the continuity equation and the particle transport in a
slightly different way. Instead of the slip velocity, uslip, they define a particle flux, Js (kg
/(m2·s)), and write the the solid phase transport according to

∂φ s ∇⋅J
-------- + ∇ ⋅ ( φ s u ) = – --------------s- (5)
∂t ρs

By comparing Equation 5 with Equation 4, it is clear that they are equivalent if

Js
u slip = -------------------------------
φs ρs ( 1 – c s )

In this example you use the model for the particle flux, Js, as suggested by Subia and
others (Ref. 3) and Rao and others (Ref. 2), but the open and editable format of
COMSOL Multiphysics makes it possible to specify the expression arbitrarily.

Following Rao and others, the particle flux is

Js
------ = – [ φD φ ∇ ( γ· φ ) + φ γ· D η ∇ ( ln μ ) ] + f h u st φ
2
ρs

3 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

Here, ust is the settling velocity (m/s) of a single particle surrounded by fluid and Dφ and
Dη are empirically fitted parameters (m2) given by

2
D φ = 0,41a
2
D η = 0,62a

where a is the particle radius (m).

·
The shear rate tensor, γ (1/s), is given by
˜
· T
γ = ∇j + ∇j
˜
and its magnitude by

· --- ( γ· :γ· )
1
γ =
2 ˜ ˜

which for a 2-dimensional problem is

· 1 2 2 2
--- ( 4u x + 2 ( u y + v x ) + 4v y )
γ =
2

The settling velocity, ust, for a single spherical particle surrounded by pure fluid is given by

2
2 a ( ρs – ρf )
u st = --- ---------------------------- g
9 η0

For several particles in a fluid, the settling velocity is lower. To account for the surrounding
particles, the settling velocity for a single particle is multiplied by the hindering function,
fh, defined as

μ f ( 1 – φ av )
f h = ----------------------------
μ

where φ av is the average solid phase volume fraction in the suspension, μf is the dynamic
viscosity of the pure fluid (Ns/m2), and μ is the mixture viscosity (Equation 1).

The following table gives the physical properties of the solid and the liquid phases:

NAME VALUE DESCRIPTION

3
ρs 1180 kg/m Density of particles
ρf 1250 kg/m3 Density of pure fluid

4 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

NAME VALUE DESCRIPTION

a 678 μm Particle radius

μf 0.589 Pa·s Viscosity of pure fluid

BOUNDARY CONDITIONS
The suspension is placed in a Couette device, that is, between two concentric cylinders.
The inner cylinder rotates while the outer one is fixed. The radii of the two cylinders are
0.64 cm and 2.54 cm, respectively. The inner cylinder rotates at a steady rate of 55 rpm.
With the cylinder centered at (0,0), this corresponds to a velocity of

110π
( u, v ) = ------------- ( y, – x )
60

The fluid and particle motion is small along the direction of the cylinder axis. You can
therefore use a 2-dimensional model. Figure 1 shows the corresponding geometry.

Figure 1: Geometry of the Couette device. The inner cylinder rotates, the outer one is fixed.

There is no particle flux through the boundaries, and the suspension velocity satisfies no
slip conditions at all walls.

INITIAL CONDITIONS
There are two different initial particle distributions. In the first example, the particles are
evenly distributed within the device. In the second example, the particles are initially
gathered at the top of the device.

5 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

Results

CASE 1—INITIALLY EVENLY DISTRIBUTED PARTICLES

A suspension with particles lighter than the fluid is placed in a concentric Couette device.
Initially, the particles are evenly distributed with a constant volume fraction of 0.35. The
shear rate in the device varies radially across the gap and thus it is expected that particles
migrate (shear-induced migration) from regions of high shear to regions of low shear
(toward the outer wall). Because the particles are lighter than the fluid, they also rise.

Figure 2: The particle concentration φ s at different times. The particles move to regions with
lower shear rate and rise because of buoyancy.

Figure 2 shows the particle concentration φ s in the device at t = 0 s, t = 30 s, t = 100 s and

t = 1000 s. The migration of the particles toward the outer wall is apparent. As a result of
the shear induced migration and gravity, the solid phase volume fraction approaches the
value for maximum packing close to the upper right outer wall. The suspension viscosity
thus becomes high in this region. The results compare well with those presented in Ref. 2.

6 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

CASE 2—PARTICLES INITIALLY GATHERED AT THE TOP OF THE DEVICE
In this case the particles are initially gathered at the top of the device. The particle volume
fraction is initially zero in the lower part, while it is 0.59 at the top.

Figure 3: Particle concentrations for t = 0 s, 10 s, 20 s, and100 s with particles initially at the

top. Note that the same color range values are used in each of the plots.

Figure 3 shows the numerically predicted particle concentration at times 0 s, 10 s, 20 s,

and 100 s. Initially, the particle motion is dominated by inertia and the effect of the shear-
induced migration is not visible. At later times, shear-induced migration causes the
particles to move toward the outer boundary. In this case also, the results agree well with
the results in Ref. 2.

References
1. R.J. Phillips, R.C. Armstrong, R.A. Brown, A.L. Graham, and J.R. Abbot, “A
Constitutive Equation for Concentrated Suspensions that Accounts for Shear-induced
Particle Migration,” Phys. Fluids A, vol. 4, pp. 30–40, 1992.

7 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

2. R. Rao, L. Mondy, A. Sun, and S. Altobelli, “A Numerical and Experimental Study of
Batch Sedimentation and Viscous Resuspension,” Int. J. Num. Methods in Fluids,
vol. 39, pp. 465–483, 2002.

3. S.R. Subia, M.S. Ingber, L.A. Mondy, S.A. Altobelli, and A.L. Graham, “Modelling of
Concentrated Suspensions Using a Continuum Constitutive Equation,” J. Fluid Mech.,
vol. 373, pp. 193–219, 1998.

Notes About the COMSOL Implementation

To set up the model with COMSOL Multiphysics, open the Mixture Model, Laminar
Flow interface. The shear rate is discretized as an additional equation to improve accuracy
because the particle flux contains derivatives of this quantity, which in turn depend on the
derivatives of the velocity.

Application Library path: CFD_Module/Multiphase_Benchmarks/

dense_suspension

Modeling Instructions
From the File menu, choose New.

NEW
In the New window, click Model Wizard.

MODEL WIZARD
1 In the Model Wizard window, click 2D.
2 In the Select Physics tree, select Fluid Flow>Multiphase Flow>Mixture Model>
Mixture Model, Laminar Flow (mm).
3 Click Add.
4 In the Select Physics tree, select Mathematics>PDE Interfaces>General Form PDE (g).
5 Click Add.
6 In the Dependent variables table, enter the following settings:

gamma

7 Click Study.

8 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

8 In the Select Study tree, select General Studies>Time Dependent.
9 Click Done.

GLOBAL DEFINITIONS
1 In the Model Builder window, under Global Definitions click Parameters 1.
2 In the Settings window for Parameters, locate the Parameters section.
3 Click Load from File.
4 Browse to the model’s Application Libraries folder and double-click the file
dense_suspension_parameters.txt.

GEOMETRY 1

Circle 1 (c1)
1 In the Geometry toolbar, click Primitives and choose Circle.
2 In the Settings window for Circle, locate the Size and Shape section.
3 In the Radius text field, type 0.0064.
4 Right-click Circle 1 (c1) and choose Build Selected.

Circle 2 (c2)
1 In the Geometry toolbar, click Primitives and choose Circle.
2 In the Settings window for Circle, locate the Size and Shape section.
3 In the Radius text field, type 0.0254.
4 Right-click Circle 2 (c2) and choose Build Selected.

Compose 1 (co1)
1 In the Geometry toolbar, click Booleans and Partitions and choose Compose.
2 Click in the Graphics window and then press Ctrl+A to select both objects.
3 In the Settings window for Compose, locate the Compose section.
4 In the Set formula text field, type c2-c1.
5 Right-click Compose 1 (co1) and choose Build Selected.
6 Click the Zoom Extents button in the Graphics toolbar.

DEFINITIONS

Variables 1
1 In the Home toolbar, click Variables and choose Local Variables.
2 In the Settings window for Variables, locate the Variables section.

9 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

3 Click Load from File.
4 Browse to the model’s Application Libraries folder and double-click the file
dense_suspension_variables.txt.

MIXTURE MODEL, LAMINAR FLOW (MM)

Now, define the Slip model, the Mixture Properties, and the Initial Values.

1 In the Model Builder window, under Component 1 (comp1) click Mixture Model,
Laminar Flow (mm).
2 In the Settings window for Mixture Model, Laminar Flow, locate the Physical Model
section.
3 From the Slip model list, choose User defined.

Mixture Properties 1
1 In the Model Builder window, under Component 1 (comp1)>Mixture Model,
Laminar Flow (mm) click Mixture Properties 1.
2 In the Settings window for Mixture Properties, locate the Continuous Phase Properties
section.
3 From the ρc list, choose User defined. In the associated text field, type rho_f.
4 From the μc list, choose User defined. In the associated text field, type eta_f.
5 Locate the Dispersed Phase Properties section. From the ρd list, choose User defined. In
the associated text field, type rho_s.
6 Locate the Mixture Model section. Specify the uslip vector as

J1/(phid*rho_s*(1-mm.cd)) x
J2/(phid*rho_s*(1-mm.cd)) y

7 In the φmax text field, type phi_max.

Initial Values 1
1 In the Model Builder window, under Component 1 (comp1)>Mixture Model,
Laminar Flow (mm) click Initial Values 1.
2 In the Settings window for Initial Values, locate the Initial Values section.
3 In the φd text field, type phi0.

Gravity 1
1 In the Physics toolbar, click Domains and choose Gravity.
2 In the Settings window for Gravity, locate the Domain Selection section.

10 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

3 From the Selection list, choose All domains.

Wall 2
1 In the Physics toolbar, click Boundaries and choose Wall.
2 Select Boundaries 3, 4, 6, and 7 only.
3 In the Settings window for Wall, click to expand the Wall Movement section.
4 From the Translational velocity list, choose Manual.
5 Specify the utr vector as

c_vel*y x
-c_vel*x y

For the continuous phase, the default boundary condition, No slip, is correct for the
remaining boundaries. For the dispersed phase, the default condition, No dispersed
phase flux, is correct for all boundaries.

Pressure Point Constraint 1

1 In the Physics toolbar, click Points and choose Pressure Point Constraint.
2 Select Point 1 only.

GENERAL FORM PDE (G)

Specify an appropriate unit for shear rate variable gamma.

1 In the Model Builder window, under Component 1 (comp1) click General Form PDE (g).
2 In the Settings window for General Form PDE, locate the Units section.
3 Click Define Dependent Variable Unit.
4 In the Dependent variable quantity table, enter the following settings:

Dependent variable quantity Unit

Custom unit 1/s

5 In the Source term quantity table, enter the following settings:

Source term quantity Unit

Custom unit 1/s

General Form PDE 1

1 In the Model Builder window, under Component 1 (comp1)>General Form PDE (g) click
General Form PDE 1.

11 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

2 In the Settings window for General Form PDE, locate the Conservative Flux section.
3 Specify the Γ vector as

0 x
0 y

4 Locate the Source Term section. In the f text field, type gamma-sqrt(0.5*(4*jux^2+
2*(juy+jvx)^2+4*jvy^2)+eps).
5 Locate the Damping or Mass Coefficient section. In the da text field, type 0.
6 In the Model Builder window, click General Form PDE (g).

Flux/Source 1
1 In the Physics toolbar, click Boundaries and choose Flux/Source.
2 In the Settings window for Flux/Source, locate the Boundary Selection section.
3 From the Selection list, choose All boundaries.

MESH 1

Size
1 In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Free Triangular.
2 In the Settings window for Size, locate the Element Size section.
3 From the Calibrate for list, choose Fluid dynamics.
4 From the Predefined list, choose Finer.
5 Click the Custom button.
6 Locate the Element Size Parameters section. In the Maximum element size text field, type
0.0012.

7 Click Build All.

STUDY 1

Step 1: Time Dependent

1 In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Times text field, type 0 30 100 1000.
4 In the Study toolbar, click Show Default Solver.

12 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

Solution 1 (sol1)
1 In the Model Builder window, expand the Solution 1 (sol1) node.
2 In the Model Builder window, under Study 1>Solver Configurations>Solution 1 (sol1) click
Time-Dependent Solver 1.
3 In the Settings window for Time-Dependent Solver, click to expand the Time Stepping
section.
4 In the Study toolbar, click Compute.

RESULTS

Mixture (mm)
To visualize the volume fraction of the Dispersed Phase as a surface plot along with the
mixture velocity field as an arrow surface plot, follow the steps given below.

Dispersed Phase Volume Fraction

1 In the Model Builder window, expand the Dispersed Phase (mm) node, then click
Dispersed Phase Volume Fraction.
2 In the Settings window for Surface, click to expand the Range section.
3 Select the Manual color range check box.
Fix the color range so that the solutions for different time values use the same color
legend. The following values cover both cases.
4 In the Minimum text field, type 0.
5 In the Maximum text field, type 0.45.

Arrow Surface 1
1 In the Model Builder window, under Results right-click Dispersed Phase (mm) and choose
Arrow Surface.
2 In the Dispersed Phase (mm) toolbar, click Plot.

Dispersed Phase (mm)

1 In the Model Builder window, under Results click Dispersed Phase (mm).
2 In the Settings window for 2D Plot Group, locate the Data section.
3 From the Time (s) list, choose 0.
4 In the Dispersed Phase (mm) toolbar, click Plot.
5 From the Time (s) list, choose 30.
6 In the Dispersed Phase (mm) toolbar, click Plot.
7 From the Time (s) list, choose 100.

13 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

8 In the Dispersed Phase (mm) toolbar, click Plot.
9 From the Time (s) list, choose 1000.
10 In the Dispersed Phase (mm) toolbar, click Plot.
This completes Case 1. Now model the case with the particles initially gathered at the
top.

DEFINITIONS
Define the particle concentration distribution using the Step function.

Step 1 (step1)
1 In the Home toolbar, click Functions and choose Local>Step.
2 In the Settings window for Step, click to expand the Smoothing section.
3 In the Size of transition zone text field, type 2*2.

Variables 1
1 In the Model Builder window, under Component 1 (comp1)>Definitions click Variables 1.
2 In the Settings window for Variables, locate the Variables section.
3 In the table, enter the following settings:

Name Expression Unit Description

phi0_2 step1(y[1/mm]-8)*0.59 Initial concentration

MIXTURE MODEL, LAMINAR FLOW (MM)

In the Physics toolbar, click General Form PDE (g) and choose Mixture Model,
Laminar Flow (mm).

In the Model Builder window, under Component 1 (comp1) click Mixture Model,
Laminar Flow (mm).

Initial Values 2
1 In the Physics toolbar, click Domains and choose Initial Values.
2 In the Settings window for Initial Values, locate the Domain Selection section.
3 From the Selection list, choose All domains.
4 Locate the Initial Values section. In the φd text field, type phi0_2.
In order to avoid dividing by zero in the region where the particle volume fraction is
initially zero, modify the expressions for the slip velocity.

14 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

Mixture Properties 1
1 In the Model Builder window, under Component 1 (comp1)>Mixture Model,
Laminar Flow (mm) click Mixture Properties 1.
2 In the Settings window for Mixture Properties, locate the Mixture Model section.
3 Specify the uslip vector as

J1/(phid*rho_s*(1-mm.cd)+eps) x
J2/(phid*rho_s*(1-mm.cd)+eps) y

COMPONENT 1 (COMP1)
Create a finer mesh compared to the one used in the previous case.

1 In the Mesh toolbar, click Add Mesh.

MESH 2

Size
1 In the Model Builder window, under Component 1 (comp1)>Meshes right-click Mesh 2 and
choose Free Triangular.
2 In the Settings window for Size, locate the Element Size section.
3 From the Predefined list, choose Extra fine.
4 Click the Custom button.
5 Locate the Element Size Parameters section. In the Maximum element size text field, type
0.0006.

6 Click Build All.

Next, add a Time Dependent study to compute the particle distribution.

ADD STUDY
1 In the Home toolbar, click Add Study to open the Add Study window.
2 Go to the Add Study window.
3 Find the Studies subsection. In the Select Study tree, select General Studies>
Time Dependent.
4 Click Add Study in the window toolbar.
5 In the Home toolbar, click Add Study to close the Add Study window.

15 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

STUDY 2

Step 1: Time Dependent

1 In the Settings window for Time Dependent, locate the Study Settings section.
2 In the Times text field, type range(0,10,100).
3 Click to expand the Mesh Selection section. In the table, enter the following settings:

Geometry Mesh
Geometry 1 Mesh 2

Use manual Scaling of the variables to improve the convergence.

4 In the Study toolbar, click Show Default Solver.

Solution 2 (sol2)
1 In the Model Builder window, expand the Solution 2 (sol2) node.
2 In the Model Builder window, expand the Study 2>Solver Configurations>
Solution 2 (sol2)>Dependent Variables 1 node, then click Velocity field, mixture (comp1.j).
3 In the Settings window for Field, locate the Scaling section.
4 From the Method list, choose Manual.
5 In the Scale text field, type 0.1.
6 In the Model Builder window, under Study 2>Solver Configurations>Solution 2 (sol2)>
Dependent Variables 1 click Pressure (comp1.p).
7 In the Settings window for Field, locate the Scaling section.
8 From the Method list, choose Manual.
9 In the Scale text field, type 100.
The default absolute tolerance for phid is set assuming that phid is nowhere near the
packing limit and hence, the tolerance becomes too stringent in this case.
10 In the Model Builder window, under Study 2>Solver Configurations>Solution 2 (sol2) click
Time-Dependent Solver 1.
11 In the Settings window for Time-Dependent Solver, click to expand the Absolute Tolerance
section.
12 In the Variables list, select Volume fraction, dispersed phase (comp1.phid).
13 In the Tolerance text field, type 1e-4.
To ensure that the volume fraction of particles has a positive value, do the following:
14 In the Model Builder window, expand the Study 2>Solver Configurations>
Solution 2 (sol2)>Time-Dependent Solver 1>Segregated 1 node.

16 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

15 Right-click Study 2>Solver Configurations>Solution 2 (sol2)>Time-Dependent Solver 1>
Segregated 1 and choose Lower Limit.
16 In the Settings window for Lower Limit, locate the Lower Limit section.
17 In the Lower limits (field variables) text field, type comp1.phid 0.
18 In the Study toolbar, click Compute.

RESULTS
To visualize the volume fraction of the dispersed phase as a surface plot along with the
arrow plot of the mixture velocity, follow the steps given below.

Dispersed Phase (mm) 1

1 In the Model Builder window, under Results click Dispersed Phase (mm) 1.
2 In the Settings window for 2D Plot Group, locate the Data section.
3 From the Time (s) list, choose 0.

Dispersed Phase Volume Fraction

1 In the Model Builder window, expand the Dispersed Phase (mm) 1 node, then click
Dispersed Phase Volume Fraction.
2 In the Settings window for Surface, locate the Range section.
3 Select the Manual color range check box.
4 In the Maximum text field, type 0.67.

Arrow Surface 1
1 In the Model Builder window, under Results right-click Dispersed Phase (mm) 1 and
choose Arrow Surface.
2 In the Dispersed Phase (mm) 1 toolbar, click Plot.

Dispersed Phase (mm) 1

1 In the Model Builder window, under Results click Dispersed Phase (mm) 1.
2 In the Settings window for 2D Plot Group, locate the Data section.
3 From the Time (s) list, choose 10.
4 In the Dispersed Phase (mm) 1 toolbar, click Plot.
5 From the Time (s) list, choose 20.
6 In the Dispersed Phase (mm) 1 toolbar, click Plot.
7 From the Time (s) list, choose 100.
8 In the Dispersed Phase (mm) 1 toolbar, click Plot.

17 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION

18 | TWO-PHASE FLOW MODELING OF A DENSE SUSPENSION