A ternary phase‐field model for wetting of soft elastic structures

Soft wetting, that is, the interaction of liquid–fluid interfaces with deformable elastic structures, provides a rich variety of physical phenomena (stick‐slip motion, durotaxis, Shuttleworth effect, etc.) which are not yet understood. We propose a novel phase‐field approach to study such problems. The method uses two phase‐fields to describe the three domains (solid, liquid, ambient fluid) by the ternary Navier–Stokes Cahn–Hilliard equations. We use a recent phase‐field approach for fluid‐structure interaction to equip one phase with elastic properties. The resulting method is significantly simpler than all previous approaches for soft wetting, since no grids are used to represent the geometries. Accordingly, the method avoids grid‐based computation of interface curvature, mappings between different grids, complicated mesh generation, and retriangulation, while the elastic object can move freely through the computational domain. As a special feature, the elastic structure is for the first time represented as a neo‐Hookean Kelvin–Voigt–Maxwell material which makes it very versatile. The accuracy of the method is shown in a benchmark simulation of a droplet on an elastic substrate. The geometrical flexibility is illustrated by simulating a rotating solid object within a fluidic interface. Finally, we provide the first 3D simulations of soft wetting including solid surface tensions.

printing or additive manufacturing. Still, our understanding of the dynamics of such soft wetting processes lags by far behind of what is known about wetting of rigid structures. 7 Starting only in 2016, the increased technical importance of soft wetting has led to a couple of numerical simulation models. [8][9][10][11][12][13] These references comprise basically the whole available literature on numerical simulations of soft wetting. The mere number of these papers is vanishingly small if compared with publications on computational progress in wetting of rigid structures or ordinary fluid-structure interaction (involving only a single fluid). One reason for this discrepancy is certainly the challenging nature of the problem which involves a strong coupling of hydrodynamics and structure dynamics near the contact line. 13 Another reason for the lack of computational results in this field lies in the diversity of the required numerical techniques: The simulation of soft wetting requires to marry two-phase flow with fluid-structure interaction (FSI). While the simulation of two-phase flow typically employs interface capturing techniques (like the phase-field method 14,15 ), FSI is typically modeled by interface tracking techniques, in particular the arbitrary Lagrangian-Eulerian (ALE) method. 16 The phase-field model (also diffuse-interface model) is arguably the most popular continuum model to describe wetting phenomena with moving contact lines on rigid substrates. The idea is based on an order parameter, the phase-field function, which is used to indicate the fluid phases. The phase-field function varies smoothly across the interface leading to a (thin) diffuse interface which is physically motivated by the composition of real fluid-fluid interfaces. 17 The diffusive nature of the interface regularizes the stress singularity at the contact line, making it a very natural approach for moving contact lines. 13 This makes phase-field models physically more sound than sharp interface models, like level-set or ALE, which typically need to introduce experimental fitting parameters, 7 while offering an almost similar accuracy for two-phase flow simulations. 18,19 Phase-field models can be derived directly from variation of the underlying energy. This sound thermodynamic foundation implies a consistent description of topological transitions (e.g., Reference 20), and enables energy stable discrete formulations 15 and robust time discretizations. [21][22][23] We refer to Reference 24 for a classification of phase-field models for two-phase fluids with incompressible constituents. Consistent boundary conditions for the phase-field modeling of wetting have been derived in References 14,25,26. In consequence, all available approaches for simulation of soft wetting [8][9][10][11][12][13] use phase-field methods to describe the two involved fluids. But the third phase, namely the elastic structure, is in all approaches discretized by a fundamentally different interface tracking method which describes the solid structure by a moving computational ALE grid. The reason of this inconsistency is certainly that the ALE method is the standard method for FSI. In consequence, all previous methods involve a certain complexity, coupling explicit interface tracking with moving grids (ALE) to implicit interface capturing (phase-field) methods.
While the ALE methodology provides a sound mathematical description and leads to a very accurate domain representation, it also comes with limitations on the evolution of the solid structure, which breaks down for large deformations or large translational and rotational movements. Even worse in the context of soft wetting is the high complexity involved in the computation of surface tension forces of the fluid-solid surfaces. Accordingly, some of the previous approaches choose to ignore solid surface tension, 10,11 which is a clear handicap as soft wetting is dominated by surface tension forces. All other approaches 8,9,12,13 have only presented 2D and 3D axisymmetric results, which makes it relatively simple to explicitly compute surface curvature and hence surface tension. 27 The second issue associated to solid surface tension is that its explicit computation leads to time step restrictions 28 which can be quite severe at the small length scales of typical soft wetting (surface tension force scales quadratically with length scale), and limits current soft wetting methods significantly.
In recent years, the above problems have led to the development of alternative modeling approaches for fluid structure interaction, for example fully Eulerian formulations 29 and also phase-field models. 30,31 The goal of this article is to apply such a phase-field approach to soft wetting. That is, we propose a numerical method which uses a phase-field method not only for the two fluids, but also to describe the elastic structure. In consequence the numerical approach becomes very consistent and much simpler than previous methods for soft wetting. Additionally, surface tension forces can be included for all three surfaces in a simple yet stable way omitting time step restrictions 21 and large deformation, translation, and rotation of the elastic structure is made possible without retriangulation.
The article is organized as follows. In Section 2 we present the numerical model equations. These equations are implemented in the CUDA GPU framework with details given in Section 3. Four numerical tests are provided in Section 4 to illustrate the applicability, validity, and the potential of the method. Finally, conclusions are drawn in Section 5.

MATHEMATICAL MODEL
Computational models of soft wetting describe the interaction of three distinct phases: two-fluids and one solid phase. We will refer to the fluids as a liquid and an ambient phase. The latter is typically a surrounding gas, but could also be a second liquid. All three phases constitute the computational domain Ω, with Ω ⊂ R d , d ∈ {2, 3}. An illustration is given in Figure 1 where the typical example of a drop on a substrate is shown. Several phase-field models with three components have been proposed in the literature, for examples, References 32-34, all of them consider only fluid phases. The most commonly accepted model proposed by Boyer et al. 33 can be derived from energy variation of a Cahn-Hilliard energy for three phases augmented by a sixth-order polynomial potential which makes the system well posed. In the following we will briefly describe this model and then extend it to account for the elasticity of one phase.

2.1
Three-phase flow model The three distinct phases can be described by three-order parameters (phase-fields) c 1 , c 2 , and c 3 , each representing the "concentration" of one component (i.e., c i ≈ 1 in phase i and c i ≈ 0 elsewhere). By definition this implies that All three phase-fields will be comprised in a vector c = (c 1 , c 2 , c 3 ) for shorter notation. As typical in phase-field models the interface between two distinct phases is diffuse, meaning that it has a thickness of size over which the transition from 0 to 1 takes place, such that the phase-fields are continuous functions. The interface thickness is a parameter in the governing equations and is generally chosen such that the equations approach the corresponding sharp interface equations. The evolution of the phase-fields is governed by the following set of Cahn-Hilliard equations 33 Here u is the velocity, D > 0 is a mobility constant, i is the chemical potential. The surface tensions between the respective phases i and j are denoted by ij . From these one can compute the three capillary coefficients i = ij + ik − jk , for pairwise distinct i, j, k ∈ {1, 2, 3}, and their harmonic mean T = 3∕ . These values are used to construct the ternary bulk energy 33 F I G U R E 1 Illustration of the standard soft wetting scenario.
A solid domain Ω 1 is in contact with two viscous domains Ω 2 and Ω 3 . All domains together comprise the computational domain Ω The last term is only necessary in total spreading cases (when i < 0 for one i), otherwise one can set Λ = 0. Note, that i can be negative for one i even though the surface tensions ij are obviously always positive. The system of equations is still well posed, if the other capillary coefficients are sufficiently large. We refer to Reference 33 for the derivation of the above equations and for more specific conditions on well posedness. The equations ensure phase conservation Equation (1). Consequently, one of the Cahn-Hilliard equations can be eliminated as the third phase-field can be computed from Equation (1). To couple the system with hydrodynamics, we will exploit the fact that all soft elastic materials (soft wetting occurs typically at Young's moduli of few kPa or less), are essentially gels, which are incompressible, due to high fluid content. 31 As a result the fluid and structural mechanics can be described in a unified way by a Navier-Stokes equation, Here, p is the pressure, i c i denote the phase-dependent density and viscosity, respectively, for given values i , i in phase i. The last term in Equation (4) contributes the surface tension force of all three phases. The system (4)-(5) is the extension of the classical model-H 35 to ternary viscoelastic fluids. The term EF denotes an elastic force, which we introduce in the solid domain. For EF = 0 the equations describe the well-accepted model for ternary fluids. 33

Including solid elasticity
We will now extend the previous model from Reference 33 to account for the elasticity of one phase following the approach in Reference 31. The solid domain will be identified by c 1 , the liquid phase by c 2 , and the ambient phase by c 3 . To express elastic forces it is necessary to define the strain of the solid material. Strain is typically computed by memorizing the initial coordinates of each solid material point. As this information is not available in an interface capturing model (like the phase-field model), an evolution equation has been proposed to keep track of the strain 31 Here Given the strain B, the elastic force of a neo-Hookean material is described by where G is the solid shear modulus. The presence of c 1 in the above equations ensures that the elastic stress is restricted to the solid domain. Equations (2)-(7) comprise the complete computational model for soft wetting. Comparing to previous approaches it is significantly simpler, in particular the model avoids moving meshes, mappings between different meshes and several jump conditions which need to be specified at the solid interface to ensure no-slip, balance of interfacial forces, no penetration, and correct contact angles. 13 The underlying equation for the strain tensor (6) describes the solid material as a Kelvin-Voigt material with elastic shear modulus G and viscosity 1 . By setting (c) = 1 for all c in Equation (6) the solid model can be extended to a Kelvin-Voigt-Maxwell material with Maxwell relaxation time T. Previous models for soft wetting include only purely elastic rheology [8][9][10][11][12] or Kelvin-Voigt rheology. 13 In our proposed model the richer solid rheology comes in naturally without additional effort.

Initial and boundary conditions
Suitable initial conditions for the ternary system (2)- (7) are provided by a specification of the initial phase distribution, initial velocity and initial strain tensor: While the given functions (c 0 1 , c 0 2 , c 0 3 ) are chosen depending on the considered problem, the initial conditions for velocity and strain are fixed throughout this work: u 0 = 0 and B 0 = I. The latter describes that the elastic body is initially in its reference configuration.
The ternary Cahn-Hilliard equations are equipped with the following set of boundary conditions: where n denotes the derivative in the direction normal to Ω. Other boundary conditions are possible, for example, to impose different contact angles, we refer to References 8,15 for more details. The outer domain boundary is considered as a solid wall in this article unless otherwise stated. Hence, u = 0 on Ω.

Special cases
We will now consider special cases, when the initial condition consists of only two phases. We will show that in this case the model equations (2)-(7) reduces to standard model for two-phase flow or fluid-structure interaction.

No elastic phase
Consider the case, that c 1 = 0 at t = 0, so there being no elastic phase. A useful initial condition must satisfy, Equation (1), that is, c 2 = 1 − c 3 . Using 23 = ( 1 + 2 )∕2, one can compute 1 = 0, therefore t c 1 = 0. With this c 1 will always be equal to 0. Hence, (c) = 0 and Equation (6) reduces to B = I, and the elastic force vanishes, EF = 0. Further, the equations for c 1 and c 3 become obsolete, as a single phase-field suffices to keep track of the remaining two phases. After simplification of the remaining equations, we obtain the system, Thus, the governing equations are reduced to a usual Cahn-Hilliard Navier-Stokes system for two-phase flow, the two phases being the phases 2 and 3, and the surface tension being 23 .

Only one fluid phase
In the case, that c 3 = 0 at t = 0, so there being no ambient phase, but only a liquid and elastic phase. Following the same arguments as before, c 3 remains zero for all times. From Equation (1) the equations for c 2 and c 3 become redundant. Hence, the model reduces to the phase-field model for fluid-structure interaction described in Reference 31:

IMPLEMENTATION
The discretization of the governing equations (2)- (7) can be done with any numerical method for partial differential equations. This is in contrast to all previous approaches for soft wetting which all rely on the grid-based representation of the solid substrate, and hence employ finite element discretizations. Here, we use a simple but efficient finite difference discretization to exploit this flexibility. For the Cahn-Hilliard Navier-Stokes system we follow the validated discretization in Reference 36 and use a finite difference scheme with staggered equidistant grid and a forward Euler method in time.
The additional equation for the elastic strain tensor is also evaluated in the cell-centers. Further details can be found in Reference 36. Parallelization of the finite-difference scheme is straightforward. To obtain a highly scalable parallel implementation, we use the computational power of graphic cards (GPUs) within the CUDA framework from Nvidia. This allows any Nvidia GPU to be used for general purpose computations (GPGPU). Nvidia offers ready-made libraries for frequent use cases, such as multiplication of big matrices. But because there were no such libraries for our use case, we use the bare CUDA framework to implement all computation steps. Each discretized equation was implemented as a separate CUDA function, called a kernel. Each of these kernels calculates the corresponding equation for all grid points of the domain in parallel. Because we needed to launch a lot of these kernels and each launch is associated with some overhead, we decided to utilize a new feature in CUDA 10 called graphs. This allows to build a graph of dependent kernel launches once and to invoke it multiple times without causing the overhead of each individual kernel launch. For us this meant, that we build up a graph consisting of all the kernels of the model and invoking this graph repeatedly for each time step.
The performance benefit of using a GPU for this simulation can be seen in Figure 2, where we show the speedup of the calculation of one time step over the serial calculation with one CPU core and the parallel calculation on multiple CPU cores using OpenMP * . It is to be noted though, that the CPU implementation of this comparison is derived from the GPU implementation and therefore not further optimized. Better performance might be achieved with optimized CPU code.

Validation study
We first study a popular wetting scenario: the relaxation of an initially half-spherical drop on a viscoelastic substrate, see Figure 3. In this test case the formation of a stationary wetting ridge has been reported, the shape of which can be seen as a standard benchmark quantity for soft wetting. Wetting ridge shapes of this case have been numerically computed in References 8,9,11-13 and experimentally measured in Reference 37. The experimental case is axisymmetric, which some numerical codes mimic by using the axisymmetric Navier-Stokes equations in 2D. 8,9,12,13 Here, we use a purely 2D implementation as also done in Reference 11. We can therefore expect to get an accurate description of the wetting ridge while the indentation below the drop should be underestimated (due to the lower Laplace pressure in a 2D drop). The parameters are chosen as in References 8,9,12,13, but repeated here for completeness. 12 = 36 mN∕m, 13 = 31 mN∕m, 23 = 46 mN∕m.
The shear modulus of the soft substrate was chosen as G = 1 kPa. The width of all three interfaces was set to = 5 μm. Further parameters are taken from Reference 12, 1 = 1000 kg∕m 3 , 2 = 1260 kg∕m 3 , 3 = 10 kg∕m 3 , Since there is no viscosity specified for the solid in Reference 12, we have chosen a negligible small value for 1 . At the bottom of the solid substrate we prescribe u = 0, while we set free-slip conditions on all other boundaries. To resolve the equations the time step size was set to 0.02 ns and the equidistant grid spacing was set to h = 1 μm. Note, that due to the explicit time discretization the present time step is several orders of magnitude smaller than in the time-implicit formulations. 12,13 The time step restrictions emerge mostly from the explicit treatment of the fourth-order Cahn-Hilliard equation as well as from the explicit coupling of elasticity and hydrodynamics. The vast number of computed time steps is compensated by the very low computational load per time step and the high parallel scalability of the explicit time-discretization. The constant mobility was set to D = 5 × 10 −14 m 3 /s, which was checked to be small enough to have negligible influence on the geometry evolution, but large enough to keep the interface profile in shape. The diffusion constant of the elasticity was chosen as D B = 0 as no visible oscillations occurred and it was therefore not necessary. Note, that this is due to the very mild advection of the present test case. For more pronounced motion, oscillations in B would occur and require stabilization by choosing D B > 0.
The first comparison considers the profile of the substrate in the equilibrium state, such that the choice of viscosities and densities is not relevant. The (almost) stationary state is shown in Figure 4. The surface tension between the fluids has pulled up the substrate at the triple point leading to a wetting ridge of approximately 9 μm height. Figure 5(A) shows the resulting substrate profile at t = 2 ms compared with the results in References 11,12. A comparison with both references shows a good agreement of the substrate profile along the wetting ridge, as well as under the ambient fluid.
In the region under the droplet, the substrate profiles of both references deviate significantly from each other, with our result being in between the other two. The indentation is about half as much as the result from Reference 12. This behavior can be explained by the different Laplace pressures in a two-and three-dimensional droplet. The Laplace pressure p in the drop is a result of the surface tension force contracting the droplet surface and its curvature: p = 23 (1∕R 1 + 1∕R 2 ), where R 1 , R 2 are the principal radii of curvature. One of these radii is missing (i.e., infinite) in a two-dimensional simulation, such that the Laplace pressure of a 2D droplet is only half as big as that of an equal sized spherical three-dimensional droplet. Having that the drop exerts only half the pressure onto the substrate, it is expected to indent the substrate roughly half as much as compared with Reference 12. A similar behavior is seen for the result in Reference 11, which also used a two-dimensional setup, and shows an even smaller substrate indentation. Figure 5(B) shows the evolution of the height of the wetting ridge in the first 2 ms, where most of the dynamic behavior takes place. A good agreement between the results is observable for the first 800 μs. Reference data for later times is not available. 12 One can clearly see, that the speed of growth of the wetting ridge gets smaller the closer it gets to the equilibrium state, reaching an almost stationary state at the end. Finally, let us note that this benchmark is an extreme case for our model as the benchmark quantity (height of the wetting ridge) is almost of the same order of magnitude as the interface thickness itself. While the reference data was produced with a sharp interface representation of the substrate, the present model uses a diffuse interface to model this slight elevation of the substrate. Accordingly, the obtained agreement is impressive.

Fluid substrates
The model allows to continuously tune all fluid parameters from purely viscous to purely elastic. To illustrate this we simulate the three phase flow scenario in which a liquid droplet slides along the interface of two other fluids. The initial conditions can be seen in Figure 6(A). The domain consists of a rectangle with 1200 μm length and 440 μm height. The droplet is prescribed as a half-spherical cap, again with a radius of R = 176 μm. The movement of the droplet was initiated with a constant force GF = gdc 2 added to the Navier-Stokes equation, where g = 981 m/s 2 is a constant force factor, for example gravity, d = (0, −1) T is a vector in which to apply the force. Note, that we choose a force 100 times higher than earth gravity to obtain a faster dynamics. The force is only added to the droplet, as in a Boussinesq approximation, therefore it was multiplied with the droplet phase-field c 2 . All three liquids are equipped with equal parameters: 23  The results are displayed in Figure 6, the corresponding video can be found in the supplementary material. One can see that the interfaces first deform to adjust the contact angle condition, that is, such that all three angles are approximately 120 • , Figure 6(B). In the meantime the droplet accelerates slowly, driven by the imposed force, Figure 6(C). Finally, a quasi-stationary flow profile has developed around t = 14.4 ms, Figure 6(D). In this state the droplet shape and the flow profile remain steady as they move with the droplet. One can clearly observe the contact angle hysteresis between the front and rear of the droplet, which is known from experiments with moving droplets. The droplet shape is more flattened at the front while it is more pointed at the rear. Contrary to two-phase flow, where a moving droplet usually exhibits two internal vortices, the present droplet assumes an almost homogeneous internal velocity profile. This seems to be a consequence of the additional fluid-fluid interface (green/gray interface in Figure 6). The presence of this interface induces contact angle conditions which F I G U R E 6 Simulation of three phase flow. All three liquids have the exact same parameters and surface tensions. The blue droplet is additionally propelled to the right by an external force change the droplet shape at the front and rear of the droplet in a way that the development of inner vortices is inhibited.

Droplet sliding over elastic pillars
To show the potential of the method to handle larger elastic deformations we consider a liquid droplet sliding over a viscoelastic microstructured substrate. The substrate consists of small elastic pillars (see Figure 7). Droplet motion on such structures is of high practical relevance as micropillar arrays are widely used to create hydrophobic surfaces in technical and everyday applications. The initial setup can be seen in Figure 7(A). The computational domain is a rectangle of 800 μm length and 400 μm height. A droplet is initiated as a half spherical cap of radius R = 180 μm placed upon the pillars. Each pillar has a height of 140 μm and a width of 50 μm. To initiate movement of the droplet an additional gravitational force GF = gdc 2 was added to the Navier-Stokes problem. Here, g denotes the gravity constant and d a vector indicating its direction. This force is only added to the droplet, as in a Boussinesq approximation, therefore it was multiplied with the droplet phase-field c 2 . Furthermore, we use the following parameters 23 = 46 mN∕m, 13 Figure 7, the corresponding video can be found in the supplementary material. As expected the tops of the pillars stay attached to the droplet and the pillars bend by the moving droplet until the elastic force is big enough to detach the backmost pillar from the droplet (Figure 7(B)). Shortly after its detachment this first pillar springs back toward its original position. Thereby, the pillar behaves like a damped harmonic oscillator: it oscillates slightly around its initial position until its movement is fully damped by fluid viscosity (Figure 7(C)).
We also set up the same experiment within a three-dimensional domain, the results of which can be seen in the Figure 8. The parameters we chose for this experiment were exactly the same as for the two dimensional setup except for the distance between the grid nodes, which we chose as 4 μm, the time step size of 1.5 ns and = 10 μm. A video of the 3D simulation can be found at: https://youtu.be/cKfIQtODdRI. Note, that this simulation is not only the first simulation of a droplet on a nanostructured surface, but also the first 3D simulation of soft wetting including solid surface tensions.

Water wheel
In the last scenario we want to show that the approach is capable of handling big deformations of the initial domain, which is a clear advantage with respect to previous approaches for soft wetting. To this end, we set up a test case in which a viscoelastic wheel rotates through a basin of water, see Figure 9(A). The computational domain consists of a square with 400 mm side length. A wheel is placed in the center of the domain, with four wings, each of which has a length of 100 mm and a width of 45.5 mm. The domain is filled with water up to a height of 100 mm. The wheel is propelled by prescribing a rotational velocity field in a circular region around its center (x 0 , y 0 ): Here v disc = 100 mm/s is the outer velocity of the prescribed circular disc of radius R = 50 μm.
To keep the water domain from just following the wheel movement an additional gravitational force GF = gdc 2 was added to the Navier-Stokes problem, similar to Section 4.3. In this case we chose g = 9.81 m/s 2 and d = (0, −1) T so the force is pointing downward to the bottom of the domain. The grid was set to be an equidistant grid with a spacing of h = 2 mm.
Due to the larger spatial scale of the present test case, the time step size could be increased to 3 μs and. The constant mobility was set to D = 6 × 10 −7 m 3 /s. All boundaries of the domain were set to no-slip conditions. Furthermore we chose the following parameters: 23 = 36 mN∕m, 13 = 10 mN∕m, 12 = 40 mN∕m, 1 = 2000 kg∕m 3 , 2 = 1000 kg∕m 3 , 3 = 10 kg∕m 3 , 3 = 0.0018 Pa ⋅ s, = 10 mm, G = 4000 Pa, T = 1 s Figure 9 shows snapshots of the simulation, the corresponding video can be found in the supplementary material. Due to the prescribed constant velocity in the center of the wheel and the elasticity of the solid the whole wheel starts spinning counterclockwise. As soon as one of the wings comes close to the water surface it starts bending backward due to the higher viscosity and density of the water phase compared with the ambient phase. At the same time another wing leaves the water surface while pulling up some water which sticks to it due to surface tensions. As the wing continues its movement the attached water eventually detaches due to the prescribed gravitational force.
Over the course of the simulation one can observe a slight rounding of the elastic wheel. This is caused by the solid surface tension which constantly contracts the solid surface. Contrary to sharp interface models, the solid elasticity cannot completely counteract such forces in a phase-field model. Instead small errors may accumulate over time, in particular when the elastic structure undergoes a large translation. This rounding effect is in our view the only downside of the present model.
Nevertheless, the present simulation is the first soft wetting simulation which involves any movement of the elastic domain at all. Having a phase-field to represent the domain makes such a scenario simple, since no retriangulation of the mesh is needed (as opposed to all previous approaches for soft wetting simulations). This further illustrates the potential of the method.

CONCLUSION
Soft wetting is an emerging field of research with many potential applications and a richness of open questions and interesting physical phenomena (e.g., stick-slip motion, 38 durotaxis, 37 Shuttleworth effect 39 ) which demand for numerical simulation tools. In this article we present a novel method to simulate soft wetting problems which differs largely from previous methods by using a phase-field to describe the solid phase. The result is a certain simplicity of the model: a unified Navier-Stokes equation describes the momentum and mass balance in all three phases and intrinsically includes all coupling conditions (no-slip, surface force balances, contact angle conditions, no penetration). Since no grids are used to represent the solid geometry, the method avoids grid-based computation of interface curvature, mappings between different grids, complicated mesh generation, and retriangulation. Hence, the method is so simple that it does not need a very long explanation, and it can be implemented with any numerical discretization method for partial differential equations. Here we use a simple finite-difference scheme with forward Euler method to discretize in space and time. While this discretization scheme is extremely fast per time step and offers a high parallel speed-up, it comes at the expense of strong time step restrictions. We note that the use of small time steps is not a restriction of the model itself.
In fact, by virtue of the fixed-grid character of the model it is relatively straightforward to set up a time-implicit formulation which should allow several orders of magnitude larger time steps, as in References 12,13. As a special feature, the elastic structure is for the first time represented as a neo-Hookean Kelvin-Voigt-Maxwell material which makes it very versatile.
We illustrate the accuracy of the method by using the standard benchmark test of soft wetting and find impressively good agreement despite the challenging length scale for a diffuse interface model. We further provide the first simulations of wetting on a microstructured surface and the first 3D simulations of soft wetting including solid surface tensions. Finally, we illustrate the capability of the method to freely move elastic objects through the computational domain (without retriangulation) which has not been shown with any other method so far.
The only downside of the used phase-field approach is a slight rounding of the elastic object which can accumulate over time. This effect scales with (see Reference 31), and can hence be seen as a usual numerical error, which can be controlled by refining and grid size. Nevertheless, the flexibility, efficiency and superior simplicity of the method will permit investigations of many soft wetting phenomena which cannot be done with any other numerical approach so far. We will explore these capabilities in a set of follow-up articles.