Influence of surface tension in the surfactant-driven fracture of closely-packed particulate monolayers

TL;DR

This paper presents a phase-field model to study how surface tension influences fracture patterns in closely-packed particulate monolayers, validated by experiments showing high sensitivity to surface tension differences.

Contribution

The study introduces a new phase-field model for surfactant-driven fracture in particulate monolayers, incorporating surface tension effects and validated by experimental data.

Findings

Model qualitatively matches experimental fracture patterns

Surface tension differences critically affect fracture regimes

Thresholds identified for different fracture behaviors

Abstract

A phase-field model is used to capture the surfactant-driven formation of fracture patterns in particulate monolayers. The model is intended for the regime of closely-packed systems in which the mechanical response of the monolayer can be approximated as a linearly elastic solid. The model approximates the loss in tensile strength of the monolayer as the surfactant concentration increases through the evolution of a damage field. Initial-boundary value problems are constructed and spatially discretized with finite element approximations to the displacement and surfactant damage fields. A comparison between model-based simulations and existing experimental observations indicates a qualitative match in both the fracture patterns and temporal scaling of the fracture process. The importance of surface tension differences is quantified by means of a dimensionless parameter, revealing thresholds that separate different regimes of fracture. These findings are supported by newly performed experiments that validate the model and demonstrate the strong sensitivity of the fracture pattern to differences in surface tension.

Figures

Tables

Table 1. Table 1: List of parameters in the model.

parameter	description
${\varphi }_j$	jamming packing fraction of monolayer
${\varphi }_m$	maximum packing fraction of monolayer
${\overline{r}}_p$	average particle radius
${\gamma }_f$	surface tension of liquid layer
${\gamma }_s$	surface tension of surfactant
$G_0$	fracture toughness of undamaged monolayer
$F_0$	surfactant force constant
$E$	Young’s modulus of monolayer
$\nu$	Poisson’s ratio of monolayer
${\mu }_f$	dynamic viscosity of monolayer
${\beta }_0$	kinetic parameter

Table 2. Table 2: Dimensional variables. 10 , 8

Property	Description	Value	Unit
$r_d$	radius of Petri dish	$8.5\cdot {10}^{-2}$	m
$\left[{\varphi }^0\left(x\right)\right]$	range of initial packing fraction	$(7.5$–$8.4)\cdot 10{}^{-1}$	none
${\varphi }_m$	packing fraction of jammed monolayer	$9.0\cdot {10}^{-1}$	none
$G_0$	fracture toughness of undamaged monolayer	$5.7\cdot ({10}^{-4}$–${10}^{-2})$	N/m
$F_0$	surfactant force constant	$1.0\cdot ({10}^{-1}$–${10}^2)$	1/m
${\gamma }_f$	surface tension of liquid layer (water)	$7.2\cdot {10}^{-2}$	N/m
${\gamma }_s$	surface tension of surfactant (oleic acid)	$4.0\cdot {10}^{-2}$	N/m
${\mu }_f$	range of dynamic viscosity of liquid layer	$1.0\cdot ({10}^{-3}$–${10}^{-1})$	Pa s
${\overline{r}}_p$	Mean particle radius	$1.0\cdot {10}^{-4}$	m
$E$	Young’s modulus of monolayer	$3.178\cdot {10}^3$	Pa
$\nu$	Poisson’s ratio of monolayer	$1/\sqrt{3}$	none
${\beta }_0$	kinetic parameter	$1.0$	1/m
$\lambda$	characteristic length scale	$1.0\cdot {10}^{-2}$	1/m

Table 3. Table 3: Dimensionless parameters.

Parameter	Value
$G^{\ast }$	${10}^{-3}$–${10}^{-1}$
$F^{\ast }$	${10}^{-4}$–${10}^{-1}$
$L^{\ast }$	$8.5\cdot {10}^2$
${\beta }_0^{\ast }$	$2.7\cdot ({10}^{-8}$–${10}^{-6})$
${\gamma }^{\ast }$	$5.6\cdot {10}^{-1}$
${\lambda }^{\ast }$	$1.5\cdot {10}^{-1}$

Equations

$$\mathbf{\epsilon }=\frac{1}{2}(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^{\top }).$$

$$\varphi (\mathbf{\epsilon })=\frac{{\varphi }_0}{1+\text{mbox}tr\mathbf{\epsilon }}.$$

$$W_{\pm }(\mathbf{\epsilon })=\frac{E}{2(1+\nu )}|{\mathbf{\epsilon }}_{\pm }{|}^2+\frac{\nu E}{(1-\nu {)}^2}(\text{mbox}tr{}_{\pm }\mathbf{\epsilon }{)}^2,$$

$$\text{mbox}tr{}_{+}\mathbf{\epsilon }=\{\begin{array}{ll}\text{mbox}tr\mathbf{\epsilon },& \text{mbox}if\text{mbox}tr\mathbf{\epsilon }\ge 0,\\ \text{vskip}4.0pt& \\ 0,& \text{mbox}otherwise,\end{array}$$

$$\text{mbox}tr{}_{-}\mathbf{\epsilon }=\{\begin{array}{ll}\text{mbox}tr\mathbf{\epsilon },& \text{mbox}if\text{mbox}tr\mathbf{\epsilon }\le 0,\\ \text{vskip}4.0pt& \\ 0,& \text{mbox}otherwise.\end{array}$$

$$\psi (\mathbf{\epsilon },d,\nabla d)=\underset{\text{tensile}}{\underbrace{(1-d{)}^2{W}_{+}(\mathbf{\epsilon })}}+\underset{\text{compressive}}{\underbrace{\xi (\varphi (\mathbf{\epsilon }))W_{-}(\mathbf{\epsilon })}}+\underset{\text{fracture}}{\underbrace{Q(d,\nabla d)}}+\underset{\text{surfactant}}{\underbrace{F(\mathbf{\epsilon },d)}}.$$

$$Q(d,\nabla d)=\frac{G(d)}{2}\left(\frac{d^2}{\lambda }+\lambda |\nabla d{|}^2\right)$$

$$G(d)=G_{0}\Big{(}\frac{\gamma_{s}}{\gamma_{f}}+\Big{(}1-\frac{\gamma_{s}}{\gamma_{f}}\Big{)}(1-d)\Big{)}=G_{0}\Big{(}1-\Big{(}1-\frac{\gamma_{s}}{\gamma_{f}}\Big{)}d\Big{)},$$

$$G_r=G(1)=G_0\frac{{\gamma }_s}{{\gamma }_f}$$

$$F(\mathbf{\epsilon },d)=\frac{F_0({\gamma }_f-{\gamma }_s)(\text{mbox}tr\mathbf{\epsilon })d^2}{2},$$

$$\text{div}\Big{(}\frac{\partial\psi({\bm{\varepsilon}},d,\nabla d)}{\partial{\bm{\varepsilon}}}\Big{)}={\bf 0}$$

$$\text{div}\Big{(}\frac{\partial\psi({\bm{\varepsilon}},d,\nabla d)}{\partial\nabla d}\Big{)}-\frac{\partial\psi({\bm{\varepsilon}},d,\nabla d)}{\partial d}=\begin{cases}\mskip 8.25mu\beta\dot{d},&\text{if}\quad\dot{d}>0,\\ -\pi_{r},&\text{if}\quad\dot{d}=0,\end{cases}$$

$$\beta ={\beta }_0\frac{1-{\varphi }_m}{1-\varphi },{\beta }_0=\frac{C{\mu }_f{\bar{r}}_p}{1-{\varphi }_m},$$

$$\frac{1-\phi(\bm{\varepsilon})}{\beta_{0}(1-\phi_{m})}\Big{\langle}\text{div}\Big{(}\frac{\partial\psi({\bm{\varepsilon}},d,\nabla d)}{\partial\nabla d}\Big{)}-\frac{\partial\psi(\bm{\varepsilon},d,\nabla d)}{\partial d}\Big{\rangle}=\dot{d},$$

$$⟨h⟩=\{\begin{array}{ll}0,& h\le 0,\\ h,& h>0,\end{array}$$

$$\text{div}\Big{(}(1-d)^{2}\frac{\partial W_{+}(\bm{\varepsilon})}{\partial\bm{\varepsilon}}+\frac{\partial(\xi(\phi(\bm{\varepsilon}))W_{-}(\bm{\varepsilon}))}{\partial{\bm{\varepsilon}}}+\frac{F_{0}(\gamma_{f}-\gamma_{s})d^{2}}{2}\bm{I}\Big{)}={\bf 0}$$

$$\frac{1-\phi(\bm{\varepsilon})}{\beta_{0}(1-\phi_{m})}\Big{\langle}G_{0}\lambda\Big{(}1-\Big{(}1-\frac{\gamma_{s}}{\gamma_{f}}\Big{)}d\Big{)}\Delta d-\frac{G_{0}\lambda}{2}\Big{(}1-\frac{\gamma_{s}}{\gamma_{f}}\Big{)}|\nabla d|^{2}\\ +2(1-d)W_{+}(\bm{\varepsilon})-F_{0}(\gamma_{f}-\gamma_{s})(\mbox{tr}\,\bm{\varepsilon})d\\ -\frac{G_{0}d}{\lambda}\Big{(}1-\frac{3}{2}\Big{(}1-\frac{\gamma_{s}}{\gamma_{f}}\Big{)}d\Big{)}\Big{\rangle}=\dot{d},$$

$$\mathbf{u}=0$$

$$\nabla d\cdot \mathbf{n}=0.$$

$$\varphi (x,0)={\varphi }_0(x),d(x,0)=d_0(x),$$

$${\mathbf{x}}^{\ast }=\frac{\mathbf{x}}{L}\text{and}{t}^{\ast }=\frac{t}{T}.$$

$${\nabla }^{\ast }d=L\nabla d,{\Delta }^{\ast }d=L^2\Delta d,\text{and}\frac{\partial d}{\partial t^{\ast }}=T\dot{d}$$

$$\left.\begin{array}[]{c}\displaystyle G^{*}=\frac{G_{r}}{E\bar{r}_{p}},\qquad F^{*}=\frac{F_{0}(\gamma_{f}-\gamma_{s})}{E},\qquad L^{*}=\frac{L}{\bar{r}_{p}},\cr\vskip 4.0pt\cr\displaystyle\beta_{0}^{*}=\frac{\beta_{0}}{ET},\qquad\gamma^{*}=\frac{\gamma_{s}}{\gamma_{f}},\qquad\text{and}\qquad\lambda^{*}=\frac{\lambda}{L},\end{array}\right\}$$

$$W_{\pm }^{\ast }(\mathbf{\epsilon })=\frac{1}{2(1+\nu )}|{\mathbf{\epsilon }}_{\pm }{|}^2+\frac{\nu }{(1-\nu {)}^2}(\text{mbox}tr{}_{\pm }\mathbf{\epsilon }{)}^2,$$

$$\text{div${}^{*}$}\Big{(}(1-d)^{2}\frac{\partial W^{*}_{+}(\bm{\varepsilon})}{\partial\bm{\varepsilon}}+\frac{\partial(\xi(\phi(\bm{\varepsilon}))W^{*}_{-}(\bm{\varepsilon}))}{\partial{\bm{\varepsilon}}}+\frac{F^{*}}{2}d^{2}\bm{I}\Big{)}={\bf 0}$$

$$\frac{1-\phi(\bm{\bm{\varepsilon}})}{\beta_{0}^{*}L^{*}(1-\phi_{m})}\Big{\langle}\frac{G^{*}}{\gamma*}\lambda^{*}(1-(1-\gamma^{*})d)\Delta^{*}d-\frac{G^{*}}{\gamma*}\frac{\lambda^{*}}{2}(1-\gamma^{*})|\nabla^{*}d|^{2}\\ -\frac{G^{*}}{\gamma*}\frac{d}{\lambda^{*}}(1-\frac{3}{2}(1-\gamma^{*})d)\\ +2(1-d)L^{*}W^{*}_{+}(\bm{\varepsilon})-L^{*}F^{*}(\mbox{tr}\,\bm{\varepsilon})d\Big{\rangle}=\frac{\partial d}{\partial t^{*}}$$

$$\chi =\frac{(F^{\ast }{)}^2}{G^{\ast }}.$$

$$\chi \propto \frac{(F_0{\bar{r}}_p{)}^2}{G_0}\frac{({\gamma }_f-{\gamma }_s{)}^2}{{\gamma }_s}.$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Department of Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA*††footnotetext: b* Mathematical Soft Matter Unit, OIST Graduate University, Onna-son, Okinawa, 904-0495, Japan ††footnotetext: c Collective Interactions Unit, OIST Graduate University, Onna-son, Okinawa, 904-0495, Japan *

1 Introduction

When a densely packed monolayer of hydrophobic particles is placed on the surface of a liquid, the particles interact through capillary bridges that form,1 leading to the formation of particle rafts. The macroscopic properties of these rafts reflect an interplay between fluid and solid mechanics,2, 3, 4 giving rise to novel physics. This interplay is relevant to a wide range of applications, from the synthesis of “liquid marbles” 5 to the design of drug delivery systems 6 to the stabilization of drops.7

The interest in particle rafts has driven researchers to investigate their mechanical properties.4, 8 It is now known that densely packed monolayers exhibit a two-dimensional linearly elastic solid response, and that the mechanical properties of such monolayers depend proportionally on the surface tension of the liquid layer. The introduction of a controlled amount of surfactant generates a surface tension gradient, producing Marangoni forces 9 and causing the surfactant to spread, fracturing the monolayer. Studies of surfactant-driven fracture have examined the role of viscosity and the initial packing fraction on the evolution of cracks in closely and loosely-packed systems, respectively.10, 11 Surprisingly, the potentially important role of differences in the surface tension of the surfactant and the underlying liquid has not been explored. The magnitude of this difference is interesting because it determines the Marangoni force exerted on the particulate monolayer and it is the main driving force for the fracture process. Modulating the surface tension difference also provides a way to probe mechanical properties that are difficult to measure directly, such as the critical failure stress. Finally, the surface tension difference can easily be controlled in the laboratory by modifying the composition of the surfactant or the underlying liquid.

In this article, we focus on closely-packed systems. We develop a phase-field model that takes into account a two-way coupling between the flow of surfactant and the motion of the monolayer. Through model-based simulations and accompanying experiments, we demonstrate that surface tension differences play a vital role in the overall fracture response of particle rafts.

The general setup for a surfactant injection experiment, illustrated in Figure 1, consists of a circular Petri dish containing a liquid layer onto which hydrophobic particles are deposited. The surfactant is introduced with a needle near the center of the dish. The surface tension of the liquid-vapor interface decreases where the surfactant is present. Marangoni forces then cause the surfactant to spread over the surface of the liquid and through the monolayer. The subsequent response is sensitive to the ratio of the fraction of the area of the liquid-vapor interface that is occupied by particles, which we refer to as the packing fraction and denote by $\phi$. As $\phi$ is increased, the properties of the liquid-vapor interface change from liquid to solid.11

Consistent with out interest in closely-packed systems, we consider situations in which the packing fraction is high ($0.7\leq\phi\leq 0.9$). We model the particle laden liquid-vapor interface as a continuous two-dimensional linearly elastic solid, capable of supporting both tension and compression.4, 8 The mechanical properties of such a particulate monolayer 12 are basic to understanding its response to surfactant-driven stresses. An estimate of those properties is provided by Vella et al.,8 based on experimental measurements and geometrical arguments. Surface wave experiments have been used to characterize the stretching and bending stiffnesses of particulate monolayers,13 which are found to present a soft granular character with nonclassical response under fluid-driven deformation.14

Experiments show that when such a system is stimulated by the localized introduction of a surfactant, an advancing front forms and fractures the monolayer.10 This type of surfactant-induced effect has also been observed in other systems, such as agar gels.15 The resulting fracture patterns, illustrated in Figure 1, are reminiscent of those observed in classical brittle materials but with significant differences. Bandi et al. 11 suggested that the fracture patterns can be very sensitive to variations in the initial distribution of particles. Additionally, in conventional elastic solids, crack branching is mostly associated with inertial effects and dynamic instabilities (e.g., bifurcations which occur as the crack tip velocity approaches 60% of the Rayleigh wave speed 16). In contrast, crack branching has been observed in particulate monolayers at crack speeds as low as 0.2% of the shear wave speed.10 Regarding the time scales of the fracture process, the most obvious distinction is that crack tip velocities do not appear to be influenced by the elastic properties of the monolayers. Rather, experimental observations suggest that the fracture time scale is mostly governed by variations in the viscosity of the underlying liquid 10 and the packing fraction.11 These observations raise a number of questions regarding the fluid-driven fracture of elastic media in general and in monolayers in particular.

The mechanisms underlying the surfactant-driven fracture of particulate monolayers are challenging to study experimentally. There are practical limits to the range of particle sizes that can be used. Sufficiently small particles fall below the current resolution limit of imaging equipment when a full field of view is needed. Particles that are too large have too much inertia to move significantly in response to surfactant flow. Modeling and simulation efforts can address these concerns and provide detailed insight concerning the basic mechanisms and sensitivities. However, computational studies of these systems are in the early stage of development. Previous attempts to model particulate monolayers have focused on loosely-packed systems, which admit several simplifying assumptions. For example, Bandi et al.,11 developed a discrete-element method to examine the influence of the initial packing fraction on the number of fractures in loosely-packed systems. In that work, a one-way coupling in which surfactant flow influences the motion of the particles, but not vice-versa, as assumed. While simulations based on this approach accurately reproduce the limiting number of cracks that develop, the underlying assumptions limit its applicability to situations where the packing fraction is low.

The numerical simulation of closely-packed systems is challenging due to the high number of particles and the complexity of the physics involved in fracture. In this work, we propose a model based on a phase-field method which makes it possible to smoothly represent transitions between the damaged and undamaged zones of the monolayer as the surfactant advances. Phase-field models are suitable for this kind of problem since they avoid the need to track the propagating front, while allowing for a simple and powerful way to introduce the essential physics. We take as a starting point the work of Miehe et al. 17 and Borden et al. 18 for phase-field regularizations of the Griffith model for fracture. Based on these approaches, we present a new model that includes several features that are important to characterize the fracture of particulate monolayers. In particular, our model incorporates the distribution of particles, the force on the monolayer due to the presence of surfactant, and the viscosity of the underlying liquid. Our objective is to develop the simplest model capable of capturing the salient features of this system, namely the sensitivity of fracture patterns to differences in surface tension and the temporal scaling of the fracture process.

This paper is organized as follows. In Section 2, we present a phase-field approach for modeling the fracture mechanics of particulate monolayers, placing emphasis on the fundamental nature of the terms used to describe each contribution to the physics, after which we propose a governing free energy per unit area from which the model is derived. In Section 3, we describe the materials and experimental methodology used to explore the fidelity of the computational model. In Section 4, we present numerical results which demonstrate the capability of the model to reproduce different cases of surfactant-driven fracture. We then present a phase diagram which delineates different fracture regimes as a function of the surface tension difference and fracture resistance of the monolayer. Finally, in Section 5, we summarize our main results and propose directions for further research.

2 Phase-field model

We model the surfactant-driven fracture of particulate monolayers by adapting recently developed phase-field approaches to fracture.17, 18, 19 The model is designed to capture how a drop of surfactant introduced at the center of the domain can spread to form cracks or fissures in the particulate monolayer (Figure 1). Accordingly, our model incorporates a number of postulates that are based on experimental observations, as detailed below. Central among these is the assumption that the fractured zones are completely filled with surfactant, without appreciable penetration of the surfactant beyond the contours of the cracks. On this basis, we use a single “surfactant damage” field as an indicator function for the surfactant concentration and for the damage to the monolayer.

2.1 General considerations

Consider a monolayer of particles floating on the surface of a liquid layer. Let the two-dimensional region occupied by that monolayer be denoted by $\Omega$. In Figure 2(a), the dark portion of $\Omega$ represents the intact portion of the monolayer and the light portion of $\Omega$ represents the portion of the monolayer damaged by the surfactant. The state of the monolayer is described by two independent variables, its vector displacement field $\bm{u}$ and a scalar surfactant damage field $d$. As indicated in Figure 2(a), $d$ takes values in $[0,1]$, with $d=1$ representing completely damaged material. For closely-packed monolayers, experimental observations suggest that fracture occurs even for small values of the (infinitesimal) strain tensor

[TABLE]

The microstructure of the monolayer is described by a scalar packing fraction $\phi$ that takes values in $[0,1]$ and depends on position $\bm{x}$ in $\Omega$, as illustrated in Figure 2(b). This field characterizes the ratio of the subset of the surface of the liquid layer that is covered by particles to the total area of that surface and is thought to influence the overall fracture pattern, crack kinking, and branching.11 We assume that the particles are rigid, so that any local contraction or expansion of the monolayer is accommodated solely by variations of the surface area not covered by particles. Thus, given an initial packing fraction $\phi_{0}$, the packing fraction $\phi$ depends on the trace, $\mbox{tr}\,\bm{\varepsilon}$, of the strain $\bm{\varepsilon}$ through

[TABLE]

We further assume that surfactant damage acts only to degrade the tensile resilience of the monolayer and that crack propagation is prohibited under compression. This is achieved by employing a spectral decomposition 20, 17 of $\bm{\varepsilon}$ into positive and negative components $\bm{\varepsilon}_{+}$ and $\bm{\varepsilon}_{-}$. It is then possible to define tensile and compressive strain-energy densities $W_{+}$ and $W_{-}$ through

[TABLE]

where $E>0$ and $0<\nu<1$ are the Young’s modulus and Poisson’s ratio of the undamaged monolayer and the positive and negative trace operations $\mbox{tr}\,_{+}$ and $\mbox{tr}\,_{-}$ are defined in accord with

[TABLE]

and

[TABLE]

We assume that the free-energy density of the monolayer is a function $\psi$ of $\bm{\varepsilon}$, $d$, and $\nabla d$ with the form

[TABLE]

As in conventional phase-field models for brittle fracture,17, 18 the tensile contribution to $\psi$ decays quadratically with the surfactant damage $d$. The remaining contributions of the energy are, however, nonstandard and therefore require further discussion.

The compressive contribution to $\psi$ accounts for increases in compressive energy that accompany increases in the packing fraction through the jamming factor $\xi$. This factor penalizes compression beyond a jamming threshold $\phi_{j}$ and prevents packing fractions from exceeding a maximum value $\phi_{m}$. The variation of $\xi$ with $\phi$ — and, thus, with reference to (2), the initial packing fraction $\phi_{0}$ and strain $\bm{\varepsilon}$ — is illustrated in Figure 3. Its value starts at $\xi=1$ for $0\leq\phi\leq\phi_{j}$ and increases monotonically for $\phi_{j}\leq\phi<\phi_{m}$, exhibiting a vertical asymptote as $\phi\to\phi_{m}$. The thresholds $\phi_{j}=0.84$ and $\phi_{m}=0.9$ correspond to a random close-packing in two space dimensions21 and the maximum packing density for two-dimensional discs, respectively. Numerically, it is more convenient to use expressions that effectively raise the compressive contribution without becoming unbounded. In this work, we use an expression with the form $1.0+f(1.0+\text{tanh}((\phi-0.5(\phi_{m}-\phi_{j}))/l)$, with $f$ being the factor of amplification and $l$ being the length scale controlling the width of the regularization.

For the contribution to $\psi$ associated with fracture, we choose a modified version

[TABLE]

of the expression proposed by Miehe et al. 17, with $G$ defined by

[TABLE]

where $\gamma_{f}$ and $\gamma_{s}$ are the surface tension of the liquid layer and the surfactant, respectively, $G_{0}>0$ is a constant that represents the fracture toughness of the undamaged monolayer, and $\lambda>0$ is a constant proportional to the characteristic thickness of a layer between damaged and undamaged material. This modification accounts for a reduction in the fracture toughness with increasing surfactant concentration. Accordingly, it is convenient to introduce

[TABLE]

as a measure of the reduced fracture toughness of the monolayer.

Finally, the surfactant contribution to $\psi$ accounts for the interplay between the monolayer and the surfactant (Figure 4) and is assumed to be of the form

[TABLE]

where $F_{0}>0$ is a constant.

2.2 Governing equations

Following Silva et al.19 (but neglecting inertia), the governing equations of the phase-field model consist of the macroforce balance

[TABLE]

and the microforce balance

[TABLE]

where $\beta\geq 0$ represents the kinetic modulus that controls the rate at which cracks can propagate through the monolayer. The alternative on the right-hand side of (12) embodies the requirement that, consistent with experimental observations, cracks that form in the monolayer never heal. This requirement takes the form of the constraint $\dot{d}\geq 0$, and a reactive microforce $\pi_{r}$ is needed to ensure satisfaction of that constraint. In particular, $\pi_{r}$ vanishes for $\dot{d}>0$ and is determined by the left-hand side of (12) for $\dot{d}=0$.

The kinetic modulus $\beta$ accounts for two effects. First, it captures the capacity of the surfactant to penetrate the particulate monolayer and thereby generate damage, which is influenced by $\phi$ and the radius of the particles. As $\phi$ and the particle size increase, the higher density of particles and the increased tortuosity impede surfactant spreading. Second, it captures the resistance of the underlying liquid layer to rearrangements of the particles, a resistance which is directly related to the viscosity $\mu_{f}$ of the liquid comprising that layer. For the foregoing reasons, we assume that $\beta$ increases with $\phi$, $\mu_{f}$, and the mean particle radius $\bar{r}_{p}>0$ in accord with the relation

[TABLE]

where $C>0$ is a constant and $\beta_{0}$ represents the kinetic modulus for the case of maximum packing $\phi=\phi_{m}$.

Since the kinetic modulus $\beta$ defined in (13) is positive, damage only increases when the left-hand side of (12) is positive. Accordingly, the evolution equation (12) for $d$ becomes

[TABLE]

where, given a scalar-valued quantity $h$,

[TABLE]

denotes its Macaulay bracket.

With reference to the right-hand side of the definition (6) of the free-energy density $\psi$ and using (7)–(10) and (13), we find the governing equations (11) and (14) can be written as

[TABLE]

and

[TABLE]

respectively, and where $\bm{I}$ denotes the (two-dimensional) identity tensor. We consider (16) and (17) subject to the zero displacement condition

[TABLE]

on the boundary $\partial\Omega$ of $\Omega$ and, letting $\bm{n}$ denote a unit normal on $\partial\Omega$, the natural boundary condition

[TABLE]

Additionally, we impose the initial conditions

[TABLE]

with $\phi_{0}$ representing the initial distribution of particles and $d_{0}$ the initial damage.

2.3 Model characterization: dimensionless number $\chi$

Introducing characteristic measures $L$ and $T$ of length and time, we define dimensionless counterparts $\bm{x}^{*}$ and $t^{*}$ of $\bm{x}$ and $t$ by

[TABLE]

Thus, defining dimensionless versions

[TABLE]

of the gradient, Laplacian, and partial time-derivative, dimensionless parameters,

[TABLE]

and dimensionless strain-energy densities

[TABLE]

we arrive at dimensionless counterparts

[TABLE]

and

[TABLE]

of the governing equations (16) and (17).

While we have not undertaken an exhaustive suite of simulations spanning the complete range of parameter space associated with the six dimensionless variables in (23), the extensive testing that we have conducted indicates that the fracture pattern is largely governed by the dimensionless driving force $F^{*}$ and the dimensionless fracture toughness $G^{*}$ of the monolayer. In particular, the number of fractures in the final configuration appears to be dictated by the ratio

[TABLE]

In Section 4.2, we discuss the threshold levels of the dimensionless parameter $\chi$ that delineate different regimes, ranging from no fractures at all to a configuration with multiple branches.

The expression (27) for $\chi$ can be further simplified by invoking the approximate scaling of Young’s modulus with the surface tension as derived by Vella et al.8 (i.e., $E\propto\gamma_{f}/\bar{r}_{p}$). The expression for $\chi$ can therefore be rewritten to isolate the influence of the surface tensions, yielding

[TABLE]

Daniels et al.22 proposed a similar relationship between the elastic and Marangoni energies to explain fracture patterns in agar gels. In that work, the number of crack branches was found to scale linearly with the difference in surface tension. As described in Section 4, our experimental observations and model-based simulations suggest the system considered here exhibits a much stronger sensitivity to surface tension differences.

3 Experiments

3.1 Materials

All videos from experimental work were acquired using a Phantom V641 high-speed camera equipped with an AF Nikkor 50 mm f/1.8 D lens. Videos were saved using PCC software provided by Phantom. A 5 W white LED board used to illuminate the mixture was supplied by VANCO (series #33342). Dispersion of silica microballoons (SIG MFG) was performed using the Active Motif Q120AM (120 W, 20 kHz) ultrasonicator equipped with a CL18 3.2 mm probe. Oleic acid (surfactant) was supplied by WAKO Chemicals and the acetone used as a cosolvent was supplied by Nacalai Tesque. The surfactant was dropped using Terumo 2.5 mL syringes (SS-02SZ) equipped with a 25 G (0.50 x 16 mm) needle. Surface tension measurements of mixtures were obtained using the KSV NIMA LB Small trough (33473).

3.2 Methodology

Figure 5 shows a schematic of the experimental setup. A clean low-form cylindrical glass vessel 17 cm in diameter and 9 cm in height filled to a height of 2.5 cm with milli-Q water (or mixed with a cosolvent) was placed atop an LED light panel. Microballoons were carefully weighed to either 0.150 g or 0.250 g, introduced at the air-water interface, and dispersed by means of ultrasonication at 25% power to form a particulate monolayer. Using a 2.5 mL syringe fitted with a 25 G metal tip, a single drop of surfactant (oleic acid with or without solvent) was introduced to the approximate center of the vessel at time t = 0 s. The water immiscible surfactant was observed to spread and push the microballoon particles radially outwards resulting in compaction of the particulate monolayers. Observations were captured using a Phantom high-speed camera at different frame rates. Data analysis was performed by first converting raw data files (.cine) to multipage image files (.tif) and these files were subsequently analyzed using ImageJ. In cases where the surfactant or the bulk water phase was mixed with a solvent to reduce surface tension, acetone was added carefully and aliquots of the liquids were quickly transferred to an LB trough for surface tension measurements.

4 Results and discussion

4.1 Pattern comparison

We discretize the model described in Section 3 by recasting the evolution equations in variational form and approximating the displacement and surfactant damage fields with finite element basis functions. Initial-boundary value problems are constructed over circular domains representative of the petri dishes used in the experiments. The experiments have a stochastic aspect corresponding to the initial particle locations. This effect is modeled at the continuum scale by using initial packing fields with spatial variability. Specifically, we numerically construct initial packing fraction fields over the circular domain that have a mean value of $\bar{\phi}=0.795$ with a 6% variation, corresponding to uniform random distributions with minimum and maximum packing fractions of 0.75 and 0.84, respectively. The initial drop of surfactant is approximated by a sufficiently small (on the order of $\lambda$), empty circular region ($\phi_{0}=0$) concentric with the center of the dish. Table 2 provides the reference values of the parameters used for the simulations in this section. All simulations were performed using MOOSE (Multiphysics Object-Oriented Simulation Environment), a finite-element framework primarily developed by Idaho National Laboratory.23

A first qualitative assessment of the model is shown in Figure 6. In this example, the parameter $\chi$ is adjusted to reproduce a three-branch pattern, which is a commonly observed configuration in the experiments for our given set of parameters. The result indicates that the model captures the general fracture process and reproduce the final three-branch configuration. Additionally, the cracks are sensitive to the packing fraction structure, which gives rise to more subtle features. These features include the particular angle in which the three-branch pattern is generated, the kinking and bending of cracks when gradients in the particle density are encountered, and the bifurcation of a crack when it encounters a high density region directly in line with its propagation. In this work, stochastic features of the experiments are approximated through random variations in the initial packing fraction field. The specific details of the numerically generated fracture patterns thus stem from the initial conditions and are sensitive to the particular random seed used in the simulations. As a result, our simulation results are only designed to capture the generic features of the experiments. Nevertheless, the model reproduces the number of cracks observed in the experiments and, at a qualitative level, the features of crack kinking and the formation of secondary branches.

We now provide a detailed description of a representative simulation in which a three-branch fracture pattern arises. Figure 7 shows the results from a simulation of this common pattern, obtained with the parameters given in Table 2 and $\chi=0.0125$ (see Table 3 for the correspondent dimensionless quantities). The columns in Figure 7 show the evolution of three fields, namely the surfactant damage, the tensile energy, and the packing fraction, along with a comparison with snapshots of a matching experiment at comparable moments in time. The distribution of $d$ closely matches the fractured portion of the particle raft in the experiment. The tensile energy field $W^{+}$, which is representative of the processes that underpin crack propagation, transitions from an initially circular tensile distribution due to the surfactant force to one with broken symmetry depending on the available energy in the system. The tensile energy drives the three-branches until they eventually arrest due to an increase in the compressive energy. The field $\phi$ depicts the evolution of the initial packing fraction. As the fractured region expands, the monolayer contracts (to approximately $\phi=0.87$) and the mobility of the surfactant in the fractured region ($\phi=0.3$) is increased.

4.2 Phase diagram for fracture pattern characterization

We now study whether the dimensionless parameter $\chi$ can be used to construct a phase diagram for the fracture patterns that are produced. In so doing, we eliminate the random aspect of the previous simulations and consider initial packing fractions that are spatially uniform with $\phi_{0}=0.80$. As illustrated in Figure 8, the results allow us to identify threshold levels of $\chi$ that delineate regions of a phase diagram for the fracture patterns obtained through simulations. The threshold $\chi=0.005$ defines a limit below which fracture is not observed. As $\chi$ is increased, the fracture pattern starts to develop as a single crack, and then into three-branch shapes with further increases. We note that the three-branch configuration is quite ubiquitous in nature, as the presence of acute angles and the compression generated by the three cracks stabilizes the system. For larger values of $\chi$, a fourth and a fifth branch are gradually realized. Larger values of $\chi$ are required to produce larger number of crack branches, resulting in unstable fracture patterns that cannot be clearly delineated in well defined regions of our phase diagram.

To validate the zones and thresholds predicted by the simulations, we first calibrated our model to experiments using oleic acid spreading over a pure water layer. The resulting fracture pattern is shown in the bottom of the center column of Figure 9, which presents images of the final configurations for three different experiments, along with simulation results for comparable values of $\chi$. Numerical simulations indicate that such a pattern can be obtained using $\chi=0.018$, which we label as $\chi_{ref}$. As indicated by expression (28), changing the surface tension difference corresponds to modulating the value of $\chi$ in the model. Accordingly, we conducted two additional experiments. In the first of these experiments, the surface tension of the liquid layer was reduced from $72$ to $56$ mN/m (colsolvent). In the second, the surface tension of the surfactant was reduced from $40$ to $28$ mN/m (acetone). These changes in surface tension correspond approximately to modifying factors of $1/3$ and $2$, respectively, on the reference $\chi_{\textit{ref}}=0.018$. We denote the corresponding levels of $\chi$ as $\chi_{\textit{down}}$ and $\chi_{\textit{up}}$.

As illustrated in Figure 9, simulations based on the values of $\chi_{\textit{up}}$ and $\chi_{\textit{down}}$ yield fracture patterns that are remarkably consistent with the corresponding experiments. Simulations based on doubling $\chi$ from $\chi_{\textit{ref}}$ to $\chi_{\textit{up}}$ result in an increase in the number of fracture crack branches, from three to five. Conversely, simulations based on a reduction of $\chi$ from $\chi_{\textit{ref}}$ to $\chi_{\textit{down}}$ result in a decrease in the number of crack branches, from three to two. In each case, a qualitative match with the corresponding experiment is clear. As $\chi$ is increased, the model presents a gradual change between patterns, from the formation of a single crack to a configuration with multiple branches. The model nicely reproduces the central damage zone and multiple crack branches observed in experiments for the range of $\chi$ in this study ($10^{-4}$–$10^{-1}$).

This result is particularly important because the conventional view has been that the packing fraction distribution is the primary variable governing the morphological features of the crack pattern. Our results suggest that, although the randomness of the initial packing fraction can underpin features such as crack deflections and arrest, the branching and general patterns observed in the experiments correspond to a much more fundamental principle related to the balance between the surfactant force and the fracture toughness. In other words, the difference in surface tension between the surfactant and underlying liquid are central to determining the crack patterns that emerge.

4.3 Temporal scaling

We now investigate the extent to which the model captures the speed of crack propagation. We examine the overall time required to fracture the particulate monolayer as well as the temporal scaling of the process. The snapshots in Figure 10 show two fracture patterns that evolve at very different rates, but for which the final configurations are visually indistinguishable. For comparison purposes, we show our results using a dimensional version of the model (i.e., viscosity in Pa s). Previous experimental observations have indicated that, for the same initial packing fractions, the geometry of the fracture pattern is insensitive to the viscosity of the underlying liquid layer.10 In the experiments, the liquid viscosity can be adjusted, for example, by increasing the ratio of glycerol in the glycerol-water mix based liquid layer, raising the viscosity from $10^{-3}$ Pa s to $10^{-1}$ Pa s. This is consistent with work showing that the rate at which the process occurs decreases by the same order of magnitude.10 We can capture this effect by adjusting the mobility defined in (13), which is inversely proportional to the viscosity $\mu_{f}$ of the liquid layer.

Additionally, experimental measurements indicate that the fracture proceeds in two stages. Crack growth rates that scale with the $3/4$ and higher powers of $t$ are observed in the early stages of the fracture process, followed by a gradual flattening as the crack advances and eventually arrests.10 Arrest may occur in response to the increase in packing around the perimeter of the domain, which gradually makes it more difficult to displace the particles of the monolayer. As shown in Figure 11, our model qualitatively reproduces this behavior with time, in terms of the normalized fractured area $A^{*}$ (namely the ratio of the area of the surface of the liquid layer that is exposed to the area of $\Omega$). As argued by Bandi et al.,11 if the radially receding front of $R(t)$ of a particle raft scales according to $R(t)\sim t^{3/4}$, then the particle free area proceeds as $R(t)^{2}\sim(t^{3/4})^{2}\sim t^{3/2}$. The growth rates in the early stages are slightly below $t^{3/2}$, which suggests the presence of the monolayer slightly retards the expansion of the surfactant from what would be expected over a pure liquid surface.24 Finally, we recognize a shift between the results from experiment and simulation with the initial expansion process being noticeably faster in the experiment. This discrepancy at early times could be related to the rapid expansion in the drop of surfactant that is not captured by the model. The subtleties of the physics at early times are likely sensitive to the volume and height of the drop of surfactant at the time of injection. Questions related to these subtleties are a subject for future work.

Following the rapid initial expansion, the dynamics slow down considerably once the particulate rafts jam together. As shown in the inset to Figure 11, the late stage dynamics exhibit approximately logarithmically slow relaxation over a decade in time. This kind of dynamics, also known as “creep” or “aging”, is observed in a variety of phenomena including granular compaction,25 electron glasses,26 polymer relaxation,27 and superconductors,28 and is now considered a hallmark of non-exponential, slow relaxation in amorphous media.29 Although the observation times in our model and experiments do not extend over several decades, they are consistent with the expected behavior of slow relaxation in frustrated, amorphous media.

5 Conclusions

We have presented a phase-field model that captures the main features of the surfactant-driven fracture of closely-packed particulate monolayers. The model and experimental results indicate that there is a competition between the spreading force of the surfactant and the fracture toughness of the monolayer that determines the number of fractures. The model gives rise to a dimensionless parameter that can be written in terms of the surface tensions, and which allows for a straightforward comparison with experimental conditions. Experiments were conducted to validate the model and delineate the regimes of fracture pattern as a function of that parameter.

Our model rests on a number of simplifying assumptions. In particular, we used the damage to the particulate monolayer as an indicator function for the surfactant concentration. As such, our model precludes investigations of the extent to which the surfactant can penetrate into the particulate monolayer network, ahead of (or behind) the crack front. We have also assumed that the strains exhibited by the monolayer are sufficiently small to justify modeling it as a linearly elastic solid. This assumption ceases to be reasonable if the surface tension difference is increased above some critical threshold. Moreover, its applicability is limited to systems in which the initial packing is low enough to allow the particulate layer to exhibit both loosely-packed and jammed behaviors as the surfactant spreads.11 In either case, transitions between linear and nonlinear response, or between fluid and solid behavior, are worthy of further study. Finally, the interplay between liquid and solid constituents in the monolayer may endow it with a viscoelastic behavior that could also be explored. As the rate of surfactant transport is a function of the Marangoni stress, modulating the surface tension difference might also be used to explore any rate dependence in the monolayer. Future work will focus on enhancing the model to incorporate these effects and to allow for detailed exploration of their consequences and significance.

6 Acknowledgments

Christian Peco, Yingjie Liu, and John E. Dolbow are grateful for the support of NSF grant CMMI-1537306, to Duke University. Wei Chen, M. M. Bandi, and Eliot Fried gratefully acknowledge support from the Okinawa Institute of Science and Technology Graduate University with subsidy funding from the Cabinet Office, Government of Japan. John E. Dolbow would also like to acknowledge support from the Okinawa Institute of Science and Technology during a sabbatical stay in late 2014, during which time this work began.