Simulating liquid–gas interfaces and moving contact lines with the immersed boundary method

Physical systems involving the coupling between a fluid and evolving immersed structures are usually impossible to describe analytically. The immersed boundary (IB) method is a numerical approach for solving such problems. Traditionally, this method is used with the assumption that the immersed boundaries are massless.¹ The method describes the fluid in Eulerian coordinates and the immersed boundaries with arrays of linked Lagrangian markers. The fluid advects the markers and the markers exert forces onto the fluid. The massless-boundary assumption is suitable for describing thin elastic membranes common in biological applications, for example, the interaction between blood flow and heart valves.¹ Similarly, the massless-boundary assumption is appropriate for liquid–gas interfaces. Thus, with a surface tension model, the IB method is capable of modeling multi-phase fluid flow.^2–5 This introduces a new challenge, which is the moving contact line (MCL) problem that emerges when a solid boundary meets with a liquid–gas interface. To avoid the force singularity at the contact line, Lai et al.⁶ proposed an MCL model that simulates Navier slip with the IB method, but it is limited to fixed solid boundaries.

In this work, we describe a static-dynamic friction condition between a moving liquid–gas interface and an evolving immersed boundary. Two types of immersed boundaries (liquid–gas interface, solid surface) coexist and are governed by the same IB numerical method. The slip condition used in this paper is informed by recent molecular dynamics (MD) simulations.⁷ The resulting numerical scheme, unlike related works that either prescribe the contact angle^8,9 or the slip velocity,^6,10 is capable of reaching the dynamic contact angle in an endogenous way, i.e., the dynamic contact angle is the result of droplet size, wall friction, gravity, and other local stresses.

In addition to the MCL model, we propose a surface tension model within the IB method. Using the IB method to simulate surface tension is often classified as a front-tracking scheme. Most front-tracking methods, including many IB variants, calculate the magnitude of surface tension by estimating the local curvature of the interface using three to six markers.^2,5 This approach is susceptible to errors and instability issues. Tryggvason et al.¹¹ computed the surface tension using tangent vector subtraction in every Eulerian cell, thus omitting the need to estimate curvature. Popinet and Zaleski¹² used tangent vector subtraction for every link between adjacent Lagrangian markers, exploiting the IB method's Lagrangian representation of the structure. Our method uses tangent vector summation for every Lagrangian marker (instead of every link). The subtle shift in perspective brings a side benefit that the component of the Young's force due to liquid–gas surface tension¹³ at the contact point is readily computed by the two markers at both ends of a liquid–gas interface. In addition, our approach uses a stepwise re-sampling procedure to ensure a near-uniform distribution of Lagrangian markers in the liquid–gas interface. Finally, we employ a stepwise interface splicing method to implement interface topological changes, including droplet coalescence and separation.

This study is motivated by a real-life industrial problem, presented by W.L. Gore & Associates in 2018 during the Mathematical Problems in Industry (MPI), which is a workshop that attracts leading applied mathematicians and scientists from universities, industry, and national laboratories.¹⁴ Catalysts are an integral part of many chemical processes. They are usually made of a dense but porous material, such as activated carbon or zeolites, which provide a large surface area. Liquids that are produced as a by-product of a gas reaction at the catalyst site are transported to the surface of the porous material, slowing down transport of the gaseous reactants to the catalyst site. One example of this is in a sulfur dioxide filter, which converts gaseous sulfur dioxide to liquid sulfuric acid.^15,16 Such filters are used in power plants to remove the harmful sulfur dioxide that would otherwise contribute to acid rain. Understanding the dynamics of liquid droplets in the gas channel of a device is critical in order to maintain performance and durability of the catalyst assembly. Among several other tests, the methods presented in this paper are applied to simulate 2D droplets moving on a vertical wall, which corresponds to this application. In this paper, we also empirically study the spatial and temporal convergence of our methods, compare simulation results at equilibrium against an analytical solution, and benchmark our methods with several standard test cases.

The paper is organized as follows: in Sec. II, we describe the IB method, the boundary conditions, the surface tension, and Young's force formulation, the moving contact point friction, the static-dynamic friction boundary condition, the stepwise interface re-sampling procedure, and the interface splicing procedure. An approach for simulating variable fluid density is also discussed. We then introduce the discretization of the governing equations in Sec. III. In Sec. IV, we verify the accuracy of our simulation results. In Sec. V, we apply our methods to several test cases. Finally, we conclude in Sec. VI with a discussion of our model and results, along with some insight into real-world applications.

In this section, we describe the IB method that we use in this paper. IB has proven to be an extremely versatile method for fluid–structure interaction problems since its development almost forty years ago.^1,17–23 The method represents the immersed structure with Lagrangian coordinates and the fluid with Eulerian or “lab frame” coordinates. The interactions between the two coordinate frames are communicated with integral transforms involving Dirac delta function kernels. For the solid surfaces and the variable density considered in this paper, we adopt the penalty immersed boundary (pIB) method.^24,25 For example, the vertical wall introduces a set of tether points that are fixed in space and represent the desired shape and location of the structure. The boundary markers are connected to the tether points via springs, and only the boundary markers interact with the fluid. Stiff springs approximate a rigid boundary. These stiff springs typically result in some numerical stiffness but simplify the implementation. In practice, accurate results can be achieved with an appropriate choice of the spring constant, spatial resolution, time stepping, and other numerical parameters.^23–25 The penalty immersed boundary method has been used to simulate a wide variety of problems that involve immersed structures with mass or variable density fluids. For example, Kim and Peskin used this approach to study oscillations of a cylinder in fluid as well as performed three-dimensional simulations of windsocks and flags.²⁶ The simulations of the oscillating cylinder were in reasonable agreement with experimental data.²⁶ For variable density problems, this method has been used to simulate bubble dynamics and Rayleigh–Taylor instabilities.²⁴ A disadvantage of this approach, as with many types of immersed boundary methods, is the regularization of the interface between low and high density fluid regions. The interface smoothing occurs because of the representation of high density regions using a discrete set of Lagrangian markers and the approximation of the Dirac delta function kernel, to be described below.

The equations of motion for the coupled fluid-structure system are

Here, ρ_g is the density of the gas, μ is the dynamic viscosity of the gas and the liquid, and t denotes the time variable. The function $u(x,t)$ is the velocity, and δ is the 2 D Dirac delta function. The vectors x and $Xi$ denote the Eulerian fluid coordinates and the Lagrangian structure coordinates. In this paper, we consider three types of immersed structures: a 1D wall boundary (⁠$X1$⁠), a 1D fluid–gas interface boundary (⁠$X2$⁠), and a 2D variable density area (⁠$X3$⁠). The vectors $X4$ and $X5$ are ephemeral markers used to impart the Young's force and the contact point friction force. For the 1D structures, s is the arc length that specifies a point on the boundary (i.e., a parametrization of the 1D boundary with respect to arc length). For the 2D structures, (r, s) is a 2D equidistant parametrization of the area. The functions $fi$ and $Fi$ are the Eulerian and the Lagrangian force densities corresponding to the immersed structures. The Lagrangian force densities $F1$ and $F2$ are 1D force densities, $F3$ is a 2D force density, and $F4$ and $F5$ are 0D force densities. Subsections II A–II G will introduce each type of immersed structure corresponding to $Xi$ and define its Lagrangian force density $Fi$⁠. Equations (1) are the Navier–Stokes equations for incompressible flow of a viscous fluid. In the regime, we study in the paper, there is little compressibility since Mach number is low; therefore, we simplify the simulation by assuming that all fluids are incompressible. Equation (2) calculates the force imparted from the immersed structures onto the fluid by converting Lagrangian force densities to Eulerian force densities using the Dirac delta function kernel. Equation (3) also uses the Dirac delta function to enforce the condition that the velocity of the immersed structure is equal to the fluid velocity, corresponding to a no-slip and no-penetration condition between the fluid and structure.

Our code is adapted from the Matlab implementation of the IB method by Peskin,²⁷ which uses periodic boundary conditions for the computational domain. In most of our simulations, we fold the domain in half to obtain a symmetry boundary condition on the left and right boundaries. This section describes our implementation of symmetry boundary conditions within Peskin's implementation. In short, symmetry is prescribed in the Eulerian force density field, then the Eulerian fluid treatment solves for both symmetric halves simultaneously on a periodic domain. Specifically, if the computational domain exposed to the fluid solver is a square with side length L, as shown in Fig. 1, then we select the left half of it to be the problem domain. According to Fig. 1, L = 2 is chosen. During each time step, once the Eulerian force densities (given by (2)) are imparted onto the grid, the simulation mirrors the Eulerian force densities matrix around $x=L/2$ (i.e., the dashed line in Fig. 1) and adds them back to the original Eulerian force field. In this way, the fluid velocity field is always symmetric around $x=L/2$⁠. Analytically, this symmetry treatment prevents structures from crossing the axis of symmetry.

The main case that motivates our study is the simulation of liquid droplets moving on a solid vertical wall. We simulate the wall as an immersed boundary. In this specific case, since the wall is static, it is arguably easier to treat the wall as a boundary condition. However, we want our MCL method to generalize to non-static solid surfaces (see, e.g., Sec. IV C 4), so we treat the wall as an immersed boundary to maximize generality. Each Lagrangian marker of the wall is tethered to its ground-truth location, thus ensuring no-penetration and no-slip. In this case, the Lagrangian force density corresponding to the 1D wall is

Equation (4) describes the tether-point construction and the associated spring force on the boundary to model a no-slip and no-penetration wall. The parameter k₁ is the spring constant (taken to be $5000 g/(s2 cm)$ in simulations) associated with the wall, which is a trade-off parameter between accuracy and numerical stability. If the spring constant is too small, the wall markers near the contact points will significantly deform under stress. If it is too large, numerical instabilities will likely arise. Since our domain is 2D, we define each solid–liquid–gas intersection to be a contact point. The vector Z is the ground-truth location for the markers. At the start of the simulation, we usually initialize $X1(s,0)=Z(s)$ (with Sec. IV B as an exception). The wall behavior will be later modified in the MCL model, as described in Sec. II E.

We use the integral formulation as described by Popinet⁴ to derive a model for surface tension. In this approach, each interface marker is pulled by its two neighbors at a constant force magnitude. In contrast to the work by Tryggvason et al.,¹¹ our method is entirely in the Lagrangian frame and does not estimate a tangent vector with a four-point polynomial fit. Instead, we use a piecewise linear fit to approximate the unit tangent vector $T̂$⁠, which is defined as

Then, the Lagrangian surface tension force density is

where σ is the surface tension coefficient. Therefore, for a segment AB of the liquid–gas interface, as shown in Fig. 2, we have

where Ω is the set of interface points s on segment AB. $T̂A$ and $T̂B$ are the unit tangent vectors at A and B, respectively. With this approach, the sum of surface tension forces in a closed loop is zero.⁴ In the discrete implementation, as shown in Fig. 2, we first compute the tangent vectors pointing from one marker to its adjacent markers and then apply the corresponding force of magnitude σ onto the marker.

A consequence of this implementation is that when the gas–liquid interface meets a solid surface, the Young's force,¹³ denoted by $FY$⁠, arises directly from the tension forces. As shown in Fig. 3, the magnitude of the Young's force is the tangential component of the sum of three tension forces

where $T̂$ is the unit vector tangent to the liquid–gas interface, $t̂$ is the upward unit vector tangent to the wall, and $t̃$ is the unit vector tangent to the wall and pointing toward the liquid phase. $σ,σsl$⁠, and $σsg$ are the liquid–gas, solid–liquid, and solid–gas tension coefficients, respectively. The Young's force is present at any three-phase contact point and is responsible for phenomena like the capillary action. The Young's force also explains why slipping behaviors are often attributed to the contact angle. The static contact angle θ_s is defined as the contact angle that will make $FY=0$⁠. Therefore, the static contact angle can be calculated by the following equation:

This also gives an alternate equation for the Young's force that is based on the static contact angle:

Here, θ_d is the dynamic contact angle. As shown in Fig. 3 and also discussed earlier, at a contact point, there is one interface marker that only has one neighbor. Therefore, (6) already supplies the first term of the three-force sum in (8). We compute the other two terms and add them to each contact point, thus completing the implementation of the Young's force by

Note that the points $X4$ are the locations where IB imparts $F4$⁠, i.e., the set of contact points.

A frictional force is applied to each moving contact point. We assume that the friction is proportional to the slip velocity and depends on the dynamic contact angle. In this section, we first assume that the wall is free-slip, contrary to the no-slip setup described in Sec. II B. Toward the end of this section, we will relax this free-slip assumption.

Specifically, for each three-phase contact point, we apply a 0D force density at that point corresponding to a frictional force as

and $t̂$ is the upward unit vector tangent to the wall (so $v=u·t̂$ is the slip velocity) and θ_d is the dynamic contact angle. The function η defines a friction coefficient (of units $g/s$⁠) in terms of the dynamic contact angle θ_d and is given by

The parameter values in (13) are based on measured results from the simulation of water-silica contact line movement by Johansson and Hess.⁷ Linear interpolation of measured data points is used to obtain $η(θd)$ as given in (13). Note that our MCL method is agnostic to the specific implementation of η; therefore, future usages are free to modify η.

The linear response of $F5$ to the slip velocity [see (12)] is supported by molecular simulations by Johansson and Hess.⁷ Such a linear relationship is related to the general Navier boundary condition (GNBC), where the slip velocity is proportional to the sum of the Young's force and the fluid viscous stress:²⁸ 

where n is the normal vector to the boundary. Similar to the Young's force, which is also a 0D force density, $X5$ is the location where IB imparts $F5$⁠. We simply set $X5$ to be the contact point locations. We can apply $F5$ to each contact point on an (otherwise) free-slip wall, but that will fail to recreate these two phenomena:

• A sufficiently small droplet should be able to hang statically on the wall without sliding down;

• The non-contact point region of the solid–liquid interface should be no-slip.

To represent the above phenomena, we apply $F5$ to a wall that otherwise follows the static-dynamic friction boundary condition, which we describe in Sec. II E.

The so-called moving contact line (MCL) problem refers to the apparent contradiction that contact lines can move on a no-slip wall (in the case of our 2D methods, it is rather a “moving contact point” problem; however, we still use the “MCL” abbreviation to align with the literature). The mechanism of an MCL and how to simulate it had been largely a mystery until 1979.²⁹ Since then, researchers have developed many numerical models of MCLs.^9,30 A recent study³¹ showed that hydrogen bonds facilitate the no-slip behavior of water on hydrophilic surfaces. These hydrogen bonds are orders of magnitudes stronger compared to the fluid's internal viscosity, so hydrophilic surfaces usually behave as if they were no-slip. This justifies the traditional no-slip assumption between water and hydrophilic surfaces. An MCL contradicts the no-slip assumption since the Young's force at the contact point is strong enough to break hydrogen bonds, therefore, allowing slip. In this light, continuum models involving MCLs should have a mechanism to allow for slip. For example, Lai et al.⁶ use the IB method with a Navier boundary condition at the wall. As opposed to the no-slip condition, the Navier boundary condition prescribes the slip velocity to be proportional to the viscous stress. However, their work has two limitations: (i) their boundary condition is static and cannot be imposed on a solid structure that changes shape; (ii) the entire boundary allows for slip, including segments far away from the MCLs. This section proposes a static-dynamic friction model for the wall, which (i) works with the IB method, (ii) can be prescribed onto an immersed structure that changes shape over time, and (iii) only allows slip near the MCLs.

Figure 4 shows that the friction force should sometimes balance with other forces while other times be too weak to balance with other forces. Panels (a) and (b) show the statically and dynamically balanced cases, respectively, where the friction balances with other forces, such as the Young's force and the internal viscous stresses. In addition, the slip velocity is zero in (a) and a constant nonzero value in (b). Panel (c) depicts the situation in which friction is not strong enough to balance with other forces, implying the slip velocity changes with time. We implement this conditional behavior using what we call length-limited springs.

The method still treats the wall as an immersed structure but changes the penalty treatment by overriding the marker advection Eq. (3) for i = 1. Let $X1̃$ denote the positions that (3) would have advected the wall markers to, i.e.,

where dt is a small difference in time. Section III will define dt more formally in the context of the numerical discretization of these equations.

Then, Table I and Fig. 5 show the boundary conditions that can be achieved by specifying how to update the tangential component of wall marker positions. In the no-slip condition [Fig. 5(a)], we set $t̂·X1(t+dt)=t̂·X1̃(t+dt)$⁠. This corresponds to the original form of (3), so the wall markers strictly follow the fluid flow. The friction always statically balances with other forces such as Young's and the internal viscous stresses. In the free-slip condition [Fig. 5(b)], we set $t̂·X1(t+dt)=t̂·Z$⁠, where Z is the ground-truth location. The wall markers are not tangentially advected by the fluid flow, and their vertical position remains at the vertical position of their ground-truth location. The springs are at rest (tangentially), and thus, no friction force is applied onto the fluid. In the static-dynamic friction condition [Figs. 5(c) and 5(d)], the length of the springs (when projected onto the wall) is limited by an upper bound q. In this case, we set $t̂·X1(t+dt)=max(t̂·Z−q,t̂·X1̃(t+dt))$⁠. The friction force is at most $q k1$ and otherwise statically balances with other forces. Effectively, if local tangential forces are strong, e.g., containing a strong Young's force from a contact point, then the spring lengths are set to q, the friction provided is $q k1$⁠, and local tangential slipping occurs. Figure 5(c) shows this slipping case. On the other hand, if local tangential forces are weak, the spring lengths do not exceed q and static friction is sufficient to produce a no-slip result. Figure 5(d) shows this no-slip case. The slip is represented by the difference between $t̂·X1̃(t+dt)$ and $t̂·Z−q$⁠. Dividing that by dt gives the slip velocity. The heat dissipation rate is quantified by multiplying the friction force with the slip velocity

The free-slip and no-slip cases can be obtained from the static-dynamic case for particular choices of q. When q goes to infinity, we obtain the no-slip condition. The length limit is effectively removed and $X1=X1̃$⁠. When q is set to zero, we obtain the free-slip condition. In this light, the static-dynamic friction condition interpolates between the no-slip and the free-slip condition by controlling the value of q.

In practice, we set $q=5Newtons/k1$⁠. This parameter value was chosen in order to obtain realistic results: The contact points can slip while the non-MCL region of the wall is no-slip, and small droplets can hang statically on the wall without sliding. To make q physically accurate, further molecular simulations or physical experiments are needed.

The static-dynamic friction condition bears similarity to the rigid-body static-dynamic friction model: The friction is a reactive force that responds to the sum of all active forces. If the sum of active forces is lower than a threshold, the rigid body is stationary and the frictional force balances with the other forces. If the sum of active forces exceeds a threshold, the rigid body slides and the dynamic friction equals the maximum static friction.

The immersed markers, as parametrized by arc length s, are initialized to be uniform. During a simulation, the markers are advected by the fluid flow, which may not conserve the arc length. The surface tension is always normal to the interface, so the fluid flow easily disrupts the distribution of the interface markers. To equi-distribute the markers on the structure, Lai et al.³² used grid redistribution, while Hou et al.,³³ and Lai et al.⁶ applied artificial tangential velocity. Our method is similar to grid redistribution in the sense that we add markers to wide gaps and remove them from tight spaces at every time step. Specifically, at each time step, the program iterates over all the interface markers. When the distance between any pair of markers exceeds $2$ times their initial distance, the program inserts a new marker between them. When the distance between the outer two of any three adjacent markers become smaller than $2/2$ times their initial distance, the program removes the inner marker of the three. Note that the value $2$ is the geometric mean of 1 and 2. Given the conditions for adding and removing markers, the average distance between adjacent markers will remain in the range $[2/2,2]$ times their initial distance.

Whenever a marker is removed, the two adjacent markers are joined together and their coordinates are left untouched. However, a sharpness check compares $cos α$ (where α is the angle formed by the three markers) with a threshold (0.9, in our case) to prevent the program from removing a marker that represents a relatively sharp corner and breaking conservation of area and tension energy. The sharpness check becomes irrelevant as the spatial resolution N approaches infinity but improves simulation accuracy for finite values of N.

On the other hand, when a marker is inserted (see Fig. 6), the program calculates two locations: the midpoint (MP) between the two markers and the intersection (IS) of the extrapolated lines from four adjacent markers. A weighted average between MP and IS is used as the location for the newly inserted marker. The averaging is weighted by a re-sampling amendment factor. If there is too much weight on the midpoint, the removal of markers will erroneously decrease tension energy while insertion will not increase tension energy. By having an optimal weight on the extrapolated intersection, we can keep the bias of the energy error to 0 (although the variance will still increase with time). However, finding the optimal amendment factor depends on an initial guess for the distribution of interface curvatures, which is out of the scope of this work. On the other hand, too much weight on the extrapolated intersection leads to alternating jagged edges and eventually to instabilities. In our numerical experiments, we set the re-sampling amendment factor to be 0.5. The amendment factor becomes irrelevant as N approaches infinity but improves simulation accuracy for finite values of N.

Our re-sampling routine constantly computes and records energy errors due to the artificial removal and insertion of markers. At the end of a simulation, these errors can be plotted for evaluation. See Sec. V B for the apparent improvement that this procedure brings. This re-sampling procedure, however, is very difficult to extend to 3D with a surface mesh.

Figure 7(a) shows six possible scenarios in which the liquid–gas interface changes topology and the chain of interface markers needs to be spliced: two droplets merging, one droplet splitting into two, a droplet attaching to the wall, a droplet detaching from the wall, a droplet splitting on the wall, and two droplets merging on the wall. Those six cases form three pairs of reversible processes as well as three pairs of scenarios whose implementations are exactly the same. The same-implementation cases are locally indistinguishable, since our splicing implementation is agnostic to the phase (liquid or gas) on each side.

To implement splicing, we store the interface markers in circular doubly linked lists. The linking direction preserves polarity information; if one follows the links in the positive direction, the liquid will always be on the right. Figure 7(b) illustrates a splicing event. When the distance between two interfaces is smaller than a threshold, we check whether the two interfaces are approaching using the sign of the dot product of their velocities. If both conditions hold true, the interfaces are spliced together. The method has two different distance thresholds, one for interface-interface events (h), and the other for interface-wall events (⁠$2.3h$⁠), where h is the meshwidth (i.e., diameter of one Eulerian grid cell). These specific parameters result in stable splicing events according to our numerical tests. Finally, a splicing event forces the simulation to skip the splicing subroutine for two subsequent time steps. This rule makes splicing events more atomic so that there will not be multiple splicing events competing to accomplish the same macroscopic effect.

Our model uses the technique proposed by Kim and Peskin²⁴ to describe the different densities for the liquid and gas within and outside of the droplet, respectively. In this approach, a uniform 2D grid of Lagrangian markers is constructed in the liquid phase to represent the density difference. This difference is achieved by allowing the markers to have an effective mass by tethering them to their massive, Newtonian counterparts using the penalty immersed boundary (pIB) method. More formally, the Newtonian particles have velocity $v$⁠,

where Y is the location of the Newtonian particles. The forces on the Newtonian particles are the tether force and gravity

where $Δρ=ρl−ρg$ is the 2D density difference between the liquid and gas, G is the gravitational constant, $ĵ$ is the vertical unit vector, and k₃ is the spring constant [taken to be $1000 g/(s2 cm2)$ in our simulations] associated with the variable density method. If k₃ is too high, numerical instability occurs and if it is too low, the allowed distance will be violated (explained below). A tether force is applied as a force density onto the fluid,

to couple the particles $X3$ with the Newtonian particles Y. Our implementation uses the variable density method,²⁴ with several modifications. Kim and Peskin²⁴ used a five-step update method for time stepping, while we use a coarser, midpoint method (see Sec. III). We also set the allowed distance³⁴ to be one tenth the mesh-width, i.e., $h/10$⁠.

For the spatial discretization, we use a global Navier–Stokes solver based on harmonics in the velocity field. This solver uses the Fast Fourier Transform to accelerate the computation, but at the cost of not being able to handle variable density/viscosity out-of-the-box (which is why we need the variable density method (Sec. II H) to distinguish the two phases). The time step is denoted by dt with units of seconds. The spatial resolution parameter N is the number of Eulerian cells in the problem domain along the y axis. The length of cell side in the spatial discretization, i.e., the meshwidth, is denoted by h with units of centimeters. Therefore, we always have $N h=L$⁠, where L is the height of the square computational domain as shown in Fig. 1.

The time step index is denoted by $n∈{0,1,2,…}$ and the physical time at time step n is given by $tn=n dt$⁠. The discrete Eulerian space coordinates that label the cells in the Cartesian mesh are $(i,j)∈{0,1,2,…,N−1}2$⁠. A physical cell position in the problem domain is given by $xij=(i h,j h)$⁠. The fluid velocity $u(xij,tn)$ at position x_ij and time tⁿ is approximated by the discrete fluid velocity, denoted by $uijn$⁠. The array of discrete velocities at time step n is denoted by $un=(uijn)$⁠. The discrete spatial parametrization variables for the immersed structures are $ℓ∈{0,1,2,…,ℓmax}$ and $m∈{0,1,2,…,mmax}$⁠, where $ℓmax$ and m_max can change during simulation due to the interface re-sampling. A point within the parametrization of the immersed structure is given by $(rℓ,sm)=(ℓ Δr,m Δs)$⁠, where $Δr$ and $Δs$ control the discrete spatial resolution of the structures. In our numerical experiments, we initialize $ℓmax$ and m_max so that $Δr=Δs=h/2$⁠. The spatial resolution of the structures corresponding to the set of Lagrangian markers has implications on discrete mass conservation and correspondingly the computed pressures. If the resolution of the structures is too coarse, for example, the droplet might “leak” through the interface and begin to nonphysically shrink in volume. The Lagrangian resolution of $h/2$⁠, which is explicitly finer than the fluid resolution, is chosen to avoid these numerical artifacts. The discrete 2D structure positions corresponding to the variable density area $X3$ are indexed as $X3,ℓmn≈X3(rℓ,sm,tn)$⁠, and we collect them in a single array denoted by $X3n=(X3,ℓmn)$⁠. We use the same notation for the discretizations of the 1D wall boundary $X1$ as well as the 1D fluid–gas interface boundary $X2$⁠, although in these cases, there is only a single parametrization variable. We represent the array of discrete structure positions at time step n by $Xn=(X1n,X2n,X3n)$⁠. $X4$ and $X5$ correspond to 0D points, so there is not a need for their discretization.

We use a midpoint method for evolving the system in time that relies on quantities at intermediate time steps ${12,32,…}$⁠. This method improves stability and constrains numerical errors. Following the approach from works by Peskin,^27,35 the procedure below describes the time discretization for the equations of motion (1), (2), and (3):

1. Given the marker positions $Xn$ and the fluid velocity $un$⁠, compute the first-order approximation of the marker positions at the midpoint step $Xn+1/2$ using Eq. (3).

2. Given the marker positions at the midpoint step $Xn+1/2$⁠, compute the force densities at the midpoint, denoted $fn+1/2$⁠, using Eq. (2).³⁶ 

3. Given the fluid velocity field $un$ and the force densities at the midpoint step $fn+1/2$⁠, compute the fluid velocity field at both the midpoint $un+1/2$ and the next time step $un+1$⁠, using Eq. (1).

4. Given the marker positions $Xn$ and the fluid velocity at the midpoint step $un+1/2$⁠, advect the markers for time dt to obtain the marker positions at the next time step $Xn+1$ using Eq. (3).

This section focuses on benchmark and convergence tests for our methods. In Sec. IV A, we simulate a hanging droplet on a wall and compare our steady state results with an analytical solution.¹⁴ In Sec. IV B, we compare our simulation results with an established variable density benchmark in which a bubble rises in a liquid column.³⁷ Finally, in Sec. IV C, we describe results from nontrivial experiments to empirically study the convergence of our methods.

In this test, we compare our simulation results with an analytical solution for the interface shape of a droplet hanging on a vertical wall.¹⁴ The model describes the hydro-static equilibrium shape of a droplet hanging on a vertical wall as a partial differential equation expressed in Cartesian coordinates. The solution to this problem implies a linear relationship between the local curvature of the interface and the altitude. Details regarding derivation of the analytical solution are given in the work by Cimpeanu et al.¹⁴ 

For these simulations, the time step is dt = 0.0001 s, the spatial resolution parameter is N = 96, the height of the computational domain is L = 2 cm, the surface tension coefficient is σ = 100 g cm/s², the gas density is $ρg=0.1$ g/cm², the liquid density is $ρl=1$ g/cm², the viscosity is μ = 1 g/s, and the static contact angle is $θs=π/2$⁠. We run several simulations and vary the magnitude of gravity coefficient G to compare curvatures in the equilibrium state with those from the analytical solution given in the work by Cimpeanu et al..¹⁴ To estimate local curvature at each marker, we compute the circle that contains the three interface markers and calculate the inverse of the circle's signed radius. The simulations are terminated at t = 0.2 s, at which time all three cases have reached a steady state. Results are shown in Fig. 8. From left to right, the gravitational constant G (cm²/s) is gradually increased to alter the droplet shape. The upper panel depicts the equilibrium droplet shape computed by our simulations. A linear relation is observed between the 2D curvature and the altitude as shown in the lower panel. In addition, there is excellent agreement between our numerical results (blue scattered dots) and the analytical solution (orange line) given in the work by Cimpeanu et al..¹⁴ In both panels, the green line represents the location where curvature is zero, i.e., the inflection point in the droplet shape.

In this section, we describe an application of our methods to a recently proposed transient multi-phase flow benchmark which describes a gas bubble rising in a liquid column.³⁷ This benchmark provides a reasonable experiment for the surface tension model, making it suitable for testing multi-phase flow methods. Specifically, we compare our simulation results to those calculated from Featflow, an open-source CFD package by Turek and Mierka.³⁷ Details of the benchmark setup can be found in “test case 1” of the work by Turek and Mierka,³⁷ except that instead of a 1:10 viscosity ratio, our simulations use a 10:10 viscosity ratio. Turek and Mierka³⁷ kindly provided us with the 10:10 viscosity ratio simulation results to benchmark against.

In this test, the time step is dt = 0.002 s, the spatial resolution parameter is N = 128, the height of the computational domain is L = 2 cm, the surface tension coefficient is $σ=24.5$ g cm/s², the gravitational acceleration is G = 0.98 cm/s², the gas density is ρ_g = 100 g/cm², the liquid density is ρ_l = 1000 g/cm², the viscosity is μ = 10 g/s for both fluids, and the Reynolds number is Re $=35$⁠.³⁷ We use a symmetry boundary condition on the left and right boundaries that corresponds to a no-penetration condition. We treat the top wall (see Fig. 9) as a penalty immersed boundary with a forcing scheme similar to Sec. II B to make it a no-penetration and no-slip boundary. Note that the altitude of the top wall is initialized so that the spring force is already in equilibrium with the gravitational force on the fluid. In other words, the initial location of the top wall is slightly below the top boundary. Otherwise, the fluid column would undergo damped oscillations at the beginning of the simulation.

Figure 10 depicts the center of mass and the circularity index calculated from our simulation results compared against those from Featflow.³⁷ To keep track of the center of mass of the bubble, the simulation initializes the problem domain with $X6$⁠, an additional 2D equidistant grid of markers inside the bubble. $X6$ is solely for accounting purposes and does not exert any forces. The mean of the current positions of the markers $X6$ corresponds to the bubble's center of mass. The circularity index is defined in the work by Turek and Mierka,³⁷ which is the perimeter of an area-equivalent circle divided by the perimeter of the bubble. Our simulation does not calculate an area-equivalent circle but instead uses the initial circle circumference as a proxy, since in this test, the area is well conserved. The center of mass and circularity index calculated from our simulations show good agreement with the Featflow benchmark results. Figure 11(a) shows excellent agreement in the position and shape of our simulated bubble vs the benchmark results. Figure 11(b) (Multimedia view) provides more details on the rising bubble.

In addition, Table II shows the order of accuracy with regard to N for this test. The finest simulation (N = 162, dt = 0.0016) is used as the ground truth, since in this case, we do not have an analytic solution. The fluid velocity (u) error is computed in the L² norm. The table shows that the velocity field converges at a second order rate and the center of mass position converges with a rate of approximately 3/2.

In this section, we describe the results of several tests in which we vary the time step and the spatial resolution to empirically verify convergence of our methods. For all tests in this section, the height of the computational domain is L = 2 cm, the surface tension coefficient is σ = 50 g cm/s², the gas density is $ρg=0.1$ g/cm², the liquid density is $ρl=1$ g/cm², and the viscosity is $μ=0.01$ g/s. The size of the computational domain is 2 × 2 cm², and the size of the problem domain is 1 × 2 cm². Therefore, $x∈[0,L/2], y∈[0,L]$⁠, where L = 2. Periodic boundary conditions are used for the top and bottom boundaries, and a symmetry boundary condition is used for the left and right boundaries (see Sec. II A).

For this test, we simulate a droplet sliding down a vertical wall as shown in Fig. 12 (Multimedia view). The Reynolds number is Re $=200$⁠, where the typical length scale is the diameter of the droplet, and the characteristic velocity is selected as the downward velocity of the droplet at its maximum toward the end of the simulation.

There is a vertical wall at x = 0. The static contact angle is $θs=0.3 π$⁠. There is a uniform upward flow of gas, which corresponds to a prescribed vertical component of the velocity at y = 0. More precisely, the boundary condition at y = 0 (i.e., top and bottom boundary) for the vertical velocity component is $30tanh(4x/L)$ cm/s for $x∈[0,L/2]$⁠.

We vary the time step dt and the spatial resolution N. The sliding droplet shape at time t = 0.09 s is shown. The results appear to converge if we increase the spatial resolution N to 256.

In test 2, we simulate two droplets coalescing without gravity. Each droplet is given a uniform initial velocity field, one upward and one downward, to simulate their coalescence. Results are shown in Fig. 13(a). The coalescence of two droplets is an extremely challenging test case since the moment of droplet merging amplifies previous numerical perturbations. We increase the spatial and temporal resolutions (holding other parameters constant), and the simulation results show good spatial and temporal convergence with our surface tension model and interface splicing procedures. Such observations can be more vividly seen in Fig. 13(b) (Multimedia view).

Figure 14(a) shows our simulation results for coalescing six droplets for three levels of spatial and temporal resolutions. An artificial radially symmetric, inward pointing body force for gravity (G = 900 cm/s²) pulls the droplets together. This merging test is difficult since the variance of the system's response to the first coalescence event is amplified by the successive merging events. More details are given in Fig. 14(b) (Multimedia view) for this test.

Here we initialize liquid droplets in the shape of the letters “IB” and let them fall into an elastic pouch. Results are shown in Fig. 15(a). The Reynolds number Re $=8000$⁠, where the typical length scale is 2 cm, the width of the letter “I.” At the beginning of the simulation, before the droplets interact with the pouch, they slightly change shape because of surface tension. As the droplets hit the pouch, the pouch deforms and the letter “B” changes topologically as the gas bubbles exit the droplet. This test demonstrates the capabilities of our method in an all-in-one setting. Note that here the solid membrane is dynamic, rendering static boundary conditions unusable. This shows that our methods are capable of simulating a moving contact point between a liquid–gas interface and a changing-shape solid surface, simulating fluid dynamics, surface tension, Young's force, moving contact line, and solid elasticity all at the same time. See Fig. 15(b) (Multimedia view) for a video of this test.

In this section, we describe several applications of our methods. In Sec. V A, we consider the effect of droplet size on its dynamics as it slides down a wall. The importance of our interface re-sampling procedure is explored in Sec. V B. Finally, in Sec. V C, we showcase some additional tests in which interfaces are spliced to represent topological changes. All tests in Sec. V share the setup for the computational and problem domains found at the beginning of Sec. IV C 1, as well as the following parameters. The surface tension coefficient is σ = 50 g cm/s², and the gravitational acceleration is G = 980 cm/s². The gas density is $ρg=0.1$ g/cm², and the liquid density is $ρl=1$ g/cm². The viscosity for both the liquid and gas phases is $μ=0.01$ g/s. The static contact angle is $θs=0.3 π$⁠.

In this section, we consider droplets sliding down a wall, where gravity competes against an upward air flow.

In Fig. 16(a), we analyze the dynamics of various sized droplets. The smallest droplet corresponding to the top panel of Fig. 16(a) does not slide; rather, it obeys the “no-slip” rule. This is consistent with reality and is implemented via the static-dynamic friction condition. Our method tracks the wall markers that are vertically displaced by the MCL model at every time step. We observe that vertical displacements only occur near the contact point (usually within 6–8 cells). This means that although our MCL model applies to the entire wall, the majority of the wall that interacts with the droplet behaves as a no-slip boundary, and only the contact region slips. This is consistent with previous literature, which suggests that the fluid's internal viscous stress is not enough to incur apparent slip in subsonic flows with hydrophilic surfaces.³⁸ The Young's force is the dominating tangential force near the contact point. Such observations can be more vividly seen in Fig. 16(b) (Multimedia view).

In Fig. 17(a), we initialize two droplets of different sizes on the wall. The larger droplet slides down faster than the smaller one. Therefore, it catches the smaller droplet and coalesces with it. This test case shows the system's response to two different droplet sizes in one simulation and, at the same time, demonstrates the capabilities of our interface splicing capabilities at the wall. More details are given in Fig. 17(b) (Multimedia view) for this experiment.

Figure 18(a) shows a comparison study on the effect of our stepwise re-sampling technique. In the upper row, we do not re-sample, resulting in badly distributed Lagrangian markers and, eventually, topological errors. In the lower row, the re-sampling procedure is used and resolves this issue. See Fig. 18(b) (Multimedia view) for a side-by-side comparison video.

Figure 19(a) shows results from a simulation in which a single droplet separates into two smaller ones. After the splicing moment, surface tension quickly brings the sharp edge into the body of the droplet. See Fig. 19(b) (Multimedia view) for a video of this experiment. In Fig. 20(a), two droplets coalesce on the wall. This is one of the six basic cases of interface splicing introduced earlier in Sec. II G. See Fig. 20(b) (Multimedia view) for a video of this experiment. Additional simulations can be found in Figs. 21(a)–21(c) (Multimedia views). In Fig. 21(a), a droplet on the wall separates into two droplets as a results of a collision with the wall. A droplet attaches to the wall in Fig. 21(b). Finally, in Fig. 21(c), an air bubble collides with a droplet, showing the effect of variable density.

In this paper, we use the immersed boundary method to simulate droplets and their interaction with a wall. The static-dynamic friction is modeled with an immersed boundary. The surface tension and the Young's force were readily integrated with the IB method by computing local unit vectors tangent to the interface. We showed that interface marker re-sampling is necessary to obtain accurate results. Our interface splicing method handles droplet coalescence, separation and other topological changes with stable splicing events. We also demonstrated empirical convergence of our methods on non-trivial examples, and we applied them to several benchmark cases, including a slipping droplet on a wall and a rising bubble.

The success of our methods demonstrates the flexibility and extensibility of the IB method. We show that complex dynamics, like a moving contact point on an evolving solid surface, can be successfully modeled with the IB method. In the formulation of the static-dynamic friction condition, we choose to implement contact point slip via length-limiting the tether springs so that IB imparts the correct friction forces onto the fluid. Conceivably there are other ways to impart this frictional force, which we did not investigate in this study. For example, one could formulate the friction as an additional external force that is imparted directly onto the fluid without intermediate markers, or one could prescribe a boundary condition for the fluid velocity according to viscous stress.

Our work is limited in several ways. The methods presented here assume that the fluids are incompressible. This is a result of the immersed boundary approach used here in which the entire domain is modeled as an incompressible fluid. We expect the interface splicing method to be fragile in the compressible-flow setting because inertia and compressibility will likely make splicing events much less atomic. In addition, although our methods can handle variable density, they do not handle variable viscosity. Finally, the methods are limited to 2D. In 3D, the interface surface becomes a surface mesh (most commonly, a triangle mesh). Care needs to be taken in computing surface tension in this case. Furthermore, the stepwise interface re-sampling procedure, in 3D, needs to also repair local topology of the mesh while preserving accuracy and numerical stability. Therefore, extending our methods to 3D will be a nontrivial task.

An important application of our work is a real-life industrial problem corresponding to the removal of gaseous reactants from power plant exhaust by converting them into liquid at the catalyst site. Future work will be geared toward understanding the dynamics of liquid droplets for this application, in order to maintain performance and durability of the catalyst assembly. More specifically, these types of models will provide details for classifying droplet movement and coalescence in the gas channel based on the filter surface and gas flow properties, which will be used to guide operating conditions for sulfur dioxide removal devices. Our goal is to develop such a mathematical model, which will help us to understand how different surface properties and operating conditions affect the dynamics of the droplet or the formation of a liquid film.

C.P. and P.S gratefully acknowledge support from the National Science Foundation (NSF) under Grants No. RTG/DMS-1646339. P.S. is also supported by an Institutional Support of Research and Creativity (ISRC) grant provided by New York Institute of Technology and financial support from NSF under Grants No. DMS-2108161. We thank Otto Mierka and Stefan Turek for providing the 10:10 viscosity benchmark data for rising bubble test case 1. Several very helpful conversations with Guanhua Sun are gratefully acknowledged.

The authors have no conflicts to disclose.

M.Y.L. and D.C. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request. All simulations in this study can be reproduced with the source code openly available at repository https://github.com/daniel-chin/droplet. The repository also contains the scripts to plot all the figures.