The effect of lubrication forces on the collision statistics of cloud droplets in homogeneous isotropic turbulence
Institute of Meteorology and Water Management, National Research Institute (IMGW–PIB, Warsaw)
Institute of Geophysics, University of Warsaw (IGF UW, Warsaw)
April 30th, 2021
Overview
Collision-coalescence of inertial particles in turbulent flows
Motivation of study:
To compute kinematics and collision statistics of water droplets typical of atmospheric clouds
Focus of study:
To investigate the effect of lubrication forces on the dynamics of aerodynamically interacting droplets
• Turbulent flows and particle transport are very common phenomena in nature. The processes occur continuously, with different intensity and at different scales.
• Particle transport by a turbulent fluid phase is also an important mechanism for many technological processes (combustion of pulverized coal in boilers, pneumatic transport in pipelines, combustion of the fuel in engines or spraying of fertilizers and plant protection products)
• In this study, the general focus is on modeling of cloud microphysical processes.
“Warm rain processes account for about 31% of the total rain fall and 72% of the total rain area in tropics. This process is active in most climate zones in all seasons.”
Motion of particles in homogeneous isotropic turbulence
Turbulent background flow U(x, t)
Particle motion Y(k) (t), V(k) (t)
Contours of vorticity (Color: magnitude)
Motion of particles (Cone: velocity magnitude
and direction)
w v u w v u
i j k
y z z x x y
=
= − + − + −
ω U
Problem definition
Simulation of the homogeneous isotropic turbulence (HIT)
• Governing equations
• Discretization (Eulerian)
3D Cartesian mesh of N equally spaced grid points in each spatial direction
• Basic assumptions
3D incompressible homogeneous isotropic turbulent flow with periodic boundary conditions (also for droplets) on a cube of size 2π
• Solution Method
Pseudo-spectral method for direct numerical simulation (DNS)
The momentum equation is solved (integrated) in the spectral space using discrete Fourier transform as
( ) ( )
( )( )
( ) ( )
1 2 3
1 2 3
ˆ ˆ
, , ,
where 0, 1,
s 2,...,
2 ˆ ˆ
and , , .Theinverse Fourier transformi
i k x k y k z i
k k k
i
t t e t e
k N
t i t
+ +
=
=
=
k xk
U x U k U k
ω k k U k
( ) ( )
( 1 2 3 )
3
1 1 1
3
0 0 0
2 2 2
ˆ , 1 ,
1 2 2 2
, , ,
i
N N N
j m n
j m n
i jk N mk N nk N
t t e
N
j m n
x y z t
N N N N
e
−
− − −
= = =
− + +
=
= = = =
k x x
U k U x
U
2 2
1 ( , )
2 0
P t
t
= − + + +
=
U U ω U U f x
U
where vorticity
Algorithm for simulating the HIT
Solution
Assuming we have the solution in the following form U kˆ
(
,dt) (
; U kˆ , 2dt)
; ... ;U kˆ(
, ndt)
⎯⎯⎯⎯⎯To compute→U kˆ(
, n 1(
+)
dt)
Step 1. Apply FFT-1 to obtain velocity in physical space
Algorithm
( )
FFT-1( )
ˆ ,t ⎯⎯⎯→ ,t
U k U x
Step 2. Obtain vorticity in spectral space and then apply
FFT-1 to obtain vorticity in physical space ω kˆ
( )
,t =ik × U kˆ( )
,t ⎯⎯⎯→FFT-1 ω x( )
,tStep 3. Calculate N1(x, t) ≡ U × ω (the first nonlinear term RHS of N–S) in physical space and FFT to spectral space
(
, ndt)
(
, ndt)
⎯⎯⎯FFT→ ˆ1(
, ndt)
U x ω x N k
Step 4. Evolve the velocity in time Note: N2(x, t) ≡ ∇(P/ρ + U2/2)
( )
( ) ( ) ( ) ( ( ) )
( )
( )
( ) ( ( ) )
( )
( )
2 2
3 1
ˆ , n 1 ˆ , n ˆ , n ˆ , n-1
2 2
ˆ , n 1
1 ˆ , n ˆ , n 1
2 ˆ , n
dt dt dt dt dt dt
i dt dt
k dt dt dt
dt dt
+ − = −
− +
− + +
+
1 1
U k U k N k N k
k N k
U k U k f k
( )
1( )
2 2
ˆ ,
in which ˆ , i t ; Because 0
t k
= − k N k =
N k U
Adams–Bashforth (AB2)
Crank–Nicolson
Euler → Deterministic
Problem definition
Evolution of the dispersed phase and aerodynamic interactions
( ) ( ) ( )
( )
( )
( )
d ( ) ( ) ( ( ), )
d d ( )
d ( )
k k k
k p k
k
t t t t
t
t t
t
= − − +
=
V V U Y
g
Y V
For k-th particle: V(k) (t) is the particle velocity, τp(k)is its Stokes inertial response time, and U(Y(k) (t), t) is the undisturbed fluid velocity, U(x, t), at the location of the particle: Y(k)(t)
• Governing equations (typical):
where u(k)is the disturbance velocity felt at the location of k-th particle, and U(k)= U(Y(k)(t), t).
• Aerodynamic interactions (AIs):
Two droplets sedimenting under gravity in still air
g ↓
The motion of each droplet is generating a perturbation field that induces a disturbance velocity at the location of the other droplet
( )
( ) ( ) ( )
( )
( )
( )
( )
d ( ) ( ) d d ( )
d ( )
k k k
k
k p k
k
t t t
t t
t
− +
= − +
=
V U u
V g
Y V
• Equations of motion considering AIs:
+ AIs
Limitation: using this set of equations, the effect of interaction among particles would be overlooked
(PR: Wang, 2017)
How to determine?
x1 x2 x3
a1
a2
V2 V1
Y1(x1, x2, x3) Y2(x1, x2, x3)
Perturbation field for: Mathematical description Position vectors*
A single droplet r = x – Y
Two droplets using
(original) superposition method
r1 = x – Y1 r2 = x – Y2 Two droplets using
improved superposition method (ISM)
( ) St( ) 3 2 ( ) 3 ( )
3 3 3 1
, ; ,
4 4 4 4
a a a a
t a
r r r r r
= − + + = +
u x = u r V r r V V r r V V
( ),t St( 1; ,a1 1 − St( 1− 2; ,a2 2)) St( 2; ,a2 2 − St( 2 − 1; ,a1 1))
u x = u r V u Y Y V + u r V u Y Y V
( ) ( ) ( )
( ) ( )
St 1 1 1 St 1 2 2 2 2 St 2 2 2 St 2 1 1 1 1
1 2
St 1 1 1 1 St 2 2 2 2
, ; , ; , ; , ; ,
; , ; ,
t a a a a
a a
− − − + − − −
= − + −
u x = u r V u Y Y V u u r V u Y Y V u
u u
u r V u u r V u
Mathematical description of disturbance fields induced by particles using the solution to Stokes equation
* The position relative to the center of the particle/droplet
x1 x2 x3
Y (x1, x2, x3)
V x
r
Y
arbitrary:
x (x1, x2, x3) A single droplet Two droplets
n droplets
2 droplets N-body ISM?
Aerodynamic interactions
Extension of ISM to an arbitrary number of droplets (Np) and hybrid DNS (HDNS)
Number of droplets Disturbance velocity (perturbation felt at the location of droplet) Flow 2
Still air Np
Np Turbulent flow field
( )
( )
1 St 1 2 2 2 2
2 St 2 1 1 1 1
; ,
; , a a
= − −
= − −
u u Y Y V u
u u Y Y V u
( )k mNp1 St
(
(k m, )*; ( )m , ( )m ( )m)
, ** 1,2, , pm k
a k N
=
= − =
u u r V u
( )k mNp1 St
(
(k m, ); ( )m , ( )m(
( )m ( )m) )
, 1,2, , pm k
a k N
=
= − + =
u u r V U u
Using ISM to compute disturbance velocities
*r(k, m)= Y(k)– Y(m)
**For an arbitrary droplet k, all the other droplets are indexed m
( ) ( )
(
( ) ( ) ( ))
( ) ( )
(
( ) ( ))
( ) ( )( ) ( )
( ) ( )
, , ,
St St
, , , ,
, ,
,
; ,
k m k m k m m m
k m k m k m m k m m
k m m i j j
k m m
a
v
= +
=
u u r v
r r v v
α v
( ) ( )
( )
( )
( )
( )
( ) ( )
( )
( )
( )
3
, , 2
, ,
3 ,
, ,
3 3
4 4
3 1
4 4
m m
k m k m
k m k m
m m
k m
k m k m
a a
r r r
a a
r r
−
+
Now we define: then:
In the last equation let’s move the unknowns, u(m), to the LHS: ( ) p St
(
( , ) ( ) ( ))
p St(
( , ) ( ) ( ) ( ))
p1 1
; , ; , , 1, 2, , .
N N
k k m m m k m m m m
m m
m k m k
a a k N
= =
+ = − =
u u r u u r V U
(,k m, )
( )
(k m, ) (,k m, ) ( )i j r ri j ij j i m k
= + =
such that α(k, m) is a 3×3 symmetric matrix with components:
( , ) ( , ) ( ) St
k m k m m
u = α v
( ) p St
(
( , ) ( ) ( ))
p St(
( , ) ( ) ( ) ( ))
p1 1
; , ; , , 1, 2, , .
N N
k k m m m k m m m m
m m
m k m k
a a k N
= =
+ = − =
u u r u u r V U
Superposed Stokes disturbance velocities of Np – 1 droplets, u(k), felt at the location of an arbitrary droplet k can be described as:
It was shown:
That is:
( ) ( ) ( ) ( ) ( ) ( ( ) ( )) ( )( ( ) ( )) ( )( ( ) ( ))
( ) ( ) ( ) ( ) ( ) ( )( ( ) ( )) ( ( ) ( )) ( )( ( ) ( ))
( ) ( ) ( ) ( ) ( ) ( )( ( ) ( )) ( )( ( ) ( )) ( ( ) ( ))
p p p p p
p p p p p
p p p p p p p
1, 1,
1 1,2 2 1 1 1,2 2 2
3 3
2, 2,
2,1 1 2 2,1 1 1 2 2
3 3
,1 1 ,2 2 ,1 1 1 ,2 2 2
3 3
N N N N N
N N N N N
N N N N N N N
+ + + = − + − + + −
+ + + = − + − + + −
+ + + = − + − + + −
I u α u α u 0 V U α V U α V U
α u I u α u α V U 0 V U α V U
α u α u I u α V U α V U 0 V U
where 3 3
1 0 0 0 0 0
0 1 0 and 0 0 0
0 0 1 0 0 0
= =
I 0
and each row consists of 3 equations along 3 spatial directions. Thus, this set of 3Npequations in a compact notation can be described as:
( ) ( )
( ) ( ) ( )
( ) ( )
( )
( )
( )
( ) ( )
( ) ( ) ( )
( ) ( )
( )
p p
p p
p p p p p
1, 1,
1, 1 1, 1
3 3
, ,
,1 , ,1 ,
,1 , ,1 ,
3 3
N N
m m
k N m k N
k k m k k m
N N m N N N m
=
I α α u 0 α α V
α α α u α α α
α α I u α α 0
A x B
( )
( ) ( )
( ) ( )
( ) p
p p
1
m m , so
N
N N
k m
−
− = +
−
U
V U A B I
V U
y
( )
( )( ) ( ) ( )
( )
( ) ( )
( )
d ( ) ( ) d
k k k k
k
k k
p k
t t
t
→
− +
= − + → →
→
V U u V
V g U
u
Eqs. of motion DNS
N-body ISM (?)
→ HDNS HDNS
Parallel GMRes solver
Aerodynamic interactions
Inaccurate force representation of ISM
Normalized magnitude of the drag force acting on two same-size spherical particles, whether approaching or receding, along their line of centers as a function of nondimensional gap between their surfaces.
Wang et al. (2005):
“While our improved formulations still perform better than the original formulation, all the formulations based on the superposition method fail to predict the lubrication effect.”
Focus of our study
Lubrication effect:
s – 2 → 0 then F → ∞ (singular) so strong dependency to Δt is
expected s – 2
a V V a
No rm al iz ed dr ag: F /6 πμ aV
Normalized gap s – 2
r
s = r / a Two orders of magnitude larger for gap 0.001a →
a1 = a2 (monodisperse) Two blue lines
Computing AIs from the analytical solutions of Jeffrey & Onishi (1984)
Ignore rotational motion
• No rotation considered in ISM (HDNS)
• Considering rotational motion JO solutions is very time-consuming
Aerodynamic interactions
Computing AIs from the analytical solutions of Jeffrey & Onishi (1984)
Coupled implementation of HDNS and analytical solutions of JO84
50a
3a
The interaction regions (spheres) around each droplet:
• red: particle of consideration;
• blue: particles with r < 3a to the particle of consideration;
• black: particles with distances 3a < r < 50a;
• grey: distant particles r > 50a.
For each particle look for an interacting nearby
particle: r < 50a
Compute distance Compute constants
for N-body ISM and fill the matrices
Finished checking
all particles? Compute forces
Compute X (s, λ) &
Y (s, λ) for defined values of s
Interpolate X (s, λ) &
Y (s, λ) based on the distance s
Solve the set of equations of N-Body
ISM
Yes No
Add short- and long- range interaction field to the equations
of motion:
u= uHDNS+ ulub r > 3a
Compute factors Pnpq, Vnpq, Qnpq
JO84
HDNS Common
Compute functions fk(λ)
r < 3a
Collision statistics
RDF, RRV, kinematic and dynamic collision kernels: gr(r/R), ⟨|wr(r/R)|⟩, ГK and ГD
( ) ( )
( )
( ) ( )
( )
( ) ( )
( )
pair
1 2
pair shell
p p box
1
pair
K 2
D c
p p box
2
; ;
1 1
2
; ;
2
1 1
2
r
n
k n
n r
r r
R a a a
n r R t V r R g r t
R N N V
r r
w t
R n r R t
R w r R g r R
N
N N V
=
= + =
=
−
−
=
= = =
=
−
V V r
→ Measure of clustering
gr = 1: uniform distribution gr > 1: clustering
Radial distribution function (RDF):
Radial relative velocity (RRV):
Kinematic collision kernel:
Dynamic collision kernel:
→ Volumetric rate of collision Collision sphere radius:
1. Introduce the particles
2. Evolve the system for 10Te (statistically stationary) 3. Start collecting samples from the following:
→ Rate of collision
Sensitivity of droplet collision statistics to time step size
Two-point collision statistics at contact computed using different time steps:
(a) normalized radial relative velocity and (b) radial distribution function.
saturate
lubrication better represented (prevents collision) low-inertia droplets (more critical)
Results
Sensitivity of droplet collision statistics to time step size
RRV computed using different time step sizes for the normalized separation distance in the range:
(a) 1 < r / R < 10 and (b) 1 < r / R < 1.5 with the set of droplets ap = 40 µm and Np = 50 000.
The black rectangle in (a) marks the region that is enlarged and shown in (b).
AI Not much change
No changes at longer separation distances
Sensitivity of droplet collision statistics to time step size
RDF computed using different time step sizes for the normalized separation distance in the range:
(a) 1 < r / R < 10 and (b) 1 < r / R < 2 with the set of droplets ap = 40 µm and Np = 50 000.
The black rectangle in (a) marks the region that is enlarged and shown in (b).
AI + relocation
Not much change
accuracy and simulation time because Δt↓ then evolution time↓ hence uncertainty↑
No changes at longer separation distances
Results
Sensitivity of droplet collision statistics to time step size
The effect of time step size on the dynamic collision kernel for
(a) the same number of droplets: Np = 50 000 and (b) the same mass loading: φ = 1.22×10−2, considering various sets of same-size droplets with different radii.
Almost no change in CK RRV↓ + RDF↑
Setting an optimal location for the matching point δ
Variations in the radial relative velocity for the sets of droplets
(a) ap = 40 µm; Np = 50 000, and (b) ap = 30 µm; Np = 118 500 when the lubrication forces are considered within different ranges of interaction δ.
(δ = 3)
50a
3a
HDNS handles many-body
interaction
JO solutions capture lubrication effects
Results
Setting an optimal location for the matching point δ
Variations in the radial relative velocity for the sets of droplets
(a) ap = 40 µm; Np = 50 000, and (b) ap = 30 µm; Np = 118 500 when the lubrication forces are considered within different ranges of interaction δ.
lubrication effect lubrication effect
losing many-body
interaction effects losing many-body
interaction effects
Setting an optimal location for the matching point δ
Variations in the radial distribution function for the sets of droplets
(c) ap = 40 µm; Np = 50 000, and (d) ap = 30 µm; Np = 118 500 when the lubrication forces are considered within different ranges of interaction δ.
lubrication effect lubrication effect
losing many-body interaction effects
losing many-body interaction effects
Results
Setting an optimal location for the matching point δ
Variations in the (a) at-contact radial relative velocity and (b) at-contact radial distribution function for three sets of droplets with the same mass loading, φ = 1.22×10−2, when the lubrication forces are considered within different ranges of
interaction.
HDNS HDNS
losing the accuracy and efficiency of many-
body interactions
Setting an optimal location for the matching point δ
Variations in the dynamic collision kernel for three sets of droplets with the same mass loading, φ = 1.22×10−2, when the lubrication forces are considered within different ranges of interaction.
HDNS
less change in DCK
Results
Effects of AIs on kinematic and dynamic collision statistics
The effect of turbulent energy dissipation rate on the at-contact RRV for different (a) droplet sizes and (b) Stokes numbers.
In panel (b) the theoretical model of Saffman & Turner (1956) is shown.
The DNS results of Rosa et al. (2013, figures 8a and 13a there) are additionally included in panels (b).
lubrication effect (though not influencing collision statistics)
Effects of AIs on kinematic and dynamic collision statistics
The effect of turbulent energy dissipation rate on the at-contact RDF for different (c) droplet sizes and (d) Stokes numbers.
The DNS results of Rosa et al. (2013, figures 8a and 13a there) are additionally included in panels (d).
for low inertial droplets not many pairs at contact (high
RRV is not important)
AI + relocation
Results
Effects of AIs on kinematic and dynamic collision statistics
The dynamic collision kernel as function of (a) droplet size in flows with different energy dissipation rates for Np = 50 000, and (b) the normalized dynamic collision kernel for corresponding values of the Stokes number.
The results of Rosa et al. (2013) are added to panel (b) for comparison.
higher inertia (curves)
higher inertia higher inertia
(each curve)
mainly RDF
Comparison of the new lubrication-included HDNS with standard HDNS and simulations without AIs
Changes in the at-contact (a) RRVs and (b) RDFs when the effect of lubrication forces is taken into consideration compared with the standard HDNS without lubrication effects as well as the case without aerodynamic interaction for Np = 50 000
droplets in a flow at a low dissipation rate ε = 50 cm2/s3. lubrication effect (though not
influencing collision statistics)
lubrication effect
lubrication effect
Results
Comparison of the new lubrication-included HDNS with standard HDNS and simulations without AIs
Comparison of at-contact RRV in case when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of droplets,
Np = 50 000 and (b) the same mass loading, φ = 1.22×10−2 at the dissipation rate ε = 400 cm2/s3. results for
higher inertia + gravity both more pronounced!
lubrication effect
lubrication effect
Comparison of the new lubrication-included HDNS with standard HDNS and simulations without AIs
Comparison of at-contact RDF in case when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of droplets,
Np = 50 000 and (b) the same mass loading, φ = 1.22×10−2 at the dissipation rate ε = 400 cm2/s3.
The empirical model of Ayala, Rosa & Wang (2008, (84) and (85) for Rλ = 76.86) is also included for comparison.
results for
higher inertia + gravity
lubrication
effect lubrication effect
clustering
Results
Comparison of the new lubrication-included HDNS with standard HDNS and simulations without AIs
Comparison of the dynamic collision kernel when there is aerodynamic interaction, both including and excluding lubrication forces, and when there is no aerodynamic interaction, all with and without gravity, for (a) the same number of
droplets, Np = 50 000 and for (c) the same mass loading, φ = 1.22×10−2 at the dissipation rate ε = 400 cm2/s3.
Not much change (g)
slight reduction
slight reduction (g) slight reduction
HDNS + Flub vs HDNS (slight decrease)
Corrections to kinematic formulations due to non-overlapping droplets condition
a2 a1
Collision sphere R = a1 + a2
≡ Cw
≡ Cg
Results
Corrections to kinematic formulations due to non-overlapping droplets condition
a2 a1
Collision sphere R = a1 + a2
0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
1 10
C g&C w
r / R δsh = 5%R
Validation of corrected kinematic formulations
Comparisons between the collision kernel obtained from the dynamic formulation with that computed using the kinematic formulation, both before (UC) and after (C) applying corrections,
when (a) gravity is not considered and (b) when there is gravity.
matched after correction
matched after correction
ΓK = ΓD no relocation after collision
Non-overlapping droplets:
relocation after collision
Overlapping (no relocation)
Results
Effect of particles number density and mass loading
Changes in the at-contact (a) RRV, (b) at-contact RDF and (c) the dynamic collision kernel for the sets of droplets with different radii and numbers, i.e. different mass loadings.
Kinematics and dynamics are more sensitive to φ with gravity and inertia (decorrelation)
• Due to the non-linearity of lubrication forces, at-contact kinematics strongly depend on the chosen time step size, while for larger separation distances there is not much change in kinematics by various choices.
• For moderate-to-large inertia droplets, considering lubrication decreases RRV and increases RDF, whereas for small inertia droplets the RRV increases which does not affect collision statistics due to very low clustering.
• By additionally considering lubrication effects, there is a slight decrease in the rates of collision.
• The effect of lubrication forces is more pronounced in systems with larger energy dissipation rates, especially if gravitational settling is considered.
• The corrections (due to relocation) applied to kinematics of droplet statistics can accurately recover values corresponding to dynamic formulation.
• When the mass loading grows, the kinematic collision statistics reveal an opposite trend, namely the RRV increases with the droplet number density, while the RDF monotonically decreases.
Thank you
