A fully Eulerian hybrid Immersed Boundary-Phase Field Model for contact line dynamics on complex geometries

We present a fully Eulerian hybrid immersed-boundary/phase-field model to simulate wetting and contact line motion over any arbitrary geometry. The solid wall is described with a volume-penalisation ghost-cell immersed boundary whereas the interface between the two fluids by a diffuse-interface method. The contact line motion on the complex wall is prescribed via slip velocity in the momentum equation and static/dynamic contact angle condition for the order parameter of the Cahn-Hilliard model. This combination requires accurate computations of the normal and tangential gradients of the scalar order parameter and of the components of the velocity. However, the present algorithm requires the computation of averaging weights and other geometrical variables as a preprocessing step. Several validation tests are reported in the manuscript, together with 2D simulations of a droplet spreading over a sinusoidal wall with different contact angles and slip length and a spherical droplet spreading over a sphere, showing that the proposed algorithm is capable to deal with the three-phase contact line motion over any complex wall. The Eulerian feature of the algorithm facilitates the implementation and provides a straight-forward and potentially highly scalable parallelisation. The employed parallelisation of the underlying Navier-Stokes solver can be efficiently used for the multiphase part as well. The procedure proposed here can be directly employed to impose any types of boundary conditions (Neumann, Dirichlet and mixed) for any field variable evolving over a complex geometry, modelled with an immersed-boundary approach (for instance, modelling deformable biological membranes, red blood cells, solidification, evaporation and boiling, to name a few)

as well. The procedure proposed here can be directly employed to impose any types of boundary conditions (Neumann, Dirichlet and mixed) for any field variable evolving over a complex geometry, modelled with an immersed-boundary approach (for instance, modelling deformable biological membranes, red blood cells, solidification, evaporation and boiling, to name a few).

Symbols
Definitions units in (SI)

AP1
Averaging points for the first interpolation point −

AP2
Averaging points for the second interpolation point −

Introduction
Motion of a three-phase contact line occurs in a variety of industrial fields from coating to energy conversion processes, nucleate boiling, droplet dynamics, two-phase flow in porous media, and microelectronics cooling, to name a few (Sui and Spelt (2013a,b); Yarin (2006)). Despite numerous studies have been performed on the contact line motion, the underlying physics is still a matter of debate. The difficulty in studying the contact line movement originates in the so-called "contact line singularity" which was first discussed by Moffatt (1964) and Huh and Scriven (1971). These authors showed that the fluid flow, close to the contact line, is in the Stokes regime and exhibits singularities in both the shear stress and the pressure (Krechetnikov, 2019). In general, three main solutions have been proposed to remove the singularities close to the moving contact line. As the first solution, a slip velocity (Navier boundary condition) can be applied at the surface near the contact line (Dussan, 1979). Modelling the dynamic contact angle and the formation of a precursor film are the two other well-known solutions (Sui and Spelt (2013a)).
In order to model the moving contact line problem, different methods have been implemented for tracking the interface and reconstructing the contact line. (Renardy et al., 2001) used a Volume of Fluid (VOF) method to model the contact line problem. A piecewise linear interface construction scheme was used to reconstruct the interface based on an indicator function. Afkhami et al. (2009) presented a mesh-dependent contact angle model to remove the stress singularity at the contact line. Mukherjee and Kandlikar (2007) proposed a Level-Set approach to study bubble growth during an ebullition cycle. Spelt (2005) used an extended Level-set method to simulate multiple moving contact lines. The model accounts for flow inertia, contact line slip velocity and contact-line hysteresis. A front-tracking method was used by Muradoglu and Tasoglu (2010) to model the impact and spreading of viscous droplets on solid walls. Izbassarov and Muradoglu (2016) studied the effects of viscoelasticity on drop impact and spreading on a flat solid surface using a front-tracking method together with finitely extensible nonlinear elastic-Chilcott-Rallison model for fluid elasticity.
They showed that during the spreading phase, viscoelastic effects increase the spreading. ? studied the role of surface wetting on interface instability and penetration modes in a porous medium consisting of two immiscible fluids. Their results suggest that the displacement patterns depend on both the capillary number and the surface wetting properties. They also examined the well-known Haines jumps (sudden interface jumps from one site to another) by analysing the characteristic time and length scales of the jumps.
Diffuse interface models have also been used extensively to study the contact line dynamics. Among different diffuse interface models, the phase-field method has drawn more attentions during the last decades. In a diffuse interface model, an order parameter (concentration) is defined to distinguish between different phases. The interface has a finite thickness within which the order parameter and the fluid properties vary smoothly (but significantly) from one phase to another one. This assumption allows tracking the interface by solving an advection-diffusion equation for the order parameter. In this context, the Allen-Cahn model (Allen and Cahn, 1979) is a reaction-diffusion phase-field equation which has been used to study phase separation in multi-component systems.
Among others, Ben Said et al. (2014) studied the equilibrium wetting of a system of immiscible fluids on a flat substrate using the Allen-Cahn model. As concerns dynamic wetting, however, the Cahn-Hilliard model (see Cahn, 1961) has been shown to be more reliable to remove the contact line singularity: imposing a dynamic contact angle boundary condition together with slip velocity at the wall is straight-forward in the Cahn-Hilliard formulation. Moreover, the  (Peskin, 2002), a powerful tool to model fluid flows over complex geometries. This approach has been extensively used to simulate fluid-solid interaction problems, mainly for a single phase fluid. The idea of the immersed boundary method is to solve the system of equations on a cartesian numerical mesh (regardless of the solid geometry), and imposing the boundary conditions by adding forces at specific grid points close to the boundary.
During the last decades, immersed boundary methods have been adopted with different interface tracking approaches, front-tracking (Deen et al., 2009), volume of fluid (Patel et al., 2017), level set (Wang and Desjardins, 2018), phase field (Liu and Ding, 2015;Nishida et al., 2018), with the aim to simulate the interaction between a multiphase fluid flow and a solid boundary. However, to the best of our knowledge, a fully Eulerian numerical approach for modelling the dynamic motion of a three phase contact line together with contact line slip velocity over any arbitrary geometry has not been reported yet. Since both the immersed boundary method and the phase field model proposed here are fully Eulerian, parallelisation of the algorithm is straight-forward and potentially highly scalable. The existing parallelisation of the underlying Navier-Stokes solver can be used for the multiphase part as well. increases the stability of the method, and reduces the numerical error (Dong and Shen, 2012;Shen et al., 2015;Yu and Yang, 2017;Huang et al., 2020). Here, we report numerical validations to show that the proposed algorithm provides accurate results for various test cases, for both choices of time integration (explicit and semi-implicit).
The outline of the manuscript is as follows. In section 2, we explain the main concept behind the phase field model together with the corresponding boundary conditions. In section 3, we summarise the implemented numerical schemes, whereas we elaborate on the proposed hybrid phase field-immersed boundary model in section 4. To validate and test our implementation, we report results from several numerical tests in section 5. We conclude our work in section 6, and finally present a more accurate semi-implicit version of the algorithm in Appendix A.

Phase Field Model
Eulerian interface tracking approaches can be divided into two main groups, namely, sharp interface methods and diffuse interface methods. In diffuse interface methods, the interface is assumed to have a finite thickness. Although the interface is much thicker than the real physical one, this assumption pro-vides resolvable properties which vary continuously within the interface. Such a continuum model avoids any requirement for jump conditions at the interface or interface reconstruction. Moreover, fluid properties are conserved within the interface. During the last decades, the Phase Field Model (PFM) has become more and more popular in the multiphase flow community for these properties.
It originates from Van der Waals model for free energy (van der Waals, 1979) where the bulk free energy and the interfacial free energy are added to give the total free energy per unit volume, f , of a system of two immiscible fluids as follows: where C is the concentration (order parameter) which distinguishes different phases; it varies from −1 in one phase to +1 in the other. The variation of C through the interface is smooth but rapid. ψ is a double-well function, ψ = (C 2 − 1) 2 /4, with two minima for each stable phase. σ and denote the surface tension coefficient between the two phases and the interface thickness, respectively. The first term in equation (1) represents the contribution of the interfacial energy, whereas the second one models the bulk free energy density (Jacqmin, 1999).
Considering the requirement of minimum free energy in the equilibrium state, and defining the chemical potential φ as the variation of the total free energy (F = f dV ) with respect to the concentration, i.e., φ = ∂F /∂C, Cahn and Hilliard proposed an equation for the evolution of the concentration where the motion of a diffuse interface within a binary fluid is modelled by the so-called Cahn-Hilliard equation (Cahn and Hilliard, 1958;Cahn, 1961): where u i , and M represent the fluid velocity vector and the mobility coefficient.
The Cahn-Hilliard equation is an advection-diffusion equation which, in the limit of zero diffusivity, reduces to a sharp interface model. Plenty of studies have considered the sharp interface limit, obtaining an appropriate range of values for the mobility coefficient based on the interface thickness (Magaletti et al., 2013;Xu et al., 2018).
The difference in the chemical potential between the two phases at the interface is the mechanism driving the interface motion (besides the advection of the interface by the mean flow). Theoretically, the chemical potential is the variation of the free energy with respect to the concentration and can be calculated with the following equation: To couple the Cahn-Hilliard equation (2) with the fluid flow, a term is added to the right hand side of the Navier-Stokes equation which accommodates for the surface tension forces at the interface (Jacqmin, 1999).
where ρ and µ are the density and dynamic viscosity of the fluid, varying from ρ 1 and µ 1 in one phase to ρ 2 and µ 2 in the other one, defined as and µ(C) = [(C + 1)µ 2 − (C − 1)µ 1 ] /2; P is the pressure and φ∂C/∂x i represents the surface tension force at the interface. The second equation represents mass conservation for incompressible fluids.
In the presence of a solid substrate, a third term is added to eq. (1), the contribution of the solid substrate to the total free energy of the system (Carlson et al., 2011): where σ sg and σ sl are the surface tension coefficients between solid-gas and solid-liquid, respectively.
!" # $! $" Figure 1: Sketch of contact angle. According to Young's equation (Young, 1805), the equilibrium contact angle can be expressed as function of the surface tension coefficients between the different phases.
According to Young's equation (Young, 1805) , the equilibrium contact angle θ e depends on the surface tension coefficient between each pair of the three phases, σ, σ sl , and σ sg (see figure 1) as The boundary condition for the concentration necessary to impose a prescribed contact angle is therefore obtained by minimising the energy at the solid wall (Jacqmin, 1999(Jacqmin, , 2000, where µ f , θ eq , and n i are the contact line friction coefficient, the equilibrium contact angle, and the vector normal to the solid surface. The first boundary condition in eq. (6) models the dynamics of the contact line motion, where the contact line friction coefficient is inversely proportional to the time needed by the contact line to relax to its prescribed static contact angle (Carlson et al., 2012;Xu et al., 2018). In this equation, g(C) = (2 + 3C − C 3 )/4 is a function which varies smoothly between zero and 1 from one stable phase to the other.
The second boundary condition in equation (6) guarantees impermeability at the wall.
Finally, in the case of non-zero wall slip velocity, the following equation can be solved together with the other boundary conditions to obtain the slip velocity at the wall (Carlson et al., 2012): where u s , l s , and t j are the slip velocity, slip length, and the unit vector tangent to the surface.

Numerical method
The full system of equations introduced above can be summarised as follows: where f i represents the immersed boundary force used to account for the complex wall geometry, explained in section 3.2, and f bi indicates all the other body forces (e.g. the gravitational force). The complete boundary conditions at the wall are As mentioned above, in the following we will introduce the numerical algorithm assuming a fully explicit approach. However, this set of equations can also be solved semi-implicitly as discussed in Appendix A. Note that we use both the algorithms alternatively in the numerical tests discussed later on, the differences between the results being negligible once the step is chosen correctly. However, the semi-implicit algorithm allows, on average, a 10 times larger time step.

Time integration and spatial discretisation
We solve the system of equations (8 and 9) on a Cartesian mesh with a staggered arrangement, where the velocity components are defined at the faces and the pressure, the chemical potential and the order parameter are defined at the cell centres. The second-order finite difference scheme is used for spatial discretisation and the different terms are advanced in time explicitly using the Adams-Bashforth scheme. Finally, the fractional-step method for incompressible two-fluid systems is implemented as in Dodd and Ferrante (2014). The baseline solver has been extensively validated in the previous works (see among others Rosti et al., 2019;Francesco De Vita, 2020;.

Immersed boundary method
During the last decades, a variety of immersed boundary methods (IBM) have been used for modelling fluid-solid interactions with moving and fixed bodies (Mittal and Iaccarino (2005)). In most of the IBM formulations, the solid boundary is represented by a set of Lagrangian points whose locations are tracked by solving an extra set of equations. An auxiliary force is added to the mesh cells surrounding each of the Lagrangian points to impose the no-slip and no-penetration conditions at the solid boundaries. However, there are also fully Eulerian IBM formulations where the immersed boundary forces are computed directly on the numerical grid points, especially for the case of fixed objects.
In this paper, we use a simple Eulerian IBM formulation, namely, the volume penalisation model proposed by Kajishima et al. (2001) where u * * ijk is the second prediction velocity and u sol is the solid wall velocity within the corresponding grid cell. In the case of stationary wall (u sol = 0), equation (10) reduces to the following simpler form,

The hybrid PFM-IBM Algorithm
To impose the boundary conditions for the order parameter and slip velocity at the wall, we need to calculate the values of the different quantities appearing in equations (9) (the order parameter, velocity components, derivatives of the order parameter, etc.) at the wall, separating tangential and normal components with respect to the wall surface. To accomplish this, we need to pre-compute some different quantities, e.g. solid volume fraction, normal vector, averaging weights, etc. Having performed this initialisation step, it is possible to integrate in time the governing equations with general boundary conditions at the wall.
We now proceed to provide a numerical recipe to solve the system at hand with corresponding boundary conditions for slip velocity and order parameter at the wall for a two-dimensional system. The method can be easily extended to three-dimensional problems, as shown in the result section.

Preliminary computations
In this section, we elaborate the details of the preliminary calculations needed for the proposed algorithm (green dashed box in figure 5).

Volume fraction and normal vectors
In order to impose all the boundary conditions introduced in equations 6 and 7 together with the volume penalisation immersed boundary method (equation 11), we should first calculate the liquid and solid volume fraction of each grid cell, α ijk . This needs to be done for all the four numerical cells on a staggeredgrid, namely, one cell-centred cell and three face-centred cells for each grid point ijk.
We propose a simple approach for computing the volume fractions. Let us consider a cell-centred cell at the grid point (i, j, k) and, first, divide the cell into a sufficient number of subgrid points in both directions as shown in figure 2 (our tests suggest that 100 points in each direction are enough). Under the assumption that the coordinates of each subgrid point in the cell are known, we can determine whether the point is inside or outside the solid wall. Hence, the volume fraction α is simply the ratio of the number of subgrid points outside the solid to the total number in the cell. Having calculated the volume fraction on each cell, the solid wall normal vectors are approximated with the gradient of the volume fraction. To this aim, we first compute the derivative of the volume fraction at the four cell corners (Ii et al., 2012), e.g., The value of the derivatives can be calculated in the three other corners in a similar fashion and then be used altogether to compute the normal vector at the cell center as where 0 is a very small positive number used to avoid division by zero and

Ghost points, intersections, mirror and interpolation points
Next, we need to identify four groups of points: the ghost, intersect, mirror, and interpolation points. In this section we use the term temporary array to refer to lists that are defined and used only in the initialisation steps and can be deallocated later on. On the other hand, the term permanent array refers to lists that contain variables required during the whole simulation. The coordinates are saved in a three dimensional temporary array whose index Finally, we identify two additional sets of points, the interpolation points denoted IP 1 and IP 2 . These points will be used to extrapolate the magnitude of any quantity at the mirror point as discussed later. The interpolation points are located on the same straight line as the mirror, intersect, and ghost points.
Note that the distance between the two interpolation points and between the first interpolation point and the mirror one is the same (i.e. d 2 in figure 3).
The magnitude of d 2 is chosen smaller that d 1 , as function of the Reynolds number of the problem under investigation. As explained later, the magnitude of any arbitrary parameter at the mirror point is linearly extrapolated from the magnitudes at the corresponding interpolation points. This assumptions is valid if the interpolation points and the mirror points are located inside the inner-part of a boundary layer where a linear profile can be safely assumed. In general, it is well known that as the Reynolds number increases, the boundary layer becomes thinner; thus the numerical grid should be finer to resolve the boundary layer properly (regardless of the IBM algorithm). Nevertheless, it is important to check that the interpolation points and the mirror points are located inside the inner part of the boundary layer (and if not, tune the value of d 2 ). Note, finally, that we also need to identify the cells in which the interpolation points are located and store their indices in a permanent array.

Coordinate transformation
Although we solve the governing equations using a simple Cartesian coordinate system (X 1 , X 2 for two-dimensional problems), the phase field boundary conditions are expressed in a coordinate system following the solid boundary, with normal and tangential vectors n j and t j . Therefore, we need to transform between the two systems using the calculated normal vector, n i . For an arbitrary variable Ω, subject to a coordinate transformation, we have where the suffixes n and t indicate the components in the frame following the boundary.

Inverse distance weighting averaging and linear extrapolation
Let us consider any arbitrary variable, say Ω, whose approximated value is needed at point M. As first step, we average the value of Ω at the 2 interpolation points (IP 1 and IP 2 in figure 3). To do so, we identify from the coordinates of each of the interpolation points the grid points surrounding each of them (AP 1 to AP 5 in figure 4). We then perform the averaging through an inverse distance weighting: for a two-dimensional problem, the average uses five points (AP 1 , AP 2 , AP 3 , AP 4 and AP 5 ), with weights equal to the inverse of the squared distance between the averaging points (AP i ) and the interpolation point (IP ), denoted here 1/h 2 i . The interpolated value can thus be calculated as follows: where w(i) is the weight of the i th averaging point and q is the sum of all the weights. Depending on the geometry, the interpolation point can overlap with one of the averaging points, in which case h i goes to zero and the interpolated value can be taken as the value at the averaging point. To carry out the computation of the boudnary condition, we therefore define two other permanent arrays, W 1(i, j, k, AP ) and W 2(i, j, k, AP ): W 1(i, j, k, AP ) contains the weight of each averaging point (AP 1 to AP 5) corresponding to the first interpolation point IP 1 of any ghost point G ijk , and similarly for W 2(i, j, k, AP ) containing the weights for the second interpolation point IP 2.
At the end of this preliminary phase, we are ready to use, during the simulation, the weight arrays to compute the averaged value of the variable Ω at each of the interpolation points and then use these to find the value at the mirror point by linear extrapolation

Solution algorithm
In this section we describe the solution of the system of equations and the algorithm used to impose the boundary conditions on a fixed wall of arbitrary shape.

Solving the Cahn-Hilliard equation
We first solve the Cahn-Hilliard equation to update the order parameter from time C n to C n+1 . Note that irrespective of the time integration method (explicit or semi-implicit), by solving the Cahn-Hilliard equation, we update the correct value of the order parameter in all the numerical grid points (even inside the solid) except for the ghost points. The details of the semi-implicit algorithm adopted for the Cahn-Hilliard equation are presented in appendix A. For the explicit algorithm, we use the second-order central finite difference scheme for the spatial discretisation and second order Adam-Bashforth for the temporal discretisation.

Imposing the boundary conditions for the order parameter
Of relevance here, we impose the contact angle and the no mass penetration boundary conditions at the ghost points using the IBM algorithm. According to equation (6) To impose the boundary conditions, we proceed as follows. Let us define γ as the ratio of the distance between the ghost and the intersect points to the distance between the intersect point and the mirror point (γ = d 1 /dn). The value of the order parameter, the derivatives of the order parameter, and the velocity components at the intersect point can be calculated as As already mentioned, we do not update the magnitude of the order parameter at the ghost point when solving the Cahn-Hilliard equation. Therefore, we use hereĈ G which is an estimation of the order parameter at the ghost point at time n + 1, defined asĈ G = 2C n G − C n−1 G . The different quantities are then projected from X 1 , X 2 to n, t Finally, we can compute ∂C/∂t at the wall as By integrating in time, we find the updated value of the order parameter at the wall C (n+1) I which is used to update the order parameter at the corresponding ghost point as We can now update the chemical potential and impose the corresponding boundary condition (∂φ/∂n = 0) in a similar way. Finally, the density and viscosity are updated using the new values of C.

Calculating the first and the second prediction velocities
As illustrated in figure 5, the next step is to solve the Navier-Stokes equations. Similar to the Cahn-Hilliard equation, we calculate the first and the second prediction velocities using either an explicit or semi-implicit algorithm.
The details of the semi-implicit implementation are provided in appendix A.
For the explicit algorithm, we use second order central finite differences for the spatial discretisation and we integrate all the terms in time using the second order Adam-Bashforth scheme. Note that we solve the Navier-Stokes equation for all the numerical grid points except at the ghost points, which we will use to impose the boundary conditions.
We recall that f bi is the summation of all the external body forces (such as gravity). The last equation, the step between u * i and u * * i , is the IBM penalisation discussed above, which imposes zero velocity inside the solid. In the cases with slip, the boundary condition at the wall is modified using the ghost point, so that the fluid has a slip velocity. This is detailed in the next section.

Enforcing the velocity boundary conditions
We impose the slip velocity boundary condition at the ghost point using the IBM algorithm. To impose the velocity boundary conditions, we first interpolate the calculated second prediction velocity u * * at the cell centres. Next, we use the interpolation scheme introduced in section 4.1.4 to calculate the magnitude of any cell-centred velocity component at the mirror points. Due to the no penetration boundary condition at the wall, the normal component of the velocity at the wall is equal to zero. Therefore, the normal component of the velocity at the ghost point can be updated as The boundary conditions for the tangential component of the velocity at the wall can be obtained from equation (7) µ l s V * * ts = µ After discretizing the equation, we can obtain the tangential velocity at the ghost point by solving the following equation: By transforming back from boundary-fitted (n, t) to the cartesian coordinates (X 1 , X 2 ), the updated value ofũ * * i at the ghost points can be found Finally, the velocity components at the corresponding faces are found by interpolating the cell-centred values.

Correction step
For the current implementation, we follow the approach in Dodd and Ferrante (2014) to correct the calculated second prediction velocity and satisfy the divergence free condition when the density is not uniform. First, we update the pressure field by solving the following equation: where ρ 0 = min(ρ 1 , ρ 2 ) andP = 2P n − P n− Having updated the pressure field, we correct the second prediction velocity and calculate the divergence-free velocity as follows: Details of the algorithm can be found in the above reference. The algorithm proposed here, with the different steps, is summarised in figure 5.

Numerical tests
To

Droplet spreading on a flat plate
A two-dimensional circular droplet is initially placed just above a flat wall so that the interface is almost tangent to the wall. Nakamura et al. (2013) assumed that the droplet is so small that the gravitational forces are negligible compared to the surface tension forces. The Reynolds (Re), capillary (Ca), and Cahn numbers (Cn) are defined with reference to the initial droplet radius (R), the reference velocity (U ref = σ/µ), the density of the liquid phase (ρ L ), the surface tension coefficient (σ), and the interface thickness ( ). The density and viscosity ratio between the two phases are equal to 100 and Note that throughout the paper we use the same definition for all the nondimensional parameters as in equation 30. In this simulation, the wall friction coefficient µ f is set to zero. The computational domain is [0, 10R] × [0, 10R] and the number of grid points per droplet diameter is equal to 320. Figure 6 shows the evolution of the normalised wetting radius versus the non-dimensional time (t * = t/(ρR 3 /σ) (1/2) ) for two different contact angles (θ eq = 45 • and 135 • ) and As shown in the figure, the results of our simulations are in good agreement with those by Nakamura et al. (2013).

Two phase Couette flow
As the second validation for the PFM module, a two phase Couette flow is simulated and the results compared with those by Bao et al. (2012). The simulation domain is [0, 1] × [0, 0.25]. Two interfaces are initially located at x = 0.25 and x = 0.75 with a tangent hyperbolic variation of the order parameter from one phase to the other according to: In this simulation, the wall velocity is u w = 0.2, Re = 5 (defined as above We report in figure 7 the contour of the order parameter at t * = 1.875 and, for an easier comparison, the digitised interface location from the simulation by Bao et al. (2012) in a second panel. The blue lines represent the results of our simulation and the red dots those by Bao et al. (2012), with a good match between the two data sets. C(X 1 , X 2 , t = 0) = 0.5 + 0.12 cos(2πX 1 ) cos(2πX 2 ) + 0.2 cos(πX 1 ) cos(3πX 2 ) Two cases are simulated: in the first case, the phase separation occurs in the absence of any velocity field, while in the second one, the initial velocity field is initialised as u 1 (X 1 , X 2 , t = 0) = −sin 2 (πX 1 ) sin(2πX 2 ), u 2 (X 1 , X 2 , t = 0) = sin 2 (πX 2 ) sin(2πX 1 ).
(34) The total energy of the system is calculated and compared with the one obtained by Nishida et al. (2018). This is defined as the summation of the kinetic and interfacial energy The evolution of the total energy is depicted in figure 9 for

Droplet spreading on a rotated wall without gravity
To validate our hybrid code for different slip lengths, we reproduce once more the results of a droplet spreading on a flat wall by Nakamura et al. (2013) but   x 3 s z z M w L E i k M u u 6 3 s 7 K 6 t r 6 x W d g q b u / s 7 u 2 X D g 6 b J k 4 1 4 w 0 W y 1 i 3 A 2 q 4 F I o 3 U K D k 7 U R z G g W S t 4 L R z d R v P X J t R K w e c J x w P 6 I D J U L B K F r p / q l X 7 Z X K b s W d g S w T L y d l y F H v l b 6 6 / Z i l E V f I J D W m 4 7 k J + h n V K J j k k 2 I 3 N T y h As shown in figure 13, the results of the hybrid PFM-IBM code are in good agreement with the reference data.

Droplet spreading over sinusoidal surfaces
To examine the capability of the developed model to simulate more complex geometries, we simulate a droplet spreading over two sinusoidal surfaces with same wave length but a phase shift of π/2 with respect to the initial droplet position. The wall geometry is defined by  In this case, we define an equivalent spreading radius as the horizontal distance between two contact lines, reported in figure 15 normalised with the initial radius of the droplet (for m = 0 in fig 15a and m = 1 in fig 15b). The red and the blue colours indicate the cases with slip and no-slip velocity boundary conditions, respectively. The solid line denotes the results for the hydrophilic wall and the dashed-line those for the hydrophobic wall. Figures 14 and 15 show that the wall geometry makes a significant difference in the droplet spreading. When the droplet is initially placed on a cavity (see fig 14c and fig 14d), the trapped gas inside the cavity cannot leave it; hence, to conserve the mass of the trapped gas, the droplet wets the surface by forming an arc shape. As a consequence the static contact angle is reached faster and the spreading is limited. Therefore, a static configuration is achieved faster without significant spreading. On the other hand, in the absence of any trapped gas (see fig 14a and fig 14b), the droplet fully wets the wall and spreads over the surface. The spreading continues until the equilibrium contact angle is attained; however, given the computational cost, we stopped the simulations at t * ≈ 47.
In addition, we note that a slip velocity remarkably speeds up the spreading when the droplet fully wets the wall (compare fig 14a and fig 14b). Finally, as shown in figure 15, the spreading radius on a hydrophobic wall is much less than that on its hydrophilic counterpart. To examine the convergence of the algorithm, we performed the same simulations as in figure 14b for three additional grid sizes, namely, 2560 × 1280, 1600 × 800, and 960 × 480. For all of the simulations, we calculate the evolution of the wetting radius. We consider the results of the finest grid size (2560×1280) as the reference values (R ref ) and calculate the normalised error obtained with the coarser numerical grids, see gigure 17a. Note that to reduce the computational costs, we performed the simulations up to t * = 2. According to figure   17a, for the grid sizes equal to or smaller than 122/640, the results of the simulations are almost independent of the grid size. It is worth to mention that due to the combinations of the different parameters affecting the necessary grid size (geometry, velocity boundary condition, contact angle boundary condition, etc.), a grid study should be performed for each specific problem under study.
Mass leakage is a well-known problem of any phase-field model (Huang et al., 2020). Different numerical methods have been proposed for solving the Cahn-Hilliard equation and to reduce the droplet shrinkage. In figure 17b, we report the evolution of the relative change in the volume of the droplet for a case with slip velocity boundary condition (l s = 0.25R) at the sinusoidal wall without any phase shift (m = 0) obtained with 1280 × 640 grid points. Our result illustrates that up to t * = 12, the mass loss is less than 0.3% which shows that the proposed IBM algorithm does not introduce additional mass leakage.

Three-dimensional droplet spreading
As mentioned in section 4, the algorithm presented above can be extended to three-dimensional formulations. In this section, we present the results of a simulation performed to model the spreading of a three-dimensional initially spherical droplet over a three-dimensional surface. The simulation is performed using the semi-implicit algorithm with a time step equal to 10 −4 .
A droplet of radius R is initially placed at [X 1 , X 2 , The solid wall is generated as a portion of a sphere with radius R w = 9R and centre located at [X 1 , X 2 ,

Conclusion
We have presented a fully Eulerian hybrid immersed-boundary phase-field model for simulating contact line dynamics on any arbitrary fixed solid wall. The algorithm consists of two independent modules, namely, a phase field and an immersed boundary module. The Navier-Stokes and Cahn-Hilliard equations are solved on a cartesian numerical mesh yet imposing boundary conditions on any complex wall geometries using a volume-penalisation and ghost-cell immersed boundary method. The proposed algorithm is capable modelling both static and dynamic-contact angle boundary conditions with possibly slip velocity at the wall. The proposed algorithm has the following novel properties: • The fully Eulerian approach facilitates an efficient implementation and, in particular, parallelisation and accelerated architectures.
• It consists of an initialisation step during which all the auxiliary quantities necessary for the immersed boundary treatment of the complex wall are calculated and stored (coordinates of ghost points, normal vectors at the wall, coordinates of averaging points, and averaging weights). Hence, the IBM module does not add significant computational cost to the base solver for the system formed by the Navier-Stokes and Cahn-Hilliard equations.
• It suits both two-and three-dimensional simulations, without additional complexities in three-dimensions.
• Due to the modular feature of the algorithm, the phase field formulation, particularly the free energy of the system, can be modified with no effect on the solution algorithm, Thus the same strategy presented here can be employed to model different near-wall physical phenomena (for instance solidification over a complex wall geometry).
The numerical tests reported in this manuscript validate the algorithm against different results from the literature on wetting and two-fluid systems.Hence, the proposed algorithm can be see as an efficient and powerful method to study multiphase flows near solid boundaries using free-energy formulations, in particular contact line dynamics over complex geometries. However, the IBM module presented here, with the accurate and efficient calculations of normal and tangential vectors, interpolation and extrapolation through an immersed boundary, can be used to impose arbitrary mixed boundary conditions for flow problems over complex geometries, not only for single phase flows but also multiphase flow simulations using other Eulerian approaches to track an interface, e.g. volume-of-fluid and level-set methods. Examples of simulations where the approach proposed here could be of help are large-eddy and RANS simulations where wall models are necessary (Roman et al., 2009;Bhattacharya et al., 2008) and heat-transfer problems where boundary conditions involve both the temperature and concentration field as well as their gradients (Lupo et al., 2019).

Appendix A. Semi-implicit algorithm
The semi-implicit algorithm follows the same steps as the fully explicit one but with different numerical procedure for the Cahn-Hilliard and Navier-Stokes equations. In this appendix the semi-implicit algorithms are presented as these turn out to be more stable and to allow for a longer time step, so they are preferable for an efficient implementation when dealing with larger problems.

(A.2)
S is the stabilisation parameter and is chosen such that S ≥ 2 6 λM ∆t . The coefficient α s of the second Helmhotz equation is also computed following Dong and Shen (2012), i.e. and update the value of the order parameter at time n + 1.
Here, we solve both equations A.1 and A.4 by taking Fourier transforms.
Note that the boundary conditions for the order parameter and the velocity at the wall are imposed through the IBM algorithm. Therefore, in the wall-normal direction X 2 direction, we simply consider Neumann boundary conditions for both Γ and C together with the no-slip boundary condition for the velocity, whereas periodic boundary conditions are considered in the flow direction, X 1 , for all the variables.