A Robust Method for Wetting Phenomena Within Smoothed Particle Hydrodynamics

This work is dedicated to improved modelling of wetting phenomena and contact angles by means of meshless Smoothed Particle Hydrodynamics (SPH) method. For this purpose we modify the surface tension model in the neighbourhood of the triple line. For higher accuracy we also alter the general method of computing the interface normals, by applying a specific smoothing to the field of vectors, and not the phase indicator function. The necessity of proper boundary conditions in SPH is concisely discussed with multiphase systems in mind. The implemented model yields accurate results for a set of static cases, even for very low values of the imposed contact angle, both in the absence and presence of gravity. The results of the calculations of droplets moving along the wall show a great potential of the proposed approach in future applications.


Introduction
Wetting phenomenon is one of the common features of multiphase flows observed in daily life. For sure anyone has encountered it in the form of rain droplets lazily sliding down the windowpanes or freshly waxed car. From industrial point of view, wetting is of importance when it comes to coating of surfaces to obtain desired hydrophilic or hydrophobic properties. For example, icing of vehicles, especially airplanes, can be remedied by using hydrophobic surfaces [1] and the hydrophilic ones can be used for the sake of anti-fogging [2]. Possible area of applications is extremely vast [3,4], hence the necessity of researching and investigating the wetting phenomena. Despite its superficial simplicity the physics of wetting is complex. The macroscopic effect observed by human eye is the result of microscopic forces resulting from the triumvirate of solid, gas and liquid phases. The understanding and control of this phenomenon has been an object of research for years; a comprehensive overview can be found in [5].
Experimental studies are an obvious source of information on the physics of this peculiar multiphase flow. Data in the form of correlations are of importance in designing and researching new materials and their properties [6][7][8]. On the other hand, numerical techniques are available as a cheaper and more flexible alternative in terms of setting physical parameters, etc. The grid based Eulerian approaches in computational fluid dynamics (CFD) are considered to be the state of the art tools and are widely applied to the modelling of interfacial flows [9]. The models of wetting were successfully applied in the most popular approaches, i.e. Volume of Fluid (VOF) [10], Level-Set [11] and Phase Field Model [12]. In order to impose desired properties of the solid walls, often characterised by the contact angle, the basic numerical schemes require proper alteration of the surface tension model and/or modification of the vectors normal to the interface at the triple point. These methods are quite mature and can be successfully applied to very challenging (from numerics point of view) setups like wetting of the rough surface [13].
In the last 20 years the Lagrangian meshless methods appeared as an alternative to the Eulerian ones, with Smoothed Particle Hydrodynamics (SPH) as the most popular and relatively mature at the moment [14]. Why even consider "alternative" approaches when well-known, established methods are at hand? No need of grid generation is very convenient and helpful when dealing with problems involving significant deformations and topology changes, with interfacial flows being a perfect example of such a case [15,16]. Furthermore, SPH deals easily with high density and viscosity ratios and handles interfaces in a fairly natural and straightforward way. Since the continuum is discretised as a set of particles that can be thought of as Lagrangian fluid elements, the separate phases are represented by separate sets of particles. The position of the interface is defined by particles' positions and no reconstruction or tracking algorithm is required. This makes SPH an attractive approach when it comes to modelling of interfacial flows. What is more, SPH is a very flexible tooldifferent physical models can be easily incorporated and coupled, e.g. heat transfer [17], two-fluid formulation for dispersed flows [18,19] or rheological model for granular flows [20].
Of course SPH is not all roses-it has some drawbacks that make CFD community a tad sceptical about the method as a general-purpose tool and its usefulness beyond specific applications. Due to particular interpolation method, SPH has a considerable computational cost when compared to the grid based methods. In SPH interpolation nodes are moving, hence there is necessity for updating the list of particles interacting with each other every time-step. The algorithms searching for neighbouring particles are a challenge in itself [21,22], while in Eulerian approaches this problem does not exist at all. Furthermore, the SPH interpolation requires calculation of many more interactions than grid-based methods. This is particularly problematic due to rather poor convergence rate of the method, as shown in [23]. To make things worse, SPH still lacks a decent and reliable method of adaptive resolution refinement. Some works were done in this matter [24,25], however, SPH has a long way to go before reaching the level of grid-based methods, even in standard techniques like simple remeshing of near-wall regions. To overcome these problems, promote and utilise its full potential, the SPHERIC group-ERCOFTAC Special Interest Group no. 40 [26], was formed. Wide area of possible applications and ability to go beyond limitations of grid-based methods result in growing interest in SPH.
The present study is dedicated to the application of SPH to computations accounting for the wetting phenomena. For this purpose we use a modification of the Continuum Surface Force (CSF) model of surface tension, with a specific ghost particles approach for solid walls. Although some works in SPH have been done in this matter, the present study is believed to be a step forward. Our approach bears a similarity to the one of Breinlinger et al. [27], however as proved by a later study, that model requires artificial scaling factors [28]. It is worth to mention the relevant work of Huber et al. [29] where CSF was supplemented with a contact line force model. The proposed method, although based on a physically sound modelling of the triple point dynamics, does not seem offer much room for improvement concerning superhydrophilic and -hydrophobic cases.
In a recent work of Nair and Pöschel [30] the contact angle model was successfully implemented within free-surface, incompressible SPH formulation. However, such a framework can not be applied to the simulations of the droplet interacting with an external air flow. Furthermore, that work uses a microscopic surface tension model, also applied in [31], which is not straightforward to couple with other physical models. The microscopic formulation has sound thermodynamical foundations and excels in problems where small (meso) scales are of importance, like the pore scale modelling of flow in porous medium [32]. Also, that approach can be applied to the computations of the droplet interacting with a rough surface, whose geometry is explicitly modelled in its complex entirety [33]. On the other hand, using the microscopic model for problems where macroscales are of interest, can be problematic due to orders of magnitude differences in, e.g. lengths considered.
From the mechanical engineering point of view the macroscopic description of wetting phenomena remains useful and is the most often applied for reasons of computational efficiency. Our idea is to improve the existing formulation for the CSF, as a reliable tool for modelling capillary effects. We also intend to do it for a sharp interface formulation, not a diffused one like in the diffused interface approach [34]. The method proposed in the following can be summarised as: (i) simple and robust, (ii) free of scaling parameters, (iii) based purely on CSF modification, (iv) maintaining the natural interface treatment of SPH.
The paper is organised as follows: Section 2 introduces fundamentals of the SPH method for multiphase flows; in Section 3 the modified surface tension model and the wall boundary conditions (b.c.) are discussed; results of the proposed model are presented in Section 4, and finally conclusions are drawn in Section 5, together with the outlook for future work.

Basics of the SPH
Interpolation theory is the very core of SPH. We now briefly recall the essentials required for a better understanding of this specific method. Let us consider any field A. The integral formula where δ(r) is the Dirac delta function, can be used to express the field value at the point r in space Ω. To obtain the SPH approximation, we first replace δ(r) with the weighting kernel function W (r, h) which should be normalised, symmetrical and converge to δ(r) with h → 0 [35]. Argument h is the so-called smoothing length and it determines the interpolation range. In our work we use the quintic kernel proposed by [36] W (r, h) = C 1 − q 2 4 (2q + 1) , for q < 2 0, otherwise, where q = |r|/h and C is the normalisation constant (C = 7/4πh 2 in 2D and C = 21/16πh 3 in 3D), as guaranteeing good stability of computations [37,38]. The second step consists in the discretisation of space into a set of particles of volume Ω b = m b / b , where m b is the mass and b is the density of b-th particle. As a result, the integral from Eq. 1 is approximated by a sum, i.e.
( 3 ) In the shorthand notation, the SPH approximation A a of field A at any point a is defined as where Thanks to the properties of W (r, h), namely compactness of the support and symmetry, differentiation can be shifted from the field to the kernel function yielding Further derivatives can be obtained in a similar way, however, for an improved accuracy (habitually impaired by irregular distribution of the particles) more sophisticated schemes are used [39].
Using the above method various kinds of differential equations can be rewritten into the SPH formalism and solved by calculating interactions between particles, hence its wide application. The readers interested in the detailed information on the derivation of SPH and basics of the method are referred to the monograph by Violeau [40], where the approach is deeply discussed in the context of fluid mechanics.

Governing equations
The set of governing equations for the Newtonian viscous incompressible flow consists of the Navier-Stokes (momentum) equation (6) and the continuity equation where u denotes the velocity vector, is the fluid density, p is the pressure, μ is the dynamic viscosity and g is the mass force. In Eq. 6, f st stands for the interfacial force, detailed in Section 3.1. Due to the Lagrangian nature of SPH we also include the advection equation Depending on the purpose and assumptions, different SPH formalisms for the fluid flow can be obtained by using Eqs. 4 and 5. In the present study we use a formulation proposed in [41]. To the best of our knowledge, this approach is best suited for modelling multiphase flows with large density ratios [42]. The pressure term in Eq. 6 will become where θ is the inverse of particle volume. The viscous term of Eq. 6, obtained by combining SPH formalism and Finite Difference Method, takes the form where The key feature of the approach proposed in [41] is the treatment of the continuity equation. Instead of rewriting Eq. 7 into the SPH language, density is found from This allows the density field to be represented only by the spatial distribution of particles and not by their masses. Thanks to this, densities of particles near the interface are not affected by the other fluid, which is important in multiphase flow modelling.
In this work we use the Weakly Compressible SPH approach (WCSPH). The set of equations is closed with the equation of state where s is the artificial speed of sound, 0 is the reference (initial) density and γ is a numerical parameter. The values of c and γ are chosen to ensure density fluctuations at the level of 1% or below. In multiphase flow modelling it is a common practice to treat liquid as incompressible with γ = 7, and gas as compressible with γ = 1.4 [43]. We follow this approach in this study.

Surface tension model
The surface tension force acting on the interface is given as where σ is the surface tension coefficient, κ is the local curvature of the interface andn is the unit vector normal to the interface. The most popular model of the influence of surface tension is the classic Continuum Surface Force method (CSF), originally proposed by [44], with SPH implementation due to [45]. In this approach the surface tension forces (per unit area) are converted into a force per volume where δ s is a suitably chosen surface delta function. Using the phase indicator or the socalled colour function c (say, c = 0 for the first phase and c = 1 for the second one),n can be calculated using the formulan In SPH the vector n is obtained from wherec stands for the smoothed colour function, i.e.
The above operation is performed to diminish the so-called parasitic currents, i.e. spurious velocity field arising from the unavoidable inaccuracies due to discrete formulation. Interestingly, the smoothing of phase indicator function with the SPH kernel was first used in VOF computations, also as a remedy for parasitic currents [46]. The local curvature is obtained from which in SPH language takes the form Note that in Eqs. 16 and 19 an alternative form of SPH gradient is used to ensure zero value of the derivatives for the constant field. Assuming that δ s = |n|, the influence of surface tension can be included in Eq. 6 as This pretty straightforward translation of well-acknowledged CSF model into the SPH language yields accurate results in various cases [15,16]. However, for the modelling of wetting some modifications are in order. In such cases, as will be shown later, the accurate calculation of ∇c is of the essence, as it is the main factor affecting the movement of the interface near the wall. During initial tests of the model we noticed some problems with estimation of ∇c field in hydrophilic cases where droplet may become extremely thin. Due to this fact the tip of the elongated droplet is comprised of only a few particles of the phase. This will always happen, independently of the resolution. Increasing it will only result in smaller area being represented by insufficient number of particles. The tip containing the contact line, however, will always be underresolved. Due to the nature of SPH interpolation, this will result in underestimated length of the ∇c vectors. To counter this issue we decided to employ the so-called Shepard's kernel [47]. The purpose of this technique is to correct under-and over-estimations arising from unavoidable irregular distribution of particles, or from a deficient number of particles under the kernel hat as originally done in free-surface flows. Shepard's kernel for a-th particle is defined as The approximation of any field A with Shepard's kernel, i.e.
results in renormalisation of the considered field, but also its additional smoothing. The first outcome is always most welcomed, the second one-not necessarily. Our idea is to apply Eq. 22 to the field of ∇c. This operation will allow to calculate ∇c more accurately for fine structures, however, additional smoothing would make the numerical representation of the interface two-times wider. This is very undesired since it would increase twofold the distance of premature coalescence, a well known issue among CSF users. To avoid this problem, but still use Shepard's kernel, we propose the following scheme for the surface tension computations: 1) calculate ∇c from the sharp colour function, i.e.
2) then use Shepard's kernel on the resulting field of n a but without any boundary conditions, see Section 3.2; the superscript Sh is skipped in the following. We note that Huber et al. [29] recently applied the Shepard kernel and its gradient for all the governing equations. Unlike them, we deem it necessary only for the field of normals, as other quantities are sufficiently well estimated. 3) finally proceed further as in the implementation of Morris.
This way lengths of n are better estimated, and their field is as smooth as in the basic implementation of Morris, Eqs. 17 and 16.
To illustrate the problem we performed calculations of the ∇c field for a droplet placed on a superhydrophobic surface (for the contact angle θ D = 15 • ). As an initial condition we used steady-state solution of SPH simulation with model described in Section 3.3. Only field of ∇c was computed, no flow equations were solved. The top panel of Fig. 1 shows raw particles coloured by their respective phase indicator. The middle panel demonstrates the field of |∇c| (projected onto a uniform grid of size r, which denotes particles' initial spacing) obtained with Morris implementation [45]. As we can see, the biggest magnitude is observed close to the interface, which is not surprising, however it is underestimated close to the triple line. The improved scheme, as shown in the bottom panel, offers the same accuracy in the part of the droplet away from contact line, but with increased maximum value of |∇c| in its closest neighbourhood. In Fig. 2 we compared the |n| profile along the wall computed with the Morris approach and the new scheme. The maximum value of |n| can be interpreted as an indicator of the interface position, hence for the profile at y wall , it marks the triple point. The profile at the wall distance y = 2h served as reference, since it was exactly the same for both methods. The profile for y = 2h was shifted for easier   [45] implementation and the new scheme; also shown the reference solution at the distance 2h from wall comparison. It is clear that using Eqs. 16 and 17 results in underestimation of ∼25%, while Shepard's kernel yields overestimation at the level below 10%.
Following [45] we also exclude from calculations of surface tension effects the particles for which |n| is below the threshold of 0.01/h. This greatly improves accuracy in curvature estimation on the fringes of the interfacial area. To mitigate errors caused by exclusion of some particles from the summation formula in Eq. 19, the resulting curvature has to be rescaled by a correction factor given as where N = 0 if |n| < 0.01/h and N = 1 otherwise, for a and b alike. The above equation will also come in handy later on, see Eq. 29.

Solid walls model
The The biggest advantage? Easiness in implementation, robustness and appealing simplicity of the idea, but most importantly of all, the applicability in modelling complex geometries. ii) the ghost (mirror) particles [38,51]-in this method a flat wall acts as a mirror-each time step the particles close to the wall (at distance less or equal to the interpolation range of the kernel function used) are mirrored onto its opposite side. The physical properties of such ghosts can be adjusted so that the desired boundary conditions are satisfied. The biggest disadvantage of this approach is limitation only to flat walls. iii) unified semi-analytical walls (USAW) and others [52,53]-a complex and sophisticated method, in which the so-called wall renormalisation of the equations is used to account for the particles missing from the kernel support. This method is typically employed in turbulent single-phase flows, which are very challenging for SPH [54].
In the case of multiphase systems the USAW method can be used, e.g. a recent work on the SPH mixture model [55]. Since it was developed with a different purpose in mind, it does not have any specific feature which would make it tempting for interfacial flows modelling. The most popular solid walls model, in general SPH applications, is the dummy particles approach. By far, it is the easiest one in implementation and despite its simplicity can produce accurate results. However, for multiphase systems a simple question immediately arises-from which phase should the wall be built? The answer is by no means trivial. In the work of Yeganehdoust et al. [56] this approach was used for modelling the droplet-wall interaction. The colour function value for the dummies was derived from the instantaneous position of droplet, to allow for accurate calculation of the interface normals in the vicinity of the triple line. The proposed approach allowed for satisfactory results, but explicit tracking of the relative position of the droplet and the wall is problematic from the algorithmic point of view and limits versatility of the method.
In our work we use the ghost particle approach. In the case of multiphase flows the material properties of mirrored particles-density, viscosity and colour function-are exactly the same as those of parent ones. The velocities are constructed to enforce no-slip b.c., and pressure is simply the same as it ensures the homogeneous Neumann b.c. This method is in fact very convenient for interfacial flows, since: (i) the boundary conditions are being constructed on the fly and there is no need for tracking the phase distribution in the flow, (ii) all summation formulae are calculated accurately without risk of under/overestimation of density or similar error. The main drawback of the ghost particle approach is its applicability only to cuboid domains, however, for the purpose of basic research and tests of new models this is more than enough.

Modified CSF -contact angles
We are not trying to reinvent the wheel. We follow the idea proposed in the original CSF paper by Brackbill et al. [44], namely modification of the unit normalsn in the vicinity of the triple line. This approach was successfully applied within the VOF framework in a number of papers [57,58]. Since we are dealing with SPH, we had to make some adjustments in the existing formulation. We are also obliged to mention that this work is not the first attempt to do so in SPH, the approach of modifyingn was already applied by Breinlinger et al. [27], however, they used dummy particles approach and a different formulation of CSF model. They also did not address the problem of ∇c calculation for very thin droplets in hydrophilic cases. Now, let us cut to the chase. As mentioned, we will modify the unit normals to the interface in the neighbourhood of the triple line, which boils down to modifyingn for the SPH particles in this region. For the reader's convenience, please refer to Fig. 3 with a schematic illustration of the objects used in the following. First of all, how do we define particles "in the neighbourhood of the triple line"? The following straightforward conditions have to be satisfied: (i) the particle is within the distance less or equal to 2h from the solid boundary and (ii) it has a non-zero ∇c. The first requirement is dictated by the interpolation range, while the second one ensures that considered particle is, in a numerical sense, part of the interface. Conjunction of those two conditions holds true only for particles close to the other phase and the solid surface. In this triple point region, to obtain the desired contact angle θ D we imposen asn where n ⊥ and n are the unit vectors perpendicular and parallel to the wall, respectively. Initial tests showed that the sudden change from original to modified normals will result in deformed interface. The same problem was encountered by Breinlinger et al. [27]. To alleviate this issue they used a blending function working as a continuous seam between modified and primal normals. The function used is a simple, linear relation of the wall distance of the considered particle, given as α(y) = 1 − y/2h. The auxiliary normal n aux (note that it is not exactly a unit vector here) is defined as and the final form of modifiedn is of coursê The ghost particles are completely excluded from the curvature calculation, thus there is no need to define theirn mod . Correction by Eq. 25, originally employed due to omitting particles with normals length below the defined threshold, is enough to alleviate errors that could be caused by this move. The curvature is then calculated using new values ofn, hence Eq. 20 becomes Note that the actual force moving the particle acts according to the direction and magnitude of unmodified ∇c, hence the modification in calculating its field, as described in Section 3.1.
We also need to stress that within the ghost particles approach framework, the solid wall is not physically present in the calculations. Due to the formulation, the particles are not penetrating it or getting exactly to its position: y a = y wall will always hold true for any particle a. Unlike a common practice in Eulerian approaches [57], in the present simulation Fig. 3 The vicinity of the triple point: illustration of the scheme used in modification ofn we did not see necessity to add a microscopically justified slip condition for the triple point vicinity, see also [41].

Results
We present and discuss here some results obtained with the model described in previous sections. All computations were performed with our in-house SPH code. Particles were initially uniformly distributed on a cartesian grid pattern. The interface was taken as the isoline c = 0.5 of the colour function interpolated with Eq. 4 onto an auxiliary uniform grid of size r [59]. The benchmarks and their analytical solutions reported in Sections 4.1 and 4.2 are taken from [57].

Static case
The first validation case we considered was the two-dimensional liquid droplet on a solid surface, with prescribed desired contact angle θ D . The droplet was initially semicircular, with radius R 0 , placed on a bottom wall of a rectangular channel (filled with gas) with the height of 4R 0 and length 8R 0 (or 16R 0 for θ D < 30 • ). The density and viscosity ratios were respectively L / G = 1000 and μ L /μ G = 100. The Laplace number is defined as La = σ L R 0 /μ 2 L and is set as La = 10 4 . The resolution used was h = R 0 /16 and h/ r = 2, yielding 32 SPH particles per droplet's radius. Due to the imposed θ D the droplet-wall contact area will either spread out (θ D < 90 • ) or contract (θ D > 90 • ). In steady state (without gravity), the droplet should assume the shape of a spherical cap, cut from a circle of radius R at a chord (namely wall), so that the angle between the tangent at a common point of circle and chord is equal to θ D . The droplet can be more easily described by its height H and area of contact with wall w (wetting area or width in 2D). This problem can be analytically described by the following relationships of R 0 , θ D , R, H and w: w an = 2R an sin θ D . (32) We decided to use values of H an and w an as the reference data, since they are both dependent on R, which corresponds to some pressure jump according to the Young-Laplace equation.
In Fig. 4, the interfaces for θ D ranging from 30 • to 150 • are plotted. The qualitative results are satisfying in general-for hydrophilic surface we observe expected spreading,   Fig. 5, and it turns out less pleasing. While the results for θ D from 45 • to 120 • are very accurate, those for more extreme values are way off the mark. The inaccuracy becomes even more evident when we look at Fig. 6. The problem of poor results for superhydrophobic and -philic surfaces can not be explained just by the resolution applied. Increasing the resolution twofold, i.e. by using two times smaller h, while keeping h/ r = 2, availed us nothing, except increased computational time. This also bares a drawback of SPH-the lack of well-established adaptive mesh refinement and poor convergence rate.
To alleviate this issue we decided to change the blending function a bit. As mentioned earlier, the particles will never get directly to the position of the wall, hence using function α(y) = 1 − y/2h implies that no particle will ever have exactly the normals imposed by Eq. 26. Also, from Figs. 5 and 6 we see that the value of w is underestimated for hydrophilic cases and overestimated for hydrophobic ones. This points towards conclusion that the a) compariso n with analytica l solut ion b) accurac y of th e result s from Fig. 8(a)   Fig. 8 Comparison of the results obtained with α(y) (× symbols) and α + (y) (square symbols) Fig. 9 Comparison of the results obtained by using α(y) (× symbols) and α + (y) (square symbols) with analytical solution (lines) acceleration resulting from Eq. 29 is too small next to the wall. To increase it, we changed α(y) in following manner: In Fig. 7 we plotted the comparison of the interface shape for cases of θ D = 15 • and 30 • obtained with α(y) and α + (y). As can be seen, for α + (y) the droplets are spread wider. The analysis of the w and H values, shown in Fig. 8a, proves higher accuracy of the modified blending function, also fortified by Fig. 8b. Alas, this method works slightly worse for moderate and high values of θ D , see Fig. 9, hence it can not be treated as the universal approach. We need to mention that prior to applying Shepard's kernel, Eq. 24, we used Eq. 16 for the normals calculation, however, the results were far less accurate (hence we do not present them herein). Additionally, we took the relative errors in contact angle reported in [29] (see Fig. 11(a) there) and recomputed them to assess (at comparable resolution) the resulting inaccuracies in droplet's width w. Specifically, for the contact angle of 30 • , the error level was about 9% there and 4% in the present work. Then, for the angle of 120 • , the errors were of 8% and again 4%, respectively. This indicates that the contact line force model, notwithstanding its elegance, may be more problematic in SPH computations for extreme values of θ D . Despite inaccuracies for some cases, the proposed model works very well for a wide range of contact angles. Results obtained with α + (y) allow us to think that this model, after some further modifications, can also be used for cases with extreme values of θ D .

Spreading effects in the presence of gravity
The second validation case was investigation of a droplet spreading due to the effects of gravity. The resolution and initial setup are the same as in Section 4.1, however this time gravity g is introduced into the system, to test whether the model works well when taking into account the balance between surface tension forces and external forcing. Calculations were performed using α(y). To describe this case we introduce the Eötvös number, Eö = L gR 2 0 /σ , which describes the ratio of gravitational and capillary forces. We have tested the model for θ D = 50 • and 130 • , and Eö ranging from 10 −3 (capillarity dominated regime) to 50 (gravity dominated regime). In the case of Eö 1 the height of the droplet can be calculated as analytical solution when no gravity is present, i.e. by using Eqs. 30 and 31. When Eö 1, the droplet height (or thickness) is equal to capillary length defined as [57] In Fig. 10a selection of results for the case of θ D = 50 • is shown. As expected, we observe increasing flattening of the droplet with growing Eö. Quantitative results in Fig. 11 show that the model correctly predicts transition from surface tension to gravity dominated regime. Up to Eö ∼ 10 −2 the droplet has shape as in the case without gravity, then at Eö ∼

3D simulations and 2D droplet in air flow
As the first of the final tests we have applied the model to 3D computations of sessile droplet. We do not present a proper validation for those simulations, as they should be considered an eye candy showing upper limit of the model capabilities. For the 3D implementation, the problem of modelling and coding poses no additional difficulty. Since SPH is inherently multidimensional, only vector n in Eq. 26 requires a bit more of simple algebra than in 2D setup to be derived. The actual problem is extensive computational cost of SPH. Yet, this should not be thought of as a dealbreaker when considering SPH as a numerical tool. Thanks to the local character of the particle interactions, the approach is suitable for massive parallelisation and execution on Graphic Processing Units (GPUs). Nonetheless, we were using the code which is optimised for a single GPU only. For this reason we performed calculations only in a very coarse resolution to check the correctness of the implementation. We choose the static case from Section 4.1, with resolution defined by h = R 0 /16 and h/ r = 1.5625. Such a low value of h/ r is dictated by computational time only. Figures 14 and 15 present outcomes of computations for θ D = 30 • and 150 • , respectively. As can be seen, the model works just fine in 3D, which is no surprise to be honest.   16 Initial setup for the moving droplet case. Royal blue and green-liquid, revolutionary red and light blue-gas, mass force g acts from left to right More interesting is the case of the droplet movement on a solid surface, driven by air flow-the cherry on top of the present study, even if it is only in humble 2D configuration. The problem of moving contact lines introduced by the droplet motion is a challenge for grid based methods. When no-slip condition is imposed at the wall, the so-called stress singularity caused by the movement of contact line appears [60,61]. To remedy this difficulty, the Navier slip condition is usually applied in the computational cells containing triple point. This is also often accompanied by further modifications in the wetting model, for a brief overview see [61]. This seems to be not the case in the SPH modelling. To the best of our knowledge no SPH study dedicated to the wetting phenomena reported any singularity issues [27, 29-31, 34, 56], except for [62]. In our work we also did not observe problems related to the stress singularity. Increasing spatial resolution also did not uncover or strengthened any numerical artifacts-in the grid based methods the stress at the contact line is diverging when refining the mesh size. Why SPH seems to be free of this problem? A possible explanation coming to our minds is that this is the Lagrangian method. Surface tension works as a force, similar to, e.g. gravity or pressure gradient, that compels particles to move in a given direction. Moreover, in SPH computations, especially when the ghost particles approach is used for the solid boundaries like in the present study, the walls are not physically present. This means that we never model the situation exactly at the wall position, we only advect the particles according to the acceleration computed from the momentum balance equation above it. For those reasons, and due to satisfactory results, we did not attempt to implement any model for dynamic contact angles.  Now, let us serve the promised cherry. The initial setup is the same as for the static case, except that the mass force g acts parallel to the walls. We have set periodicity at the left and right boundaries of the domain, while the top and bottom ones are modelled as solid walls with no-slip b.c. Due to the presence of g, a Poiseuille-like flow will form in the channel. We also used higher resolution, i.e. h = R 0 /32 with h/ r = 2. The initial setup with colours used for marking purposes is shown in Fig. 16. The magnitude of g was chosen so that the Reynolds number estimated with channel's width would be Re 1 (creeping flow regime). The qualitative results are, well, spectacular. In Fig. 17 the movement of the droplet on hydrophilic wall (θ D = 30 • ) is shown. As expected, we can observe stretching of the droplet in the flow direction and its motion in such elongated form. For the hydrophobic wall (θ D = 150 • ), the results are in agreement with physical intuition as well, the droplet simply starts to roll to the right, see Fig. 18. The results for moderate contact angles are slightly less visually enjoyable, the droplets are assuming the final shape and just steadily moving according to the flow, for example see Fig. 19. The quantitative analysis can be performed with e.g. the Cox-Voinov law, however, for now we restrained ourselves from such investigation. Before doing so, a more thorough study should be conducted on the dynamic contact angle models and their implementation within SPH framework, which also should address the problems with spatial accuracy and convergence of the method. The very first try looks promising and encourages us to further tread the chosen path.

Conclusions and Outlook
In the present study we have presented a simple and robust, yet effective implementation of the contact angle model for SPH. The model consists in modification of the interface normals enforcing the desired contact angle. The approach was implemented within the ghost particles b.c. formulation. The model was tested against a wide range of static cases, and yielded accurate results for the contact angles θ D ∈ [45 • , 135 • ] when compared to analytical solution for droplet's height and width. For superhydrophilic surfaces, i.e. θ D < 45 • , a different blending function can be used with significant increase in accuracy. The results for superhydrophobic walls were less accurate, but still satisfactory. The model was also tested for the effects of gravity, yielding satisfactory estimation of the droplet's height and accurate prediction of the transition between capillarity-and gravity-dominated regimes. We also need to stress a small, but important, modification in the SPH handling of the CSF model. Instead of smoothing the colour function, like it is usually done, we proposed to smooth the field of its gradient with the Shepard interpolation. This operation allowed for better estimation of the normals to the interface near the triple point. The model presented in this work can be easily extended to the 3D setup, as proved in Section 4.3.
Finally the approach was tested in computations of droplet exposed to the creeping air flow in a channel. The obtained results are in agreement with physical intuition [63][64][65]: rolling and tumbling of droplet was observed for hydrophobic walls, while smearing was the effect of simulating the hydrophilic surfaces. For moderate values of θ D we observed steady movement of droplet whose shape was dictated by prescribed contact angle. The next step to take should be focused on a thorough validation of the model in the dynamic cases. The problem of stress singularity, which is observed in grid-based methods, does not seem to be present in SPH. Such a study should answer the question whether Lagrangian methods are free of this drawback. For that purpose detailed investigation of dynamical models and their implementation in SPH should be put into the test. Wall-bounded laminar flows, dominated by surface tension effects, like formation of droplets in microchannels [66] (the so-called lab-on-chip devices) can be considered as possible field of application in which SPH may excel despite its extensive computational cost.
Development of reliable, robust model for contact angles is a prerequisite for SPH to become a general purpose tool for multiphase flows computations, possibly even considered as a respectable rival for the well-established Eulerian methods. Presented results, their accuracy and promising outlook allow us to state that this moment already looms on the horizon.