RESEARCH ARTICLE | MAY 07 2025
The diffuse solid method for wetting and multiphase fluid simulations in complex geometries
Special Collection: Proceedings of the 33rd Conference on Discrete Simulation of Fluid Dynamics
Fandi Oktasendra ; Michael Rennick ; Samuel J. Avis ; Jack R. Panter ; Halim Kusumaatmaja 
Physics of Fluids 37, 052103 (2025)
https://doi.org/10.1063/5.0267171

View

Export

Online Citation

23 February 2026 16:50:51

Articles You May Be Interested In
Beyond primitive bionic surfaces: Bionic structures coupled surface with superhydrophobicity and programmable directional adhesion
Physics of Fluids (June 2025)
On interfacial deformation and resistance characteristics over superhydrophobic surfaces: A numerical study
Physics of Fluids (November 2024)
On the drag reduction of an inclined wing via microstructures with the immersed boundary-lattice
Boltzmann flux solver
Physics of Fluids (August 2023)

pubs.aip.org/aip/pof
The diffuse solid method for wetting and multiphase fluid simulations in complex geometries
Cite as: Phys. Fluids 37, 052103 (2025); doi: 10.1063/5.0267171
Submitted: 22 February 2025 . Accepted: 17 April 2025 .
Published Online: 7 May 2025
Fandi Oktasendra,1,2	Michael Rennick,1,3	Samuel J. Avis,1	Jack R. Panter,4,a) and Halim Kusumaatmaja1,3,b)
AFFILIATIONS
1Department of Physics, Durham University, Durham DH1 3LE, United Kingdom
2Department of Physics, Universitas Negeri Padang, Padang 25131, Indonesia
3Institute for Multiscale Thermofluids, School of Engineering, University of Edinburgh, Edinburgh EH9 3FD, United Kingdom 4School of Engineering, Mathematics and Physics, University of East Anglia, Norwich NR4 7TJ, United Kingdom
Note: This paper is part of the Special Topic, Proceedings of the 33rd Conference on Discrete Simulation of Fluid Dynamics. a)Electronic mail: jack.panter@uea.ac.uk
b)Author to whom correspondence should be addressed: halim.kusumaatmaja@ed.ac.uk
ABSTRACT
We develop a diffuse solid method that is versatile and accurate for modeling wetting and multiphase flows in highly complex geometries. In this scheme, we harness N þ 1-component phase field models to investigate interface shapes and flow dynamics of N fluid components, and we optimize how to constrain the evolution of the component employed as the solid phase to conform to any pre-defined geometry. Implementations for phase field energy minimization and lattice Boltzmann method are presented. Our approach does not need special treatment for the fluid–solid wetting boundary condition, which makes it simple to implement. To demonstrate its broad applicability, we employ the diffuse solid method to explore wide-ranging examples, including droplet contact angle on a flat surface, particle adsorption on a fluid– fluid interface, critical pressure on micropillars and on Salvinia leaf structures, capillary rise against gravity, Lucas-Washburn’s law for capillary filling, and droplet motion on a sinusoidally undulated surface. Our proposed approach can be beneficial to computationally study multiphase fluid interactions with textured solid surfaces that are ubiquitous in nature and engineering applications.
VC 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0267171

I.		INTRODUCTION
Inspired by the extraordinary diversity of surface structures adapted to manipulating liquids in nature,1,2 there is a rapid and increasing drive to fabricate complex physical and chemical textures with functional wetting properties.3 This is particularly prevalent in the field of surface science, which has seen a range of recent innovations across fog harvesting,4 robust liquid-repellent surfaces,5 directional wetting,6 and wettability control for sensing.7 This has been enabled by advances in fabrication methods of increasingly complex surface microstructures, such as additive manufacturing (3D print- ing),8 photolithography,9 and self-assembly.10 However, the interest in understanding wetting interactions with complex surface structures is not only limited to fundamental wetting physics but also spans a diverse range of fields, from open microfluidics in biomedicine11 to sub-milimeter manufacturing in engineering12 and carbon storage in geology.13
To complement experimental and theoretical efforts, there is a need for efficient and accurate simulation methods that are capable of modeling the statics and dynamics of capillary phenomena on highly complex solid structures. Ideal simulation methods must achieve the following aims: (1) accuracy in modeling surface tensions and contact angles, (2) generality and accuracy in describing the solid geometry, and (3) computational efficiency. In this work, to achieve these aims, we present the diffuse solid method, which is applied both to the lattice Boltzmann method (for fluid dynamics), and to a phase-field energy minimization method (for static equilibria). In this, we use a N þ 1-phase diffuse interface method that models N evolving fluid components and 1 solid component that we limit the evolution of. A principal outcome is to show how to control this partial evolution of the solid to conform to any user-defined geometry, without compromising the accuracy of the liquid contact angles. While we demonstrate the scheme for N ¼ 2, the scheme can be easily generalized to an arbitrary number of components. Importantly, the interfacial energy balance is intrinsically accounted for by the phase field free energy, and no separate scheme is required for the fluid–solid wetting boundary condition.
Previously, a number of alternative simulation methods have been developed toward achieving accurate wetting simulations on curved and complex boundaries. At the highest level, these simulation strategies are split into two classes: sharp interface models, where the fluid interface is modeled as infinitely thin and explicitly tracked during simulations;14,15 and diffuse interface models, where each fluid phase is associated with a smoothly varying order parameter, and fluid–fluid interfaces have a characteristic thickness.16,17 Sharp interface models are often more efficient computationally; however, they are cumbersome for problems that entail topological changes in the fluid interface shapes, such as due to droplet breakup and coalescence, as well as interface pinning and depinning events.17–21 Since these effects are essential to capture for the problems of interest in this work, here, we shall focus on diffuse interface models.
For diffuse interface models, there are already well-established methods for incorporating fluid–solid interactions on flat surfaces,22–25 including liquid contact angle and contact line pinning and depinning due to regions of different surface chemistry. However, current extensions to complex, and often curved, solid shapes are either cumbersome or introduce unwanted artifacts. This is a principal issue when employing diffuse interface models on methods that are based on fixed, regular grids, such as the lattice Boltzmann method. However, it is desirable to use such methods as they are highly efficient for parallel and graphical processing unit (GPU) computing and hence are becoming increasingly popular for solving large-scale simulation challenges.26–28
Within the lattice Boltzmann method, a common approach to multiphase flows on curved and complex solids is to treat the surface via the conceptually simple staircase approximation,29–31 or similarly, as in the case of the homogenized lattice Boltzmann method, by defining voxellated regions of space within the solid body to be momentum sinks.32 This, however, leads to interface pinning at the corners of the approximated shapes, and in problems where surface energies are important (such as capillary rise as we shall demonstrate), the solid– fluid interaction area is dramatically overestimated.
A key strategy to enable simulations on solid surfaces, which do not conform to the underlying (Eulerian) mesh, is the immersed boundary method, in which the solid surface is described via a superimposed non-conformal (Lagrangian) mesh. Recent advances have overcome core challenges to this method for use in multiphase flows (namely, poorly obeyed mass conservation and limited stable fluid density ratios) that now allow accurate wetting boundary condi- tions.33,34 Nonetheless, an inherent downside of the immersed boundary method is the need for additional interpolation steps between the Eulerian and Lagrangian meshes.
A contrasting approach, which entirely avoids the explicit description of the solid boundary, is fluid particle dynamics, in which the solid is modeled as a highly viscous liquid.35,36 A core advantage of this method is that no special wetting boundary conditions are required: the correct contact angles emerge implicitly from the diffuse interface interactions. Although this method is successful for, e.g., colloidal systems, where the solids are always spherical, this cannot achieve solid geometries with arbitrary user-defined shapes.
Sitting between the solid and liquid methods is the diffuse domain approach, currently employed in Cahn–Hilliard and Navier–Stokes solv- ers.37 In this, the solid domain is identified with an order parameter, exactly as for the fluid phases. In contrast to fluid particle dynamics, however, the solid order parameter is not updated in time and is instead fixed upon initialization. In order to accurately capture the contact angles on the solid phase, it is essential to have the correct solid–fluid diffuse interfacial profiles. Recently, this has been achieved through using a pre-simulation stage to relax the solid order parameter via a diffusion process.38,39 However, the diffuse domain approach suffers from an inherent limitation: relaxing the solid diffuse interface also allows the overall shape of the solid object to evolve away from the desired geometry. A principal achievement in our work is to develop a method of constraining the solid morphology, while allowing the correct diffuse interface to form. This method is essential in order to accurately model wetting on complex geometries and to allow implementation in free energy minimization and the lattice Boltzmann methods.
Following the description of our diffuse solid method in Sec. II, we rigorously demonstrate it’s broad applicability and accuracy across a range of benchmarks. In Sec. III, static phenomena are investigated, featuring droplets on flat surfaces and spherical particles at flat fluid interfaces as tests of the contact angle accuracy, equilibrium capillary rise to test the energy contribution accuracy, and critical intrusion pressures on micropillared surfaces as a test of the pressure implementation. A case study is then shown to show how our method can be applied to understand the wetting properties of nature-inspired surface structures. In Sec. IV, we advance to study dynamic phenomena via the lattice Boltzmann method: capillary filling in micro-channels and droplet motion on an undulating surface.
II.		DIFFUSE SOLID METHOD
A.		Phase field energy minimization
The diffuse solid method uses a phase field model with N þ 1 immiscible components in which the local composition of each component is represented by an order parameter Cn, where n ¼ f1;…; N þ 1g. Here, we focus on N ¼ 2 for concreteness. In this context, we designate one component, C1, as a solid phase and the rest, C2 and C3, as fluids, which shall be referred to as liquid and gas phases, respectively.
The ternary model (N þ 1 ¼ 3) used in this work is based on the work of Semprebon et al.,40 in which the thermodynamics of the system are described by a free-energy functional of the fluids that captures the immiscibility of the fluid components and the surface tensions between different fluids. The phase interface is modeled using a diffuse interface with a finite thickness, which yields a smooth transition region between two different phases.
For the phase field energy minimization approach, the total free energy of the system can be written as
WðCnÞ ¼ ðWbþWiþWcsþWcfÞdV; (1) where all of the energy densities are integrated over the system’s volume, V. The bulk energy density of the phases, Wb, is modeled using a double well potential

3
Wb¼	j2nCn2ð1
n¼1

2 nÞ :

(2)

Here, the jn ’s are tunable parameters related to the interfacial tension between different fluids. The form of potentials chosen in this model has two minima at Cn ¼ 0 and 1. We then relate Cn ¼ 0 to indicate the absence of component n and Cn ¼ 1 to indicate pure component n. We also normalize the total composition to unity so that C1 þ C2 þ C3 ¼ 1toobtainC3 ¼ 1 ^ C1 ^ C2, which reduces the number of the order parameter in the system that needs to be optimized. It is worth noting that this constraint allows C1 and C2, in principle, to have any value.
Wi is referred to as the gradient energy density, which accounts for the energy penalty for having interfaces. The gradient term takes the following form:
30
Wi ¼ j2n j$Cnj2; (3) n¼1 where j0n ’s are also tunable parameters whose values determine the contribution of the gradient terms to the total free energy.
In combination, Wb and Wi energetically penalize the formation of interfaces between each pair of phases, such that the surface tension between phases n and m, cnm, is given by
cnm ¼ a6 ðjn þ jmÞ:	(4)
a is the characteristic interface width,
a ¼ ðj0n þ j0mÞ=ðjn þ jmÞ;	(5) that describes the smooth variation of Cn across an interface. For a planar interface (e.g., spanning the z ^ y plane) between phases n and m, this profile takes the following form:
11 x
Cn ðxÞ¼2 þ 2 tanh 2a :	(6)
Here, the value of Cn varies between 0 and 1 across the interface, and we take the value of Cn ¼ 0:5 as the interfacial boundary between the two phases. We can further use the definition of a to simplify the simulation variables by relating j0fn;mg and jfn;mg via ðj0n þ j0mÞ ¼ a2ðjn þ jmÞ. In this work, we take j0n ¼ a2 jn and use a ¼ 1 (this and all other parameters are in lattice units, L.U., unless otherwise stated). While we will call a the interface width, it is noted that for a ¼ 1, the transition between Cn ¼ 0:1 and Cn ¼ 0:9 occurs over a range of 4:4a.
In order to relate the phase field model parameters j to physical interfacial characteristics (namely the surface tensions c), Eq. (4) can be inverted to yield
3
j1 ¼ ðc12 þ c13 ^ c23Þ;
a
3
j2 ¼ ðc12 ^ c13 þ c23Þ;	(7) a
3
j3 ¼ ð^c12 þ c13 þ c23Þ: a
Wcs is the constraining energy density that allows the system to preserve either the pressure difference between liquid and gas phases, DP, or the volume of the liquid phase, Vl, during the minimization. Following Oktasendra et al.,41 when the former is desired, the constraining potential takes the following form:
Wcs ¼^DPVl; (8)
whereas for the latter,
1
Wcs¼2kcsðVl^ V0Þ2:	(9)
Here, kcs > 0 is a constant and V0 is the target volume. The volume of the liquid phase, Vl, can be calculated as
Vl ¼	C2dV:
V

(10)

Thus far, the described model is similar to that we have previously published for phase field energy minimization.41–44 For the diffuse solid method, a key addition is the confining potential Wcf . This potential acts only on the component representing the solid, which throughout we choose to be phase 1. Wcf must achieve two purposes: (1) to confine phase 1 to a defined region of arbitrary geometry (which we call the “solid region,” S) and (2) to allow the interface of phase 1 to relax into the equilibrium diffuse interfacial profile. This latter point is essential to accurately define the surface tensions with the solid phase, as defined in Eq. (4). These can both be achieved by defining
> wcf if C1 > 0:5 and r 2 S;
WcfðrÞ ¼ < wcf if C1 < 0:5 and r 62 S; (11)
: 0 otherwise;
where r defines the position of each point in the system, and wcf ¼ bðC1 ^ 0:5Þ2. b is a positive variable controlling the strength of the potential. This form of Wcf energetically penalizes the solid phase from moving out of the solid region and likewise energetically penalizes fluid phases from entering the solid region, so centering the solid– fluid interface at the solid region boundary.
In the simulation, we discretize W in Eq. (1) into Lx ^ Ly ^ Lz cubic lattice points (called nodes), in which each node is associated with a value of fCngijk, where the indices i, j, andk indicate the spatial dimension in the x-, y-, and z-directions. The grid spacing, G, issetto be equal to the interface width a, which is taken to be 1. Depending on the case being investigated, we employ two types of boundary conditions: the periodic and symmetry boundary conditions. The energy minimization is then carried out using the computationally efficient L- BFGS algorithm,45,46 optimized for both memory usage and processing speed. During the minimization routine, the value of fCngijk of each node is updated in order to achieve the minimum total energy configuration.
The implementation of the diffuse solid method during the energy minimization process in the simulation is carried out in two stages, as illustrated in Fig. 1. In the first stage, energy minimization is performed to create a solid geometry with a smooth and diffuse solid interface.47 The solid is initialized as component 1 (C1) within the desired geometry. The remaining system is filled with pure gas (C3) and the liquid phase (C2) is fixed at 0 [Fig. 1(a)]. The total free energy is subsequently minimized using the L-BFGS algorithm. During this stage, the energy minimization process is truncated at a specific

FIG. 1. Schematic of implementing the diffuse solid method for a droplet on a flat surface in the energy minimization. (a) Initialization of the flat solid geometry (C1) and the gas phase (C3). (b) The diffuse solid interface with a thickness of a is formed after a number of iterations as highlighted by the inset which also shows the tanh profile of the solid component (C1) across the interface. (c) Initialization of the liquid phase (C2) using a staircase approximation. (d) The final result after the minimization process is completed.

@tCi þ $ ^ ðCivÞ ¼ $ ^ ðM$liÞ: (14)
We will assume the mobility M is a constant for simplicity, and the chemical potential li ¼ ddCW can be derived from the free energy of the system.
As with the free energy minimization scheme, there are two stages of computations when employing the diffuse interface method. The first stage is to create the smooth and diffuse solid interface. This can be done by employing an energy minimization scheme as discussed in Sec. II A. Alternatively, we can also use lattice Boltzmann approach to relax the C1 ^ C3 (solid–gas) interface to its tanh profile. To do this, we can write the free energy from Eq. (1) using the constraint C1 þ C3 ¼ 1,

W¼

ðj1 þ2 j3Þ C12ð1 ^ C1Þ2 þ Wcf ðrÞþ a2ðj12þ j3Þ þ XC12HðC1Þ þ Xð1 ^ C1Þ2Hð1 ^ C1Þ dV;

ð$C1Þ2

(15)

number of iterations (typically, niter ¼ 10) to allow the solid–gas diffuse interface to form without strongly perturbing the intended solid geometry’s shape [Fig. 1(b)]. Once this stage is complete, the liquid phase (C2) is initialized in the system according to a predefined configuration [Fig. 1(c)]. The solid phase (C1) is now fixed, while the liquid (C2) and gas (C3) phases are allowed to vary. The total free energy is now minimized under these constraints [Fig. 1(d)]. Convergence of the minimization process is determined based on an rms gradient condition: convergence is achieved when j$Wj < 10^5.
The determination of niter, a, b, and the domain dimensions are chosen to balance the accuracy of the solid geometry, accuracy of the wetting properties described in Sec. III, and computational time. A systematic investigation of the impact of each of these parameters is performed in the supplementary material, Sec. A.
B.		Lattice Boltzmann method
Next, we will describe how the diffuse solid method can be integrated with lattice Boltzmann method to model wetting dynamics involving solid boundaries of complex shapes. Generally, the continuum equations of motion solved by the lattice Boltzmann method comprise of three equations. First, the continuity equation
@tq þ $ ^ðqvÞ¼0;	(12) which describes the conservation of mass in a fluid with density q and velocity v.
Second, the incompressible Navier–Stokes equation
@tðqvÞþ$ ^ ðqvvÞ ¼ ^$ ^ P þ F þ gDv;	(13) which describes the conservation of momentum. Here, P denotes the pressure, g is the dynamic viscosity, and F ¼ Fs þ Fg is the body force applied to the fluid. Fs is the fluid–solid and fluid–fluid surface tension force and Fg is the external body force (e.g., gravitational force).
Third, we need another equation to enforce the phase separation of the two fluid components and evolve their interface shapes. The Cahn–Hilliard equation can be written in terms of the relevant order parameter Ci ,
where the constraining term, Wcs, is ignored because the lattice Boltzmann method conserves total volumes. The additional final terms using the Heaviside step function, H, and parametrized by X48 are used to prevent the order parameter deviating from the range C1 2ð0; 1Þ, which can improve simulation stability during the solid initialization when the value of b in Wcf ðrÞ is large. For this initialization stage, for simplicity, we set the density, viscosity, and mobility to be uniform everywhere with q ¼ 1, g ¼ q3 ðs ^ 0:5Þ ¼ 16 obtained by setting s ¼ 1, and M ¼ ðsg ^30:5Þ ¼ 16 obtained by setting sg ¼ 1. As is common in lattice Boltzmann method,49 s and sg are relaxation parameters that are related to the fluid viscosity and mobility parameter.
Given the free energy model in Eq. (15), the surface tension force and chemical potential are given by
Fs ¼ ^C1$l;	(16)
dW	2	dWcf ðrÞ
l ¼ ¼ðj þ j ÞC ð1 ^ 3C þ 2C Þþ
dC1	1	3 1	1	1	dC1
^ a2ðj1 þj3Þr2C1 þ2XC1HðC1Þ^2Xð1 ^ C1ÞHð1 ^ C1Þ;

dWcf ðrÞ
dC1

2bðC1 ^	0:5Þ	if	C1	>	0:5	and	r 2	S;
2bðC1 ^	0:5Þ	if	C1	<	0:5	and	r 62	S;
0 otherwise:

(17)
(18)

If the solid has sharp features or corners, to avoid numerical instabilities, we can carry out the diffuse solid initialization in two steps. First, we perform a short lattice Boltzmann run, typically with larger b and X in Eq. (15), and then redefine the solid locations as the positions where C1 > 0:5. Typically, this step is only necessary where the solid geometry has corners with a radius of curvature that is comparable to 4:4a. Following this, a further initialization run is performed, where the diffuse solid is evolved using the updated solid locations. After this point, the order parameter C1 is not evolved anymore.
The second stage of the simulations is to evolve the fluid order parameters (C2 and C3) for a fixed diffused solid (C1 order parameter distribution) using the continuity, Navier–Stokes and Cahn–Hilliard equations. Starting from Eq. (1), we can leave out the (now constant) terms in C1 and apply the constraint that C3 ¼ 1 ^ C1 ^ C2 to get

W¼	j22C22ð1	^	C2Þ2	þ j23	ð1 ^C1	^C2Þ2ðC1	þ	C2Þ2
þ a2 j22 ð$C2 Þ2 þ a2 j23 ð$C1 þ $C2Þ2idV:	(19)
The benefit of this approach is that the surface tension between fluids C2 and C3 and the contact angle of the C2 component with the diffuse solid are given from the following relations (see supplementary material, Sec. B):
aðj2 þ j3Þ	j3 ^ j2
c23 ¼	; cos h ¼	:	(20)
6	j2 þ j3
From the definition of the free energy in Eq. (19), we are able to derive the surface tension force for the Navier–Stokes equation and chemical potential for the Cahn–Hilliard equation, given by
Fs ¼ ^C2$l;	(21)
with
dW
l¼
dC2
¼ j2C2 ð1 ^ 3C2 þ 2C22Þ
þ j3ð1 þ C2ÞðC12 þ 2ð^1 þ C2ÞC2 þ C1 ð^1 þ 3C2ÞÞ
^ a2j2r2C2 ^ a2j3r2C1:	(22)
accounts for the advective terms in the Navier–Stokes equation. The form of the force Fi, following Refs. 48, 50, and 51, is given by
1	ci ^ F
Fi ¼ 1 ^ xi ðci ^ uÞ ^ $qsiðuÞ þ : 2sf	cs2
(28)
Although, in principle, incompressibility is not vital for this method, it has been shown to reduce spurious velocities around the diffuse interface50 and this form of f eq can aid with stability when a density ratio between fluid components is included. This ensures that our model can be compatible with problems where inertia is important. Directional gradients $ are given by a second order isotropic central difference scheme. Moments of the distribution function provide the hydrodynamic pressure and velocity
8	dt
P¼	fiþdtu^$q;	(29)
i¼0	2
8	Dt
qu ¼	cifi þ l$C:	(30)
i¼0	2
Now, we will give specific details for the lattice Boltzmann49 implementation in this study, which is used to discretise equations (12) and (13). The evolution of the macroscopic variables is expressed using distribution functions fi, which are stored over a discrete regular grid. Distributions can travel to adjacent grid points in a discrete number of velocity directions ci with weightings xi in each direction. In this study, we use the D2Q9 stencil with
01^10	0	1^1
c 0001^11^1
and
41111	1	1
x¼
9 9 9 9 9 36 36
1
^1
1
36
^1
1;
(23)
(24)
For simplicity, we choose the BGK collision operator, although extension to other collision models is possible. This choice leads to the following lattice Boltzmann equation:
fiðx þ cidt; t þ dt Þ ^ fiðx; tÞ¼^1s ^fiðx; tÞ ^ fieqðx; tÞ^
þ dtFiðx; tÞ:	(25)
It is customary in lattice Boltzmann to take dt ¼ 1.
Given that we solve the incompressible Navier–Stokes equation (13), fieq takes the following form:48,50
fieq ¼ xi^P2 þ qsiðuÞ^;
cs
where cs2 ¼ 13 is the sound speed and
(26)
siðuÞ ¼
u ^ c	ðu ^ ciÞ2	u ^ u
cs2iþ	2cs4i	þ 2cs2
(27)
A consequence of choosing fieq from Eq. (26) is that the fluid density no longer arises from the distribution functions fi. Instead, it is a function of the fluid concentration, which evolves through the Cahn– Hilliard equation (14). The lattice Boltzmann equation for the Cahn– Hilliard equation is given by
giðxþ cidt; t þ dtÞ ^ giðx; tÞ¼^1 ^giðx; tÞ ^ gieqðx; tÞ^; sg	i
with the equilibrium distribution gieq given by Ref. 52,
eq ¼ ð1 ^ x0Þl þ Cð1 þ x0si ðuÞÞ; if i ¼ 0;
gi	xil þ xiCsiðuÞ;	otherwise:
The zeroth moment of the distribution function gives the concentration
8
C ¼	gi :
(31)
(32) fluid
(33)
The mobility appearing in the Cahn–Hilliard equation is given by
M ¼ cs2ðsg ^ 0:5Þdt:	(34)
In our simulations, we relate the fluid density and viscosity as functions of the concentration as given by
C2
q¼q3þðq2^q3Þ	;	(35)
C2 þ C3
s ¼ s3 þðs2 ^ s3ÞHðC2 ^ C3Þ;	(36)
g¼q3ðs^0:5Þ:	(37)
Moreover, in our lattice Boltzmann scheme, to implement no-slip boundary condition at the solid boundary, we apply an interpolated bounce back boundary condition.53 The interface distance for the interpolation is approximated by assuming the fluid concentration follows a tanh profile normal to the interface:
11	x
C ¼ þ tanh :
n 2 2	2a
(38)

D¼

x
n^ ^ ^d

(39)

The distance to the interface along an arbitrary direction d is then given by where n^ ¼ $Cn=j$Cn j is the normal direction to the interface. Pseudocode to describe the full simulation procedure is provided in the supplementary material, Section C.
III.		PHASE FIELD ENERGY MINIMIZATION
A.		Liquid contact angle
We first benchmark the diffuse solid method by measuring the contact angle of liquid wetting on both flat and curved solid surfaces within two-dimensional (2D) and three-dimensional (3D) frameworks. For wetting on flat surfaces, we simulate a liquid droplet deposited on a chemically homogeneous solid surface. The simulation domain for the 2D and 3D cases are 300 ^ 100 and 300 ^ 300 ^ 80, respectively. A sharp-interface solid surface of thickness 10 is initialized at the base and allowed to evolve into a diffuse interface during the minimization process. A sessile liquid drop of radius 30 is then placed on this diffuse-interface solid, and the system is further minimized until equilibrium is achieved, yielding the final configuration shown in Figs. 2(a) and 2(b). Periodic boundary conditions are applied along the horizontal directions, while symmetry boundary conditions are enforced in the vertical direction [Figs. 2(a) and 2(b)].
For wetting on curved solid surfaces, we position a spherical solid particle of radius 30 at the liquid-gas interface within simulation domains of 100 ^ 100 (2D) and 60 ^ 60 ^ 100 (3D). The solid particle is initialized using a staircase approximation and undergoes relaxation steps to form a smooth surface. The liquid is then introduced into the system. In the 2D case, periodic boundary conditions are applied along the horizontal directions, while symmetry boundary conditions are enforced in the vertical direction to ensure a flat liquid–gas interface, as shown in Fig. 2(c). For the 3D case, computational efficiency is optimized by simulating only a quarter of the sphere, with symmetry boundary conditions imposed in all directions [Fig. 2(d)]. ffiffiffiffiffiffiffi
The parameters c23 ¼ 8=9, c12 ¼ 1, b ¼ 10, and niter ¼ 10 are used for these benchmark simulations. The contact angles are varied between 10^ and 150^ to show a broad applicability of our method. The measurement of the contact angle is done by fitting the liquid–gas interface line, for the droplet case, or the solid interface line, for the spherical particle case, to a circle, both for the 2D and 3D cases. Note that, as shown in Figs. 2(a) and 2(c), there is small, localized distortion at the three phase contact line. However, this deformation has negligible effect as the droplet or the particle size is always much larger and hence is ignored in the contact angle measurement.

FIG. 2. Equilibrium shapes of liquids wetting (a) and (b) flat and (c) and (d) curved solid surfaces in 2D and 3D cases, respectively. Images also show boundary conditions used in simulations for each case. (e) Comparison of the input (hanalytical) and measured (hmeasured) contact angles for wetting cases on flat and curved surfaces in 2D and 3D. (f), (g) Pinning investigation for a drop on an inclined plane using a binary phase field model with a staircase approximation to the surface and using the diffuse solid method (DSM), respectively. Images show a 40^ surface inclination angle with a drop with hanalytical ¼30^. (h) Comparison between the staircase approximation and DSM mean contact angle and range of contact angles exhibited as the volume of a drop on an inclined plane is gradually increased and reduced. Mean angles are indicated by a dashed black line (staircase) or solid black line (diffuse solid method). The range of contact angles is indicated by the gray-shaded (staircase) or blue-shaded (diffuse solid method) regions.

FIG. 3. (a) A snapshot of the liquid–gas interface at equilibrium in a capillary rise simulation within a cylindrical tube. The tube radius rc, the contact angle h, the rise height hc and the meniscus height hm are shown in the enlarged view. (b) The measured rise height hc (data points) plotted against contact angle h with rc ¼ 26 (main panel) and inner tube radii rc with h ¼ 60^ (inset), and compared to the numerical solution given by Eq. (43) (solid line). (c) The change in the gravitational potential energy DEg and the surface energy DEsurface obtained from simulation and numerical solution plotted against contact angle h, where rc ¼ 26.

The comparison between the input [hanalytical from Eq. (20)] and measured (hmeasured) contact angles is shown in Fig. 2(e). From this figure, it can be seen that the simulation results are in good agreement with the analytical prediction for a wide range of contact angles. We note that for low contact angles (h ^ 10^) in the spherical case, the three-phase contact line meets near the north pole of the sphere, and due to the diffuse interface, the accuracy of the contact angle measurement reduces, yielding a relatively large error. Otherwise, we find that the errors of the presented data are typically less than 2^ . Similar accuracy for the contact angle is also obtained for the lattice Boltzmann method, as illustrated in the supplementary material, Section D.
At this stage, it is instructive to compare the accuracy of the contact angles when using the diffuse solid method, with the widely used staircase approximation. To achieve this, we perform a contact angle hysteresis test, whereby the volume of a sessile droplet is incrementally increased and then reduced, with the contact angle measured at each step. This procedure is performed on a planar surface at varying inclination angles relative to the underlying mesh, from 0^ to 45^ (all other angles being equivalent to this range). Furthermore, this test is carried out at four representative analytical contact angles, hanalytical ¼ 30^, 60^,90^,and 120^.
In Figs. 2(f) and 2(g), an example is shown of a 2D drop with hanalytical ¼ 30^, on an surface of inclination angle 40^, for the staircase approximation and diffuse solid method, respectively. Magnifications of the surfaces show the roughness introduced by the staircase approximation, compared to the smooth surface when using Diffuse Solid Method. In Fig. 2(h), the results of the contact angle hysteresis test are shown. For the staircase approximation, the mean measured contact angle hmeasured, indicated by the black dashed line, deviates substantially from hanalytical by up to 25^. The gray area indicates the range of hmeasured observed, which has a maximum value of 20^ . Both of these effects are caused by the underlying surface roughness, by increasing the contact area of the drop compared to a flat plane, and introducing artificial pinning. In comparison, for the diffuse solid method, the mean hmeasured, indicated by the solid blue line, deviates by a maximum of 0:4^, with a maximum range (light blue area) of 2^ . Thus, the diffuse interface method is effective at creating smooth, off-lattice surface geometries.
B.		Capillary rise
We next consider the fundamental wetting phenomenon of capillary rise of a liquid in a smooth cylindrical tube.54,55 The equilibrium rise height of liquid within the tube can be understood from an energetic perspective. On changing the height of liquid in the tube, the total energy change of the system, DE, can be described as a sum of interfacial (DEsurface) and gravitational (DEg) contributions,
DE ¼ DEsurface þ DEg: (40)
DEsurface can be broken into a sum of individual interfacial energy contributions,
DEsurface ¼ DAlgclg þ DAlscls þ DAgscgs; (41) where DAlg, DAls, and DAgs are the area changes of the liquid–gas, liquid–solid, and gas–solid interfaces, respectively, and clg, cls, and cgs, are the surface tensions of the same respective interfaces. If we assume a dry capillary tube as an initial state, then we have DAgs ¼ ^DAls ¼^2prcðhc þ hmÞ, where, as shown in Fig. 3(a), rc is the tube radius,

h¼

2clg cos h

qgrc

(42)

(43)

Wg ¼	^fgC2z^dV;
V

(44)

FIG. 4. (a) Schematic of a droplet undergoing a transition from a suspended state to a collapsed state on a two-dimensional periodically patterned surface with a rectangular micropillar geometry having width a, pillar spacing s, pillar height h, and periodicity b. h is the equilibrium contact angle on the corresponding flat surface. Pl and Pg are the pressure in the liquid and gas (surrounding air), respectively. The enlarged view shows the liquid–gas interface underneath the droplet hanging between two adjacent pillars in a suspended state and is used as the system setup in simulations. (b) Measured critical pressure plotted against contact angle, where s ¼ 81, and pillar spacing (inset), where h ¼ 112^, and compared with the analytical prediction to Eq. (45). Snapshots show liquid–gas interface shapes at the critical pressure for h ¼ 60^ and h ¼ 112^ .

hc is the capillary rise height (here, defined from the bottom of the tube to the bottom of the meniscus) and the meniscus height hm . Famously, Jurin’s Law describes the equilibrium rise height in such a system, but it has only recently been shown that the precise rise height that Jurin’s law describes is the mean rise height h ,56
The rise heights hc and hm cannot in general be analytically determined, and likewise neither can the energetic contributions in Eq. (40). However, as shown previously, these quantities can be determined numerically by solving the Young-Laplace equation for this axisymmetric system,56
z00 þ	z0	¼ z
ð1 þ z02Þ32 xð1 þ z02 Þ12	kc2 ;
where kc is the capillary length (kc ¼ clg =qg). Here, as in Ref. 56, the shooting method is employed in MATLAB to solve the meniscus profile, rise heights, and energetic contributions.
To evaluate our simulations against this established numerical method, we begin by initializing a cylindrical tube with an inner radius of rc using a staircase approximation of a hollow cylinder with a sharp interface. We set b ¼ 10 in the confining potential and niter ¼ 10 steps to smooth out the solid interface. The gravitational force is applied as a body force per unit volume, Fb ¼ qg, in the vertical direction after the liquid phase is initialized in the system. This gravitational force enters the total free energy functional as an additional external energy term as follows:
where fg is the gravitational force density. We use water as the rising liquid with properties: q ¼ 1000 kg m^3 and clg ¼ 0:0728 N m^1. The gravitational acceleration is g ¼ 9:81 m s^2. All of these physical parameters are then converted into simulation units via the gravitational force density: fg ¼ qgðc23P02Þ=ðclgP2Þ, where P0 and P are the length scale of the capillary tube in physical and simulation units, respectively, and clg and c23 are the liquid–gas interfacial tensions in physical and simulation units, respectively. We use these values in simulations: c23 ¼ pffi8ffiffi=ffiffi9ffiffi, P ¼ 80, and P0 ¼ 10^2 m. To reduce the computational cost, we only simulate half of the cylinder and apply symmetry boundary conditions along all directions. The typical simulated capillary rise in a smooth cylinder tube is shown in Fig. 3(a).
Figure 3(b) presents the measured rise height of liquid in a smooth cylinder tube for varying wetting contact angles and tube radii (see the inset). The simulations accurately capture the dependence of the rise height on h and rc, demonstrating a strong agreement with the numerical solution of Eq. (43). Additionally, we analyzed the variations in interfacial energy, DEsurface, and gravitational potential energy, DEg, obtained from simulations for different contact angles, comparing them with the theoretical prediction of Eq. (40), as shown in Fig. 3(c). The change in interfacial energy is computed from the bulk and squared gradient terms within the free energy model in Eq. (1). The change in gravitational potential energy is derived from the gravitational force term in Eq. (44). The energy calculations show excellent agreement with theoretical predictions.
C.		Critical pressure on micropillared surfaces
We study the wetting transition on superhydrophobic surfaces, also known as the Cassie-Baxter to Wenzel transition. In this, a liquid droplet transitions from a state suspended on top of the surface texture to a collapsed state where the droplet fully wets the structure. This is associated with a reduction in the apparent contact angle of the droplet and a loss of superhydrophobicity. The wetting transition occurs if the pressure difference between the liquid and gas is larger than a critical pressure given by the surface textures, which depends upon the equilibrium contact angle. At higher equilibrium contact angles, it is harder for the liquid to enter the texture, resulting in a higher critical pressure.
We consider a 2D system where the solid is textured with rectangular micropillars, as shown in Fig. 4(a). We simulate a liquid layer suspended on top of two pillars separated by a gap s instead of simulating a whole droplet on a patterned surface [see enlarged view in Fig. 4(a)]. We make this approximation assuming that the droplet size is much larger than the pillar size and spacing. In our simulations,
these values are used: a ¼ 40, h ¼ 170, b ¼ 160, b ¼ 10, and niter ¼ 10. The system is fitted into a simulation box of 160 ^ 200. We incorporate pressure into the system by applying a pressure difference between the liquid and gas phase, DP, in the free energy functional in Eq. (8). In order to obtain an accurate critical pressure, we apply a binary search algorithm57,58 around the critical pressure value with a tolerance level of 10^4.
Theoretically, for the 2D model system used here, the critical 59–61
pressure is given by59
2clg cos h
;

DPc ¼^

(45)

s where clg is the liquid–gas surface tension; h is the equilibrium contact angle of the droplet on the corresponding flat surface, and s ¼ b ^ a is the wall-to-wall spacing between two adjacent pillars.
The equilibrium shapes of the liquid–gas interface at the critical pressure are shown in the snapshots of Fig. 4(b), in which concave and convex menisci occur due to the liquid–solid surface tension for h ¼ 60^ and h ¼ 112^, respectively. At the transition (DP > DPc), the suspended liquid–gas interface abruptly collapses, enabling the liquid phase to displace the trapped gas and completely fill the gap between the pillars. The critical pressure dependency on the equilibrium contact angle, h, and the pillar spacing, s, is shown Fig. 4(b). As predicted by Eq. (45), the critical pressure is linearly dependent on cos h. For h < 90^, the critical pressure is negative, indicating that there is no energy barrier and the transition to the collapsed state is energetically favorable. On the contrary, the positive critical pressures for h > 90^ suggest the system has an energy barrier, making the transition not energetically favorable to occur unless some external perturbations from the liquid pressure are added to the system. As the separation increases, the critical pressure reduces making it more favorable for the liquid to transition to a collapsed state. The simulated critical pressures are in good agreement with the prediction.
D.		The Salvinia Paradox
Having validated the diffuse solid method, we now demonstrate its utility in representing real structures of complex geometry through a case study of the so-called “Salvinia Paradox.” First investigated by Barthlott et al.,62 the Salvina Paradox concerns the wetting properties of the hair-like structures that coat the leaves of the Salvina water fern. In the species Salvina molesta, these mm-long hairs have a complex, “eggbeater” geometry, which is modeled in Fig. 5(a). The habitat of Salvina molesta is on the surface of rivers and lakes, and plants are susceptible to being submerged. However, the plants also require an air layer close to the leaf (plastron) to maintain gaseous exchange. The eggbeater hairs are an adaptation that allow the plastron to be maintained when submerged. The four arms are coated in a hydrophobic wax and through their physical structure maintain a superhydrophobic suspended state when submerged. This resists the intrusion of water into the plastron. The paradox arises because the uppermost tip of the hair is hydrophilic, a feature which conventionally is assumed to be detrimental to superhydrophobicity. However, it has been shown experimentally that whereas the hydrophobic structured base resists the intrusion of water, the hydrophilic tip resists the exit of air as bubbles from the plastron. These two different wetting behaviors are suggested to work synergistically to maintain the plastron, even in flowing water.62
To measure the effectiveness of these structures at resisting breakdown of the plastron, it is useful to consider the extrusion and intrusion pressures, the pressures at which the liquid–gas interface either detaches upwards from the tip of the structure, or moves spontaneously down the structure toward the base, respectively. To date, these properties have largely been investigated experimentally on plant samples, with reported extrusion pressures of (^6.1 6 0.9) kPa63 and intrusion pressures of (12 6 2) kPa.64 However, there remains a lack of fundamental understanding about the relationship between the structure, surface wettability, and the extrusion/intrusion pressures. Toward achieving this understanding, some progress has been made via computational studies. The extrusion and intrusion characteristics on a 2D simplified geometry showed how the interaction between wettability patterning and structure enhanced the difference between extrusion and intrusion pressures.65 Meanwhile more recently, the apparent contact angle of the suspended state was investigated, alongside some aspects of the intrusion characteristics using the proprietary software Ansys Fluent;66 however, wettability patterning was not considered.67
Using the diffuse solid method, the extrusion and intrusion of water on the complex geometry with patterned wettability is investigated. To begin with, a model of the geometry was designed in Blender68 (although any suitable modeling or computer-aided design (CAD) software is appropriate). The triangulated mesh of this model is shown in Fig. 5(a). Within the computational domain of the phase field model, the face and vertex data of the triangulate mesh was then used to distinguish the solid interior and exterior, and so set up the confining potential Wcf described in Eq. (18). Following the solid relaxation step, the resultant solid geometry is shown in Fig. 5(b), in which the hydrophobic base of contact angle hbase is shown in gray, and the hydrophilic tip of contact angle htip is highlighted in red.
When choosing the dimensions of the model hair and domain size, it is recognized that real plant specimens show substantial variability in the hair geometry and arrangement.69 Key geometric parameters commonly examined are the hair height Hs, outermost diameter Ds, and center-to-center distance Cs, indicated in Fig. 5(b). Here, for this demonstration of the diffuse solid method, we opted for a typical hair geometry characterized by Ds=Cs ¼ 0.77, and Hs=Cs ¼ 0.99. Periodic boundary conditions were imposed in the x and y directions, with symmetry boundary conditions in the z directions.
To examine the extrusion pressure, DPExt, a liquid layer was initialized in contact with the tip of the hair, and the energy of the system minimized. The pressure in the liquid phase was reduced relative to the gas phase and the system energy re-minimized. This step was repeated with variable steps in the pressure changes to locate the pressure at which the interface detached from the top of the structure: DPExt. To explore the impact of the base and tip wettability on DPExt, this process was repeated across a range of hbase and htip from 20^ to 160^, capturing the range of biologically relevant contact angles.
The variation of DPExt with hbase and htip is shown in Fig. 5(c), with the color scale located at the bottom of the figure. Note that DPE0 xt is reported, which is a nondemensionalised form of DPExt: DPE0 xt ¼ DPExt=ðclg=LxÞ, where clg is the liquid-gas interfacial tension (clg ¼ 1:0) and Lx is the number of cells in the computational domain in the x-direction (Lx ¼ 120). Two different modes of interface detachment are observed, labeled base-mediated and tip-mediated, which are separated by the thick, blue line. In base-mediated detachment, at the point of detachment, the interface is pinned to the solid

FIG. 5. (a) The triangulated 3D mesh that models the “eggbeater” shaped leaf hair of the water-fern, Salvinia molesta. (b) The diffuse solid conforming to the inputted triangulated mesh, showing the grid dimensions and different contact angles imposed on the base, hbase, and the tip, htip.Key hair parameters are indicated: the hair height Hs, outer- 0
most diameter Ds, and center-to-center distance Cs. (c)–(e) Surveys of the nondimensionalised extrusion pressure DPExt, intrusion pressure DPInt, and difference between these pressures, respectively. Boundaries between the two collapse modes (tip-mediated and base-mediated) are indicated by a thick blue line (extrusion) and thick red line (intrusion). (f) A 2D slice through the Salvinia hair, showing how the liquid–gas interface changes across the pressure range between the intrusion and extrusion pressures. Interfaces are shown in blue at equal increments of DP0 ¼ 0.12. The hydrophilic tip is shown in red and the hydrophobic base in gray.

base, and so htip does not effect DPExt. Instead, the smaller the value of hbase, the more negative DPExt becomes. Conversely, in tip-mediated detachment, the interface is pinned to the tip at the point of detachment and so DPExt becomes independent of hbase. The smaller the value of htip, the more negative DPExt becomes. However, the small perimeter of the three-phase contact line in tip-mediated detachment limits the maximum magnitude of DPExt. This is compared to basemediated detachment where the relatively large three-phase contact line perimeter enables larger magnitude extrusion pressures to be achieved. Thus, if preventing extrusion was the sole purpose of the Salvinia hairs, then this would most effectively be achieved by minimizing hbase.
To measure the intrusion pressures, DPInt, the pressure in the liquid relative to the gas was iteratively increased up until the point where the liquid spontaneously filled the system. We note as an aside that depending on hbase, the interior of the eggbeater structures either trap an air bubble, or fill with liquid. The variation of DPI0nt with hbase and htip is shown in Fig. 5(d), where DP0Int has been nondimensionalised in an identical manner to DPE0 xt. Here, we also see two separate detachment mechanisms, base-mediated and tip-mediated. In base-mediated detachment, the three-phase contact line pins close to the widest part of the hair geometry at the point of detachment. The eggbeater structure traps air within the space between the four arms, making it energetically expensive for the liquid-gas interface to move down the outside of the hair. Thus, the intrusion pressure is positive, even for hydrophilic values of hbase as low as 50^. In scenarios where the base is hydrophilic and the tip highly hydrophobic (the reverse of that seen in the Salvinia plant), the tip becomes responsible for resisting intrusion in tip- mediated detachment. However, due to the small contact line perimeter of the tip, only modest positive intrusion pressures can be achieved. If preventing intrusion was the sole purpose of the Salvinia hairs, then this would most effectively be achieved by maximizing hbase.
To explore the proposed function of the wettability-patterened Salvinia hair (to maintain the largest range between intrusion and extrusion pressures62) we show the difference, DPInt-DPExt, in Fig. 5(e). This shows that, in order to maximize DPInt-DPExt, hbase should be maximized, and htip minimized. Although not yet explicitly studied for salvinia, the nano-structured wax coating of the hydrophobic base is anticipated to show similar contact angles to other plants, approximately 160^.62,70 Likewise, although not yet robustly

pubs.aip.org/aip/pof

c Hcosh
L^L ¼ lg d
0	3gl

FIG. 6. (a) Visualization of simulation domain, with a reservoir of liquid on the left, which will enter the solid structure. The channel is inclined at 20^ . The boundaries on the left and right are periodic with symmetry conditions on the top and bottom.49 (b) Simulation data points and analytical curves (shown by solid lines), for non- dimensional distance and time.

confirmed, the hydrophilic tip is suggested to have a contact angle of approximately 35^ .63 At these contact angles, our simulations show that DPE0 xt ¼ ^0.81, DPI0nt ¼ 2.94, with a difference of 3.75. Choosing a typical value of the center-to-center hair spacing of 790 lm,69 these convert to DPExt ¼ ^8.85 kPa and DPInt ¼ 32.1 kPa, with a difference of 41.0 kPa. This compares sensibly to the experimental values of DPExt ¼ (^6.1 6 0.9) kPa63 and DPInt ¼ (12 6 2) kPa.64 Note that the magnitudes of these pressures are reduced compared to the simulated values. This is because in real plant samples, substantial geometric variability between hairs means that measured detachment pressures will be dominated by the local regions which are least resistant to intrusion or extrusion. In the supplementary material (Sec. E), a sensitivity analysis is performed to show that the intrusion and extrusion critical pressures are highly sensitive to Ds =Cs, and only moderately sensitive to hbase and htip. Finally, an interesting further analysis is that if the hair exhibited a uniform contact angle, the maximum extrusion–intrusion pressure difference occurs at hbase ¼ htip ¼ 160^, with a value of 32.6 kPa. Thus, for the geometry tested here, wettability patterning increased the plastron stability pressure range by 26%.
To further illuminate the Salvinia hair function, in Fig. 5(f), a 2D slice through the 3D hair tip is shown. Between the intrusion and extrusion pressures, the liquid-gas interface is shown at uniform intervals of DP0 ¼ 0.12. At zero pressure, the interface is pinned on the hydrophobic/hydrophilic boundary. If the pressure is made increasingly negative (extrusion conditions), it can be seen that the three- phase contact line remains pinned to this location up until PInt is reached. If instead the pressure is increased gradually from 0 (intrusion conditions), at DP0 ¼ 0.48, the three-phase contact line depins and slides partially down the hair. As the pressure is increased further, the contact line continuously slides down the hair until DP0 ¼ 1.80, at which point another jump in contact line location is observed. From here, the interface slides toward the widest part of the hair, up until PExt is reached.
IV.		LATTICE BOLTZMANN SIMULATIONS
A.		Capillary filling dynamics
In order to benchmark the accuracy of the diffuse solid for dynamic problems, we investigate the filling of a capillary channel.71 The setup involves a finite reservoir of liquid entering a narrow channel, which the liquid preferentially wets over the gas. The curved liquid–gas interface in the channel creates a pressure gradient between the liquid and the surrounding gas, which subsequently drives the liquid to fill the channel. In the case where the liquid is much more viscous than the gas, the distance traveled by the liquid into the channel follows a square-root dependence on time72
1
2
ðt ^ t0Þ2: (46)
Here, gl is the dynamic viscosity of the liquid, hd is the dynamic contact angle, clg is the liquid-gas interfacial tension, H is the channel width, and t is time. L0 and t0 are constants of integration that take into account the initial dynamics of the liquid as it enters the channel.
Figure 6 shows how the distance traveled by the liquid into the channel depends on time. The simulation parameters are a surface tension clg ¼ 2:5 ^ 10^3, an interface width a ¼ 1 and viscosities of 1=6 and 2=300 for the liquid and gas phases, respectively. The densities of
each component are kept identically at 1. The value of mobility M in the Cahn–Hilliard equation is 1=6. The solid is initialized by equilibrating for 5 ^ 103 iterations using b ¼ 1 ^ 10^3 and X ¼ 2 ^ 10^4.
The full simulation size is 1750 ^ 450 and the height and width of the channel are 900 and 40 respectively. We incline the channel by 20^ in order to validate the dynamics in the case that the diffuse solid does not align to the grid.
These simulations are run for equilibrium contact angles of 60^ and 75^. However, as expected from previous literature,71,73 the dynamic contact angles are larger than these values. The dynamic contact angle is measured by fitting a circle to points on the liquid-gas interface and finding the tangential angle at the point where this intercepts the solid plane. As the liquid front moves along the channel and the velocity drops, the measured contact angles approach their equilibrium values, see the supplementary material (Sec. F), providing another validation that we achieve the correct equilibrium wetting behavior from Eq. (19).
In Fig. 6(b), the time and distance are non-dimensionalised by choosing
L^ ¼ L=H; t^ ¼ tclg=glH : (47)
We fit the data to Eq. (46) using a least squares algorithm, using the instantaneous dynamic contact angle at each data point when fitting the curve. As shown by Fig. 6(b), the distance traveled by the liquid interface follows a square root dependence on time. The fit agrees well with the expected curve given by Eq. (46). These results provide
confidence that the capillary dynamics are captured accurately. The curved interface provides the correct pressure gradient to drive the capillary flow with the expected velocity for a large viscosity contrast.
B.		Droplet dynamics on sinusoidal surfaces
To highlight a complex dynamic wetting application of the diffuse solid method, we now demonstrate simulations of droplets moving on sinusoidally corrugated surfaces. The solid initialization approach enables the generation of arbitrary solid shapes. This capability is vital to simulate multicomponent flow on real geometries, which very rarely conform to a regular simulation grid. Sinusoidal surfaces are a key example for this, and have found growing interest as superhydrophobic surfaces due to their unique dynamic behavior. In contrast to traditional pillared surfaces, the smooth texture limits pinning,74 so droplets can move more easily. Additionally, these wavy surfaces more closely match biological superhydrophobic surfaces, which have shown remarkable capabilities to repel liquid droplets.75,76
We investigate a droplet moving on these surfaces in the cases where the droplet is in the suspended and collapsed states, shown in Figs. 7(a) and 7(b). We initialize the solid component following a sinusoidal profile with an amplitude of 20 and a period of 50. We then allow the interface to smooth over 2500 timesteps with b ¼ 2 ^ 10^3 and X ¼ 0:1, followed by a further 47 500 timesteps with b ¼ 2 ^ 10^3 and

FIG. 7. (a) Droplet in the suspended state at a tilt angle of 45^, raised above the corrugated texture as it moves. The equilibrium contact angle is 150^ . (b) Droplet in the collapsed state at a tilt angle of 45^ that has penetrated the texture. The equilibrium contact angle is 90^ . (c) Terminal capillary number as a function of bond num- ber/tilt angle. Insets show that the velocity is close to constant in the suspended state and exhibits stick–slip motion in the collapsed state. The droplet in the suspended state moves at a significantly smaller tilt angle of 15^ compared with 37:5^ for the collapsed state. In addition to this, the terminal capillary number is consistently larger by approximately an order of magnitude.
X ¼ 2 ^ 10^3. To achieve the suspended and collapsed states, we set the equilibrium contact angles to 150^ and 90^ respectively. The other simulation parameters are identical to the values in Sec. IVA. The full simulation size is 250 ^ 130 and each simulation was run for 2 ^ 106 timesteps. Periodic boundary conditions are used at the left and right of the domain, and the halfway bounce-back scheme49 is used at the top. At the beginning of the simulation, a circular droplet with a radius of 40 is initialized in the collapsed state. To apply gravity, we include a constant body force to the droplet (we neglect gravity on the surrounding component with the assumption that the density ratio is large). The magnitude of the body force is calculated to match the Bond number

Bo ¼

ðqwater ^ qair ÞD g
;
cwater;air

(48)

of a water droplet with a diameter D ¼ 3:37 mm. Here, g^ is the gravitational acceleration parallel to the surface.
The terminal velocity V is calculated as the average center of mass velocity parallel to the solid surface for the final 1 ^ 106 timesteps, after the variation in average velocity of the droplet over successive periods of the sinusoidal texture was below 1%.
The results of these simulations is presented in Fig. 7(c), showing the terminal capillary number
Ca ¼ gl V : (49) clg
These results demonstrate the abilities of sinusoidal surfaces to repel liquid droplets. If the underlying substrate is hydrophobic, the droplet will enter the suspended state, and easily roll over the top of the corrugated texture with a small tilting angle of 15^. Even while the surface is neutrally wetting (h ¼ 90^), the collapsed state still begins to move at a moderate tilt angle of 37:5^, due to the lack of sharp corners on the sinusoidal surface, unlike other superhydrophobic surfaces patterned with posts.
Insets in Fig. 7(c) display the typical velocity variation with time for the suspended and collapsed states. For the suspended state, the velocity stabilizes quickly, with small oscillations as the droplet passes the sinusoidal profile. In contrast, for the collapsed state, the droplet must move vertically to overcome the peaks of the sinusoidal profile, which leads to a much smaller velocity with large oscillations as the droplet moves from one sinusoidal period to the next. The alternating pattern in the peak velocity of these oscillations arises from the different curvatures of the front and rear of the droplet in Fig. 7(b). Consequently, this leads to a slightly larger velocity as the front of the droplet overcomes the texture, in contrast with the rear. From these results, we also observe that the diffuse solid approach allows these surfaces to be fully resolved without a staircase-like approximation. The latter would create artificial pinning locations. Moreover, the solid initialization method allows for the correct interface profile to be achieved when analytical smoothing is impossible, which is essential for wetting on arbitrary geometries.
V.		CONCLUSION
In this work, we have presented a diffuse solid method that allows us to study quasi-static and dynamic wetting and multiphase flow phenomena in complex geometries. The diffuse solid approach is simple to implement and can be embedded in different numerical solvers, as shown here for phase field energy minimization and lattice Boltzmann method. Moreover, we have obtained excellent agreement against analytical results for all the benchmark scenarios; equilibrium contact angles, capillary rise height, critical pressure on pillared surfaces, and dynamic capillary rise. In turn, the diffuse solid model allowed us to easily study phenomena involving highly non-trivial geometries, highlighted here using the Salvinia leaf structure with its eggbeater shape and variation in contact angle at the base and the tip as well as droplet motion on sinusoidal surfaces.
We believe the proposed diffuse solid method will open up new opportunities to investigate wide ranging interfacial phenomena. For instance, it is useful for rationalizing wetting dynamics on bio-inspired surfaces, ranging from Lotus leaves and Rose petals to springtails and butterfly wings,77–80 which have been shown to exhibit a rich display of variations in shapes, sizes, and chemical variations. It can also be harnessed as a powerful design tool, especially when coupled with CAD software (as has been done here for the Salvinia leaf study), and potentially 3D printing to experimentally realize the computationally optimized surface structures.
SUPPLEMENTARY MATERIAL
See the supplementary material for additional information on (A) a parameter sensitivity study to obtain specific solid geometries, (B) the derivation of interfacial tensions and contact angle in the diffuse solid model, (C) pseudocode for the lattice Boltzmann algorithm employed in this work, (D) equilibrium contact angle validation for the lattice Boltzmann method simulations, (E) parameter sensitivity survey for the Salvinia paradox example, and (F) dynamic contact angle data for the lattice Boltzmann simulations shown in Sec. IV B.
ACKNOWLEDGMENTS
F.O. acknowledges a BPPLN scholarship from the Directorate General for Science Technology and Higher Education, Republic of Indonesia. H.K. acknowledges funding from the Leverhulme Trust (Research Project Grant No. RPG-2022-140) and UKRI Engineering and Physical Sciences Research Council (No. EP/ V034154/2). H.K. and S.J.A. acknowledge funding from EPSRC Impact Acceleration Account and ExxonMobil.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Fandi Oktasendra: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Michael Rennick: Investigation (equal); Methodology (equal); Writing – original draft (equal). Samuel J. Avis: Methodology (equal); Writing – review & editing (equal). Jack R. Panter: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). Halim Kusumaatmaja: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.
REFERENCES
1W. Barthlott, M. Mail, B. Bhushan, and K. Koch, “Plant surfaces: Structures and functions for biomimetic innovations,” Nano Micro Lett. 9, 1–40 (2017).
2P. Comanns, “Passive water collection with the integument: Mechanisms and their biomimetic potential,” J. Exp. Biol. 221, jeb153130 (2018).
3E. Arzt, H. Quan, R. M. McMeeking, and R. Hensel, “Functional surface microstructures inspired by nature – from adhesion and wetting principles to sustainable new devices,” Prog. Matter. Sci. 120, 100823 (2021).
4Z. Yu, T. Zhu, J. Zhang, M. Ge, S. Fu, and Y. Lai, “Fog harvesting devices inspired from single to multiple creatures: Current progress and future perspective,” Adv. Funct. Mater. 32, 2200359 (2022).
5F. Chen, Y. Wang, Y. Tian, D. Zhang, J. Song, C. R. Crick, C. J. Carmalt, I. P. Parkin, and Y. Lu, “Robust and durable liquid-repellent surfaces,” Chem. Soc. Rev. 51, 8476–8583 (2022).
6L. Hou, X. Liu, X. Ge, R. Hu, Z. Cui, N. Wang, and Y. Zhao, “Designing of anisotropic gradient surfaces for directional liquid transport: Fundamentals, construction, and applications,” The Innovation 4, 100508 (2023).
7P. Vineeth, A. Peethan, and S. D. George, “Special wettability for sensing: Drawing inspiration from nature,” Chem. Eng. J. 459, 141615 (2023).
8C. Yan, P. Jiang, X. Jia, and X. Wang, “3D printing of bioinspired textured surfaces with superamphiphobicity,” Nanoscale 12, 2924–2938 (2020).
9Y. Gao, Z. Peng, T. Shi, X. Tan, D. Zhang, Q. Huang, C. Zou, and G. Liao, “Bio-inspired fabrication of complex hierarchical structure in silicon,” J. Nanosci. Nanotechnol. 15, 5918–5923 (2015).
10K. Brassat and J. K. N. Lindner, “Nanoscale block copolymer self-assembly and microscale polymer film dewetting: Progress in understanding the role of interfacial energies in the formation of hierarchical nanostructures,” Adv. Mater. Interfaces 7, 1901565 (2020).
11E. Berthier, A. M. Dostie, U. N. Lee, J. Berthier, and A. B. Theberge, “Open microfluidic capillary systems,” Anal. Chem. 91, 8739–8750 (2019).
12K. S. Kwok, Q. Huang, M. Mastrangeli, and D. H. Gracias, “Self-folding using capillary forces,” Adv. Mater. Interfaces 7, 1901677 (2020).
13J. Liu, K. Duan, Q. Zhang, Y. Zheng, H. Cao, and Y. Zhang, “Pore-scale insights into CO2-water two-phase flow and implications for benefits of geological carbon storage,” Adv. Water Resour. 191, 104780 (2024).
14N. T. Chamakos, D. G. Sema, and A. G. Papathanasiou, “Progress in modeling wetting phenomena on structured substrates,” Arch. Comput. Methods Eng. 28, 1647–1666 (2021).
15P. Yue, C. Zhou, and J. J. Feng, “Sharp-interface limit of the Cahn–Hilliard model for moving contact lines,” J. Fluid Mech. 645, 279–294 (2010).
16D. Jacqmin, “Calculation of two-phase Navier–Stokes flows using phase-field modeling,” J. Comput. Phys. 155, 96–127 (1999).
17D. Anderson, G. McFadden, and A. Wheeler, “Diffuse-interface methods in fluid mechanics,” Annu. Rev. Fluid Mech. 30, 139–165 (1998).
18H.-G. Lee, J. Lowengrub, and J. Goodman, “Modeling pinchoff and reconnection in a Hele-Shaw cell. I. The models and their calibration,” Phys. Fluids 14, 492–513 (2002a).
19H.-G. Lee, J. Lowengrub, and J. Goodman, “Modeling pinchoff and reconnection in a Hele-Shaw cell. II. Analysis and simulation in the nonlinear regime,” Phys. Fluids 14, 514–545 (2002b).
20V. E. Badalassi, H. D. Ceniceros, and S. Banerjee, “Computation of multiphase systems with phase field models,” J. Comput. Phys. 190, 371–397 (2003).
21Q. Li, P. Zhou, and H. Yan, “Pinning–depinning mechanism of the contact line during evaporation on chemically patterned surfaces: A lattice Boltzmann study,” Langmuir 32, 9389–9396 (2016).
22R. Ledesma-Aguilar, D. Vella, and J. M. Yeomans, “Lattice-Boltzmann simulations of droplet evaporation,” Soft Matter 10, 8267–8275 (2014).
23K. A. Raman, R. K. Jaiman, T.-S. Lee, and H.-T. Low, “Lattice Boltzmann simulations of droplet impact onto surfaces with varying wettabilities,” Int. J. Heat Mass Transfer 95, 336–354 (2016).
24B. Dong, F. Wang, X. Zhang, and X. Jiang, “3D lattice Boltzmann simulation of droplet evaporation on patterned surfaces: Study of pinning–depinning mechanism,” Int. J. Multiphase Flow 125, 103218 (2020).
25S. Gong and P. Cheng, “Numerical simulation of pool boiling heat transfer on smooth surfaces with mixed wettability by lattice Boltzmann method,” Int. J. Heat Mass Transfer 80, 206–216 (2015).
26X. Shan and H. Chen, “Lattice Boltzmann model for simulating flows with multiple phases and components,” Phys. Rev. E 47, 1815–1819 (1993).
27Y. Ye, K. Li, Y. Wang, and T. Deng, “Parallel computation of entropic lattice Boltzmann method on hybrid CPU-GPU accelerated system,” Comput. Fluids 110, 114–121 (2015).
28J. T€olke, M. Krafczyk, M. Schulz, and E. Rank, “Lattice Boltzmann simulations of binary fluid flow through porous media,” Philos. Trans. R. Soc. London, Ser. A 360, 535–545 (2002).
29B. Stahl, B. Chopard, and J. Latt, “Measurements of wall shear stress with the lattice Boltzmann method and staircase approximation of boundaries,” Comput. Fluids 39, 1625–1633 (2010).
30J. Lu, H. Lei, and C. Dai, “A simple difference method for lattice Boltzmann algorithm to simulate conjugate heat transfer,” Int. J. Heat Mass Transfer 114, 268–276 (2017).
31Y. Yu, Q. Li, Z. Wen, and R. Huang, “Investigation on boundary schemes in lattice Boltzmann simulations of boiling heat transfer involving curved surfaces,” Phys. Fluids 32, 063305 (2020).
32M. P. Lautenschlaeger, J. Weinmiller, B. Kellers, T. Danner, and A. Latz, “Homogenized lattice Boltzmann model for simulating multi-phase flows in heterogeneous porous media,” Adv. Water Resour. 170, 104320 (2022).
33Y. Yao, Y. Liu, X. Zhong, and B. Wen, “Multiphase curved boundary condition in lattice Boltzmann method,” Phys. Rev. E 106, 015307 (2022).
34C. S. Li, R. R. Cai, and L. Z. Zhang, “An improved immersed moving boundary-based curved wetting condition method for the pseudopotential multiphase lattice Boltzmann model,” Phys. Fluids 36, 112119 (2024).
35H. Tanaka and T. Araki, “Simulation method of colloidal suspensions with hydrodynamic interactions: Fluid particle dynamics,” Phys. Rev. Lett. 85, 1338– 1341 (2000).
36A. Furukawa, M. Tateno, and H. Tanaka, “Physical foundation of the fluid particle dynamics method for colloid dynamics simulation,” Soft Matter 14, 3738– 3747 (2018).
37S. Aland, J. Lowengrub, and A. Voigt, “Two-phase flow in complex geometries: A diffuse domain approach,” Comput. Model. Eng. Sci. 57, 77–106 (2010).
38J. Yang, Y. Li, and J. Kim, “Modified multi-phase diffuse-interface model for compound droplets in contact with solid,” J. Comput. Phys. 491, 112345 (2023).
39F. Zhu, K. Sun, G. Zhang, and J. Yang, “Diffuse-interface modeling and energystable numerical framework for the heat transfer-coupled two-phase fluids in contact with solids,” J. Comput. Phys. 524, 113699 (2025).
40C. Semprebon, T. Kr€uger, and H. Kusumaatmaja, “Ternary free-energy lattice Boltzmann model with tunable surface tensions and contact angles,” Phys. Rev. E 93, 1–11 (2016).
41F. Oktasendra, A. Jusufi, A. R. Konicek, M. S. Yeganeh, J. R. Panter, and H. Kusumaatmaja, “Phase field simulation of liquid filling on grooved surfaces for complete, partial, and pseudo-partial wetting cases,” J. Chem. Phys. 158, 204501 (2023).
42J. R. Panter and H. Kusumaatmaja, “The impact of surface geometry, cavitation, and condensation on wetting transitions: Posts and reentrant structures,” J. Phys. Condens. Matter 29, 084001 (2017).
43J. Panter, Y. Gizaw, and H. Kusumaatmaja, “Multifaceted design optimization for superomniphobic surfaces,” Sci. Adv. 5, eaav7328 (2019).
44M. S. Sadullah, J. R. Panter, and H. Kusumaatmaja, “Factors controlling the pinning force of liquid droplets on liquid infused surfaces,” Soft Matter 16, 8114–8121 (2020).
45J. Nocedal, “Updating quasi Newton matrices with limited storage,” Math. Comp. 35, 773–958 (1980).
46D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Program. 45, 503–528 (1989).
47J. R. Panter, A. R. Konicek, M. A. King, A. Jusufi, M. S. Yeganeh, and H. Kusumaatmaja, “Rough capillary rise,” Commun. Phys. 6, 44 (2023).
48T. Lee and L. Liu, “Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces,” J. Comput. Phys. 229, 8045–8063 (2010).
49T. Krueger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. Viggen, The Lattice Boltzmann Method: Principles and Practice, Graduate Texts in Physics (Springer, 2016).
50T. Lee, “Effects of incompressibility on the elimination of parasitic currents in the lattice Boltzmann equation method for binary fluids,” Comput. Math. Appl. 58, 987–994 (2009).
51H. Liang, J. Xu, J. Chen, H. Wang, Z. Chai, and B. Shi, “Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows,” Phys. Rev. E 97, 033309 (2018).
52L. Zheng, S. Zheng, and Q. Zhai, “Lattice Boltzmann equation method for the Cahn-Hilliard equation,” Phys. Rev. E 91, 013309 (2015).
53M. Bouzidi, M. Firdaouss, and P. Lallemand, “Momentum transfer of a Boltzmann-lattice fluid with boundaries,” Phys. Fluids 13, 3452–3459 (2001).
54T. Young, “III. An essay on the cohesion of fluids,” Philos. Trans. R. Soc. 95, 65–87 (1805).
55P.-G. De Gennes, “Wetting: Statics and dynamics,” Rev. Mod. Phys. 57, 827 (1985).
56S. Liu, S. Li, and J. Liu, “Jurin’s law revisited: Exact meniscus shape and column height,” Eur. Phys. J. E41, 46 (2018).
57L. F. Williams, Jr., “A modification to the half-interval search (binary search) method,” in Proceedings of the 14th Annual Southeast Regional Conference (ACM, New York, NY, 1976), pp. 95–101.
58D. E. Knuth, The Art of Computer Programming: Fundamental Algorithms, Volume 1 (Addison-Wesley Professional, 1997).
59S. Mitra, N. S. K. Gunda, and S. K. Mitra, “Wetting characteristics of underwater micro-patterned surfaces,” RSC Adv. 7, 9064–9072 (2017).
60N. A. Patankar, “Consolidation of hydrophobic transition criteria by using an approximate energy minimization approach,” Langmuir 26, 8941–8945 (2010).
61A. Giacomello, M. Chinappi, S. Meloni, and C. M. Casciola, “Metastable wetting on superhydrophobic surfaces: Continuum and atomistic views of the Cassie- Baxter–Wenzel transition,” Phys. Rev. Lett. 109, 226102 (2012).
62W. Barthlott, T. Schimmel, S. Wiersch, K. Koch, M. Brede, M. Barczewski, S. Walheim, A. Weis, A. Kaltenmaier, A. Leder, and H. F. Bohn, “The Salvinia Paradox: Superhydrophobic surfaces with hydrophilic pins for air retention under water,” Adv. Mater. 22, 2325–2328 (2010).
63D. Gandyra, S. Walheim, S. Gorb, P. Ditsche, W. Barthlott, and T. Schimmel, “Air retention under water by the floating fern Salvinia: The crucial role of a trapped air layer as a pneumatic spring,” Small 16, 2003425 (2020).
64M. J. Mayser and W. Barthlott, “Layers of air in the water beneath the floating fern Salvinia are exposed to fluctuations in pressure,” Integr. Comp. Biol. 54, 1001–1007 (2014).
65M. Amabili, A. Giacomello, S. Meloni, and C. M. Casciola, “Unraveling the Salvinia Paradox: Design principles for submerged superhydrophobicity,” Adv. Mater. Interfaces 2, 1500248 (2015).
66A. Inc., see http://www.ansys.com/products/fluids/ansys-fluent for “Ansys fluent” (2024).
67S. Y. Zhang and J. G. Guo, “Superhydrophobicity and compressive stability of Salvinia surface: Mechanical analysis and structural design,” Appl. Surf. Sci. 627, 157274 (2023).
68B. O. Community, Blender - a 3D Modelling and Rendering Package (Blender Foundation, Stichting Blender Foundation, Amsterdam, 2024).
69P. Ditsche, E. Gorb, M. Mayser, S. Gorb, T. Schimmel, and W. Barthlott, “Elasticity of the hair cover in air-retaining Salvinia surfaces,” Appl. Phys. A 121, 505–511 (2015).
70W. Barthlott and C. Neinhuis, “Purity of the sacred lotus, or escape from contamination in biological surfaces,” Planta 202, 1–8 (1997).
71C. Pooley, H. Kusumaatmaja, and J. Yeomans, “Modelling capillary filling dynamics using lattice Boltzmann simulations,” Eur. Phys. J. Spec. Top. 171, 63 (2009).
72E. W. Washburn, “The dynamics of capillary flow,” Phys. Rev. 17, 273–283 (1921).
73N. Bala, M. Pepona, I. Karlin, H. Kusumaatmaja, and C. Semprebon, “Wetting boundaries for a ternary high-density-ratio lattice Boltzmann method,” Phys. Rev. E 100, 013308 (2019).
74G. Wells, E. Ruiz Gutierrez, Y. Le Lirzin, A. Nourry, B. Orme, M. Pradas, and R. Ledesma-Aguilar, “Snap evaporation of droplets on smooth topographies,” Nat. Commun. 9, 1380 (2018).
75Y. Cheng, D. Rodak, C. Wong, and C. Hayden, “Effects of micro- and nanostructures on the self-cleaning behaviour of lotus leaves,” Nanotechnology 17, 1359–1362 (2006).
76G. Carbone and L. Mangialardi, “Hydrophobic properties of a wavy rough substrate,” Eur. Phys. J.E 16, 67–76 (2005).
77Y. Bami Chatenet and S. Valette, “Elucidating the lotus and rose-petal effects on hierarchical surfaces: Study of the effect of topographical scales on the contact angle hysteresis,” J. Colloid Interface Sci. 676, 355–367 (2024).
78B. Mehrjou, Y. Wu, P. Liu, G. Wang, and P. K. Chu, “Design and properties of antimicrobial biomaterials surfaces,” Adv. Healthcare Mater. 12, 2202073 (2023).
79S. V. Oopath, A. Baji, M. Abtahi, T. Q. Luu, K. Vasilev, and V. K. Truong, “Nature-inspired biomimetic surfaces for controlling bacterial attachment and biofilm development,” Adv. Mater. Interfaces 10, 2201425 (2023).
80B. Jin and J. He, “Bio-inspired nanoarchitectonics of Re-entrant geometries toward transparent omniphobic surfaces,” ACS Appl. Opt. Mater. 1, 1620– 1626 (2023).