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

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 that of a linearly elastic solid. The model approximates the loss in tensile strength of the monolayer with increasing surfactant concentration 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.

---

1 Introduction

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

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. The interaction of the particles through the capillary bridges allows the monolayer to support both tension and compression. Moreover, 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⁹ 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 potential importance 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 force driving 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 Fig. 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 ϕ. As ϕ is increased, the properties of the liquid–vapor interface change from liquid to solid.¹¹

Fig. 1 Left, center: States before and after the introduction of surfactant. The surfactant is injected at the center, decreasing the surface tension and generating an advancing front at which fracture initiates. Right: Micrograph showing the microstructure of the closely-packed particulate monolayer.

Consistent with out interest in closely-packed systems, we consider situations in which the packing fraction is high (0.7 ≤ ϕ ≤ 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¹² are basic to understanding its response to surfactant-driven stresses. An estimate of those properties is provided by Vella et al.,⁸ based on experimental measurements and geometrical arguments. Surface wave experiments have been used to characterize the stretching and bending stiffnesses of particulate monolayers,¹³ which are found to present a soft granular character with nonclassical response under fluid-driven deformation.¹⁴

Experiments show that when such a system is stimulated by the localized introduction of a surfactant, an advancing front forms and fractures the monolayer.¹⁰ This type of surfactant-induced effect has also been observed in other systems, such as agar gels.¹⁵ The resulting fracture patterns, illustrated in Fig. 1, are reminiscent of those observed in classical brittle materials but with significant differences. Bandi et al.¹¹ 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¹⁶). In contrast, crack branching has been observed in particulate monolayers at crack speeds as low as 0.2% of the shear wave speed.¹⁰ 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 monolayer. Rather, experimental observations suggest that the fracture time scale is mostly governed by variations in the viscosity of the underlying liquid¹⁰ and the packing fraction.¹¹ 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.,¹¹ 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 large 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.¹⁷ and Borden et al.¹⁸ 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 needed 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 a 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–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 (Fig. 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 Ω. In Fig. 2(a), the dark portion of Ω represents the intact portion of the monolayer and the light portion of Ω 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 u and a scalar surfactant damage field d. As indicated in Fig. 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

Fig. 2 (a) Schematic of a fractured particulate monolayer. The fractured zone is represented as a smooth transition zone across which the damage field d takes values between zero and unity. (b) The microstructure of the elastic domain Ω is captured by an additional packing fraction field ϕ representing the solid to fluid ratio of the medium.

The microstructure of the monolayer is described by a scalar packing fraction ϕ that takes values in [0,1] and depends on position x in Ω, as illustrated in Fig. 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.¹¹ 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 ϕ₀, the packing fraction ϕ depends on the trace, trε, of the strain ε through

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^17,20 of ε into positive and negative components ε₊ and ε₋. It is then possible to define tensile and compressive strain-energy densities W₊ and W₋ through

where E > 0 and 0 < ν < 1 are the Young's modulus and Poisson's ratio of the undamaged monolayer and the positive and negative trace operations tr₊ and tr₋ are defined in accord withand

We assume that the free-energy density of the monolayer is a function ψ of ε, d, and ∇d with the form

As in conventional phase-field models for brittle fracture,^17,18 the tensile contribution to ψ 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 ψ accounts for increases in compressive energy that accompany increases in the packing fraction through the jamming factor ξ. This factor penalizes compression beyond a jamming threshold ϕ_j and prevents packing fractions from exceeding a maximum value ϕ_m. The variation of ξ with ϕ—and, thus, with reference to (2), the initial packing fraction ϕ₀ and strain ε—is illustrated in Fig. 3. Its value starts at ξ = 1 for 0 ≤ ϕ ≤ ϕ_j and increases monotonically for ϕ_j ≤ ϕ < ϕ_m, exhibiting a vertical asymptote as ϕ → ϕ_m. The thresholds ϕ_j = 0.84 and ϕ_m = 0.9 correspond to a random close-packing in two space dimensions²¹ and the maximum packing density for two-dimensional discs, respectively. Numerically, it is more convenient to use expressions that monotonically increase the compressive contribution without becoming unbounded. In this work, we use an expression with the form 1.0 + f(1.0 + tanh((ϕ − 0.5(ϕ_m − ϕ_j))/l)), with f being the factor of amplification and l being the length scale controlling the width of the regularization.

Fig. 3 Variation in the jamming factor ξ with the packing fraction ϕ. The jamming factor tends towards ∞ as ϕ surpasses the jamming threshold ϕ_j = 0.84 and approaches the maximum packing fraction ϕ_m = 0.9.

For the contribution to ψ associated with fracture, we choose a modified version

of the expression proposed by Miehe et al.,¹⁷ with G defined bywhere γ_f and γ_s are the surface tension of the liquid layer and the surfactant, respectively, G₀ > 0 is a constant that represents the fracture toughness of the undamaged monolayer, and λ > 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 introduceas a measure of the reduced fracture toughness of the monolayer.

Finally, the surfactant contribution to ψ accounts for the interplay between the monolayer and the surfactant (Fig. 4) and is assumed to be of the form

where F₀ > 0 is a constant.

Fig. 4 Representation of the surfactant force driving the system. A collection of particles, initially in equilibrium at a constant surface tension γ_f, are reached by the surfactant. The gradient between its lower surface tension γ_s and the surrounding media generates a Marangoni flow which pushes the particulate monolayer outward.

2.2 Governing equations

Following Silva et al.¹⁹ (but neglecting inertia), the governing equations of the phase-field model consist of the macroforce balanceand the microforce balancewhere β ≥ 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 ḋ ≥ 0, and a reactive microforce π_r is needed to ensure satisfaction of that constraint. In particular, π_r vanishes for ḋ > 0 and is determined by the left-hand side of (12) for ḋ = 0.

The kinetic modulus β incorporates two effects. First, it captures the capacity of the surfactant to penetrate the particulate monolayer and thereby generate damage, which is influenced by ϕ. As ϕ increase, the higher density of particles 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 μ_f of the liquid comprising that layer. For the foregoing reasons, we assume that β increases with ϕ and μ_f in accord with the relation

where C > 0 is a dimensionless constant and β₀ represents the kinetic modulus for maximal packing ϕ = ϕ_m.

Since the particular choice of the kinetic modulus β defined in (13) is positive, damage only increases when the left-hand side of (12) is positive. Accordingly, the evolution eqn (12) for d becomes

where, given a scalar-valued quantity h,denotes its Macaulay bracket.

With reference to the right-hand side of the definition (6) of the free-energy density ψ and using (7)–(10) and (13) with reference to the parameters listed in Table 1, we find the governing eqn (11) and (14) can be written as

andrespectively, and where I denotes the (two-dimensional) identity tensor. We consider (16) and (17) subject to the vanishing displacement condition

u = 0 (18)

on the boundary ∂Ω of Ω and, letting n denote a unit normal on ∂Ω, the natural boundary condition

∇d·n = 0. (19)

Additionally, we impose the initial conditions

ϕ(x,0) = ϕ₀(x), d(x,0) = d₀(x), (20)

with ϕ₀ representing the initial distribution of particles and d₀ the initial damage.

Table 1 List of parameters in the model

Parameter	Description
ϕ j	Jamming packing fraction of monolayer
ϕ m	Maximum packing fraction of monolayer
r̄ p	Average particle radius
γ f	Surface tension of liquid layer
γ s	Surface tension of surfactant
G 0	Fracture toughness of undamaged monolayer
F 0	Surfactant force constant
E	Young's modulus of monolayer
ν	Poisson's ratio of monolayer
μ f	Dynamic viscosity of monolayer
β 0	Kinetic parameter

---

2.3 Model characterization: dimensionless number χ

Introducing characteristic measures L and T of length and time, we define dimensionless counterparts x* and t* of x and t byThus, defining dimensionless versionsof the gradient, Laplacian, and partial time-derivative, dimensionless parameters,and dimensionless strain-energy densitieswe arrive at dimensionless counterpartsandof the governing eqn (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

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

The expression (27) for χ can be further simplified by invoking the approximate scaling of Young's modulus with the surface tension as derived by Vella et al.⁸ (i.e., E ∝ γ_f/r̄_p). The expression for χ can therefore be rewritten to isolate the influence of the surface tensions, yielding

Daniels et al.²² 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 × 16 mm) needle. Surface tension measurements of mixtures were obtained using the KSV NIMA LB Small trough (33473).

3.2 Methodology

Fig. 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.

Fig. 5 Setup: a clean low-form cylindrical vessel containing distilled water was placed on an LED tablet. A camera suspended vertically above records particulate monolayer compaction dynamics at the air–water interface when surfactant was introduced using a 25 G steel needle.

4 Results and discussion

We discretize the model described in Section 2 by recasting the evolution equations in a variational form and approximating the displacement and surfactant damage fields with finite element basis functions. Temporal integration is performed with a second-order accurate, unconditionally stable backward Euler algorithm (BDF2) with adaptive time stepping. The coupled system of nonlinear equations are solved with the Newton method in a staggered manner (i.e., alternating the updates for the displacement and damage fields).

4.1 Pattern comparison

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 [small phi, Greek, macron] = 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 λ), empty circular region (ϕ₀ = 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 developed primarily by Idaho National Laboratory.²³

Table 2 Dimensional variables^8,10

Property	Description	Value	Unit
r d	Radius of Petri dish	8.5 × 10^−2	m
ϕ 0	Range of initial packing fraction	(7.5–8.4) × 10^−1	None
ϕ m	Packing fraction of jammed monolayer	9.0 × 10^−1	None
G 0	Fracture toughness of undamaged monolayer	5.7 × (10^−4–10^−2)	N m^−1
F 0	Surfactant force constant	1.0 × (10^−1–10^2)	m^−1
γ f	Surface tension of liquid layer (water)	7.2 × 10^−2	N m^−1
γ s	Surface tension of surfactant (oleic acid)	4.0 × 10^−2	N m^−1
μ f	Range of dynamic viscosity of liquid layer	1.0 × (10^−3–10^−1)	Pa s
r̄ p	Mean particle radius	1.0 × 10^−4	m
E	Young's modulus of monolayer	3.178 × 10^3	Pa
ν	Poisson's ratio of monolayer		None
β 0	Kinetic parameter	1.0	Pa s
λ	Characteristic length scale	1.0 × 10^−2	m

---

A first qualitative assessment of the model is shown in Fig. 6. In this example, the parameter χ 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 reproduces 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 are designed to capture only 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.

Fig. 6 Comparison between the result of a simulation (left) and an experimental image (right) for a three-branch fracture pattern. The model captures the main traits of the final configuration, as represented by the number of fractures, and more detailed features corresponding to variations in the distribution of particles, such as secondary fracturing, crack bending, and branching. The particular orientation of the fracture pattern shown in the simulation result is arbitrary, as it is sensitive to the random seed used to construct the initial packing field.

We now provide a detailed description of a representative simulation in which a three-branch fracture pattern arises. Fig. 7 shows the results from a simulation of this common pattern, obtained with the parameters given in Table 2 and χ = 0.0125 (see Table 3 for the correspondent dimensionless quantities). The columns in Fig. 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 strain-energy density 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 strain-energy density drives the three-branches until they eventually arrest due to an increase in the compressive energy. The field ϕ depicts the evolution of the initial packing fraction. As the fractured region expands, the monolayer contracts (to approximately ϕ = 0.87) and the mobility of the surfactant in the fractured region (ϕ = 0.3) increases.

Fig. 7 Numerical fields obtained for a simulation of a three-branch fracture pattern. From left to right, each column shows the initial (top) and final (bottom) states of surfactant damage d, tensile strain-energy density W₊ and ϕ fields. The rightmost column shows the reference experiment, obtained by adding oleic acid to a 0.25 g mass monolayer of particles laying on a 18% salt saturation water solution.

Table 3 Dimensionless parameters

Parameter	Value
G*	10^−3–10^−1
F*	10^−4–10^−1
L*	8.5 × 10^2
	2.7 × (10^−8–10^−6)
γ*	5.6 × 10^−1
λ*	1.5 × 10^−1

---

4.2 Phase diagram for fracture pattern characterization

We now study whether the dimensionless parameter χ 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 ϕ₀ = 0.80. As illustrated in Fig. 8, the results allow us to identify threshold levels of χ that delineate regions of a phase diagram for the fracture patterns obtained through simulations. The threshold χ = 0.005 defines a limit below which fracture is not observed. As χ 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 χ, a fourth and a fifth branch are gradually realized. Larger values of χ 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.

Fig. 8 Phase diagram for fracture pattern characterization. Threshold levels of the dimensionless number χ = (F*)²/G* delineate various patterns. In the region below χ = 0.005 (zone i), the initially fractured region does not grow with time. As χ increases, the fracture patterns gradually change from the single branch configuration (zone ii) to the three (zone iii), four (zone iv) and five (zone v) branch configurations.

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 Fig. 9, which presents images of the final configurations for three different experiments, along with simulation results for comparable values of χ. Numerical simulations indicate that such a pattern can be obtained using χ = 0.018, which we label as χ_ref. As indicated by expression (28), changing the surface tension difference corresponds to modulating the value of χ 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 χ_ref = 0018. We denote the corresponding levels of χ as χ_down and χ_up.

Fig. 9 Comparison between fracture patterns from experiments (top) and simulations (bottom). In all cases, the final configuration is shown. Simulation results were obtained for values of the dimensionless parameter χ (χ_down = 0.006, χ_ref = 0018, χ_up = 0.035), achieved by modulating the surface tension of the liquid layer and the surfactant used in each of the experiments, consistent with (28).

As illustrated in Fig. 9, simulations based on the values of χ_up and χ_down yield fracture patterns that are remarkably consistent with the corresponding experiments. Simulations based on doubling χ from χ_ref to χ_up result in an increase in the number of crack branches, from three to five. Conversely, simulations based on a reduction of χ from χ_ref to χ_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 χ 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 (10⁻⁴–10⁻¹) of χ in this study.

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 interplay 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 Fig. 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.¹⁰ 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⁻³ Pa s to 10⁻¹ Pa s. This is consistent with work showing that the rate at which the process occurs decreases by the same order of magnitude.¹⁰ We can capture this effect by adjusting the mobility implicitly defined in (13), which is proportional to the viscosity μ_f of the liquid layer.

Fig. 10 Identical fracture patterns are produced at different speeds due to changes in the viscosity μ_f of the underlying 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.¹⁰ 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 Fig. 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 Ω). As argued by Bandi et al.,¹¹ if the radially receding front of R(t) of a particle raft scales according to R(t) ∼ t^3/4, then the particle free area proceeds as R(t)² ∼ (t^3/4)² ∼ 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.²⁴ 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, a process 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.

Fig. 11 Comparison between simulation and experiment for the temporal scaling for the normalized fractured area A* (the ratio of the portion of the area of the liquid layer that is exposed to the area of Ω), indicating two regimes of crack evolution. For reference, a t^3/2 curve is provided (red dashed line). The inset presents a zoom of the latter stages using semi-log axes, revealing a slow logarithmic relaxation of the fracture pattern.

Following the rapid initial expansion, the dynamics slow down considerably once the particulate rafts jam together. As shown in the inset to Fig. 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,²⁵ electron glasses,²⁶ polymer relaxation,²⁷ and superconductors,²⁸ and is now considered a hallmark of non-exponential, slow relaxation in amorphous media.²⁹ 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 numerical simulations 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.¹¹ 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.

Acknowledgements

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.

References

1. A. Varshney, A. Sane, S. Ghosh and S. Bhattacharya, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 31402.
2. S. U. Pickering, J. Chem. Soc., Trans., 1907, 91, 2001–2021.
3. B. P. Binks and P. D. I. Fletcher, Langmuir, 2001, 17, 4708–4710.
4. P. Cicuta, E. J. Stancik and G. G. Fuller, Phys. Rev. Lett., 2003, 90, 236101.
5. P. Aussillous and D. Quéré, Nature, 2001, 411, 924–927.
6. N. Tsapis, D. Bennett, B. Jackson, D. A. Weitz and D. A. Edwards, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12001–12005.
7. A. B. Subramaniam, M. Abkarian and H. A. Stone, Nat. Mater., 2005, 4, 553–556.
8. D. Vella, P. Aussillous and L. Mahadevan, Europhys. Lett., 2004, 68, 212–218.
9. L. Scriven and C. Sterling, Nature, 1960, 187, 186–188.
10. D. Vella, H. Y. Kim, P. Aussillous and L. Mahadevan, Phys. Rev. Lett., 2006, 96, 1–4.
11. M. M. Bandi, T. Tallinen and L. Mahadevan, Europhys. Lett., 2011, 96, 36008.
12. P. Cicuta and D. Vella, Phys. Rev. Lett., 2009, 102, 138302.
13. C. Planchette, E. Lorenceau and A.-L. Biance, Soft Matter, 2012, 8, 2444–2451.
14. C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Phys. Rev. X, 2015, 5, 11020.
15. C. Spandagos, T. B. Goudoulas, P. F. Luckham and O. K. Matar, Langmuir, 2012, 28, 7197–7211.
16. E. H. Yoffe, Philos. Mag., 1951, 42, 739–750.
17. C. Miehe, M. Hofacker and F. Welschinger, Comput. Methods Appl. Mech. Eng., 2010, 199, 2765–2778.
18. M. J. Borden, C. V. Verhoosel, M. A. Scott, T. J. R. Hughes and C. M. Landis, Comput. Methods Appl. Mech. Eng., 2012, 217-220, 77–95.
19. M. N. da Silva Jr., F. P. Duda and E. Fried, J. Mech. Phys. Solids, 2013, 61, 2178–2195.
20. C. Miehe, Comput. Struct., 1998, 66, 37–43.
21. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 11306.
22. K. E. Daniels, S. Mukhopadhyay, P. J. Houseworth and R. P. Behringer, Phys. Rev. Lett., 2007, 99, 124501.
23. D. Gaston, C. Newman, G. Hansen and D. Lebrun-Grandié, Nucl. Eng. Des., 2009, 239, 1768–1778.
24. O. E. Jensen, J. Fluid Mech., 1995, 293, 349–378.
25. E. Ben-Naim, J. B. Knight, E. R. Nowak, H. M. Jaeger and S. R. Nagel, Phys. D, 1998, 123, 380–385.
26. A. Amir, Y. Oreg and Y. Imry, Phys. Rev. Lett., 2009, 103, 126403.
27. J. D. Ferry, Viscoelastic Properties of Polymers, John Wiley and Sons, 1980.
28. A. Gurevich and H. Küpfer, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 6477–6487.
29. A. Amir, Y. Oreg and Y. Imry, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 1850–1855.

---