Experimental investigations on the nonequilibrium dynamics of pattern   formation in fluid and granular systems

Xiaolei Ma

arXiv:1905.12201·cond-mat.soft·May 30, 2019

Experimental investigations on the nonequilibrium dynamics of pattern formation in fluid and granular systems

Xiaolei Ma

PDF

Open Access

TL;DR

This paper presents experimental studies on pattern formation in fluid and granular systems, revealing hydrodynamic and desiccation mechanisms underlying diverse multiscale patterns.

Contribution

It uncovers the hydrodynamic origin of star-shaped oscillations in Leidenfrost drops and identifies distinct desiccation crack mechanisms with universal scaling laws.

Findings

01

Hydrodynamic coupling excites star-shaped Leidenfrost oscillations

02

Multi-scale polygonal crack patterns follow a universal power law

03

Effective particle repulsion influences sedimentation patterns

Abstract

Patterns are quotidian in nature. Distinct multiscale patterns are generally a consequence of nonequilibrium dynamical processes associated with mechanical or hydrodynamic instabilities. In this thesis, I report experimental investigations on pattern formation in a few examples of fluid and granular systems, and uncover the underlying mechanisms that give rise to those patterns. Leidenfrost drops are known to experience star-shaped oscillations with little damping. However, the underlying mechanism remains unclear. Here I report that the hydrodynamic coupling between the rapid evaporated vapor flow and vapor-liquid interface excites the star-shaped oscillations, suggesting a purely hydrodynamic origin. Polygonal desiccation crack patterns are commonly observed in natural systems. However, it is unclear whether similar crack patterns spanning multiple length scales share the same…

Tables3

Table 1. Table 1: Physical properties of different liquids at the boiling point T b subscript 𝑇 𝑏 T_{b} (K). Units are as follows: γ 𝛾 \gamma (mN/m), ρ l subscript 𝜌 𝑙 \rho_{l} (kg/m 3 ), η l subscript 𝜂 𝑙 \eta_{l} (mPa s), l c subscript 𝑙 𝑐 l_{c} (mm). Data was taken from Ref. [ lemmon2011nist ] . The last three columns indicate the range of substrate temperatures, T s subscript 𝑇 𝑠 T_{s} (K), used in the experiments, the observed mode numbers ( n 𝑛 n ), and the Reynolds number computed in Section 5 .

liquid	$T_{b}$	$γ$	$ρ_{l}$	$η_{l}$	$ł_{c}$	$T_{s}$	Modes	$R e_{l}$
water	373	59.0	958	0.282	2.5	523-773	2-13	1340
liquid N₂	77	8.90	807	0.162	1.1	298	3-5,7	539
acetone	329	18.2	727	0.242	1.6	523	5-10	601
methanol	338	18.9	748	0.295	1.6	523	6-10	511
ethanol	352	18.6	750	0.420	1.6	523	7-11	355
isopropanol	356	15.7	723	0.460	1.5	523	9,10	283

Table 2. Table 2: Physical properties of water at the boiling points [ lemmon2011nist ] .

Symbol	Quantity	Value
$ν_{l}$	liquid kinematic viscosity	2.94 $\times$ 10^-7 m²/s
$D_{l}$	liquid thermal diffusivity	1.743 $\times$ 10^-7 m²/s
$α$	thermal expansion coefficient	7.52 $\times$ 10^-4 /K
$κ_{l}$	liquid thermal conductivity	3.2 $\times$ 10^-2 W/m/K
$h$	drop thickness	0.005 m
$δ T$	temperature difference	$T_{hot} - T_{cold}$
$β$	heat transfer coefficient	depends on $δ T$
$\partial γ / \partial T$	surface tension gradient	1.46 $\times$ 10^-4 N/m/K
$T_{1}$	surface temperature	$>$ 367 K

Table 3. Table 1: Physical properties of the parameters in Eq. 1 .

Symbol	Meaning	Value
$G$	particle shear modulus	4 GPa (cornstarch)
$G$	particle shear modulus	32 GPa (CaCO₃)
$R$	particle radius	5 $μ$ m (cornstarch)
$R$	particle radius	1 $μ$ m (CaCO₃)
$γ$	surface tension	72 mN/m (water)
		23 mN/m (IPA)
		16 mN/m (silicone oil)
$M$	coordination number	5
$ϕ$	random close packing	0.67
$h_{c}$	critical thickness
$σ_{c}$	critical stress

Equations129

\frac{\partial u}{\partial t} + (u \cdot \nabla) u = g - \frac{1}{ρ} \nabla p + ν \nabla^{2} u,

\frac{\partial u}{\partial t} + (u \cdot \nabla) u = g - \frac{1}{ρ} \nabla p + ν \nabla^{2} u,

\frac{\partial ρ}{\partial t} = 0,

\frac{\partial ρ}{\partial t} = 0,

\nabla ρ = 0,

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = \frac{\partial ρ}{\partial t} + \nabla ρ \cdot u + ρ (\nabla \cdot u) = 0,

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = \frac{\partial ρ}{\partial t} + \nabla ρ \cdot u + ρ (\nabla \cdot u) = 0,

\frac{D ρ}{D t} = \frac{\partial ρ}{\partial t} + u \cdot \nabla ρ = - ρ (\nabla \cdot u) = 0,

\frac{D ρ}{D t} = \frac{\partial ρ}{\partial t} + u \cdot \nabla ρ = - ρ (\nabla \cdot u) = 0,

\nabla \cdot u = 0,

\nabla \cdot u = 0,

\nabla p = η \nabla^{2} u,

\nabla p = η \nabla^{2} u,

6 π η a U = \frac{4}{3} (ρ_{p} - ρ_{f}) π a^{3} g,

6 π η a U = \frac{4}{3} (ρ_{p} - ρ_{f}) π a^{3} g,

U_{s t o k es} = \frac{2}{9} \frac{a ^{2}}{η} (ρ_{p} - ρ_{f}) g .

U_{s t o k es} = \frac{2}{9} \frac{a ^{2}}{η} (ρ_{p} - ρ_{f}) g .

u_{r} = (1 - \frac{3}{2} \frac{a}{r} + \frac{1}{2} \frac{a ^{3}}{r ^{3}}) U cos θ,

u_{r} = (1 - \frac{3}{2} \frac{a}{r} + \frac{1}{2} \frac{a ^{3}}{r ^{3}}) U cos θ,

u_{θ} = - (1 - \frac{3}{4} \frac{a}{r} - \frac{1}{4} \frac{a ^{3}}{r ^{3}}) U sin θ,

u_{ϕ} = 0,

p = - \frac{3}{2} η \frac{a}{r ^{2}} U cos θ,

u \sim 1 - \frac{1}{r} ≃ \frac{1}{r},

u \sim 1 - \frac{1}{r} ≃ \frac{1}{r},

p \sim \frac{1}{r ^{2}} .

p \sim \frac{1}{r ^{2}} .

S (r) = e + \frac{r ^{2}}{2 R _{s}} + ϵ cos (\frac{2 π r}{λ _{c}}),

S (r) = e + \frac{r ^{2}}{2 R _{s}} + ϵ cos (\frac{2 π r}{λ _{c}}),

P_{0} \approx γ κ ∣_{r \to 0}

P_{0} \approx γ κ ∣_{r \to 0}

P_{1} \approx γ κ ∣_{r \to λ_{c} /2} - ρ_{l} g Δ z

λ_{c} = 2 (2 ϵ R_{s} - \frac{2 ρ _{l} g ϵ R _{s} ( - 3 π ^{2} γ + 2 ρ _{l} g ϵ R _{s} )}{ρ _{l} g})^{1/2} .

λ_{c} = 2 (2 ϵ R_{s} - \frac{2 ρ _{l} g ϵ R _{s} ( - 3 π ^{2} γ + 2 ρ _{l} g ϵ R _{s} )}{ρ _{l} g})^{1/2} .

\frac{λ _{c}}{l _{c}} = 2 (2 χ - 2 χ (2 χ - 3 π^{2}))^{1/2},

\frac{λ _{c}}{l _{c}} = 2 (2 χ - 2 χ (2 χ - 3 π^{2}))^{1/2},

v = \frac{R ^{'} ( t )}{R ( t )} (r, 0, - 2 z) .

v = \frac{R ^{'} ( t )}{R ( t )} (r, 0, - 2 z) .

T =

T =

T =

T =

U = \frac{1}{2} ρ_{l} V g h + 2 π γ R (h + R) = \frac{g V ^{2} ρ _{l} + 4 π γ R ( V + π R ^{3} )}{2 π R ^{2}} .

U = \frac{1}{2} ρ_{l} V g h + 2 π γ R (h + R) = \frac{g V ^{2} ρ _{l} + 4 π γ R ( V + π R ^{3} )}{2 π R ^{2}} .

\frac{d U}{d R} = 0 = 4 π γ R_{0} - \frac{2 V γ}{R _{0}^{2}} - \frac{g ρ _{l} V ^{2}}{π R _{0}^{3}} .

\frac{d U}{d R} = 0 = 4 π γ R_{0} - \frac{2 V γ}{R _{0}^{2}} - \frac{g ρ _{l} V ^{2}}{π R _{0}^{3}} .

h \approx 2 \frac{γ}{ρ _{l} g}, R_{0} \to \infty,

h \approx 2 \frac{γ}{ρ _{l} g}, R_{0} \to \infty,

R^{''} = \frac{6 ρ _{l} V ^{2} g π R ^{4} + 12 γ V π ^{2} R ^{5} - 24 γ π ^{3} R ^{8} + 24 ρ _{l} V ^{3} R ^{'2}}{8 ρ _{l} V ^{3} R + 3 ρ _{l} V π ^{2} R ^{7}} .

R^{''} = \frac{6 ρ _{l} V ^{2} g π R ^{4} + 12 γ V π ^{2} R ^{5} - 24 γ π ^{3} R ^{8} + 24 ρ _{l} V ^{3} R ^{'2}}{8 ρ _{l} V ^{3} R + 3 ρ _{l} V π ^{2} R ^{7}} .

ω^{2} = \frac{576 γ ^{2} ( γ + Γ ) + 36 ρ _{l}^{2} g ^{2} R _{0}^{4} ( γ + 2Γ ) + 384 γ ρ _{l} g R _{0}^{2} ( 5 γ + 2Γ )}{ρ _{l} R _{0}^{3} ( 1120 γ ^{2} + 9 ρ _{l}^{2} g ^{2} R _{0}^{4} + 192 γ ρ _{l} g R _{0}^{2} )},

ω^{2} = \frac{576 γ ^{2} ( γ + Γ ) + 36 ρ _{l}^{2} g ^{2} R _{0}^{4} ( γ + 2Γ ) + 384 γ ρ _{l} g R _{0}^{2} ( 5 γ + 2Γ )}{ρ _{l} R _{0}^{3} ( 1120 γ ^{2} + 9 ρ _{l}^{2} g ^{2} R _{0}^{4} + 192 γ ρ _{l} g R _{0}^{2} )},

ω \approx (\frac{36 γ}{35 ρ _{l} R _{0}^{3}})^{1/2}, R_{0} \to 0

ω \approx (\frac{36 γ}{35 ρ _{l} R _{0}^{3}})^{1/2}, R_{0} \to 0

ω \approx \frac{4}{R _{0}} (\frac{g γ}{ρ _{l}})^{1/4}, R_{0} \to \infty.

f_{n} = \frac{1}{2 π} \frac{γ n ( n - 1 ) ( n + 2 )}{ρ _{l} R ^{3}} .

f_{n} = \frac{1}{2 π} \frac{γ n ( n - 1 ) ( n + 2 )}{ρ _{l} R ^{3}} .

f_{n} = \frac{1}{2 π} \frac{γ n ( n ^{2} - 1 )}{ρ _{l} R ^{3}} \frac{1}{1 + ( 2 - \frac{π}{2} + \frac{n - 3}{4} ) \frac{l _{c}}{R}} .

f_{n} = \frac{1}{2 π} \frac{γ n ( n ^{2} - 1 )}{ρ _{l} R ^{3}} \frac{1}{1 + ( 2 - \frac{π}{2} + \frac{n - 3}{4} ) \frac{l _{c}}{R}} .

\frac{d ^{2} r}{d t ^{2}} + ω^{2} r = 0.

\frac{d ^{2} r}{d t ^{2}} + ω^{2} r = 0.

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Taxonomy

TopicsFluid Dynamics and Heat Transfer · Fluid Dynamics and Thin Films · Nanomaterials and Printing Technologies

Full text

\dsp

See pages 1,2 of Special_Pages

\degreemonth\degreeyear

2019 \degreeDoctor of Philosophy \fieldPhysics \departmentPhysics \advisorJustin Clifford Burton

Experimental investigations on the nonequilibrium dynamics of pattern formation in fluid and granular systems

Xiaolei Ma

Abstract

Patterns are quotidian in nature and occur on multiple length scales. Distinct multiscale patterns are generally a consequence of nonequilibrium dynamical processes associated with particular mechanical and hydrodynamic instabilities, which play a vital role in shaping the pattern geometry. In this thesis, I report experimental investigations on the pattern formation in a few examples of fluid and granular systems, and uncover the underlying mechanisms that give rise to those patterns.

Leidenfrost drops are known to experience star-shaped oscillations with little damping. However, the underlying mechanism remains unclear. Here I report that the hydrodynamic coupling between the rapid evaporated vapor flow and vapor-liquid interface excites the star-shaped oscillations, suggesting a purely hydrodynamic origin. In addition, I also give an analytical explanation for an oscillatory “breathing mode” found in small Leidenfrost drops.

Polygonal desiccation crack patterns are commonly observed in natural systems. However, it is unclear whether similar crack patterns spanning multiple length scales share the same underlying physics. I also report experimental investigation on polygonal cracks in drying suspensions of micron-sized particles. In cornstarch-water mixtures, multi-scale crack patterns were observed due to two distinct desiccation mechanisms. In addition, we find that the characteristic area of the polygonal cracks ( $A_{p}$ ), and film thickness ( $h$ ) obey a universal power law, $A_{p}=\alpha h^{4/3}$ . Thus we provide a robust framework for understanding multiscale polygonal crack patterns.

Finally, I report experimental results on sedimentation of non-Brownian particles in viscous fluids, which is crucial in both nature and industrial processes. We observed an effective repulsion between particles with nonuniform density in both two-body and many-body systems, in contrast to particles with uniform density. This trend holds true in two and three dimensions. In addition, we also characterize the statistical properties of the sedimentation patterns of particles in three dimensions. Our results also shed light on the potential for controlling the uniformity of particle layers after sedimentation.

The patterns I report in this thesis represent typical examples in fluid and granular systems that are driven by nonequilibrium dynamics, and the underlying mechanisms we uncover are expected to enhance our understanding of how these seemingly simple patterns can arise in natural systems.

\makecoverpage

0 Introduction
1 Pattern formation in nature
2 Patterns in fluids
3 Patterns in colloidal films
4 Patterns driven by sedimentation
1 Oscillations of Leidenfrost drops
1 Introduction
2 Experimental setup
3 Results and discussion
1 The geometry of Leidenfrost drops on curved surfaces
2 Breathing mode of small Leidenfrost drops
3 Star-shaped oscillations of Leidenfrost water drops
4 Pressure oscillations in the vapor layer
5 Star-shaped oscillations of different liquids
6 Origin of pressure oscillations in the vapor layer
7 Thermal effects
4 Summary and outlook
2 Polygonal desiccation crack patterns
1 Introduction
2 Experimental setup
3 Results and discussion
1 Multiscale cracks in cornstarch-water suspensions
2 Critical condition for primary cracks
3 Cornstarch particle deswelling drives secondary cracks
4 Primary cracks in different particle suspensions
5 Cornstarch-water suspensions in thin chambers
6 Universal scaling of multiscale polygonal cracks
7 Effective film thickness for cracks in thick cornstarch-water suspensions
4 Summary and outlook
3 Sedimentation of non-Brownian particles
1 Introduction
2 Experimental setup
3 Results and discussion
1 Effect of density distribution on sedimentation
2 Sedimentation of multiple particles in quasi two dimensions
3 Sedimentation of multiple particles in three dimensions
4 Summary and outlook
4 Summary

List of Figures

1 Patterns in nature
2 Patterns driven by hydrodynamic instabilities
3 Star-shaped oscillation modes of Leidenfrost water drops
4 Crack patterns by drying colloidal suspensions
5 Patterns driven by sedimentation
6 Sketch of the sedimentation of a spherical particle in a viscous fluid
7 Sketch of Stokes flow around a sphere in a Cartesian coordinate system
1 Leidenfrost effect and the geometry of Leidenfrost drops
2 Experimental setup for star-shaped oscillations of Leidenfrost drops
3 Schematic of the suppression of Rayleigh-Taylor instability by surface curvature
4 Threshold of the perturbation amplitude for the suppression of Rayleigh-Tayor instability
5 Drop radius $R$ versus volume $V$ for star-shaped Leidenfrost water drops
6 Breathing oscillation mode of small Leidenfrost drops
7 Drop oscillation wavelength $\lambda_{d}$ and frequency $f_{d}$ versus drop radius $R$ for star-shaped Leidenfrost water drops
8 Initiation of star-shaped oscillation of a Leidenfrost water drop
9 Pressure oscillation frequency $f_{p}$ and drop oscillation frequency $f_{d}$ versus drop radius $R$
10 Images of star-shaped oscillation modes with $n=2$ to 13 of six different liquids when the lobes are at their maximum displacement.
11 Parametric oscillations of star-shaped oscillations of Leidenfrost drops
12 Threshold for star-shaped oscillations of Leidenfrost drops
13 Sketch of the geometry of a large Leidenfrost drop
14 Variations of vapor film thickness $e$ , radial flow velocity $u$ , vertical flow velocity $v$ , and the Reynolds number of the vapor $Re_{v}$ with respect to temperature difference $\Delta T$
15 Sketch of the experimental setup for imaging capillary waves underneath Leidenfrost drops
16 Fourier power spectra of the capillary waves underneath star-shaped Leidenfrost drops
17 Temperature profile of a large Leidenfrost water drop
18 Fourier power spectra of the pressure variations in the vapor layer at different environment temperature $T_{1}$
1 Multiscale polygonal crack patterns
2 Experimental setup for drying cornstarch-water suspensions in a petri dish
3 Effect of stress boundary condition on crack formation
4 Measurement of critical film thickness for crack formation
5 Deswelling of cornstarch particles during drying
6 Gallery of polygonal desiccation crack patterns of different particulate suspensions
7 Packing behavior of cornstarch particles dispersed in different liquids
8 Packing behavior of cornstarch, CaCO3 and glass beads particles dispersed in water
9 Effect of suspending fluids on drying behavior of particulate suspensions
10 Schematic for drying cornstarch-water suspensions in quasi-two-dimensional chambers
11 Polygonal crack patterns of fully-desiccated cornstarch-water suspensions in thin chambers
12 Scaling behavior of polygon crack area $A_{p}$ vs. film thickness $h$ for cornstarch-water suspensions with various initial volume fractions $\phi_{i}$ after drying in thin quasi-2D chambers
13 Characteristic area ( $A_{p}$ ) vs. the film thickness ( $h$ ) of the multiscale polygonal cracks observed in both thin chambers and petri dishes
14 Schematic for cracking in a material driven by tensile stress at the substrate
15 The variation of $\kappa$ with respect to $N$
16 Modulus measurement of particulate suspensions at different drying stages
17 Modulus measurement of fully desiccated particulate suspensions
18 Rheology measurements of particulate suspensions with different initial volume fractions $\phi_{i}$
19 Schematic for the effective film thickness during drying thick samples of cornstarch-water suspensions
1 Response to fluid disturbance for particles with different density distributions
2 Sketch of the sedimentation of a rod-like object
3 Sketch of quasi-2D chamber for sedimentation
4 Sketch of 3D chamber for sedimentation
5 Sedimentation of particles with different density distributions
6 Sedimentation of two Al-St particles in a quasi-2D chamber
7 Sedimentation of two Al-Al particles in a quasi-2D chamber
8 Sedimentation of many particles in a quasi-2D chamber
9 Probability density function (PDF) of particle spacing for sedimentation of many particles in a quasi-2D chamber.
10 Images of the sedimentation pattern of Al-Al particles in 3D
11 Images of the sedimentation pattern of Al-St particles in 3D
12 Probability density function (PDF) of particle clusters after sedimentation in 3D
13 Theoretical prediction of the effect of geometry and density distribution of particles on the sedimentation dynamics

List of Tables

1 Physical properties of liquids at the boiling point $T_{b}$ (K). Units are as follows: $\gamma$ (mN/m), $\rho_{l}$ (kg/m3), $\eta_{l}$ (mPa s), $l_{c}$ (mm). The last three columns indicate the range of substrate temperatures, $T_{s}$ (K), used in the experiments, the observed mode numbers ( $n$ ), and the Reynolds number computed in Section 5
2 The values of the physical quantities of liquids at the boiling points
1 Physical properties of the parameters in Eq. 1

{citations}

Chapter 2 contains research from the following 3 publications:

$\bullet$ “The many faces of a Leidenfrost drop”, Xiaolei Ma, Juan-José Liétor-Santos, Justin C. Burton, Physics of Fluids, 27, 091109, (2015)

$\bullet$ “Star-shaped oscillations of Leidenfrost drops”, Xiaolei Ma, Juan-José Liétor-Santos, Justin C. Burton, Physical Review Fluids, 2, 031602(R), (2017)

$\bullet$ “Self-organized oscillations of Leidenfrost drops”, Xiaolei Ma, Justin C. Burton, Journal of Fluid Mechanics 846, 263 (2018)

Chapter 3 contains research from the following publication:

$\bullet$ “Universal scaling of polygonal desiccation crack patterns”, Xiaolei Ma, Janna Lowensohn, Justin C. Burton, Physical Review E, 99, 012802 (2019)

Chapter 4 contains research from the following paper:

$\bullet$ “Interactions between non-Brownian particles with nonuniform density during sedimentation in viscous fluids”, Xiaolei Ma, Justin C. Burton, to be submitted to Physical Review Fluids

Acknowledgements.

With the submission of this thesis, I would like to acknowledge the people who have played a key role on my way to pursing a Ph.D. degree in physics. First of all, I wish to thank Prof. Justin Burton, my thesis advisor and role model, for the wonderful research training under his supervision in the last five years. Discussing science with Justin is a really enjoyable experience during which I not only learned specific technical tools for research but more importantly got the spirit of thinking about problems as a physicist and to be an independent researcher. Besides the role as an academic advisor, he is also a great friend in my ordinary life, the activities he organized outside Emory campus are so much fun and made my life here really wonderful. It is his patience, bright personality, research passion, immense knowledge and financial support that lead to this thesis, and the times with him at Emory will always be appreciated in the rest of my life. I am especially grateful to Prof. Eric Weeks, one of my thesis committee board, both for the constructive suggestions he made for my research activities, and the permission and convenience he provided to do experiment in his lab. In addition, I also want to thank the rest of my thesis committee board: Prof. Alberto Fernández-Nieves, Prof. Hayk Harutyunyan and Prof. Connie Roth for their insightful suggestions on my research. Besides, I wish to acknowledge my colleague Dr. Juan-José Liétor-Santos for his mentoring and extensive discussions with me when I initially joined the Burton Lab, and my colleagues Cong Cao and Carlos Orellana for many discussions about rheology. I also want to thank the current and past Burton Lab members: Nicholas Cuccia, Stephen Frazier, Yannic Gagnon, Guram Gogia, Dana Harvey, Janna Lowensohn, Joshua Méndez, Asher Mouat, Justin Pye, Dominic Robe, Toyin Thompson, Clay Wood, Rui Wu, and Jiaqi Zheng for the awesome lab environment they created and useful discussions. I really thank you all for the laughters you brought to the lab, which made my times in the Burton Lab unforgettable. Special thanks go to my friends outside Emory, especially Prof. Guangyin Jing who has brought me to the field of soft matter. I really appreciate the inspiring and enlightening discussions with him, which also shed light for my future career. I also want to thank Hao Luo and Mo Zhou for the helpful discussions about drying crack patterns. And also my friends Xin Tian, Mingyuan He, Bo Jiang, Jiang Zhang, Long Yue, Zhan Zhao, Di Dong, Weixiao Hou, Fuyou Liao, Xiaohe Yan, Zhihui Yang, thanks for the joyful times with all of you, and the funny stories you shared with me. Finally and most importantly, I wish to express my heartfelt gratitude to my family, this thesis would certainly not have been possible without their unrequited love, unwavering support, and faithful encouragement.

Xiaolei Ma

Emory University, Atlanta, USA

Spring 2019

\dedication

\hsp

*Dedicated to my parents

\copyrightpage

\startarabicpagination

Chapter 0 Introduction

1 Pattern formation in nature

The nature is composed of diverse patterns existing from nanoscopic to geophysical scales, such as the examples shown in Fig. 1. The cause of patterns can be very complicated depending on particular physical mechanisms. For example, the striking six-fold symmetry of a snowflake, as shown in Fig. 1a, is basically a particular signature of crystallization of water molecules. During this process, the water molecules align themselves in order to localize in the lowest energy state. Based on this principle, snowflakes develop hexagonal symmetry, and in reality there are different types of snowflakes depending on the surrounding temperature and humidity [libbrecht2005physics]. Figure 1b shows the sand ripples that can be easily observed in a desert or beach. The formation of these nearly periodic patterns originate from the interactions between the fluid (liquid or wind) flow and bed topography. The rough process can be as follows: when the fluid is moving fast, the sand grains will be lifted, gaining translational and rotational energy. In a higher energy state, the grains on the top of a dune are unstable and prone to flow down to dissipate energy when perturbed. As a consequence, the sand grains display such an amazing pattern with the wavelength depending on the fluid flow velocity and depth of the granular material [charru2013sand, cross2009pattern]. A spectacular geologic column joint pattern found in cooled lava flows is shown in Fig. 1c. Some well-known examples of those hexagonal fracture patterns can be found for instance in Giant’s Causeway, Northern Ireland, Devils Tower, Wyoming, USA, and Devils Postpile, California, USA. As the hot lava flows are exposed to air, cooling takes place from outside in, generating enormous stress inside the rock. Fractures then form at the surface and grow into the material, which usually takes decades. Depending on the specific environmental settings, the number of sides of the joints can vary and joint structures with the number of sides ranging from 4 to 7 have been observed in reality, nevertheless, hexagonal symmetry is the most common scenario, which is likely due to a particular ordering process [goehring2014cracking, hofmann2015hexagonal].

In addition to the amazing patterns observed in non-living systems, biological organisms display a tremendous diversity of patterns, such as the branching pattern of veins in a leaf, the fractal structure in a Romanesco broccoli, the spiral pattern in a nautilus shell, the Turing pattern in a zebra’s skin, and the fractures in a crocodile’s skin. The origins for the development of the diverse and fascinating patterns observed in biological organisms are quite complicated due to the highly nonequilibrium nature in biological systems which usually involve cell differentiation, cell growth and death, and interaction with environmental conditions, however, many models have to some extent successfully explained the origin of patterns observed in biological systems at the molecular level or based on mechanical principles on the mesoscopic scale, although oversimplified in many cases [koch1994biological, ouyang1991transition, sachs2005pattern, nakamasu2009interactions, boettiger2009neural, fujita2006pattern, milinkovitch2013crocodile].

Despite the difficulty in fully understanding the formation of patterns in diverse systems in nature, there have been a number of widely applicable theories and models that are now employed as frameworks to describe the generic dynamics of the nonequilibrium growth of patterns in certain systems, for instance, the diffusion-limited aggregation of Brownian particles [witten1981diffusion, garcia2012growth] and the spinodal decomposition of a mixture of liquids [cahn1961spinodal, cahn1965phase]. Although our understanding of the fascinating patterns observed in nature is pretty much in its infancy, there are growing interests in uncovering the underlying mechanisms behind these nonequilibrium patterns physically and mathematically [langer1980instabilities, cross1993pattern, aranson2006patterns, hoyle2006pattern, cross2009pattern]. In Sections 2, 3 and 4, I will present some fundamental background for pattern formation in particular systems of fluids, colloidal suspensions, and sedimentation, respectively.

2 Patterns in fluids

The world of fluids is rich in intriguing patterns, which can essentially be explained by hydrodynamic theories [lamb1993hydrodynamics, drazin1982hydrodynamic, charru2011hydrodynamic]. Generally, the origin of patterns observed in incompressible fluid systems can be understood by the Navier-Stokes equation:

[TABLE]

In most cases, the fluid of interest is assumed to be incompressible, i.e.,

[TABLE]

and given the condition of mass conservation:

[TABLE]

we can simply arrive:

[TABLE]

which yields the continuity equation:

[TABLE]

where u is the flow velocity, $\rho$ is the fluid density, $\nu=\eta/\rho$ is the kinematic viscosity wherein $\eta$ is the dynamic viscosity, g is acceleration field of the volume forces, and $p$ is the pressure. In certain regimes, the viscous force is dominant over inertial forces and gravity, and the nonlinear, advective inertial terms in Eq. 1 can be safely ignored, leading to:

[TABLE]

which is known as the Stokes equation. Equation 7 is much easier to deal with both mathematically and numerically than Eq. 1, and most hydrodynamic phenomena in this regime can be analytically understood due to the linearity of Eq. 7.

Figure 2a shows the hexagonal convection cells which can be observed when heating a liquid from below with the top surface exposed to air. Due to the perturbation, the warm fluid will move upward, thus bringing local perturbation to the temperature of the free surface and leads to a surface tension gradient at the interface resultantly. The surface tension gradient induces a Marangoni stress, which drives the fluid to spread out thus forming the hexagonal pattern. In this systems, the resisting forces are caused by the liquid viscosity and thermal diffusion, which tend to stabilize the temperature variations in the fluid and make the fluid more homogeneous. Usually, a dimensionless Marangoni number defined as the ratio of surface tension effect to viscous and thermal diffusive effects is used to characterize to the relative importance, and when the Marangoni number reaches a critical value, hexagonal cell patterns should occur between the bottom layer and the top free surface with the size of the cell depending on the thickness of the fluid layer, which is the so-called Bénard-Marangoni thermocapillary instability [benard1901tourbillons, rayleigh1916lix, marangoni1871ausbreitung, schatz2001experiments, maroto2007introductory, guyon2001physical]. Another example is the viscous fingering patterns that usually occur when a less viscous liquid is injected into a more viscous liquid in a rectangular configuration, for instance, a Hele-Shaw cell, as shown in Fig. 2b. Once the critical injecting velocity of the less viscous liquid is reached, an instability will be excited and lends to the front to be unstable for growth of fingers, which consequently initiates a fingering pattern known as the Saffman-Taylor instability [saffman1958penetration, homsy1987viscous]. Figure 2c shows the wave pattern in a cloud layer, which is driven by a shear flow at the interface composed of two fluids with different densities (heavy fluid on the bottom and light fluid on the top). Basically, a shear flow at the interface will cause a pressure gradient in the vicinity of the interface, then the interface will be unstable and heavy fluid will be lifted up whereas the light fluid be pushed down, thus forming vortex-like wave pattern in the could layer. This is known as the Kelvin-Helmholtz instability [charru2011hydrodynamic]. Besides the pattern shown in Fig. 2c, the Kelvin-Helmholtz instability also manifests in other well-known examples, such as Jupiter’s Red Spot, the Sun’s corona, and in ocean waves [marcus1988numerical, foullon2011magnetic, smyth2012ocean].

In addition to the gorgeous patterns displayed in Fig. 2, isolated drops are also able to show fascinating patterns. Leidenfrost drops are known for the star-shaped oscillations that spontaneously form, as shown in Fig. 3. The Leidenfrost effect occurs when depositing a liquid drop onto a hot surface well above the boiling point of the liquid. Then the drop will levitate on a vapor cushion fed by the drop evaporation, consequently the Leidenfrost drops will undergo a long lifetime as the vapor insulates the liquid from the hot solid substrate [biance2003leidenfrost, quere2013leidenfrost].

The capillary length is a length scale that characterizes the relative importance of gravitational force and surface tension. Assuming a spherical liquid drop with a size, $l$ , surface tension, $\gamma$ , and density, $\rho_{l}$ , then the balance of surface tension force, $\gamma l$ , and gravitational force, $\rho_{l}gl^{3}$ , leads to a characteristic length scale $l_{c}=\sqrt{\gamma/\rho_{l}g}$ , which is known as the capillary length. For drop with a size larger than $l_{c}$ the gravitational force will dominate over surface tension forces and the drops will be flatten by gravity. More interestingly, for drops with a radius much larger than $l_{c}$ , the thickness of the drops, $h$ , will become constant, i.e., $h\approx 2l_{c}$ , which can be simply derived by balancing the surface tension force, $\gamma/(h/2)$ , and hydrostatic force, $\rho_{l}gh/2$ , per unit length near the edge of the drop. Large Leidenfrost drops easily form large-amplitude, star-shaped oscillations as shown in Fig. 3, which is generally due to the complex coupling between the solid, liquid and vapor phases. As will be shown in Chapter 1 in more detail, the star-shaped oscillations of Leidenfrost drops are parametrically driven by the pressure oscillations in the vapor beneath the drop. The pressure oscillations are initiated by the capillary waves with a characteristic wavelength (frequency) traveling from the drop center to the edge and the capillary waves are excited by a strong shear due to the rapidly flowing evaporated vapor. This suggests that the star-shaped oscillations of Leidenfrost drops are purely hydrodynamic in origin [ma2017star, ma2018self].

3 Patterns in colloidal films

Crack patterns by drying are ubiquitous in nature, artistic works and industrial processes [routh2013drying]. Typical examples include the complex crack network in a dried blood drop [brutin2011pattern], craquelures in old paintings [giorgiutti2015striped, goehring2017drying], T/Y-shaped cracks in dried mud [goehring2010evolution], and polygonal terrain cracks [goehring2014cracking, Maarry2010mars, el2014potential]. Practically, a broad range of applications, such as thin film coating, forensics, and controllable surface patterning rely on the knowledge of the physical processes that determine crack patterns [prosser2012avoiding, hatton2010assembly, liu2016surface, nam2012patterning, zeid2013influence]. To understand the underlying mechanisms behind those versatile, multiscale cracks, drying colloidal suspensions has been widely employed as a model system in the laboratory. Colloid generally refers to a type of multiphase system wherein particles (solid, liquid or gas) with a size between $10^{-9}$ m to $10^{-6}$ m are dispersed in a liquid. Prototypical colloidal systems include milk, blood and ink in our daily lives. One of the most important properties of colloidal particles is that they are susceptible to thermal fluctuations, which leads to Brownian motion due to the random collision with liquid molecules [russel1991colloidal]. However, many crack patterns in natural systems contain particles with sizes from 10 $\mu$ m to 1 mm, so larger particles are also used to investigate crack pattern formation. In contrast to colloidal particles, the major interactions between between these granular particles are from friction and collision, and when the particles size is close to 1 $\mu$ m, other factors such as van der Waals forces, humidity are also likely to play a role in the interaction between granular particles [de1999granular, andreotti2013granular, herminghaus2013wet, franklin2016handbook].

Heretofore, a large number of types of crack patterns have been observed by drying colloidal suspensions under various conditions, and those results indicate that many factors can influence the selection of crack patterns [routh2013drying, goehring2015desiccation]. Figure 4 shows a few typical crack patterns observed in drying colloidal suspensions in recent years. Different substrate wettability can lead to distinct types of crack patterns as indicated by Figs. 4a and 4b which show circular and radial cracks by drying colloidal drops on polystyrene (hydrophobic) and glass (hydrophilic) substrates, respectively [jing2012formation]. The film thickness can also play a role in selecting the crack pattern as shown in Fig. 4c in which spiral crack patterns are observed in drying a latex-water suspension in a circular container with the top surface exposed to air with a film thickness of $\approx$ 200 $\mu$ m, in contrast, only polygonal cracks can be observed in films with thickness of $\approx$ 20 $\mu$ m in the same experimental setup [lazarus2011craquelures]. In addition to the open system shown in Figs. 4a, 4b and 4c, when the colloidal suspensions are confined in a quasi-two-dimensional space, an array of roughly periodically spaced cracks are usually observed with drying as shown in Fig. 4d [dufresne2003flow].

Besides the crack patterns shown in Fig. 4, a large variety of crack patterns have been reported in dried suspensions [allain1995regular, inasawa2012self, goehring2011wavy, nandakishore2016crack, kiatkirakajorn2015formation, neda2002spiral, vermorel2010radial, darwich2012highly, jing2012formation, xu2010imaging, giorgiutti2015striped, giorgiutti2016painting, bohn2005hierarchicalI, bohn2005hierarchicalII, pauchard2003morphologies, boulogne2013annular, allen1987desiccation, gauthier2010shrinkage]. The diversity in observed patterns depends on numerous factors such as film geometry [lazarus2011craquelures, nandakishore2016crack], particle mechanics [nawaz2008effects], liquid additives [pauchard1999influence, liu2014tuning], preparation history [nakahara2006transition, nakahara2006imprinting], solvent volatility [boulogne2012effect, giorgiutti2014elapsed], and external fields [khatun2012electric, pauchard2008crack]. The difficulty in understanding pattern selection and propagation lies in the combinations of elastic instabilities, plastic deformations and multiphase interactions in the such systems [goehring2013plasticity, goehring2015desiccation, goehring2017drying].

Despite this broad range of crack patterns, we know surprisingly little about what controls the size and hierarchy of commonly observed polygonal cracks, which are visible on both the micro- and macro-scales. Although multiple experimental, numerical and theoretical studies have reported specific scalings between the crack spacing and film thickness considering the mechanical, thermodynamic, hydrodynamic, and statistical properties of cracking, however they often lead to contradictory results. Given the importance of understanding crack patterns, in Chapter 2 I report experimental results on the mechanisms underlying polygonal crack patterns by drying micro-sized particulate suspensions. We found distinct mechanisms lead to multiscale crack patterns in drying cornstarch-water suspensions, and the crack area and film thickness of the multiscale cracks we observed in our experiment obey a universal power law that can be derived from a balance of stored elastic energy in the film and the surface energy released upon cracking [ma2018universal].

4 Patterns driven by sedimentation

Sedimentation of particles in fluids plays a key role in shaping natural patterns [allen1970physical, vanoni2006sedimentation], for example the formation of river deltas (Fig. 5a) and geological deposits of sandstone and siltstone (Fig. 5b). In industrial applications, sedimentation is also widely used in the processes of decontamination, neutralization, particle separation and biological purification. Given its crucial role, sedimentation has been extensively studied with a particular focus on the settling velocity fluctuations of mono/poly-disperse, Brownian/non-Brownian particles, density fluctuations, the effects of particle geometry, system size, and many-body hydrodynamics [batchelor1972sedimentation, davis1985sedimentation, guazzelli2011fluctuations, piazza2014settled, mucha2004model], however, the sedimentation is still far from being fully understood due to the complex long-range hydrodynamic, many-body interactions between particles.

Practically, one of the important issues is to precisely control particle clustering and the uniformity of particulate films after sedimentation, which relies on one’s ability to control particle interactions. Recently, it has been theoretically suggested that controlling these interactions can be done by tuning the particle geometry and mass distribution of non-Brownian particles in Stokes flows [goldfriend2017screening]. The key idea is that as the particle geometry center is different from the center of mass, the particles will response anistropically to the fluid disturbance generated by the density fluctuations. This will eventually generate an effective repulsion between particles, suppressing the density fluctuations. In contrast, for symmetric particles (e.g., sphere) with uniform density distributions there is no such repulsion between particles, which is due to the isotropic response to the fluid disturbance endowed by the highly symmetrical geometry. In Chapter 3, I will show some experimental evidence related to those theoretical predictions.

In typical experimental, theoretical and computational studies, viscous fluids are used for sedimentation in order to simplify the analysis, in which the inertial forces can be ignored, and the sedimentation dynamics of particles are merely determined by Stokes flow (i.e., Eq. 7). For a spherical particle settling down in a very viscous fluid as illustrated in Fig. 6, the particle moves very slowly and the motions can be considered steady state. Then balancing the drag, gravitational, and buoyancy forces leads to:

[TABLE]

and the steady state sedimentation velocity (Stokes velocity) is:

[TABLE]

One of the most interesting properties of sedimentation in viscous fluids is the long-range hydrodynamic nature of the interactions between particles. This can be easily derived by considering Eqs. 6 and 7 with appropriate boundary conditions, and the disturbance velocity and pressure fields around a sphere in a spherical coordinate system ( $r$ , $\theta$ , $\phi$ ), the solution is [lautrup2011, guazzelli2011physical]:

[TABLE]

where $U$ is the uniform stream velocity at $r\rightarrow\infty$ .

Figure 7 shows the two dimensional view of the Stokes flow around a solid sphere in a Cartesian coordinate system based on Eqs. 10 to 13. What is different here from the pattern of potential flow around a sphere is that the perturbation created by a translating sphere decays slower and extends further due to the term $O(1/r)$ , which is what we mean by the long-range nature of the interactions.

To leading orders, Eqs. 10 and 11 are reduced to:

[TABLE]

whereas Eq. 13 is simplified to:

[TABLE]

This long-range hydrodynamic characteristic, as indicated by Eq. 14, can result in interesting and complex phenomena in many-body systems, such as collective motions, large-scale velocity and structure fluctuations in the sense of the dynamic coupling between hydrodynamics and the positions and orientations of particles, and considerable progress has been made regarding those problems [jeffery1922motion, xue1992diffusion, segre1997long, ramaswamy2001issues, guazzelli2011fluctuations, piazza2014settled, padding2004hydrodynamic, jung2006periodic, batchelor1972sedimentation, hinch1977averaged, koch1991screening, caflisch1985variance, brenner1999screening, hunter2012physics, chajwa2018kepler]. More recent attention has been paid to controlling the particle interactions during sedimentation so as to control the uniformity of particle cluster patterns after sedimentation.

One of the promising pathways for realizing controllable sedimentation is to tune the density distributions as theoretically suggested by Goldfriend et al. [goldfriend2017screening]. Using linear stability analysis, the authors predict that the particles (non-Brownian) with different density distributions respond very differently to the fluid disturbance, and there should be an effective repulsion between elongated particles with nonuniform density distributions, in contrast to particles with uniform density distributions, for example, a sphere or a rod will align with the fluid flow (see Fig. 1 for comparison). This will create a lateral motion of the particles because of the greater downward gravitational torque of the heavier part compared to the lighter adjacent part. This torque pulls the particle from the preferred alignment direction, consequently pulling the particles to move towards the regions of lower settling velocities and lower densities. This lateral motion is theoretically promising in suppressing density fluctuations so as to realize controllable uniformity of particles after sedimentation.

In Chapter 3, I report our we experimental investigation on the effect of particle density distributions on the sedimentation dynamics. In the experiment, we used equal-sized aluminum and steel balls to construct particles with different density distributions. We observed an effective repulsion between the particles with nonuniform density distributions in both two-body and many-body systems, in contrast to the counterpart of particles with uniform density distributions. Our results suggest that the theoretical predictions by Goldfriend et al. [goldfriend2017screening] are basically correct despite that the density fluctuations are only restricted to the $x$ - $y$ plane (see Fig. 6). It is expected that our results could open a novel pathway for tunable route towards the uniformity of particulate films via sedimentation.

The examples of pattern formation I report in this thesis represent the typical scenarios in fluid and granular systems driven by nonequilibrium dynamics, though in the real world the systems are quite complex and diverse. The underlying mechanisms giving rise to the pattern formation in fluid and granular systems we uncover in laboratory experiments (Chapters 1, 2 and 3) can enhance our understanding of how nonequilibrium dynamics shape the nature, and are also expected to serve as directive tools for manually controlling pattern formation for multiple purposes.

Chapter 1 Oscillations of Leidenfrost drops

1 Introduction

When a millimeter-scale volatile liquid drop (e.g., water) is deposited on a sufficiently hot solid, it can survive for minutes due to the presence of a thermally-insulating layer of evaporated vapor beneath the drop. This phenomenon is known as the Leidenfrost effect for example the Leidenfrost drop shown in Fig. 1a, and the geometry of the Leidenfrost drop is illustrated in Fig. 1b. In this levitated state, commonly known as the Leidenfrost regime [leidenfrost1756aquae], the supporting vapor layer is maintained by the sustained evaporation of the liquid, and individual drops are free to undergo frictionless motion due to the absence of liquid-solid contact.

The Leidenfrost effect can be easily observed by placing a water drop onto a hot pan over a cook stove in the kitchen, and has been the subject of numerous fundamental and applied studies due to the complex and rich interactions between the solid, liquid, and vapor phases [quere2013leidenfrost]. Examples include the evaporation dynamics and geometry of the drop [burton2012geometry, biance2003leidenfrost, myers2009mathematical, pomeau2012leidenfrost, xu2013hydrodynamics, sobac2014leidenfrost, hidalgo2016leidenfrost, maquet2016leidenfrost, wong2017non], the stability of the vapor-liquid interface [duchemin2005static, lister2008shape, snoeijer2009maximum, bouwhuis2013oscillating, trinh2014curvature, raux2015successive, maquet2015leidenfrost, van2019asymptotic], hydrodynamic drag-reduction [vakarelski2011drag, vakarelski2012stabilization, vakarelski2014leidenfrost], self-propulsion of droplets [linke2006self, dupeux2011viscous, dupeux2011trapping, lagubeau2011leidenfrost, cousins2012ratchet, li2016directional, soto2016surfing, sobac2017self, bouillant2018leidenfrost], impact dynamics [biance2006elasticity, tran2012drop, castanet2015drop, shirota2016dynamic], green nanofabrication [abdelaziz2013green], nanoparticle painting [elbahri2007anti], particle self-organization [maquet2014organization], explosion of Leidenfrost liquid mixtures [moreau2019explosive], chemical reactions [bain2016accelerated], fuel combustion [kadota2007microexplosion], quenching process in metallurgy [bernardin1999leidenfrost], heat transfer [talari2018leidenfrost, shahriari2014heat], directional transport [agapov2014length, li2016directional], soft heat engines [waitukaitis2017coupling, wells2015sublimation] and thermal control of nuclear reactors [van1992physics]. The rich phenomena behind Leidenfrost effect essentially stem from the complex coupling between the liquid-vapor interface and the rapid vapor flow beneath the drop, which is however far from being well understood [quere2013leidenfrost]. Here I report our experimental investigation on the coupling dynamics through the study of self-sustained oscillations of Leidenfrost drops.

In many of these examples, transient and self-sustained capillary oscillations play an important role in the dynamics. The interplay between gravity, the flow of vapor beneath the drop, and the liquid surface tension can lead to both small- and large-amplitude oscillations with very little damping. An understanding of these detailed interactions is crucial for the stability of the vapor layer, the failure of which can lead to explosive boiling upon contact with the hot surface. However, the excitation mechanism of these oscillations is complicated by the presence of both hydrodynamic and thermal effects, for example, rapidly-flowing vapor can cause a strong shear stress at the liquid-vapor interface, and temperature gradients in the liquid can lead to convective and Marangoni forces. Here we focus on capillary oscillations in individual Leidenfrost drops, where the shape is mostly determined by the competition between the gravity and surface tension, as measured by the relative size of the drop with respect to the capillary length, $l_{c}\equiv\sqrt{\gamma/\rho_{l}g}$ , where $\gamma$ and $\rho_{l}$ denote the surface tension and density of the liquid, and $g$ is the acceleration due to gravity. For drops with radius $R<l_{c}$ , surface tension forces are dominant, and the drop shape is essentially spherical except for a vanishingly small flat region near the solid surface [burton2012geometry]. Caswell [caswell2014dynamics] identified a planar vibrational mode (“breathing” mode) in the neck region of the drop closest to the solid substrate. The oscillation frequencies were found to obey a power law that is not consistent with a general three-dimensional dispersion relation for capillary waves [rayleigh1879capillary]. Here we provide an analytical expression for the breathing mode using a simple model based on gravity and surface tension which shows excellent agreement with the experimental data.

Large Leidenfrost drops form liquid puddles whose thickness is approximately 2 $l_{c}$ (see Section 2). These puddles are known to spontaneously form large-amplitude, star-shaped oscillations [holter1952vibrations, adachi1984vibration, strier2000nitrogen, snezhko2008pulsating, strier2000nitrogen, ma2015many, Maleidenfrost2016]. Similar oscillations have been observed in large drops on periodically-shaken, hydrophobic surfaces [noblin2005triplon, noblin2009vibrations], drops levitated by an underlying airflow [bouwhuis2013oscillating], and drops excited by an external acoustic or electric field [shen2010parametric, shen2010parametrically, mampallil2013electrowetting]. In studies where the frequency of external forcing is prescribed, the oscillations are excited by a parametric mechanism [brunet2011star]. The external forcing leads to variations in the drop radius with time. Since the drop radius appears in the dispersion relation for azimuthal, star-shaped oscillations, the evolution of the oscillation amplitude obeys an equation similar to the Mathieu equation. For Leidenfrost drops, however, the mechanism is less clear since there is no prescribed frequency, and the star oscillations are excited and sustained through the heat input and resulting evaporation of the liquid. It has been suggested that the star oscillations may result from modulations of the surface tension of the liquid due to temperature variations [adachi1984vibration, takaki1985vibration, tokugawa1994mechanism], or perhaps due to convective patterns [snezhko2008pulsating, strier2000nitrogen]. However, [bouwhuis2013oscillating] observed star-shaped oscillations in drops which are supported by an external, steady air flow, suggest that a hydrodynamic coupling between the gas flow and liquid interface initiates the oscillations.

Given the importance of the Leidenfrost effect in basic fluid and thermal transport, or the numerous practical applications (e.g., self-propulsion and drag-reduction), we know surprisingly little about the coupling between the evaporated vapor flow and vapor-liquid interface that lead to rich dynamical phenomena. Here we explore this coupling by investigating the self-organized, star-shaped oscillations of Leidenfrost drops using six different liquids: water, liquid N2, ethanol, methanol, acetone and isopropanol. The liquid drops were levitated on curved surfaces in order to suppress the Rayleigh-Taylor instability, and star-shaped oscillation modes with $n$ = 2 to 13 lobes along the drop periphery were observed. The number of observed modes depended sensitively on the liquid viscosity, whereas the oscillation frequency (wavelength) depended only on the capillary length but not the mode number, substrate temperature, or drop size. Accompanying pressure measurements in the center of the vapor layer indicate that the pressure variation frequency was approximately twice that the drop oscillation frequency for all of the observed modes, consistent with a parametric forcing mechanism. We show that the pressure oscillations are driven by capillary waves of a characteristic wavelength beneath the drop traveling from the drop center to the edge, and such capillary waves can be generated by a strong shear stress at the liquid-vapor interface. Additionally, we find that although thermal convection is expected to be quite strong, the robust frequency (wavelength) of star oscillations is only weakly affected by varying either the substrate or environmental temperature, suggesting that star-shaped oscillations of Leidenfrost drops are hydrodynamic in origin.

2 Experimental setup

In the experiment, blocks of engineering 6061 aluminum alloy with dimensions of 7.6 cm $\times$ 7.6 cm $\times$ 2.5 cm were used as substrates. Resistive heaters were embedded with high-temperature cement into the aluminum to control the substrate temperature, $T_{s}$ . Six different liquids were used as Leidenfrost drops: deionized water, liquid nitrogen (liquid N2), ethanol, methanol, acetone, and isopropyl alcohol (isopropanol). The relevant liquid properties at the boiling point, $T_{b}$ , such as the surface tension $\gamma$ , density $\rho_{l}$ , dynamic viscosity $\eta_{l}$ , and capillary length $l_{c}$ , are listed in Table 1. The substrates were heated to different temperatures based on $T_{b}$ for each liquid. For water, the temperature of the substrate was set from 523 K to 773 K, for ethanol, methanol, acetone, and isopropanol the temperature of the substrate was set to 523 K, whereas the substrate for liquid N2 was not heated due to its extremely low $T_{b}$ .

The upper surfaces of the substrates were machined into a concave, spherical shape in order to suppress the buoyancy-driven Rayleigh-Taylor instability at the vapor-liquid interface and keep the drops stationary [snoeijer2009maximum, quere2013leidenfrost, trinh2014curvature]. After machining, the roughness of the surface was inspected using optical microscopy with 50 $\times$ magnification. By changing the focus, we determined that over a 250 $\mu$ m $\times$ 250 $\mu$ m surface area, the peak-to-peak roughness was less than 10 $\mu$ m, and often much smaller than this value. The curved surfaces for different liquids were designed to satisfy $l_{c}/R_{s}$ = 0.03, where $R_{s}$ is the radius of curvature of the surface whose top and cross-sectional views are schematically shown in Figs. 2a and 2b, respectively. Following this principle, we fabricated three types of curved substrates, i.e. one for water, one for ethanol, methanol, acetone, and isopropanol, and one for liquid N2, considering their respective capillary lengths ( $l_{c}$ ) listed in Table 1. For some experiments, a plano-concave, fused silica lens (focal length = 250 mm) was used as the heated substrate in order to allow for optical imaging of the capillary waves underneath the Leidenfrost drop. A T-type thermocouple (maximum measurable temperature: 473 K, tip diameter $\approx$ 0.08 cm, HYP-2, Omega Engineering) was used to measure the internal temperature profile of Leidenfrost drops.

For most of the substrates, a pressure sensor (GEMS Sensors, response time: 5 ms, sensitivity: 2 mV/Pa) was connected to a hole (diameter = 1 mm) at the center of the curved substrate in order to measure the pressure variations in the vapor layer at sample rates of 500-1000 Hz during quiescent and oscillatory phases of drops as illustrated in Fig. 2b. We used a high-speed digital camera (Phantom V7.11, Vision Research) with a resolution of about 132 pixels/cm to image the motions of drops from above at frame rates of 1000 frames per second. Recorded videos were then analyzed with NIH ImageJ software to obtain the frequency and wavelength of the star-shaped oscillations (A typical star-shaped oscillation mode $n$ = 8 is schematically shown in Fig. 2c).

3 Results and discussion

1 The geometry of Leidenfrost drops on curved surfaces

The shape of a Leidenfrost drop is determined by a competition between surface tension and gravity. It has been shown that a nonwetting drop with radius $R\ll l_{c}$ will exhibit a quasi-spherical profile except for the bottom region where the drop is slightly flattened by gravity, which still holds true for small Leidenfrost drops [burton2012geometry]. In comparison, for large Leidenfrost drops, i.e. $R\gg\l_{c}$ , the shape of the drops is dominated by gravity and resembles a circular puddle with a constant thickness of approximately 2 $l_{c}$ [biance2003leidenfrost]. However, these large puddles are susceptible to the Rayleigh-Taylor instability, which is manifested by bubbles rising from the vapor layer beneath the drop [biance2003leidenfrost]. This instability can be easily understood by considering a system composed of two fluids with the heavy fluid on the top of the light fluid. Such a system is gravitationally unstable since the heavy fluid is prone to fall through the light fluid, and any fluctuation at the initial flat interface of the two fluids will cause pressure imbalance, which in turn amplifies the fluctuation and deforms the drop. In our experiment, we used curved surfaces (see Section 2), which essentially act as a restoring-force-like role, in order to suppress the gravitationally driven Rayleigh-Taylor instability and obtain large, stable Leidenfrost drops [snoeijer2009maximum, quere2013leidenfrost, trinh2014curvature].

The effects of surface curvature on the Rayleigh-Taylor instability can be seen through the following simple model, similar to the original model used by Biance et al. [biance2003leidenfrost]. Assuming the liquid-vapor interface beneath a Leidenfrost drop is perturbed with axisymmetric, sinusoidal variations, as illustrated in Fig. 3, then the shape function of the new drop interface, $S(r)$ , can be expressed in cylindrical coordinate as:

[TABLE]

where $e$ is the mean vapor layer thickness, $\epsilon$ is the perturbation amplitude, and $\lambda_{c}$ is the capillary wavelength generated at the bottom of the drop (see Fig. 13). We also assume that $R_{s}\gg e\gg\epsilon$ , implying a small perturbation to the equilibrium shape of the drop [duchemin2005static, lister2008shape]. For a stable interface, the pressure at the drop center should be smaller or equal to the pressure at $r\approx\lambda_{c}/2$ in order to drive liquid back to the center. If the pressure in the vapor layer is constant, then to leading order, the pressures at $r$ = 0, and $r$ = $\lambda_{c}/2$ are:

[TABLE]

where the height difference between the two points is $\Delta z=S(r)_{r\rightarrow 0}-S(r)_{r\rightarrow\lambda_{c}/2}=2\epsilon-\lambda_{c}^{2}/8R_{s}$ . The condition for stability can be found by equating Eqs. 2 and 3, leading to an expression for the wavelength:

[TABLE]

Since the perturbation amplitude and radius of curvature always appear as a product, Eq. 4 can also be written as:

[TABLE]

where $\chi=\epsilon R_{s}/l_{c}^{2}$ . For finite $\epsilon$ , in the limit $R_{s}\rightarrow\infty$ , Eq. 4 simplifies to $\lambda_{c}/2\approx 3.85l_{c}$ , in agreement with the prediction in Biance et al. [biance2003leidenfrost] for a flat surface. However, the addition of a curved surface couples the perturbation amplitude to the radius of curvature. The quantity $2\chi-3\pi^{2}$ must be positive in order for $\lambda_{c}$ to be real, and thus represent the condition for instability. Figure 4 shows the plots of Eq. 5, in which the red curve shows the real part, whereas the blue curve shows the imaginary part. As we can see from Fig. 4, on a curved surface, the perturbation amplitude must satisfy $\chi\gtrsim 14.8$ , i.e., $\epsilon\gtrsim 14.8l_{c}^{2}/R_{s}$ in order to lead to the Rayleigh-Taylor instability. In this sense, we can conclude that the surface curvature suppresses the Rayleigh-Taylor instability.

Although the curved surface suppresses the Rayleigh-Taylor instability, it still has an affect on the overall drop shape. We investigated the influence of surface curvature on the drop shape by depositing water drops of given volumes $V$ on a curved surface with $T_{s}$ = 623 K. We measured the radius $R$ of the drops from recorded images taken immediately after deposition prior to the onset of oscillations. As can be seen in Fig. 5, stable Leidenfrost water drops with $R\approx 11l_{c}\approx$ 28 mm are obtained on the curved surface, while on a flat substrate, the maximum radius of a stable Leidenfrost water drop is $R\approx 3.95l_{c}\approx$ 10 mm ( $l_{c}$ = 2.5 mm for Leidenfrost water drops, see Table 1) [snoeijer2009maximum, biance2003leidenfrost, burton2012geometry], suggesting that surface curvature plays a crucial role in suppressing the Rayleigh-Taylor instability. The dashed blue line in Fig. 5 shows the theoretical prediction for a drop on a flat surface by solving the Young-Laplace differential equation numerically and assuming that the contact angle of the drop is $180^{\circ}$ [burton2010experimental, burton2012geometry]. The experimental results show excellent agreement with the theoretical prediction, indicating that the shape of the Leidenfrost drops is not strongly affected by the surface curvature. More specifically, when $R<l_{c}$ (2.5 mm for water), the drops are quasi-spherical and thus $R\propto V^{1/3}$ , whereas $R\propto V^{1/2}$ for $R>l_{c}$ , which is expected for puddle-like drops with constant thickness.

However, one can also notice that most of the experimental data is slightly less than the corresponding theoretical prediction, which we attribute to two possible reasons. First, for small drops, the evaporation begins from the moment of deposition on the hot surface, which removes a small amount of water prior to imaging. Second, large drops are thicker in the center due to the underlying curved substrate, leading to a smaller apparent radius for a given volume. Thus, although surface curvature and evaporation may somewhat reduce the drop size, the effects seem to be a minor influence on the overall drop shape.

2 Breathing mode of small Leidenfrost drops

Caswell [caswell2014dynamics] experimentally characterized axisymmetric oscillations in the radius of small Leidenfrost drops using interference imaging. The drops displayed small-amplitude changes in the radius of the flat region near the surface. Due to volume conservation, an increase (decrease) in the drop radius leads to a decrease (increase) in the thickness of the drop. Caswell [caswell2014dynamics] found that the oscillation frequency of the breathing mode, $f_{b}$ , obeyed a distinct power law, $f_{b}\propto R_{0}^{-0.68\pm 0.01}$ , where $R_{0}$ is the average drop radius during the oscillation. This dependence is distinctly different than the expected three-dimensional dispersion relation for inviscid spherical drops, $f\propto R_{0}^{-3/2}$ [rayleigh1879capillary]. Here we provide an analytical model which explains this contrast and fits the experimental data with no adjustable parameters.

To leading order, and due to the axisymmetry of the breathing mode, we model a Leidenfrost drop as an incompressible liquid cylinder of volume = $V$ and a time dependent radius, $R(t)$ . We will assume that the bottom of the cylinder is fixed at $z=0$ , which is reasonable if the thickness of the vapor layer varies much less than the radius. In cylindrical coordinates, $(r,\phi,z)$ , the simplest form for the velocity which satisfies $\nabla\cdot\vec{\bf v}=0$ and the boundary conditions $\vec{\bf v}\cdot\hat{\bf r}|_{r=R(t)}=R^{\prime}(t)$ and $\vec{\bf v}\cdot\hat{\bf z}|_{z=0}=0$ is:

[TABLE]

For simplicity, we will write $R(t)=R$ , and the time derivative of $R$ as $R^{\prime}(t)=R^{\prime}$ . The total kinetic energy of the drop is then:

[TABLE]

Evaluating the integrals, and using the fact that $V=\pi R^{2}h$ , where $h$ is the time-dependent height of the cylinder, we obtain:

[TABLE]

The total potential energy is the sum of gravitational potential energy and surface energy:

[TABLE]

The equilibrium drop radius, $R_{0}$ , can be found by minimizing the potential energy with respect to $R$ :

[TABLE]

For large drops where $R\gg h$ , the first and third terms must balance. Equating these two terms and using the fact that $V=\pi R^{2}h$ , we see that:

[TABLE]

which agrees with the expected asymptotic thickness of large Leidenfrost drops [biance2003leidenfrost]. We may now define the Lagrangian of the system as $L=T-U$ , and apply the Euler-Lagrange equation to obtain a differential equation for $R$ :

[TABLE]

To proceed further, we will linearize the equation by considering only small oscillations of the radius, $R=R_{0}(1+\epsilon e^{i\omega t})$ , where $R_{0}$ is equilibrium drop radius found by solving Eq. 10, $\epsilon$ is the perturbation amplitude, and $\omega$ is the angular frequency of the oscillation. Assuming that $\epsilon\ll 1$ , to leading order Eq. 12 reduces to an expression for the angular frequency:

[TABLE]

where $\Gamma\equiv\sqrt{\gamma^{2}+4\gamma\rho_{l}gR_{0}^{2}}$ and we have substituted $V$ for the equilibrium radius $R_{0}$ using Eq. 10.

Taking typical values of the parameters in Eq. 13 for water, we can plot the oscillation frequency of the breathing mode, $f_{b}=\omega/2\pi$ , as a function of $R_{0}$ , which is shown by the solid blue curve in Fig. 6a. The theory shows excellent agreement with the data. For comparison, the power law $f_{b}\propto R_{0}^{-0.68}$ found by Caswell [caswell2014dynamics] is plotted as the dashed red line in Fig. 6a. A closer view of the comparison between the prediction by Eq. 13 and $f_{b}\propto R_{0}^{-0.68}$ can be seen in Fig. 6b. The data lie in the transition regime between two asymptotic limits, $R_{0}\ll l_{c}$ and $R_{0}\gg l_{c}$ , which may explain the anomalous power law reported by Caswell [caswell2014dynamics]. In these two limits, Eq. 13 reduces to:

[TABLE]

The scaling for small drops is independent of $g$ , which is expected based on their nearly-spherical shape. For large drops where $h\approx 2l_{c}$ (Eq. 11), the scaling is determined essentially by gravity: $\omega\sim\sqrt{2gh}/R_{0}$ .

3 Star-shaped oscillations of Leidenfrost water drops

In addition to the breathing mode, large Leidenfrost drops may also develop azimuthal, star-shaped oscillations. Similar oscillations have been observed in a variety of systems involving liquid drops [brunet2011star]. In our experiments with Leidenfrost water drops using a curved substrate, we observe star-shaped oscillation modes with $n$ = 2 to 13 lobes along the drop periphery at $T_{s}=623$ K. A typical star-shaped oscillation mode ( $n=11$ ) is shown in Fig. 7a. Figures 7b and 7c show the drop oscillation wavelength $\lambda_{d}$ and frequency $f_{d}$ for different $R$ , respectively, which are measured by analyzing the high-speed videos of the oscillating drops. Surprisingly, both $\lambda_{d}$ and $f_{d}$ remained nearly constant as we varied $R$ . Increasing $R$ thus led to an increase in the allowable number of lobes $n$ , as indicated by the red arrows. This similar trend also applies to other liquids used in our experiments, which will be discussed in Section 5.

For the free oscillations of an incompressible, axisymmetric spherical drop with infinitesimal deformations, the natural resonance frequency, $f_{n}$ , of the $n$ th-mode is given by [rayleigh1879capillary]:

[TABLE]

For a liquid puddle with $R\gg h\approx$ 2 $l_{c}$ , where $h$ is the thickness of the liquid puddle [biance2003leidenfrost], the resonance frequency $f_{n}$ takes the form [yoshiyasu1996self]:

[TABLE]

The first term under the square root is for a strictly two-dimensional drop, whereas the correction factor, $1/\sqrt{1+(2-\frac{\pi}{2}+\frac{n-3}{4})\frac{l_{c}}{R}}$ , is due to the quasi-two-dimensional nature of the puddle.

As shown in Fig. 7, for water drops where we observed $n$ = 2 modes, the average radius was $R\approx$ 4 mm, which is greater than $l_{c}$ . Hence, it is reasonable to use Eq. 16 to predict the oscillation frequency for all star-shaped modes in our experiments. The comparison between theory and experiment is shown in Fig. 7c, which indicates an excellent agreement. For smaller values of $n$ , both $\lambda_{d}$ and $f_{d}$ vary non-monotonically with $R$ . There are a few potential reasons for this. One possibility is that nonlinear effects play a more significant role for smaller values of $n$ because the ratio of the mode amplitude to the drop radius is larger [becker1991experimental, smith2010modulation]. However, the excellent agreement with the linear theory, i.e. Eq. 16, suggests that nonlinear effects may not be important.

Instead, we suggest that the behavior may be related to the mode selection mechanism. For a given radius, either the frequency or the wavelength is preferentially selected by the excitation mechanism. Once $f_{d}$ or $\lambda_{d}$ is selected, Eq. 16 will determine the other. As will be shown in Section 6, there are strong capillary waves excited at the liquid-vapor interface beneath the drop. These waves produce pressure oscillations which parametrically couple to the star-shaped modes. A determination of the expected behavior of $f_{d}$ on $R$ is thus complicated by the physics of the vapor flow beneath the drop, nevertheless, we would expect the effects to be less variable for $R\gg\lambda_{d}$ , which agrees with the data shown in Fig. 7.

4 Pressure oscillations in the vapor layer

It has been well-documented that star oscillations of large liquid drops can be initiated by a parametric forcing mechanism [brunet2011star]. For instance, when a liquid puddle is deposited onto a superhydrophobic, vertically-vibrating substrate with a prescribed frequency, then the equation of motion of the drop can be described by an equation similar to Mathieu’s equation, and star oscillations will be excited when the drop oscillation frequency is approximately half that the excitation frequency [yoshiyasu1996self, brunet2011star]. In the experiment, we observed star oscillations of Leidenfrost water drops (e.g. see Fig. 7a) in the absence of an external excitation. Therefore, we hypothesize that the star oscillations are driven by fluctuations in the vapor layer pressure that supports the drop.

In order to explore the dynamics of pressure fluctuations in the vapor layer, we utilized a pressure sensor attached to the center of the curved substrate, schematically shown in Fig. 2b. We measured the pressure variations in the vapor layer of water drops at $T_{s}$ = 623 K during the initiation of the star-shaped oscillations. We started recording the pressure immediately after the drops were placed on the substrate, and stopped after the well-defined star-shaped oscillations had been initiated for several seconds. The whole process lasted about 18.5 s.

The snapshots on the top panel of Fig. 8a show images of an $n$ = 4 drop at different stages during the oscillation. More specifically, the drop shape was nearly circular immediately upon placement on the substrate ( $t$ = 0 s), then the pressure fluctuations in the vapor layer increased as the drop shape developed well-defined lobes. The amplitude continued to grow with time ( $t$ =12 s) until the formation of a steady-state, star-shaped oscillation with large amplitude ( $t$ = 16.1 s).

Figure 8a shows the pressure variations in the vapor layer during the whole process. The mean pressure required to support the drop is $\rho_{l}gh\approx 2\rho_{l}gl_{c}=47$ Pa. However, the mean pressure we measure is slightly larger than this value because the pressure is measured at the center of the substrate where the pressure is larger in order to drive the viscous vapor to flow out to the drop edge. The pressure fluctuations around the mean increase with time until the formation of a steady star-shaped oscillation with a large amplitude (about $t$ = 14 s).

In order to gain insight into the underlying relationship between the oscillations of the drop and pressure in the vapor layer, we performed a Fast Fourier Transform on the pressure data in four different time intervals: 0-1 s, 5-6 s, 12-13 s, and 16-17 s. The results are shown in Figs. 8b-8e. Initially, the pressure remains nearly constant, and there are no sharp peaks in the power spectrum (Fig. 8b). Between $t$ = 5-6 s, the pressure fluctuations become stronger and more periodic, and several small peaks are visible (Fig. 8c). Then, as the star oscillation is further developed, a sharp peak which is located at $\approx$ 28 Hz shows up in the power spectrum (Fig. 8d). Finally, when the star oscillation is fully developed, the pressure fluctuations stop growing, and a sharp peak located at $\approx$ 28 Hz dominates the spectrum (Fig. 8e).

It is interesting to note that the location of the sharp peak is approximately twice the oscillation frequency of a fully-developed, $n$ = 4 mode Leidenfrost water drop (see Fig. 7c). This is consistent with a parametric forcing mechanism for the excitation of the star oscillations. Generally, for an oscillator with time-dependent frequency, the oscillation amplitude exponentially increases with time when the driving frequency is approximately twice that the natural frequency, which is known as parametric forcing. A well-known example is the parametric pendulum with a natural frequency of $\sqrt{g/l}/(2\pi)$ , and parametric oscillation occurs when the vertical frequency $\sqrt{k/m}/(2\pi)$ of the center of mass of the pendulum is twice that natural frequency, i.e., $\sqrt{k/m}/(2\pi)=\sqrt{g/l}/\pi$ [landau1976course], where $l$ and $m$ are the length and mass of the pendulum, respectively, whereas $k$ is the spring constant of the pendulum. Moreover, harmonics at higher frequencies are also visible as indicated by the secondary peaks in Fig. 8e. This is likely due to nonlinear effects involved in the star-shaped oscillations. In addition, the location of the peak is robust. As will be shown in Section 7, the location of the dominant peak in the power spectrum is mostly independent of the substrate temperature and the environmental temperature.

Figure 9a shows the pressure oscillation frequency $f_{p}$ in the vapor layer of Leidenfrost water drops for all of the observed modes. The data points from left to right correspond the oscillation modes from $n$ = 2 to 13, as indicated by the red arrows. By comparing this with the drop oscillation frequency $f_{d}$ (Fig. 9b), we find $f_{p}\approx 2f_{d}$ as indicated by the dashed blue and red lines in Figs. 9a and 9b, respectively. This robust relationship suggests that the star-shaped oscillations are parametrically driven by the pressure [miles1990parametrically, kumar1994parametric, yoshiyasu1996self, brunet2011star, terwagne2011tibetan]. Consider the radial position of a point on the perimeter of the drop during a star-shaped oscillation, $r(t)$ . This point will oscillate in time due to the azimuthal standing wave, and obey the following equation:

[TABLE]

To leading order, $\omega$ is the resonant frequency of the mode, Eq. 16 can be written as:

[TABLE]

For simplicity, we have ignored the correction due to the quasi-two-dimensional nature of the drop. The pressure variations will induce vertical oscillations of the drop, and in turn oscillations of the drop radius:

[TABLE]

where $R_{0}$ is the average radius of the drop, and $\epsilon$ is the amplitude of the small perturbation. Plugging Eq. 19 into Eq. 18, to leading order in $\epsilon$ we obtain $\omega^{2}=\omega_{0}^{2}(1-3\epsilon\cos\omega_{p}t)$ . Therefore, Eq. 17 now becomes:

[TABLE]

When $\omega_{p}\approx 2\omega_{0}$ , the parametric resonance is the strongest, i.e. the oscillation amplitude will exponentially increase with time until the star-shaped oscillations are fully developed. This is consistent with our measurements of the pressure oscillation frequencies in the vapor layer (Fig. 9), which are approximately twice that the drop oscillation frequencies (Fig. 7c) that we observed in our experiment. We note that it is possible that the pressure oscillations are a consequence of the star oscillations. If the star oscillations were excited by some other mechanism, then the symmetry of the star shape would necessitate a pressure extrema when the star shape reaches its maximum amplitude. However, the only energy which is available to drive the oscillation comes from the evaporation and gas flow in the vapor layer. Thus we expect that the pressure oscillations are the source of the star oscillations. Section 6 will explore the potential source of the pressure oscillations.

5 Star-shaped oscillations of different liquids

In addition to water, we performed similar experiments with five other liquids: liquid N2, acetone, methanol, ethanol, and isopropanol. For these liquids, the curved surfaces were also machined to satisfy $l_{c}/R_{s}\approx 0.03$ (see Section 2), and the relevant physical properties and substrate temperatures are listed in Table 1. Figure 10 shows snapshots of star-shaped oscillations for the six different liquids we used in the experiments, where the other five liquids show similar star-shaped patterns to water drops but with different observable oscillation modes. By analyzing the images of the star-shaped oscillations, we find that star-shaped oscillations of acetone, methanol, ethanol, and isopropanol share a similar wavelength $\lambda_{d}\approx 0.9$ cm and frequency $f_{d}\approx 14$ Hz, whereas liquid N2 shows a wavelength of $\lambda_{d}\approx 0.6$ cm and frequency $f_{d}\approx 17$ Hz. This suggests that the star oscillations may depend on the capillary length of the liquid.

Figure 11 shows the dependence of star-shaped oscillation mode $n$ , rescaled wavelength $\lambda_{d}/l_{c}$ , and rescaled frequency $f/\sqrt{g/l_{c}}$ on the rescaled drop radius $R/l_{c}$ of the six liquids used in the experiment. The mode number $n$ increases linearly with $R/l_{c}$ as indicated by the dashed black line in Fig. 11a, showing that large modes can only be observed in larger drops, as supported by Fig. 10. The rescaled wavelength $\lambda_{d}/l_{c}$ for the six liquids collapses onto a straight line as shown in Fig. 11b, which indicates that the wavelength of the star oscillations only depends on $l_{c}$ of the liquid. Figure 11c shows the oscillation frequency of both the azimuthal star oscillations, $f_{d}$ , and the corresponding pressure, $f_{p}$ , in the vapor layer. The solid symbols denote the drop oscillations, whereas open symbols with the same color represent the corresponding pressure oscillations in the vapor layer for a specific liquid. Both $f_{d}$ and $f_{p}$ collapse fairly well, and the relationship $f_{p}\approx 2f_{d}$ is consistent with the data, as indicated by the dashed black and red lines. Additionally, by comparing different values of $T_{s}$ for water, Fig. 11 shows that $\lambda_{d}$ , $f_{d}$ , and $f_{p}$ are mostly independent of substrate temperature.

The main influence of $T_{s}$ is the through the number of observable modes. For instance, at $T_{s}$ = 673 K, we did not observe oscillation modes $n$ = 12 and $n=13$ in water, which we attribute to the extremely fast evaporation rate of the liquid, as will be discussed in Section 6. Although the robust relationship between $f_{d}$ and $f_{p}$ suggests that the star oscillations of the drops are parametrically driven by the pressure in the vapor layer, there is an interesting variation in both $f_{d}$ and $f_{p}$ at small mode numbers. The data for smaller modes (e.g. $n=2,3,4$ ) deviates further from the average compared to larger modes. This behavior may be nonlinear in nature since the amplitude of the mode relative to the drop radius is much larger for small $n$ [becker1991experimental, smith2010modulation]. In addition, the quasi-two-dimensional nature of the drop is more important at small $n$ , which may also contribute to the variations in this regime. A more quantitative explanation for this dependence is left for future studies.

As shown in Table 1, different oscillation modes are observed in different liquids. In particular, some liquids failed to display smaller mode numbers. This may be due to inherent viscous damping that prevents the sustained excitation of large-amplitude oscillations. Following the analysis in Ma et al. [Maleidenfrost2016], the role of damping can be characterized by the Reynolds number associated with the liquid oscillation. For a star oscillation, the characteristic length scale and time scale are $l_{c}$ and $\sqrt{l_{c}^{3}\rho_{l}/\gamma}$ , respectively. Thus, the inertial term in the Navier-Stokes equation is estimated as $|\rho_{l}(\vec{\bf v}\cdot\nabla)\vec{\bf v}|=\gamma/l_{c}^{2}$ , and the viscous term is $|\eta_{l}\nabla^{2}\vec{\bf v}|=\eta_{l}/\sqrt{l_{c}^{5}\rho_{l}/\gamma}$ . The Reynolds number of the oscillating liquid is defined as the ratio of inertial to viscous terms: $Re_{l}=\sqrt{l_{c}\rho_{l}\gamma}/\eta_{l}$ , and the values of $Re_{l}$ for different liquids are listed in Table 1.

Figure 12 shows the dependence of the minimum mode number $n_{min}$ on $Re_{l}$ for each liquid. The dashed line represents a suggested scaling $Re^{-1}$ , indicating that the liquid viscosity damps the oscillations of smaller drops and thus sets the minimum mode number of stable star oscillations. However, the substrate temperature and thus evaporation rate will likely affect the number of observed modes as well. A higher evaporation rate will induce a stronger driving of the oscillation modes, which would reduce $n_{min}$ . This can be seen in Fig. 12 for liquid N2, which undergoes much more rapid evaporation due to the large difference between the boiling and substrate temperature, $T_{s}-T_{b}$ . In addition, as mentioned previously, the stronger evaporation may also inhibit larger modes if the pressure oscillations are less coherent. In the next section we show how the flow in the vapor layer is linked to the pressure oscillations.

6 Origin of pressure oscillations in the vapor layer

As evidenced by Fig. 11c, the star-shaped oscillations of Leidenfrost drops are parametrically driven by the pressure variations in the vapor layer. Star oscillations induced by a parametric coupling have been observed in a variety of systems where variations of the drop radius are induced by external, periodically-modulated fields [brunet2011star]. However, in our experiment, there are no obvious external fields. It is then crucial to understand the source of the pressure oscillations. In this section we show how capillary waves of a characteristic wavelength at the liquid-vapor interface induce the pressure oscillations in the vapor layer.

Figure 13 shows s sketch of a large, axisymmetric Leidenfrost drop with a maximum radius = $R$ and thickness $\approx 2l_{c}$ . The mean vertical velocity of the gas at the liquid surface is $v$ , $e$ is the mean thickness of the vapor layer, and $u$ is the outward radial velocity of the gas near $r=R$ . For such a large Leidenfrost drop, both the bottom surface of the drop and the substrate surface are assumed to be approximately flat to simplify the analysis. Following the model of Biance et al. [biance2003leidenfrost], the mass loss rate of the drop due to evaporation can be expressed as:

[TABLE]

where $\kappa_{v}$ is the thermal conductivity of the vapor, $L$ is the latent heat of the evaporation, $\Delta T=T_{s}-T_{b}$ , and $\rho_{v}$ is the density of the vapor at the liquid interface. Since $e\approx$ 100 $\mu$ m, and thus $R\gg e$ [burton2012geometry], we employ lubrication theory, leading to an expression relating the flow rate to the mean pressure, $P$ , in the vapor layer:

[TABLE]

where $\eta_{v}$ is the dynamic viscosity of the vapor. We note here that parameters such as $\kappa_{v}$ , $\rho_{v}$ , and $\eta_{v}$ will vary in the vapor layer due to the temperature gradient between the liquid and solid surfaces. For simplicity, we will assume that these represent average values of the vapor layer properties. A more detailed treatment, possibly including full numerical simulations, would be necessary to extend this simplified model.

In addition, due to mass conservation, the volumetric flow rate of gas from the bottom of the drop must be equal to the flow rate exiting the perimeter of drop:

[TABLE]

Solving Eqs. 21, 22, and 23 for $u$ , $v$ , and $e$ yields:

[TABLE]

In this model, the steady-state, linear temperature profile in the vapor layer is valid since the characteristic time scale associated with thermal diffusion, $e^{2}/D_{v}\approx$ 1 ms, is smaller than the typical residence time of the gas in the vapor layer, $R/u\approx$ 10 ms, where $D_{v}$ is the thermal diffusivity of the vapor. Let us assume a water Leidenfrost drop with $R$ = 0.01 m, and for the vapor, we use the properties of steam at 373 K and 1 atm of pressure: $\rho_{l}$ = 959 kg/m3, $\rho_{v}$ = 0.45 kg/m3, $\gamma$ = 0.059 N/m, $g$ = 9.8 m/s2, $L$ = 2.26 $\times 10^{6}$ J/kg, $\eta_{v}$ = 1.82 $\times$ 10*-5* Pa s, and $\kappa_{v}$ = 0.04 W/m/K. Although the properties of the vapor may vary somewhat in the vapor layer due to the temperature gradient, we have verified that this does not significantly affect our analysis.

We can now plot $u$ , $v$ , and $e$ with respect to $\Delta T$ based on Eqs. 24, 25 and 26. The results are shown in Figs. 14a, 14b, and 14c, respectively. The inertial term in the Navier-Stokes equation is estimated as $|\rho_{v}(\vec{\bf v}\cdot\nabla)\vec{\bf v}|=\rho_{v}v^{2}/e$ , and the viscous term as $|\eta_{v}\nabla^{2}\vec{\bf v}|=\eta_{v}v/e^{2}$ . Thus the Reynolds number in the vapor layer is $Re_{v}=\rho_{v}ve/\eta_{v}$ . Plugging in Eqs. 25 and 26, we arrive at a surprisingly simple expression for the Reynolds number in the vapor: $Re_{v}=\Delta T\kappa_{v}/L\eta_{v}$ , which is independent of the drop size. This expression is plotted in Fig. 14d. Under typical experimental conditions, $Re_{v}\approx 0.2$ , suggesting that the original lubrication flow assumption is valid, although inertial forces are non-negligible.

Based on Fig. 14, we can estimate the Bernoulli pressure in the vapor layer of Leidenfrost water drops as $\rho_{v}u^{2}/2\approx 1$ Pa, which is much smaller than the pressure variations $\approx$ 10 Pa (see Fig. 8a), suggesting that the inertial force does not account for the pressure variations. The viscous pressure is $P_{\mathrm{vis}}\sim\eta_{v}vR^{2}/e^{3}$ , thus the pressure variation due to the local variations of the film thickness is $\Delta P_{\mathrm{vis}}=\frac{\mathrm{d}P_{\mathrm{vis}}}{\mathrm{d}e}\Delta e$ . Then $\Delta P_{\mathrm{vis}}=$ 10 Pa corresponds to $\Delta e\approx$ 15 $\mu$ m considering the typical values of $v$ (Fig. 14c) and $R$ . Although this local film thickness could also be possibly induced by the vertical motion of the center of mass of the drops, the fact that small and large drops (see Fig. 11c) share a nearly constant oscillation frequency suggests this is not the case. Therefore, the pressure variations in the vapor layer are likely to be induced by the local variations of the vapor film thickness.

We propose that capillary waves with a characteristic wavelength, $\lambda_{c}$ , traveling from the center to the edge of the drop lead to local variations in vapor film thickness, and thus pressure variations in the vapor layer, as schematically shown in Fig. 13. To confirm this possibility, we used a heated, plano-concave, fused silica lens as a substrate to image the capillary waves from below, as schematically shown in Fig. 15. The results for acetone and ethanol Leidenfrost drops are shown in Fig. 16. The images were produced by averaging all frames in a given video sequence, then subtracting this background image in order to enhance contrast in the center of the drop. The white “halo” surrounding the drop is a consequence of this subtraction process. The Fourier spectra for the capillary waves in the central region of the drops were computed in both time and space using the pixel intensity as the signal. Using the Nyquist sampling theorem, the maximum frequency was limited by half the video frame rate (500-1000 Hz), and the maximum spatial frequency was set by half the camera magnification (132 px/cm).

We observed a sharp peak in the measured Fourier spectrum for all the analyzed video sequences. In Figs. 16a and 16b (8-mode, acetone), this peak is located at a capillary wave frequency $f_{c}\approx$ 24 Hz and within a range of wave numbers, $k_{c}=2\pi/\lambda_{c}$ , from 7 cm*-1* to 12 cm*-1*, where $\lambda_{c}$ is the wavelength of the capillary waves. Figures 16c and 16d show the capillary waves beneath a 11-mode Leidenfrost ethanol drop and a sharp peak at a capillary wave frequency $f_{c}\approx$ 25 Hz with $k_{c}\approx$ 7-10 cm*-1*. Similarly, for the capillary waves beneath a 12-mode Leidenfrost ethanol drop (Fig. 16e), a sharp peak exists at $f_{c}\approx$ 27 Hz and $k_{c}\approx$ 10-11 cm*-1* as shown in Fig. 16f. These frequencies show excellent agreement with the corresponding typical pressure oscillation frequencies measured in the vapor layer of star-shaped Leidenfrost ethanol and acetone drops (see Fig. 11c).

It is well-known that capillary waves can be generated at a liquid-vapor interface due to a strong shear stress in the vapor [Miles1957, Zhang1995, Paquier2015, Zeisel2008, chang2002complex]. A similar mechanism underlies the Kelvin-Helmholtz instability, however, the intermediate Reynolds number in the Leidenfrost vapor layer complicates the analysis. Nevertheless, we can estimate the strength of this shear stress. Typically, the “friction velocity” is generally used to measure the strength of shear, which is defined as $u_{*}=\sqrt{\tau/\rho_{v}}$ , where $\tau$ is the shear stress at the liquid-vapor interface. The maximum shear stress at the interface is:

[TABLE]

assuming a parabolic-flow profile in the vapor layer with mean velocity $u$ near the edge of the drop. Using water as an example, and plugging Eqs. 24 and 26 and the value of $\eta_{v}$ of water vapor into Eq. 27 yields $u_{*}\approx 2$ m/s. This friction velocity is quite strong and sufficient to lead to the growth of unstable modes with wavelengths of millimeter scale [Zhang1995, Zeisel2008].

The general dispersion relation of gravity-capillary waves with a dense upper-layer is:

[TABLE]

where $k_{c}$ = 2 $\pi$ / $\lambda_{c}$ , and $h\approx 2l_{c}$ is the thickness of the drop. For simplicity, we have assumed that the normal velocity is zero at the upper surface of the drop. For a large Leidenfrost water drop whose characteristic pressure oscillation frequency in the vapor layer is $f_{p}\approx$ 26 Hz (see Fig. 9), the corresponding capillary wavelength is calculated to be $\lambda_{c}\approx 3.03l_{c}$ , so that $k_{c}\approx$ 8.3 cm*-1* [Maleidenfrost2016]. Similarly, for the Leidenfrost acetone and ethanol drops shown in Figs. 16a, 16c, and 16e, the corresponding $k_{c}\approx 13$ cm*-1*, which is slightly larger than the positions of the peaks indicated in Figs. 16b, 16d, and 16f, where $\lambda_{c}\approx 4l_{c}$ . The agreement between the estimate for $k_{c}$ and the measurements shown in Fig. 16 is good considering the simplicity of Eq. 28, which is derived using an inviscid, semi-infinite flow in both phases.

The capillary-wave origin for the star-shaped oscillations also agrees with the minimum size of the $n$ = 2 mode. More specifically, the radius we measured in our experiments is $R$ rather than $R_{\mathrm{neck}}$ as illustrated in Fig. 13. The relationship between these two lengths is $R=R_{\mathrm{neck}}+0.53l_{c}$ [snoeijer2009maximum, burton2012geometry]. Thus, the minimum drop size required to fit one capillary wavelength beneath the drop is $2R_{\mathrm{neck}}\approx\lambda_{c}$ . Using $\lambda_{c}\approx 3.03l_{c}$ from above, we find that the radius of an $n$ = 2 mode drop should be $R\approx 2.05l_{c}$ . This is in good agreement with Fig. 11a, which shows $R/l_{c}$ is slightly less than 2 for the smallest drops. Taken together, this analysis suggests a purely hydrodynamic origin for the star oscillations based on capillary waves generated by a strong shear stress in the rapidly-flowing vapor beneath the drop.

7 Thermal effects

In the Leidenfrost effect, thermal energy is transferred from the substrate to the drop, inducing evaporation and enabling the sustained star-shaped oscillations. Here we consider the role of thermal effects, such as convection in the liquid, which may play a role in the star oscillations. At sufficiently large temperature gradients in liquid layers, convective structures can develop with well-defined length scales. A well-known example is Bénard-Marangoni convection, where hexagonal patterns are generated in a thin layer of liquid when heated from the bottom [benard1901tourbillons, rayleigh1916lix]. The convection pattern is strongly affected by the variation of surface tension with temperature [marangoni1871ausbreitung, schatz2001experiments, maroto2007introductory]. Given that our data collapses when scaled by the capillary length of the liquid (Fig. 11), this particular type of convection may be important for Leidenfrost drops since the capillary length also appears in the length scale which characterizes the size of the convective patterns.

Generally, the Marangoni number, $Ma$ , and Rayleigh number, $Ra$ , are used to characterize Bénard-Marangoni convection. These dimensionless numbers are defined as:

[TABLE]

These numbers characterize the strengths of thermocapillary and buoyancy effects, respectively. Additionally, the ratio of thermal conduction inside the liquid to the conduction at the liquid-vapor interface is characterized by the Biot number:

[TABLE]

The definitions of the symbols in $Ma$ , $Ra$ , and $Bi$ are listed in Table 2, as well as their values for water at the boiling point.

Figure 17 shows the temperature profile of a large Leidenfrost water drop levitating on its vapor layer with $T_{s}=673$ K, in which the substrate surface is defined as $z=0$ . The temperature was measured with a fine-point thermocouple, as described in Section 2. We define the temperatures of the bottom and top surfaces of the drop as $T_{\mathrm{hot}}$ and $T_{\mathrm{cold}}$ , respectively, and the temperature of the position which is slightly above the top surface of the drop is denoted as the environment temperature, $T_{1}$ , as shown in Fig. 13. From Fig. 17 we can obtain $\delta T$ = $T_{\mathrm{hot}}-T_{\mathrm{cold}}\approx$ 3.8 K, the environment temperature measured to be $T_{1}=$ 367 K, thus we can calculate $Ma\sim 2\times 10^{5}$ , $Ra\sim 7\times 10^{4}$ , and $Bi\approx$ 1.1. Both the values of $Ma$ and $Ra$ are much larger than their critical values, which are typically of order 100 for the initiation of convective instability [maroto2007introductory, leal2007advanced]. Thus it is possible that thermal convection plays a role in initiating the star-shaped oscillations. In this case, one may expect the star-shaped wavelength, $\lambda_{d}$ , to be related to the critical wavelength of the convective instability. This critical wavelength depends on the Biot number. We implemented a qualitative test for this dependence by wrapping aluminum foil around the substrate and Leidenfrost drop, which dramatically increased the environment temperature, $T_{1}$ , near the top of drop.

Figures 18a and 18d show Fourier power spectra of the pressure oscillations in the vapor layer of Leidenfrost water drops during $n$ = 4 and $n$ = 5 oscillations at $T_{1}\approx$ 367 K, respectively. Figures 18b and 18e are the spectra of pressure oscillations during $n$ = 4 and $n$ = 5 oscillations at $T_{1}\approx$ 483 K, respectively. Finally, Figs. 18c and 18f represent the same modes at $T_{1}$ = 623 K. For $T_{1}$ = 483 K and 623 K, the direction of heat transfer at the upper surface of the drop has been reversed; energy is added to the drop. Overall, dramatically increasing the surrounding temperature affects the appearance of higher harmonics in the spectra. The behavior is non-monotonic, and is likely due to the highly nonequilibrium conditions (high evaporation rate) induced by the high temperatures of the substrate ( $T_{s}$ ) and environment ( $T_{1}$ ). A more detailed understanding of the pressure oscillation spectra is left for future studies. Nevertheless, the position of the main peak in the spectra is independent of $T_{1}$ , indicating that convection and the details of thermal transport in Leidenfrost drops play a secondary role in the star-shaped oscillations.

4 Summary and outlook

Both large and small Leidenfrost drops display self-organized oscillations due to the constant input of thermal energy and continuous evaporation and flow beneath the drop. Here we have focused on radial oscillations (i.e. “breathing” mode) of small Leidenfrost drops, and the large-amplitude, star-shaped oscillations that appear in large Leidenfrost drops. We have characterized the number of observed modes for various volatile liquids, the frequency and wavelength of the oscillations, and the pressure variations in the vapor layer beneath the drops. The number of observed modes is sensitive to the properties of the liquid (see Table 1), i.e. the star-shaped oscillations of smaller Leidenfrost drops are dissipated by the liquid viscosity, which sets the minimum oscillation mode number $n_{min}$ that can be observed in experiment. The relationship between the frequency and wavelength agrees very well the quasi-two-dimensional theory proposed by Yoshiyasu et al. [yoshiyasu1996self]. The dominant frequency associated with the pressure oscillations is approximately twice the drop oscillation frequency, consistent with a parametric forcing mechanism for the star oscillations.

One of the main findings of our work is identifying the underlying cause of the pressure oscillations. By imaging the liquid-vapor interface from below the drop, and using a simplified model for the flow in the thin vapor layer, we conclude that capillary waves of a characteristic wavelength, $\lambda_{c}\approx 4l_{c}$ , lead to pressure oscillations at the experimentally measured frequency. The flow in the vapor layer is quite rapid. Near the edge of the drop, the mean radial velocity can reach 1-2 m/s or more. In the small gap ( $\sim$ 100 $\mu$ m) between the liquid and the solid surface, this flow applies a large shear stress to the liquid-vapor interface, and can easily excite capillary waves with millimeter-scale wavelengths. The dispersion relation for the capillary waves then leads to a characteristic frequency for the pressure oscillations, which in turn parametrically drive the star-shaped oscillations. Furthermore, although the vapor flow is inherently driven by evaporation and heat transfer, the substrate and surrounding temperature have little effect on the dominant frequency and wavelength of the oscillations, suggesting they are purely hydrodynamic in origin.

Although the work presented here has focused mostly on the origin of star-shaped oscillations, the coupling of the flow in the vapor layer and the liquid-vapor interface underlies a rich spectrum of dynamical phenomena observed in both Leidenfrost liquid layers and drops. In particular, of key interest is understanding the failure of the Leidenfrost vapor layer which can lead to explosive boiling. If the observed capillary waves beneath the drop act as the precipitant to vapor-layer failure, then it is possible that geometrical tailoring of the surface to be commensurate with $\lambda_{c}$ may inhibit the generation of capillary waves. In addition, patterned, ratchet-shaped substrates with wavelengths $\approx$ 1-3 mm are known to induce propulsion of small Leidenfrost drops [cousins2012ratchet, linke2006self, lagubeau2011leidenfrost], however, less is known about transport of large quantities of liquid, and the dependence on the wavelength of the surface patterns. We leave these questions open to future experiments. More generally, our results may offer insight into the direct control of oscillations in levitated drops in many other systems [haumesser2002high, duchemin2005static, paradis2005surface, ishikawa2006noncontact, lister2008shape, langstaff2013aerodynamic], for example, precise control is crucial when levitating high-temperature or harmful liquids using a gas film. We also expect our results to enhance the understanding of dynamics that couple a thin, supporting gas film, a liquid interface, and a solid surface, a scenario which occurs through forced wetting and gas entrainment in liquid coating [Xu2005, Driscoll2011, marchand2012, Kolinski2012, Liu2013, Liu2015].

Chapter 2 Polygonal desiccation crack patterns

1 Introduction

In addition to the gorgeous star-shaped oscillatory patterns in Leidenfrost drops, there also exist a large group of amazing patterns driven by drying particulate suspensions. For example, when spill a coffee drop onto a desk, the final residues of the drop after drying form a ring-shaped pattern, which is known as the coffee ring effect. The formation of the coffee ring is driven by the capillary flow drawn from the bulk of the drop to the edge as the drying rate at the drop periphery is greater than of the bulk [deegan1997capillary, deegan2000contact, deegan2000pattern]. For more concentrated particulate suspensions, the particles pack closely into a network during drying, and the particulate films are prone to form cracks if enough capillary stress is stored in the film. This capillary induced phenomena contribute largely to the diverse natural patterns, and one of the well-known category is the desiccation crack patterns [routh2013drying, goehring2015desiccation].

Desiccation crack patterns observed in natural systems span many orders of magnitude in size (Fig. 1) [ball1999self, xu2009drying, goehring2015desiccation, routh2013drying, goehring2013pattern, goehring2013evolving]. Among the large diversity of desiccation crack patterns, polygonal patterns are the most common. Typical examples include the complex crack network in dried blood [brutin2011pattern], craquelures in old paintings [giorgiutti2015striped, goehring2017drying], T/Y-shaped cracks in dried mud [goehring2010evolution], and polygonal terrain cracks [goehring2014cracking, Maarry2010mars, el2014potential]. In particulate suspensions, the formation of desiccation cracks depends on the interplay between order and disorder in granular systems, as well as mechanical instabilities initiated by local, nonequilibrium interactions between the liquid, solid, and vapor phases [routh2013drying, goehring2015desiccation]. Nevertheless, a broad range of practical applications, such as thin film coating, forensics, and controllable surface patterning rely on knowledge of the physical processes that determine crack patterns [prosser2012avoiding, hatton2010assembly, liu2016surface, nam2012patterning, zeid2013influence]. Despite numerous studies which focus on desiccation crack patterns in a diverse range of systems, a fundamental understanding of the characteristic length scales associated with polygonal crack patterns is lacking, and it is not clear if the observed patterns in both microscopic and geologic crack patterns share the same underlying mechanisms.

In the laboratory, drying particulate suspensions, both Brownian and non-Brownian, are model systems for replicating and understanding desiccation cracks in nature. For a crack to form in any material, the mechanical potential energy released during fracture must exceed the energetic cost of creating new surfaces [griffith1921vi]. In homogeneous, isotropic, elastic solids, the dynamics of fracturing have been recently characterized with exquisite detail [hutchinson1991mixed, fineberg1999instability, bouchbinder2010dynamics, bouchbinder2014dynamics, creton2016fracture]. However, the dynamics of drying-induced cracks are complicated by the lack of material homogeneity and a nonlinear relationship between stress and strain [kitsunezaki2009crack, goehring2013plasticity, goehring2015desiccation, dufresne2006dynamics, lidon2014dynamics, man2008direct, xu2010imaging, xu2013imaging]. This complexity is enhanced by the multi-phase nature of the material: liquid menisci between particles generate heterogeneous shrinkage through capillary pressure, and excess liquid-vapor surface area in the bulk of the material [dufresne2003flow, tirumkudulu2005cracking, singh2007cracking, chiu1993dryingI, chiu1993dryingII]. As a consequence, in addition to polygonal cracks, a large variety of cracks patterns have been reported in dried suspensions [allain1995regular, inasawa2012self, goehring2011wavy, nandakishore2016crack, kiatkirakajorn2015formation, neda2002spiral, vermorel2010radial, darwich2012highly, jing2012formation, xu2010imaging, giorgiutti2015striped, giorgiutti2016painting, bohn2005hierarchicalI, bohn2005hierarchicalII, pauchard2003morphologies, boulogne2013annular, allen1987desiccation, gauthier2010shrinkage]. The variability in observed patterns depends on numerous factors such as film geometry [lazarus2011craquelures, nandakishore2016crack], particle mechanics [nawaz2008effects], liquid additives [pauchard1999influence, liu2014tuning], preparation history [nakahara2006transition, nakahara2006imprinting], solvent volatility [boulogne2012effect, giorgiutti2014elapsed], and external fields [khatun2012electric, pauchard2008crack].

Despite this broad range of crack patterns, we know surprisingly little about what controls the size and hierarchy of commonly observed polygonal cracks, which are visible on both the micro- and macroscales. As shown in Fig. 1, a detailed understanding of desiccation crack patterns can lead to more accurate interpretations of planetary geomorphology, where data are limited to satellite-based imaging [Maarry2010mars, el2014potential]. Thus far, laboratory experiments and numerical models have produced contradictory results as to the mechanism and dependence of the crack spacing on desiccation conditions. For example, for regularly-spaced cracks produced by directional drying, the relationship between the crack spacing, $\lambda$ , and material thickness, $h$ , is a power law, $\lambda\propto h^{\beta}$ . Experimentally, numerous groups have reported $\beta\leq 1$ [allain1995regular, lee2004drying, smith2011effects], and various theoretical predictions give $2/3\leq\beta\leq 1$ [ma2012possible, komatsu1997pattern]. For polygonal patterns, one may expect $A_{p}\propto\lambda^{2}\propto h^{2\beta}$ . Groisman et al. [groisman1994experimental] reported $A_{p}\propto h^{2}$ in desiccated suspensions of coffee grinds, although only over a four-fold increase in $h$ . Other experimental [shorlin2000alumina] and numerical [leung2000pattern, leung2010criticality] studies have reported similar scalings. Most recently, Flores [flores2017mean] showed that $A_{p}\propto h^{4/3}$ using a model based solely on a balance of the average stress and surface energy released during cracking. Yet despite this history of investigation, there has been no systematic experimental study of the film’s thickness and mechanics on the size of polygonal crack patterns.

The pattern morphology of polygonal cracks is also of interest since it reveals information about the formation and history of the cracking process [goehring2015desiccation]. For example, the distribution of angles at crack junctions [akiba2017morphometric, shorlin2000alumina, groisman1994experimental] and statistical analysis of the correlation length in crack patterns [colina2000experimental] are commonly quantified from images of the surface. Repeated wetting and drying of the material can lead to more “Y”-shaped junctions rather than “T”-shaped junctions [goehring2010evolution, goehring2014cracking]. For some commonly used desiccation suspensions such as cornstarch-water mixtures, crack patterns with two distinct length scales can be identified, as shown in Fig. 1a. For thick samples, the smaller polygons grow into the material, resembling columnar jointing patterns often found in nature [muller1998starch, muller1998experimental, toramaru2004columnar, goehring2009drying, goehring2009nonequilibrium, akiba2017morphometric]. These smaller polygons are also known to coarsen with depth in the material [goehring2005order]. Hierarchical patterns have also been observed in the cracking glaze of ceramics [bohn2005hierarchicalI, bohn2005hierarchicalII]. It is not yet apparent why some materials display hierarchical crack patterns, and some do not.

Here we present experimental evidence which resolves many of these important, outstanding questions. Our experiments involve analysis of desiccated crack patterns in various granular materials, such as cornstarch and CaCO3, suspended in different volatile liquids: water, isopropanol (IPA), and silicone oil. We use both very thin, quasi-two-dimensional chambers, as well as open petri dishes to dry the samples. For all observed crack patterns, we find that the characteristic polygonal area is consistent with $A_{p}=\alpha h^{4/3}$ over more than three orders of magnitude in $h$ , in agreement with a recent theoretical prediction [flores2017mean] based on a balance of stress and surface energy for crack formation. This scaling is independent of the shape of the polygons, which varies considerably depending on the material-liquid combination. By characterizing the modulus of the desiccated suspension, we are also able to quantitatively predict the prefactor $\alpha$ . For all material-liquid combinations, we only observe multiscale crack patterns in cornstarch-water mixtures. We show that these cracks are due to two distinct desiccation mechanisms. Primary cracks form first due to capillary-induced shrinkage of the material. Secondary cracks form much later, and are due to deswelling of the hygroscopic cornstarch particles. Taken together, these results provide a quantitative pathway for interpreting multiscale polygonal desiccation crack patterns observed in diverse systems, from microscale colloidal films to terrestrial and extra-terrestrial planetary surfaces.

2 Experimental setup

We used commercial, polydisperse cornstarch particles from ARGO. The average radius, $R$ , of the particles is $\approx$ 5 $\mu$ m. We also used CaCO3 particles from OnlineScienceMall with $R\approx$ 1 $\mu$ m, and for some experiments, glass beads with $R\approx$ 5 $\mu$ m from Miscrospheres-Nanospheres. Particle sizes were measured using optical microscopy. Different fluids such as deionized water, low-viscosity silicone oil (0.65 cSt, ClearCo), and 99.9% isopropanol (IPA) were used as solvents to prepare particulate suspensions. For samples thicker than $h\approx$ 1 mm, we dried suspensions in polystyrene petri dishes of diameters 14 cm and 8.5 cm. Glass microscope slides were used to build thin, quasi-two-dimensional chambers as discussed in Section 5. For the thin chambers, the tunable thickness, $h$ , was set by vinyl spacers cut from a sheet. Once the vinyl spacers were placed on the edges, a sample of suspension was placed on the microscope slide, and the sample was compressed by a second glass slide and secured mechanically before gluing with optical epoxy. Similar setups have been used by previous authors [allain1995regular, dufresne2003flow, goehring2010solidification, inasawa2012self]. However, one important distinction in our experiments is that due to the relatively slow evaporation in this system, evaporation occurs nearly isotropically around the perimeter of the sample (see Figs. 10a and 10b), so drying is not uni-directional.

We used a conventional bright-field microscope to image crack patterns in the thin chambers, and a USB 3.0 digital video camera (Point Grey) connected to a macro lens to image the crack formation during drying of particulate suspensions in petri dishes from above. Recorded images were analyzed using NIH ImageJ software to obtain the area of polygonal cracks and thickness of the dried films. An electronic balance (Omega) was used to monitor the instantaneous mass of suspensions during drying. All experiments were performed at room temperature (20 ∘C) with uncontrolled relative humidity of $\approx$ 60%. Modulus measurements were obtained by slowly pressing a stainless steel ball of diameter 1.9 cm into the material using a rheometer (TA Instruments AR2000), and recording the applied normal force.

3 Results and discussion

1 Multiscale cracks in cornstarch-water suspensions

We prepared cornstarch-water suspensions and dried the samples in petri dishes. We varied the thickness of the films by controlling the initial volume of the suspensions. Above a critical thickness, $h_{c}$ , as will discussed in Section 2, we observed two distinct crack patterns that appeared at different stages of desiccation. Figure 2 shows the formation of these multiscale cracks during drying of a cornstarch-water suspension with $\phi_{i}=40\%$ . As shown in Fig. 2a, drying occurs at the air-water interface. The primary cracks (blue polygons) first appear, then as drying proceeds, secondary cracks (red polygons) appear within the larger polygons. After the suspension has dried for $t\approx$ 10.5 h, primary cracks first appear (Fig. 2c), and the number of primary cracks increases with time (Fig. 2d). At $t\approx$ 24 h, secondary cracks are visible (Fig. 2f). When $t\approx$ 42.5 h, the number of secondary cracks stops increasing though drying still proceeds (Fig. 2h). Finally, the drying contributes to the widening of the existing cracks, as shown in Fig. 2h. We note that the primary cracks penetrate completely through the sample when they form, whereas the secondary cracks grow more slowly, and their visibility increases with time.

As reported by previous authors [leung2000pattern, mizuguchi2005directional, akiba2017morphometric], the primary cracks are a result of film shrinkage induced by the Laplace pressure on the scale of the particle size. As the water evaporates, menisci form between individual particles of radius $R\approx 5$ $\mu$ m. Thus the average pressure is reduced in the suspension by $\approx\gamma/R\approx$ 15 kPa, where $\gamma=72$ mN/m is the surface tension of water. Since the suspension is partially adhered to the bottom surface of the petri dish, the cracks form almost uniformly over the sample. Without this adhesion, the suspension undergoes isotropic shrinkage, and the number of primary cracks is reduced [groisman1994experimental].

As the stress boundary condition at substrate can potentially play a role in the formation of desiccation crack patterns. We examined the effect of the boundary condition on the formation of crack patterns by drying cornstarch-water suspensions in petri dishes. We modified the bottom surfaces of the petri dishes: sandpaper sheets with grit size 120 were used to roughen the petri dish surface so as to make the surface less adhesive. This is a little bit counterintuitive since increasing surface roughness usually increases contact area and makes the surface more adhesive. However, here the boundary stress is mainly contributed by the tangential component of the adhesion, and roughening the substrate surface leads to the structures of troughs and peaks, which decrease tangential component of the adhesion such that the boundary stress decreases. A commercial coating (Rain-X) was applied in order to make the surface more hydrophobic or less adhesive. In addition, a 5-min epoxy layer was applied to the petri dish surface to enhance the adhesive properties.

We prepared cornstarch-water suspensions with the same initial volume ( $V_{i}=$ 120 ml) and the same initial volume fractions ( $\phi_{i}=40\%$ ), then deposited the suspensions into the prepared petri dishes with different boundary conditions, and the suspensions dried at room temperature. Figure 3a shows the crack patterns of the dried cornstarch-water film without any modification to the petri dish surface as a control experiment. It is evident that the increase of the surface roughness will decrease the number of large-scale cracks as indicated by Fig. 3b since increasing the surface roughness decreases the contact area between the suspension and the substrate, and the boundary stress decreases consequently, whereas the small-scale cracks are still observable as shown in the zoomed-in image of a particular region enclosed by a dashed yellow box in Fig. 3b. Similarly, as the hydrophobicity was increased (Fig. 3c), the number of large-scale cracks dramatically decreased, whereas the small-scale cracks are insensitive to the coating layer. The results shown in Figs. 3b and 3c indicate that the decrease in adhesion of the suspension to the substrate leads to dramatic decrease in the number of large-scale, primary cracks (or the increase of polygon area), which is likely due to the more isotropic shrinkage of the film during drying.

Figure 3d shows the polygonal cracks in a petri dish whose surface was coated with an epoxy resin layer. In this case, although the adhesion of the suspension to the substrate has increased, the number of large-scale cracks does not dramatically change compared to Fig. 3a. This is not surprising since shrinkage in the film will continue until the yield stress of the granular material is reached so that a crack is formed, and then this process continues. Therefore, although the boundary condition plays a key role in determining the area of the polygonal cracks, weakening the adhesion has the strongest effect, especially if the adhesive stress is smaller than the yield stress of the material. The secondary cracks are generally unaffected by the choice of boundary condition, and the origin of secondary cracks will be discussed in Section 3.

2 Critical condition for primary cracks

The appearance of primary cracks for thicker samples of cornstarch and water suspensions, as shown in Fig. 2, can be understood in terms of a well-known theory for the initiation of cracks in colloidal thin films [tirumkudulu2005cracking, singh2007cracking]. The theory assumes a no-slip boundary condition between the bottom boundary of the film and the substrate, and predicts a relationship between the critical film thickness and stress when crack should appear:

[TABLE]

where the definitions of the parameters in Eq. 1 are listed in Table LABEL:Drying_tab:1. Although the particle radius ultimately cancels from Eq. 1, it is included here so that each term is dimensionless, as in Ref. [singh2007cracking]. We have included typical values for the shear modulus of both cornstarch and CaCO3 taken from the literature [katz2013], assuming a crystalline form of CaCO3 and a Poisson’s ratio of 0.5 for cornstarch. The same values of $\phi$ and $M$ were used for all calculations.

Taking the typical values of the suspensions of cornstarch-water and CaCO3-water used in the experiment (see Table LABEL:Drying_tab:1), we can calculate the critical thickness for cracking. Since the stress driving the primary cracks is due to capillary pressure [lee2004drying], we can assume that $\sigma_{c}\approx\gamma/R$ . Solving for $h_{c}$ , we obtain $h_{c}\approx$ 1500 $\mu$ m for cornstarch, and $h_{c}\approx$ 400 $\mu$ m for CaCO3.

In order to examine whether Eq. 1 accurately predicts the critical film thickness for cracking, we dried both thin films of cornstarch-water and CaCO3-water suspensions in petri dishes, and the results are shown in Fig. 4. In Figs. 4a-4b, only small-scale, secondary cracks are observed. In Fig. 4c, the small-scale, secondary cracks are obvious and the large-scale, primary cracks are about the appear. In Fig. 4d, both large-scale, primary cracks and small-scale, secondary cracks are obvious. This indicates that the critical film thickness of cornstarch-water suspensions lies in between $h=864$ $\mu$ m (Fig. 4c) and $h=1181$ $\mu$ m (Fig. 4d). In Fig. 4g, no cracks are visible, whereas in Fig. 4h cracks appear, suggesting that the critical film thickness for CaCO3-water suspensions lies between 500 $\mu$ m and 553 $\mu$ m. Both the critical film thicknesses we observed show good agreement with the predictions by Eq. 1, though slightly different from the predicted values. However, Eq. 1 fails to explain the origin of secondary cracks in dried cornstarch-water suspensions. For the secondary cracks, there is no critical film thickness, and the cracks are visible in samples that are only a few particles thick, as will be discussed in Section 5, suggesting that the stress is not solely due to capillary pressure, as described in Ref. [lee2004drying].

3 Cornstarch particle deswelling drives secondary cracks

Multiscale crack patterns in dried cornstarch-water suspensions have been reported previously [bohn2005hierarchicalI, bohn2005hierarchicalII, akiba2017morphometric]. The secondary cracks in cornstarch-water suspensions have been used as a model system to investigate the formation of geophysical columnar joints [muller1998starch, muller1998experimental, toramaru2004columnar, goehring2008scaling, goehring2009drying, goehring2009nonequilibrium, goehring2013evolving]. However, to the best of our knowledge, the origin of the different types of cracks is not well-understood. It has been suggested that the small-scale cracks are driven by the spatial nonuniformity of the local shrinkage of the film [mizuguchi2005directional, akiba2017morphometric], which has never been confirmed. More recently, Goehring [goehring2009drying] showed that the strong separation between two distinct drying regimes dominated by liquid and vapor transport in the particle network could influence the formation of small-scale crack patterns in cornstarch-water suspensions, yet similar physics should apply in other particle networks where small-scale cracks are not observed. Consequently, the underlying mechanism for the formation of multiscale crack patterns remains unclear.

One of the main results of this work is that distinct polygonal crack patterns are due to distinct shrinkage mechanisms. The initial evaporation of the suspending liquid creates capillary stress at the interface, which shrinks the sample and induces stresses sufficient for cracking (primary cracks). For cornstarch in water, the secondary cracks are driven by a second shrinkage mechanism: the deswelling of the particles. We observed strong deswelling by drying swollen cornstarch particles dispersed in dilute suspensions, and the average deswelling ratio was $\approx$ 5-10%, as shown in Figs. 5a-5b.

In order to examine whether particle deswelling is unique for cornstarch-water suspensions, we prepared suspensions of cornstarch-water, -IPA, and -silicone oil, and suspensions of CaCO3-water and -silicone oil and compared their drying kinetics. We used an electronic balance to record the instantaneous mass of the prepared suspensions during drying, and the results are shown in Fig. 5c. The instantaneous mass of the liquid, $m_{l}(t)$ , is normalized by the initial mass of the liquid, $m_{l0}$ (vertical axis), and the drying time $t$ is normalized by the initial evaporation rate of the liquid, $m_{l0}/(\Delta m/\Delta t)$ (horizontal axis). It can be easily seen in Fig. 5c that the normalized drying dynamics of all of the suspensions follow the same curve, except for cornstarch and water. The drying dynamics are much slower at late times for cornstarch in water, and the sample takes more than twice as long to dry. This discrepancy suggests that the cornstarch particles are deswelling in the later drying stage, so that evaporation depends on diffusion of water out of the individual particles. Finally, the point where the drying rate of the cornstarch-water suspensions start to deviate from the other four suspensions is exactly when small-scale, secondary cracks show up during drying cornstarch-water suspensions in petri dishes ( $t\sim 24$ h, see Fig. 2), indicating that the secondary cracks are driven by deswelling-induced shrinkage.

Generally, the drying kinetics of granular suspensions follow a similar pattern, i.e., initially the saturated state, in which the drying rate is only determined by the evaporation at the top free surface, next drying rate almost keeps constant and then decreases followed by a saturation state as the liquid content is nearly completely gone, and these characteristics are essentially determined by the fluid-pore interactions, pore size distributions and the connectivities of pores [silverstein1997studies, lu2004unsaturated, goehring2015desiccation]. Therefore, the normalized drying rate for various particle-liquid suspensions shown in Fig. 5 furthermore suggests the unique small-scale, secondary cracks are driven by the deswelling of cornstarch particles in water given the differentiated drying behavior of cornstarch particles.

4 Primary cracks in different particle suspensions

Above the critical thickness, we explored primary crack patterns in various suspensions. We dried suspensions of cornstarch and CaCO3 in water, silicone oil, and IPA in petri dishes. We varied the initial volumes of the suspensions in order to obtain different film thicknesses and polygon areas. Figure 6a shows the multiscale cracks observed in dried cornstarch-water films, for different film thicknesses. For primary polygonal cracks, the average polygon area increases with thickness, whereas for secondary cracks, the average polygon area initially increases with film thickness and then saturates for large $h$ . For CaCO3-water suspensions (Fig. 6b), the polygonal pattern is not isotropic, and has a preferred direction. This anisotropy is well-known, and is likely due to particle chain formation induced by drying [akiba2018]. For CaCO3 particles in both silicone oil and IPA (Figs. 6c and 6d), the cracks are distinctly different than those observed in water, nevertheless, the characteristic size of the polygons increase with thickness. This dependence will be discussed in detail in Section 6.

Although the liquids used in our experiments have different surface tensions and vapor pressures, we suspect that some of the differences in patterns are mostly due to particle-liquid interactions and surface energies. In order to examine these assumptions, we have tested the packing ability of various particle-liquid combinations, and found significant differences between the same particles in different solvents (Fig. 7). We prepared cornstarch-water, cornstarch-silicone oil, and cornstarch-IPA suspensions with the same initial volume ( $V_{i}=$ 10 ml), and the same initial volume fraction ( $\phi_{i}=20\%$ ). The samples were then centrifuged at 2000 rpm for 2 mins. Longer centrifuge times did not change the result. According to particle number conservation,

[TABLE]

where $\phi_{f}$ and $V_{f}$ are the final volume fraction and the total volume of the particle suspension after being centrifuged, respectively. Figure 7a shows images of the initial cornstarch suspensions, and the images of the cornstarch suspensions after being centrifuged are shown in Fig. 7b. From Fig. 7b we find that $V_{f}\approx$ 4.6 ml, 4.4 ml, and 3.9 ml for cornstarch particles in silicone oil, water, and IPA, respectively. Therefore, we can calculate that $\phi_{f}\approx 43\%$ for cornstarch-silicone oil, $\phi_{f}\approx 45\%$ for cornstarch-water, and $\phi_{f}\approx 51\%$ for cornstarch-IPA, respectively. The results shown in Fig. 7 suggest that particles suspended in different liquids show different packing abilities, depending on the particle-liquid interactions.

Additionally, we measured the packing abilities of different particles suspended in water. To do this, we prepared the suspensions of cornstarch-water, CaCO3-water, and glass beads-water with the same initial volume ( $V_{i}=10$ ml), and the same initial volume fraction ( $\phi_{i}=20\%$ ) as shown in Fig. 8a from left to right, respectively. Figure 8b shows the images of the three suspensions after being centrifuged. From Fig. 8b, we can obtain the the final volume $V_{f}$ of different particles in water, i.e., $V_{f}\approx 4.4$ ml for cornstarch, $V_{f}\approx 3.8$ ml for CaCO3, and $V_{f}\approx 3.4$ ml for glass beads. Therefore, using Eq. 2 we can calculate the final volume fraction $\phi_{f}$ for different particles in water, i.e., $\phi_{f}\approx 45\%$ for cornstarch, $\phi_{f}\approx 53\%$ for CaCO3, and $\phi_{f}\approx 59\%$ for glass beads. The results shown in Fig. 8 indicate that the packing efficiency of different materials suspended in a given liquid can be significantly different. From Figs. 7 and 8, we can conclude that the packing efficiency of particles in liquids can potentially affect the modulus of the material and the visibility of the cracks, among other properties. It should be noted that the packing fraction of porous media is not well-defined, and here we we calculate the volume fraction by the ratio of the granule particle volume to the total volume of liquid and granule particle in a given suspension for the sake of simplicity.

The ultimate packing fraction obtained after desiccation can potentially affect the maximum strain attainable upon drying, the modulus and tensile strength of the sample, the adhesion to the underlying substrate, and the visibility of the cracks. For example, after centrifuging prepared suspensions, we found that cornstarch particles pack significantly more densely ( $\phi_{f}\approx 0.51$ ) in IPA than either water ( $\phi_{f}\approx 0.45$ ) or silicone oil ( $\phi_{f}\approx 0.43$ ), as shown in Fig. 7. In addition, both cornstarch and CaCO3 pack more loosely than spherical glass beads in water (Fig. 8).

We also observed a stark contrast in the crack patterns for cornstarch in both IPA and silicone oil as shown in Fig. 9. For IPA suspensions, we only observed very fine surface cracks in thick samples ( $h\gtrsim$ 1 cm), which did not penetrate more than $\approx$ 1 mm into the desiccated material (Fig. 9a). We suspect that this is due to high packing density of cornstarch in IPA (Fig. 8), so that only the surface layer could obtain sufficient strain to crack upon desiccation. Suspensions of cornstarch in silicone oil did not display any visible cracks for most thicknesses used in our experiments (Fig. 9b), and only small cracks for very thick samples ( $h\approx$ 2 cm). Even after full desiccation, the surface of these films looked like smooth paste (Fig. 9b), suggesting that the particles retained some sort of sticky interactions, possibly due to residual silicone oil adhered to the surface. Even weak, attractive interactions are expected to have a significant effect on granular packings and their mechanical properties for large system sizes [koeze2018sticky]. This hypothesis is consistent with the low packing density of cornstarch in silicone oil (Fig. 8), suggesting that the particles may have a strong affinity for the silicone oil.

5 Cornstarch-water suspensions in thin chambers

The small-scale secondary cracks observed in Fig. 2i are a unique feature of cornstarch-water suspensions, and did not show evidence of a critical thickness ( $h_{c}$ ), in contrast to primary cracks. In order to explore the thickness dependence of the secondary cracks, we prepared cornstarch-water suspensions with different initial volume fractions, $\phi_{i}$ , and then deposited the suspensions into the thin, quasi-two-dimensional chambers, as shown in Figs. 10a and 10b. We observed two drying stages: initially a compaction front invades throughout the film; then a second drying stage “percolates” throughout the film with a characteristic branching pattern, leading to the formation of liquid, capillary bridges between particles. Finally, the liquid bridges dried up followed by the formation of polygonal cracks after several days (see Figs. 10c-10e).

Figure 11 shows images of the polygonal cracks observed in the dried films of a cornstarch-water suspension ( $\phi_{i}=26\%$ ) in chambers with increasing $h$ . Note here the thickness of the chamber is safely taken as the thickness of the dried film since the final dried film was attached to the top and bottom surfaces of the chambers, as shrinkage mostly occurred in the plane (Fig. 10b). In very thin chambers ( $h\simeq R$ ), “dendritic” fracture patterns are observed, as shown in Figs. 11a and 11b. Since the thickness of the chamber is comparable to the particle size, these patterns were sensitive to the flow of the suspending liquid as evaporation occurred. As $h$ exceeds a critical value, the sensitivity to flow ceased, and regular, polygonal cracks appeared, as shown in Figs. 11c-11j. This trend holds true for all values of $\phi_{i}$ used in the experiments.

It should be noted here that we also prepared suspensions of glass beads and CaCO3 particles in thin chambers with water, however, no cracks appear during drying of the suspensions in these chambers. This supports the hypothesis that the secondary cracks are due to a distinct drying mechanism involving deswelling of the hygroscopic cornstarch particles (Section 3).

6 Universal scaling of multiscale polygonal cracks

Figures 6 and 11 show that the characteristic size of all observed polygonal cracks, in both petri dishes and thin chambers, increases with $h$ . Although there are many ways to characterize the polygonal patterns, such as the average number of edges, or the aspect ratio, we simply measured the average area $A_{p}$ of the polygonal cracks. The results are shown in Fig. 13. The solid symbols represent cornstarch-water suspensions desiccated in thin chambers. The data were nearly independent of the initial volume fraction deposited in the chamber (Fig. 12). The open symbols represent polygons observed in petri dishes. More specifically, the open blue circles represent the secondary cracks in dried cornstarch-water suspensions in petri dishes, which overlap with the data from the thin chambers. For very thick suspensions of cornstarch and water, the area of the small-scale, secondary polygons saturated, and did not increase further. As we will show in Section 7, the deswelling of the particles proceeds as a drying front that penetrates diffusively into the material. Thus, the effective thickness associated with the crack formation depends on the diffusion of water vapor from the film.

As mentioned previously, numerous authors have investigated the thickness dependence of the characteristic crack spacing or polygon area in desiccated suspensions [groisman1994experimental, allain1995regular, komatsu1997pattern, shorlin2000alumina, lee2004drying, ma2012possible]. However, for polygonal crack networks, most of these studies cover a very limited range in thicknesses, so that a comprehensive picture of the thickness dependence of polygonal cracks is lacking. Although there is significant variation in the data from any single set of experiments, taken together, our results in Fig. 13 strongly suggest:

[TABLE]

over a wide range of thickness, for different types of polygonal cracks, in different experimental geometries, and for different liquid-particle combinations. This scaling law is indicated by the dot-dashed blue and dashed black lines in Fig. 13. Although the prefactor, $\alpha$ , is distinct for primary and secondary cracks in different materials, the data suggest that the exponent is universal. Recently, Flores [flores2017mean] derived this simple scaling law using continuum elastic theory, and a balance of surface energy and elastic energy for the initiation of cracking. We will repeat this argument here since we will make small alterations to the expression for $\alpha$ .

Figure 14 shows a section of a thin film after the formation of cracks [flores2017mean]. The thickness of the film is $h$ , and the characteristic area of the polygons is $A_{p}$ . The main tensile stress, $P$ , acts on the bottom surface, where the film is adhered to the substrate. We assume that cracks will form when the energy cost of creating new surfaces at the sides of a polygon is equal to the elastic energy released during cracking:

[TABLE]

where $\gamma$ is the surface tension of the newly-created interfaces which have a typical area $\sqrt{A_{p}}h$ , and $\kappa$ is the ratio of the perimeter to the area of the polygons. Here, $\left<\sigma\right>$ is the volume-averaged tensile stress in the film, $E$ is Young’s modulus, and $V=A_{p}h$ is the volume of one polygon.

We can relate $\left<\sigma\right>$ to $P$ using the fact that the average stress in a volume element can be related to an integral over the forces acting on its boundaries [landau1986theory]:

[TABLE]

Since the main stress is applied at the substrate, Eq. 5 provides the approximate scaling:

[TABLE]

Although we have not rigorously evaluated the integral for a thin, adhered film, one obtains the same result, $\left<\sigma\right>\propto 1/h$ , by considering a similar problem, a thin, spherical elastic shell under uniform pressure. In this case, which can be solved exactly, the tangential, tensile stress scales in the same way as Eq. 6 [lautrup2011].

Combining Eqs. 4 and 6, we arrive at the predicted scaling relation:

[TABLE]

Given the surface tension of the new, “wet” interfaces, $\gamma$ , the prefactor, $\alpha=(2\kappa\gamma E/P^{2})^{2/3}$ , is determined by the modulus of the material when cracks form, and stress at the substrate generated by shrinkage, $P$ . When cracks form, $P$ will essentially be the yield stress, and will be smaller than $E$ [mahaut2008rheology]. Most of the polygons we observe are convex. In this case, we can estimate $\kappa$ by assuming they are regular polygons, where $\kappa$ has an analytic expression:

[TABLE]

Here, $N$ is the number of sides of the polygon. As shown in Fig. 15, for $N$ = 3, $\kappa\approx$ 4.56. As $N\rightarrow\infty$ , $\kappa\rightarrow$ 3.54. Thus, in our further discussion, we will assume $\kappa\approx$ 4 for simplicity in estimating the prefactor $\alpha$ .

In order to provide some quantitative measurement of the modulus, we prepared different mixtures at different drying stages, and used a rheometer to measure the applied normal force upon indenting the material with a stainless steel ball of radius $R_{b}\approx 9.5$ mm using the indentation load-displacement method [pharr1992generality, oliver2004measurement]. More specifically, in this method, a complete cycle of loading and unloading of the indenter into the material is performed to measure the elastic and plastic contributions to the deformation. Generally, the loading phase contains both elastic and plastic effect, whereas the initial unloading is assumed to be approximately pure elastic, such that the Young’s modulus can be calculated using the linearly-elastic Hertzian contact theory [doerner1986method, phadikar2012establishing].

The stiffness (slope) for the initial unloading is [oliver2004measurement, pharr1992generality, broitman2017indentation]:

[TABLE]

where $F_{N}$ is the normal force, $y$ is the indentation depth, $y_{m}$ is the maximum displacement when the indenter is fully loaded, $A_{\textrm{proj}}=\pi y_{m}(2R_{b}-y_{m})$ is the projected circular area of the contact for a sphere indented into a half-space by a distance $y=y_{m}$ , $E$ is Young’s modulus, and $\nu$ is the Poisson’s ratio of the film ( $\approx 0.5$ ).

Figure 16 shows the loading and unloading versus displacement for different particulate films when primary and secondary cracks form. Note for each tested film, the entire film thickness is at least 100 times greater than the indentation depth $y_{m}$ , and the maximum load was held for about 10 s before unloading started. By fitting the slope of the initial unloading data at $y=y_{m}$ using Eq. 9, the Young’s modulus of the films fall in the range from 2.6 to 20 MPa. These values are consistent with similar modulii measured in non-Brownian dense suspensions and soils [moller2007shear, mahaut2008rheology, kezdi1974, obrzud2018]. In addition, we also measured the Young’s modulus of fully-desiccated films, and the results are shown in Fig. 17, in which the modulus of films is calculated to be as follows: $3.3\times 10^{6}$ Pa (cornstarch-IPA film, see Fig. 17a), $1.0\times 10^{6}$ Pa (cornstarch-water film, see Fig. 17b), $3.3\times 10^{6}$ Pa (cornstarch-silicone oil film, see Fig. 17c), $6.9\times 10^{7}$ Pa (CaCO3-water, see Fig. 17d). The errors of the nonlinear fitting parameters using Eq. 9 range from 0.3% to 2.0%. The results of Figs. 16a and 16b are consistent with experimental observations that small-scale secondary cracks have a smaller prefactor, as shown in Fig. 13, although this alone does not explain the difference (0.02 m2/3 versus 0.5 m2/3, Fig. 13).

The value of the substrate stress, $P$ , is more difficult to measure. The maximum adhesion to the substrate clearly affects the crack pattern, as shown in Fig. 3, so that a stronger maximum stress at the substrate decreases the polygon area, as indicated by Eq. 7. Even for strong adhesion, shrinkage in the film will continue until $P$ is approximately the yield stress of the particle network.

For a random, close-packed particle network, the yield stress, $Y$ , can be estimated as [goehring2013plasticity]:

[TABLE]

where $\phi$ is the volume fraction of particles, $M$ is the coordination number, and $R$ is the particle radius (see Table LABEL:Drying_tab:1). $F_{\mathrm{max}}$ is the maximum force between particles. For the initiation of primary cracks, the suspension is still saturated with liquid, so we can assume that $F_{\mathrm{max}}/R^{2}$ can be simply estimated as the capillary pressure, $\gamma/R$ . Assuming typical values of the parameters from Table LABEL:Drying_tab:1, Eq. 11 gives:

[TABLE]

This result is not surprising for a wet sample if the inter-particle adhesion in the bulk liquid is small, i.e. the stress required to form a crack by pulling particles apart is of order the capillary forces holding them together. To confirm this, we have performed rheological measurements on cornstarch-water, CaCO3-water, and glass beads-water suspensions with different volume fractions, and the results are shown in Fig. 18. Although shear thickening is observed for larger volume fractions, the maximum shear stress is smaller than $\gamma/R$ [brown2012thickening], showing that capillary forces are larger than any inter-particle force in the bulk liquid.

When a crack forms, we can assume that the boundary stress will be of the same order as the yield stress, so that $P\sim Y$ . Plugging the values of $E$ , $\gamma$ , and $R$ for the large-scale cracks in cornstarch-water samples, and assuming $\kappa\approx 4$ , Eq. 7 yields $\alpha\approx 0.4$ m2/3, which agrees well with dot-dashed blue line in Fig. 13. Given the variability in the modulus $E$ between different suspensions, we do not currently have a way of collapsing all the data in Fig. 13 for the primary cracks. For wet samples, the modulus will likely depend on the surface tension, particle size, the modulus of the particles, the particle shape, and the inter-particle friction. Nevertheless, given this large parameter space, the good agreement with Eq. 7 suggests that the polygonal crack pattern can be quantitatively understood for a range of different particles and liquids, provided some knowledge of the modulus and yield stress of the suspension.

The small-scale, secondary polygonal cracks of cornstarch-water suspensions can be observed in both open petri dishes (Fig. 2) and thin chambers (Fig. 11). This suggests that the formation of small-scale, secondary cracks does not sensitively depend on the drying geometry, and that capillary interactions are not a dominant force, in contrast to the primary cracks. Thus, local particle adhesion likely determines the yield stress of the material. As shown in Section 3, the cornstarch particles are swollen with water, so we do not currently have a way to estimate this adhesion. In addition, the factor of $\gamma$ in the stress balance (Eq. 4) would be related to the surface energies of the particle-particle adhesion [goehring2013plasticity]. Since $F_{\mathrm{max}}\propto\gamma$ for adhesive forces, then $\alpha\propto(\gamma/F_{\mathrm{max}}^{2})^{2/3}\propto 1/\gamma^{2/3}$ . We can then conclude that this adhesion of swollen particles must be stronger than capillary interactions since the prefactor is smaller for small-scale cracks.

7 Effective film thickness for cracks in thick cornstarch-water suspensions

In Fig. 13, we showed how the area of the small-scale, secondary polygons in cornstarch-water suspensions saturated for large values of $h$ (blue open circles). Here we show that this saturation of $A_{p}$ is set by the diffusion of the water vapor in the later drying stage. The individual particles will remain swollen until their environment is sufficiently dry. In an open particle network where evaporation is occurring from above, water transport is initially limited by viscosity as the liquid is pulled through the porous network according to Darcy’s law. As evaporation proceeds, eventually the diffusion of water vapor through the top of the sample limits the transport. This diffusion-limited transport is likely to set a boundary between wet and dry layers, leading to an effective film thickness, $L$ , as illustrated in Fig. 19.

It has been suggested that the water vapor content can be described by a nonlinear, one-dimensional effective diffusion equation [goehring2009drying, goehring2009nonequilibrium]:

[TABLE]

where $\phi_{w}=\phi_{w}(z,t)$ is the spatio-temporal variation of the local volume fraction of water content, and $D(\phi_{w})$ is local diffusivity of the water vapor, and can be expressed as [pel2002analytic, goehring2009nonequilibrium]:

[TABLE]

where $0\leq z\leq h$ , assuming a no-flux boundary condition at the lower boundary $z=h$ . With Eq. 14, Goehring et al. [goehring2009nonequilibrium] measured an average $D(\phi_{w})\sim 10^{-9}$ m2/s with $\phi_{w}$ ranging from 0.1 – 0.3 g/cm3. A similar value has also been reported by Müller [muller1998starch]. In our experiment the characteristic time scale for the formation of secondary cracks in petri dishes is $T\sim$ 24 h, thus the characteristic length scale of the dry layer $L$ can be estimated as $L\sim\sqrt{TD}\sim$ 1 cm, which shows excellent agreement with the saturation thickness for small-scale cracks shown in Fig. 13. It should be noted that in our experiment, the cornstarch-water suspensions were dried under room temperature without introducing extra heat, so that $D(\phi_{w})$ should be a bit smaller than 10*-9* m2/s.

4 Summary and outlook

In summary, we experimentally investigated polygonal crack patterns in desiccated, particulate suspensions composed of various liquid and particle combinations. The thicknesses of the films ranged from $h$ = 10 $\mu$ m to 4 cm. There are two major results of this work. First, the appearance of multiple, distinct length scales associated with cracks results from distinct shrinkage mechanisms during the drying process. Whereas larger, capillary-induced crack patterns occurred in many of the liquid-particle combinations, such multiscale crack patterns only appeared in dried suspensions of cornstarch and water due to deswelling of the hygroscopic starch particles. As Fig. 1 shows, similar multiscale crack patterns can be observed in meter-scale, planetary terrain. This finding alone may help interpret geomorphological history from surface images, even though knowledge of the relevant material properties may not be known.

Second, the characteristic area of the polygons, for all observed cracks, is consistent with a power law scaling: $A_{p}=\alpha h^{4/3}$ , where the prefactor is determined by a balance of surface energy ( $\gamma$ ), film modulus ( $E$ ), and boundary stress ( $P$ ): $\alpha\sim(2\kappa\gamma E/P^{2})^{2/3}$ [flores2017mean]. The values of these parameters depend on the dominant particle-particle interaction forces at play during the initiation of cracking. By quantifying the modulus and equating $P$ with the yield stress, we are able to quantitatively predict $\alpha$ for primary cracks. We note that although this scaling law is consistent with some previous predictions [komatsu1997pattern], other authors have reported a quadratic relationship between $A_{p}$ and $h$ [groisman1994experimental, shorlin2000alumina, leung2000pattern, leung2010criticality]. However, nearly all experimental studies report less than one order of magnitude in thickness variation. In our analysis, we have assumed that for small strains, our films can be considered homogeneous, elastic materials, and it is possible that effects such as Brownian motion, particle density fluctuations, or sticky particle interactions may explain differences observed in the literature for crack spacing and polygon area. We leave this hypothesis to future studies of other materials with different particle interactions.

Although our experiments are limited to laboratory scales, the scaling law, $A_{p}\propto h^{4/3}$ reproduces reasonable values for polygon areas on larger scales. For example, if we assume that polygonal cracks commonly observed in wet mud with $R\approx$ 60 $\mu$ m are mainly due to capillary pressure during drying, and a typical modulus of 5 MPa [kezdi1974, obrzud2018], then Eqs. 7, 8, 11, and 12 give $\sqrt{A_{p}}\approx$ 1 m for a crack depth of $h\approx$ 20 cm. For polygonal crack patterns on much larger scales, such as those show in Fig. 1b, we note that the polygon area may saturate due to heterogeneity in the material properties with depth. In this case it is likely that the modulus of the material is much larger, or that the stress induced during shrinkage is much smaller, in order to produce very large polygon areas. In addition, in the scaling, the boundary stress applied by the substrate requires more detailed characterizations for substrates with diverse conditions for various particle-liquid combinations, and we leave these interesting points for our future studies.

It is expected that the scaling law we found in the laboratory experiments can enhance our general understanding of how nonequilibrium dynamics give rise to crack patterns that we encounter in our daily lives, and it could also provide a novel route for predicting the history and material properties of particular lands (e.g., Fig. 1b) by simply analyzing the images.

Chapter 3 Sedimentation of non-Brownian particles

1 Introduction

In Chapters 1 and 2, I discussed pattern formation driven by hydrodynamic and mechanical instabilities in fluid and granular systems, respectively. In addition, when solid particles are immersed in liquid, they are prone to gravitational settling in the liquid as long as the density of the particle is greater than that of the liquid. This sedimentation process gives rise to a larger number of intriguing natural structures, and is still a main research focus in fundamental fluid mechanics, mechanical, chemical, and biological engineering. In this chapter, I will report some experimental results on the dynamics of sedimentation of non-Brownian particles and the resulting pattern formation driven by such dynamics.

Sedimentation is a ubiquitous and crucial process both in nature and industry, for instance the accumulation of geological deposits [pettijohn1957sedimentary] (e.g., the structure shown in Fig. 5), the separation of particles with different sizes or densities [selim1983sedimentation], and water treatment [guo2009removal]. Driven by its importance, sedimentation has been extensively investigated and considerable progress has been made, however it is still far from being well understood [batchelor1972sedimentation, caflisch1985variance, koch1989instability, herzhaft1999experimental, brenner1999screening, davis1985sedimentation, ramaswamy2001issues, guazzelli2011fluctuations, piazza2014settled]. The difficulty in understanding sedimentation stems from the long-range, many-body hydrodynamic nature of particle interactions as mentioned in Section 4, which usually leads to problems that are difficult to solve mathematically.

Despite the complexity of sedimentation, some well-known results have been successfully proposed and experimentally verified. For instance, the settling velocity of a cloud of particles will always be faster than that of a single particle, which is due to cooperative effect of hydrodynamic interactions between particles [herzhaft1996experimental, guazzelli2011physical]. One of the vital conflicting problems in sedimentation is particle density and velocity fluctuations during sedimentation, which are likely related to the particle volume fraction, system size, side wall effect and inertia [caflisch1985variance, xue1992diffusion, nicolai1995effect, segre1997long]. In addition, the density nonuniformity also plays a key role in the upward swimming for the microorganisms with a heavy bottom [pedley2010collective, pedley1992hydrodynamic, maddock1994mechanics, pedley1990new, guasto2012fluid, pennington1990consequences, roberts1970motion, lauga2016bacterial, elgeti2015physics]. The basic principle behind this particular locomotion is that the net effect of the competing gravitational and hydrodynamic torques exerted on the microorganisms generates an upward motion while swimming [guasto2012fluid].

The dynamics of sedimentation are also strongly dependent on the geometry of the sedimenting objects [manghi2006hydrodynamic]. For example, a pair of disk-like or semi-spherical particles exhibit periodic orbits [jung2006periodic, chajwa2018kepler], and usually non-spherical particles are subject to complicated settling dynamics due the coupling between the translational and orientational degrees of freedom even for a pair of particles [kim1985sedimentation, goldfriend2015hydrodynamic], and the sedimentation dynamics for many-particle systems are even more complicated as shown in Refs. [butler2002dynamic, herzhaft1996experimental]. In what follows, I will show that symmetric rod-like particles with uniform density align with the fluid flow during sedimentation. Generally, for a rigid object moving in a fluid, the velocity $\bm{U}$ at position $\bm{r}$ can be written as:

[TABLE]

where $\bm{U^{s}}$ is the translation velocity and $\bm{\omega^{s}}$ is angular velocity.

For a particle sedimenting in a fluid, the fluid stress field will generate a hydrodynamic force $\bm{F}$ and a hydrodynamic torque $\bm{T}$ :

[TABLE]

where $\bm{\sigma}$ is the fluid stress tensor, and can be expressed as:

[TABLE]

where $\eta$ is the fluid dynamic viscosity, $\tau_{viscous}$ is the viscous stress tensor, $p$ is the absolute pressure which is defined as $p=(-1/3)\sigma_{ii}$ , $\delta_{ij}$ is the Kronecker delta function such that $\delta_{ij}=1$ if $i=j$ , otherwise $\delta_{ij}=0$ , $\varepsilon_{ij}$ is the strain rate tensor and can be expressed in a matrix form in a coordinate system with the axes of 1, 2 and 3 as:

[TABLE]

Since $\varepsilon_{ij}$ is generally symmetric, there is $\varepsilon_{ij}=\varepsilon_{ji}$ . For any material undergoing displacement, the strain rate tensor is related to velocity gradient:

[TABLE]

where $u$ is the velocity and $x$ is the direction.

In Eqs. 2 and 3, $\bm{n}$ is the outward unit vector normal to the particle surface, and $S_{p}$ denotes particle surface. In particular, for a spherical particle with radius $a$ in a uniform flow stream $\bm{U^{\infty}}$ , there is $F=6\pi\eta aU^{\infty}$ (see Eq. 8), which is Stokes drag force [stokes1851effect], and for the same particle in a rotational flow with angular velocity $\bm{\omega^{\infty}\times r}$ there is $T=8\pi\eta a^{3}\omega^{\infty}$ .

Assuming a no-slip boundary condition, the relationship between the object velocity, hydrodynamic force $\bm{F}$ and torque $\bm{T}$ exerted on the object can be expressed in a matrix form [manghi2006hydrodynamic, lauga2009hydrodynamics]:

[TABLE]

where the tensors $\bm{\mu_{t}}$ and $\bm{\mu_{r}}$ are the translation and rotation mobility tensors, respectively, and $\bm{\mu_{tr}}$ and $\bm{\mu_{rt}}$ are the coupling tensors between the translation and rotation motions. By dimensional analysis we can easily arrive $[\bm{\mu_{t}}]\sim(\eta L)^{-1}$ , $[\bm{\mu_{r}}]\sim(\eta L^{3})^{-1}$ , and $[\bm{\mu_{tr}}]\sim(\eta L^{2})^{-1}$ , where $L$ is the characteristic length scale and $\eta$ is fluid dynamic viscosity. More specifically, for symmetric objects the mobility tensors can be diagonalized, i.e., $\bm{\mu_{tr}=\mu_{rt}=0}$ , thus for symmetric rod-like objects we have [guazzelli2011physical]:

[TABLE]

where the dimensionless quantities $\mu_{\perp}$ and $\mu_{\parallel}$ represent perpendicular and parallel mobility constants, respectively, and the values of $\mu_{\perp}$ and $\mu_{\parallel}$ depend on the geometry of the particle, and for long rod-like object they obey $\mu_{\parallel}\simeq 2\mu_{\perp}$ , meaning that an rod-like object with a horizontal orientation moves about twice as fast as an identical object with a vertical orientation. Therefore, the angle $\beta$ between the symmetric axis of the rod-like object and the object velocity can be expressed as [guyon2001physical]:

[TABLE]

where $\theta$ is the angle between the symmetric axis of the object and the vertical direction as schematically shown in Fig. 2. This means for a rod-like object sedimenting in a viscous fluid with the velocity of $U_{stokes}$ will always remain initial orientation, i.e., $\theta$ is constant, which is consistent with the predictions in Refs. [oberbeck1876uber, jeffery1922motion]

Recently a theoretical study by Goldfriend et al. [goldfriend2017screening] showed that particle density distributions can strongly influence the velocity and density fluctuations, and in some certain regimes particles with nonuniform density distribution are expected to suppress velocity and density fluctuations. Figure 1 sketches that the coupling between the density distribution and the translational and rotational motion of the settling particles can lead to distinct sedimentation dynamics. Due to spherical symmetry, the spherical particles with uniform density isotropically respond to the flow disturbance generated by multiple particles, leading to the mutual cancellation of influx and outflux of the particles, and the velocity of spherical particles align with the fluid flow as indicated by the black arrows in Fig. 1a. In contrast, the prolate particles with the center of mass (red spots in Fig. 1b) different from the geometry center will respond anistropically to the flow disturbance due to the competition between the hydrodynamic and gravitational torques. Consequently, this will create a net outflux, i.e., lateral motion, of particles as illustrated in Fig. 1b.

This theoretical analysis is based on linear stability analysis, in which the authors assumed a velocity fluctuations along the horizontal direction, and this results in an imbalance between hydrodynamic and gravitational torques. For spherical particles as shown in Fig. 1a, this torque imbalance simply contributes to the rotational motion of the particles. However, surprisingly the prolate particles with nonuniform density distributions as shown in Fig. 1b are predicted to undergo a lateral motion driven by the anistropic response to the flow disturbance, which strongly depends on the orientation and geometry of the particles. The detailed derivation for differentiated responses to the flow disturbance of the particles with different density distributions and geometries is shown in Ref. [goldfriend2017screening]. It should noticed that the linear stability analysis in Ref. [goldfriend2017screening] is only restricted to fluctuations in the horizontal direction, while in what follows we show experimentally that the mechanisms shown in Fig. 1 also apply when the particles are moving in three dimensions.

As stated above, it is of vital significance to understand and control sedimentation both for the sake of fundamental curiosity (e.g., many-body physics) and practical applications (e.g., manipulation of particulate films with uniform density and homogeneities, understanding large-scale geologic sedimentation pattern). The theoretical ideas by varying the particle density distribution as shown in Fig. 1 [goldfriend2017screening] is a one way to realize that. For this purpose, we experimentally investigated the settling dynamics of non-Brownian particles with different density distributions at low-Reynolds numbers ( $10^{-3}$ ) in quasi-two-dimensional (quasi-2D) and three-dimensional (3D) space. The particles are composed of equal-sized aluminum and steel balls glued together in various configurations. In the experiment, we observed that particles with uniform and nonuniform density distributions respond differently during sedimentation. A single doublet made of equal-sized aluminum and steel balls aligns with gravity, in contrast to a doublet made of two aluminum balls, which aligns with the fluid flow. For a pair of identical particles made of equal-sized aluminum and steel balls, we observed an effective repulsion between particles during sedimentation in both two-body and many-body systems, however, the particles made of identical balls show variable settling dynamics. In addition, for three or more pairs of particles, more complex dynamics were observed that are likely due to the many-body hydrodynamic nature, yet repulsion between the particles with nonuniform densities still played a key role. It is expected that density-distribution-dependent sedimentation can provide a novel tunable route towards the control of uniformity of particulate films after sedimentation, although more work is needed to fully characterize the particle dynamics and interactions.

2 Experimental setup

In the experiment, we used equal-sized (diameter $d\approx$ 2 mm) aluminum (Al) and steel (St) balls with densities of $8\times 10^{3}$ kg/m3 and $2.8\times 10^{3}$ kg/m3, respectively. In order to prepare particles with various density distributions, we glued three different types of particle doublets, namely Al-St, Al-Al, and St-St. We used silicone oils with viscosity of 1000 cSt and 10000 cSt, respectively as fluid for sedimentation. Transparent acrylic sheets and blocks were used to construct quasi-two-dimensional (quasi-2D) and three dimensional (3D) chambers as shown in Figs. 3 and 4, respectively. For sedimentation in quasi-2D chambers, a USB 3.0 digital video camera (Point Grey) connected to a macro lens was used to image the sedimentation behavior from the front, whereas for the sedimentation in 3D chambers, the camera was held from the top to image the sedimentation pattern of particles. Some recorded images were enhanced using NIH ImageJ software for better visualization. The software Trackpy (an open source Python tracking software) was used to track down the trajectories of particles for the experiment in 2D chambers, and the particle diameter ( $\approx 2$ mm) corresponds to 24 pixels in the images.

Given the size of the particles and the viscosity of silicone oil, the Reynolds number as defined in Chapter 1 is calculated to be $10^{-3}<Re<10^{-2}$ , suggesting that inertial effects can be ignored. The Péclet number is defined as $Pe=\frac{aU}{D}$ , which characterizes the relative importance of advective effects to Brownian motion where $U$ is the characteristic velocity (i.e., $U_{stokes}$ ), $a$ is the particle radius, and $D$ is the diffusion coefficient given by $D=\frac{k_{B}T}{6\pi\eta a}$ , where $k_{B}$ is the Boltzmann’s constant, and $\eta$ as fluid dynamic viscosity. Since the particles used in the experiment are milimeter-scale, we can estimate $Pe\sim 10^{12}$ , suggesting that the particle motion was purely determined by hydrodynamics.

3 Results and discussion

1 Effect of density distribution on sedimentation

We first test the settling dynamics of single particles with uniform and nonuniform density distributions in a 3D chamber. We used single particles of Al-Al and Al-St in 3D chambers for sedimentation, and then tracked the settling dynamics of the particles.

Figures 5a and 5b show the images of a single Al-St and Al-Al particle at a certain moment during sedimentation, respectively. In the sedimentation of a single Al-St particle, the heavier part was initially pointed upwards and lighter part downwards. After releasing the particle, the heavier part will induce a torque in order to lower the center of mass, leading the particle to flip. Afterwards, the Al-St particle will align with the external force, i.e., gravity. This process is confirmed by characterizing the variations of angle ( $\theta$ ) between the vertical direction and the symmetric axis with respect to time as shown in Fig. 5c. In contrast to the settling dynamics of Al-St particle, during sedimentation the single Al-Al particle does not align with gravity but rather aligns with the fluid flow. This means that the angle $\theta$ will remain the same as the initial angle during sedimentation as indicated by Eq. 9. The variation of $\theta$ with time is plotted in Fig. 5d, in which fluctuations in $\theta$ is likely due to noise and anomalous disturbances in the fluid, however, the average value of $\theta$ almost remains constant ( $\approx\pi$ /8), and generally this angle depends on the initial conditions of the sedimentation. Note that the trajectories shown in Figs. 5c and 5d represent the variation of the angle between the geometry center axis of the doublet and the vertical direction as indicated in the inset of Fig. 5c.

2 Sedimentation of multiple particles in quasi two dimensions

We also systematically investigated the sedimentation of multiple particles with different density distributions in quasi-2D chambers. Figures 6a to 6c show the images of two identical Al-St particles during sedimentation at different times, which clearly indicate that the separation between the particles increases as sedimentation proceeds. This trend can be confirmed by tracking the trajectories of the particles during sedimentation, and it holds true for different trials as indicated in Fig. 6d. Note that the trajectories with the same color in Fig. 6c two doublet particles in the same trial. In addition, the increase of separation between particles suggests that there is an effective repulsion between particles, which agrees with the mechanism shown in Fig. 1b.

For comparison, we also investigated the sedimentation of two identical Al-Al particles. Figures 7a to 7c show the images of two identical Al-Al particles during sedimentation at different times, in which it is clearly seen that the orientations of each particle and separation between the two particles are more irregular and complicated than that of Al-St particles shown in Fig. 6. Figure 7c shows the trajectories of the particles for two different trials, we can see that the particles’ trajectories often overlap though they pass the overlapping point at different times. In addition, we also occasionally observed flipping of Al-St particles during sedimentation of two such particles. Although there are only two bodies, the complicated dynamics of the sedimentation of two Al-Al particles are likely due to the non-spherical geometry.

Next we investigated the settling dynamics of many particles with different density distributions in quasi-2D chambers. Figures 8a and 8b show the images of 29 St-St and Al-St at given times during sedimentation in a quasi-2D chamber, respectively. It is apparent that both the spacing and orientation of the St-St particles are more irregular than that of Al-St particles. This property can also be confirmed by Figs. 8c and 8d, which show the sedimentation trajectories of particles in Figs. 8a and 8b, respectively.

In order to quantity the interactions between particles during sedimentation as shown in Figs. 8a and 8b, we investigated the statistics of the particle spacing. Figure 9 shows the probability density function (PDF) of particle spacing for the sedimentation shown in Figs. 8a and 8b. The distributions shown in Fig. 9 are essentially the pair correlation function in one dimension. It is also clear that the distribution decays to zero at the particle spacing of $\approx$ 14 cm, which is simply due to the finite size of the experimental system. In other words, if the system size was infinite, the distribution will eventually saturate to a certain non-zero value, typically normalized to unity.

As shown in Fig. 9, the main difference between Al-Al particles and Al-St particles is that at a given small particle spacing ( $\approx$ 0.4 cm), the probability of finding another St-St particle (blue curve) is higher than that of Al-St particle (nearly zero, red curve). The distinct behaviors essentially suggest that there is an effective repulsion between Al-St particles during sedimentation, which is absent in St-St particles. Therefore, the results shown in Fig. 9 validate the mechanisms shown in Fig. 1 despite the many-body nature of the system. It should also be noted that the statistics shown in Fig. 9 are based on the average of five trials for both particle configurations.

3 Sedimentation of multiple particles in three dimensions

Given the distinct sedimentation statistics for Al-Al and Al-St particles in quasi-2D chambers as shown in Fig. 9, we extended our experiment to three dimensions as illustrated in Fig. 4, which is more similar to real-world situations. In the experiment, we prepared the same number of pairs of Al-Al and Al-St particles, and deposited the particles into two separate 3D chambers (see Section 2 for experimental details), then we investigated particle spatial distributions on the bottom of the chamber after sedimentation..

For the sedimentation of multiple particles in 3D, clusters of particle obstruct the view and make tracking all particle impossible. Therefore, we took images of the sedimentation pattern on the bottom of the chamber. We performed more than 100 trials for the sedimentation of Al-Al and At-St particles, respectively, and some typical images are shown in Figures 10 and 11. By simple comparison, we can conclude that the average density of Al-St particles is smaller than that of Al-Al particles.

In order to quantify the average density of the particle layer after sedimentation, we plot the probability density function (PDF) of the particle position in the radial direction, as shown in Fig. 12a, in which integral of the distribution is normalized to unity. It is apparent in Fig. 12a that the distribution curve for Al-St particles is more stretched than that of Al-Al particles, meaning that the cluster size of Al-St particles is larger than that of Al-Al particle clusters, which can also be confirmed by Figs. 12b and 12c. It should be noted that each distribution curve shown in Fig. 12a is based on 10 separate images, and the shaded area for each curve represents the standard deviation for 10 images, and the distinction between the distribution curves of Al-St and Al-Al particle clusters is pretty obvious. The larger cluster size of Al-St particles is due to the hydrodynamic repulsion interaction between the particles, which agrees with the theoretical prediction shown in Fig. 1, though the theory is formulated in 2D [goldfriend2017screening]. The tendency of Al-Al particles to form smaller clusters suggests that the particles are unstable to density fluctuations, and consequently the cluster of Al-St particles is more uniform than that of Al-Al particles (see Figs. 12b and 12c), which agrees with experimental, theoretical and numerical studies on the sedimentation of non-spherical particles in dilute suspensions [koch1989instability, herzhaft1996experimental, helzel2017kinetic, guazzelli2011fluctuations, wang2009numerical, gustavsson2009gravity, tornberg2006numerical].

4 Summary and outlook

In summary, we experimentally investigated sedimentation of particles with different density distribution in quasi-2D and 3D. For the sedimentation of a single particle with uniform density distribution, the particle aligns with the fluid flow, whereas a single particle with nonuniform density distribution aligns with gravity. In quasi-2D chambers, we observed an effective repulsion during sedimentation of a pair of doublets with nonuniform density distribution. In contrast, the sedimentation dynamics of a pair of doublets with uniform density distribution are more complex, we did not observe an repulsion between such particles and we occasional observed flipping during sedimentation. In addition, we found that the repulsion still plays a key role in the sedimentation of many particles with nonuniform density distributions, which is confirmed by the probability density function of particle spacing shown in Figure 9.

For the experiment in 3D, we quantified the uniformity of the final sedimented layer of particles, and found that the particles with nonuniform mass density form layers of lower number density when compared to particles with uniform mass density as indicated in Figs. 10, 11 and 12. This suggests that the repulsion between particles with nonuniform density distribution still plays a role in 3D scenario. Our experimental results on the sedimentation of particles with different density distributions in 2D system validate the main predictions by Goldfriend et al. [goldfriend2017screening]. Moreover, the theoretical analysis in Ref. [goldfriend2017screening] is restricted in 2D, whereas our results in 3D chambers also show that those mechanisms also apply to 3D. We expect our experimental results on the density-distribution-dependent sedimentation dynamics could provide a novel tunable route towards self-assembly and controlling the uniformity of particulate films after sedimentation.

The future direction in this project is to vary the geometry and density distribution of the particles. Figure 13 shows the phase diagram for the effective response parameter $\xi$ with respect to the aspect ratio of the particle $\kappa$ based on linear stability analysis [goldfriend2017screening]. Theoretically, when $\xi>0$ , the density fluctuation of particles will be suppressed, which is driven by the effective repulsion between particles with nonuniform density distributions or the self-alignability of particles. Whereas for the scenario of $\xi<0$ , there is no suppression of density fluctuation and the uniform suspension becomes unstable to density fluctuations, which in turn means that no effective repulsion should be observed between particles with $\kappa<1$ during sedimentation. For spherical particles, $\kappa=1$ , and the theory becomes invalid to predict the stability of the suspension. In addition, it is also interesting to investigate the effect of system size on the dynamics of sedimentation, which is still controversial hitherto [caflisch1985variance, nicolai1995effect, ladd1996hydrodynamic, segre1997long, ladd1997sedimentation], and we leave these open questions for future research.

Chapter 4 Summary

During my Ph.D. research activities, I experimentally investigated the nonequilibrium dynamics of pattern formation by studying a few laboratory-accessible examples in fluid and granular systems. I primarily focused on the formation of diverse spatial patterns in different experiments, and the mechanisms we uncover represent distinct pattern formation phenomena with different physics. This is perhaps not surprising since they are all nonequilibrium in nature. Although the nonequilibrium dynamics that give rise to patterns in real world can be rather complicated given the complex environment they are located in, we expect our findings to enhance our understanding of how nonequilibrium dynamics determine tremendous groups of patterns in nature.

In this thesis, I particularly reported three stories about spatial pattern formation under nonequilibrium conditions, and the main conclusions in each story are summarized as follows:

$\bullet$ Oscillations of Leidenfrost drops

In fluid dynamics, most interesting and intriguing patterns are driven by multiphase-interaction-induced instabilities, such as the examples shown in Fig. 2. Leidenfrost drops are known to form star-shaped oscillations, however, the complicated interplay between the solid, liquid and vapor phases as well as the highly nonequilibrium mass and heat transfer in Leidenfrost phenomenon has made the underlying mechanism for the star-shaped oscillatory pattern remain unclear. In the experiment, by focusing on the star-shaped oscillation frequency of Leidenfrost drops and pressure variations in the sustaining vapor layer, we found that the star-shaped oscillations are parametrically driven by the pressure in the vapor layer. We furthermore assumed that the pressure variations in the vapor are induced by the capillary waves underneath the drop, which was confirmed by imaging the bottom surface of the drop. The capillary waves with characteristic wavelength are generated by a large shear stress at the liquid-vapor interface due to the rapid flow of evaporated vapor, and travel form the drop center to the edge. This leads to the local variations of vapor film thickness, which generates the pressure variations in the vapor. Therefore, we conclude that hydrodynamic coupling between drop interface and vapor flow leads to the star-shaped oscillations of Leidenforst drops. We also explored potential effects of thermal convection in the liquid, and found it plays a minor role in Leidenfrost drop oscillations. In addition to the study of the star-shaped oscillations found in large Leidenfrost drops, we also gave an analytical explanation for the small-amplitude axisymmetric oscillations (breathing mode) observed in small Leidenforst drops by a simple balance of gravitational and surface tension forces. We expect the nonequilibrium dynamics we uncover behind the patterned oscillations in both large and small Leidenfrost drops are of help to understand the coupling dynamics in systems with a supporting gas layer contacting a liquid interface and a solid surface.

$\bullet$ Polygonal desiccation crack patterns

Polygonal desiccation crack patterns are commonly observed in nature, and they span multiple orders of magnitude in length scale. The formation and selection of crack patterns strongly depend on nonequilibrium multiphase interactions between fluid and solid particles as well as the drying kinetics (see Fig. 4). To date, there is still an outstanding question that whether crack patterns on different length scales share similar dynamics. To answer this question, in the experiment we dried particulate suspensions composed of various liquid and particle combinations, and found that large-scale polygonal cracks occurred in many liquid-particle combinations are driven by the shrinkage of the particulate films (capillary-induced), whereas in cornstarch-water suspensions we found that deswelling of the hygroscopic starch particles leads to hierarchical polygonal crack patterns. In addition, we found a universal power law relating polygonal crack area, film thickness, film modulus and boundary stress for multiscale polygonal crack patterns, which can help interpret geomorphological history from surface images without knowing the material property and enable to manipulate crack patterns for specific purposes. Therefore, the nonequilibrium dynamics for the formation of multiscale polygonal crack patterns we uncover in the laboratory experiment can serve as a framework for understanding polygonal crack patterns from microscopic to geologic scales.

$\bullet$ Sedimentation of non-Brownian particles

Sedimentation is a crucial process both in nature and industry, which is however not well understood due to complex coupling between the particle orientation and many-body, long-range hydrodynamic interactions. In the laboratory experiment, I investigated the sedimentation dynamics of non-Brownian particles with different density distributions. By focusing on the statistical properties of the spacial patterns of particles during sedimentation (2D) and after sedimentation (3D), we observed an effective repulsion between particles with nonuniform density distribution, in contrast to particles with uniform density distribution. These results indeed highlight the important role of particle density distribution in determining the highly nonequilibrium dynamic process of sedimentation, and can shed light into a novel pathway for large-scale self-assembly of uniform particulate films by sedimentation, and furthermore can enable to shape geomorphology manually as desired.

The underlying mechanisms we uncovered for the spatial pattern formation in the above three laboratory experimental systems represent the typical examples about how nonequilibrium dynamics give rise to pattern formation in fluid and granular systems. Although realistically the situations for pattern formation are supposed to be more complicated, our results in the limited experiment can help to model studies on the nonequilibrium dynamics that more resemble the real situations of pattern formation in nature.

Bibliography285

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] K. G. Libbrecht, The physics of snow crystals, Rep. Prog. Phys. , 68 , 855 (2005).
2[2] F. Charru, B. Andreotti, and P. Claudin, Sand ripples and dunes, Annu. Rev. Fluid Mech. , 45 , 469 (2013).
3[3] M. Cross and H. Greenside, Pattern Formation and Dynamics in Nonequilibrium Systems (Cambridge University Press) (2009).
4[4] L. Goehring and S. W. Morris, Cracking mud, freezing dirt, and breaking rocks, Phys. Today , 67 , 39 (2014).
5[5] M. Hofmann, R. Anderssohn, H.-A. Bahr, H.-J. Weiß, and J. Nellesen, Why hexagonal basalt columns?, Phys. Rev. Lett. , 115 , 154301 (2015).
6[6] K. Libbrecht, Snowflakes (Voyageur Press) (2008).
7[7] L. Goehring and S. W. Morris, Scaling of columnar joints in basalt, J. Geophys. Res.: Solid Earth , 113 , B 10203 (2008).
8[8] A. J. Koch and H. Meinhardt, Biological pattern formation: from basic mechanisms to complex structures, Rev. Mod. Phys. , 66 , 1481 (1994).

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Taxonomy

Experimental investigations on the nonequilibrium dynamics of pattern formation in fluid and granular systems

Abstract

Contents

List of Figures

List of Tables

Acknowledgements.

Chapter 0 Introduction

1 Pattern formation in nature

2 Patterns in fluids

3 Patterns in colloidal films

4 Patterns driven by sedimentation

Chapter 1 Oscillations of Leidenfrost drops

1 Introduction

2 Experimental setup

3 Results and discussion

1 The geometry of Leidenfrost drops on curved surfaces

2 Breathing mode of small Leidenfrost drops

3 Star-shaped oscillations of Leidenfrost water drops

4 Pressure oscillations in the vapor layer

5 Star-shaped oscillations of different liquids

6 Origin of pressure oscillations in the vapor layer

7 Thermal effects

4 Summary and outlook

Chapter 2 Polygonal desiccation crack patterns

1 Introduction

2 Experimental setup

3 Results and discussion

1 Multiscale cracks in cornstarch-water suspensions

2 Critical condition for primary cracks

3 Cornstarch particle deswelling drives secondary cracks

4 Primary cracks in different particle suspensions

5 Cornstarch-water suspensions in thin chambers

6 Universal scaling of multiscale polygonal cracks

7 Effective film thickness for cracks in thick cornstarch-water suspensions

4 Summary and outlook

Chapter 3 Sedimentation of non-Brownian particles

1 Introduction

2 Experimental setup

3 Results and discussion

1 Effect of density distribution on sedimentation

2 Sedimentation of multiple particles in quasi two dimensions

3 Sedimentation of multiple particles in three dimensions

4 Summary and outlook

Chapter 4 Summary