Consistent Thermo-Capillarity and Thermal Boundary Conditions for Single-Phase Smoothed Particle Hydrodynamics

A model for capillary phenomena including temperature-dependency and thermal boundary conditions is presented in the numerical framework of smoothed particle hydrodynamics (SPH). The model requires only a single fluid phase and is therefore computationally more efficient than surface tension schemes which need an explicit fluid-fluid or fluid-gas interface. The model makes use of a surface identification mechanism based on the SPH renormalization tensor. All relevant properties of the continuum surface force (CSF) based approach, i.e., the delta function, normal vector and curvature, are calculated in a consistent manner. The model is parametrized by physical material properties and is successfully validated by means of a large set of analytical test cases. The applicability of the proposed model to more complex scenarios is demonstrated.


Introduction
Modeling surface tension in the framework of smoothed particle hydrodynamics [1,2] is generally done using one of the following two approaches.
On the one hand, surface tension can be formulated as a continuum surface force. This approach introduces a surface delta function to approximate the tension, i.e., a force per area, in form of a force per volume [3]. A color function is used here to identify the involved phase and the interface between them. An early SPH implementation of the CSF approach for two phases was described and analyzed by Morris [4]. Hu and Adams [5] used an interface stress tensor to avoid the explicit calculation of the curvature at the phase interface. They also extended the formulation to more than two phase which allows modeling of contact angles. Adami et al. [6] improved the CSF scheme by an accurate curvature calculation. Their method is also able to handle large density ratios of the involved phases of up to 1000. Zhang [7] enabled surface tension modeling in a single-phase SPH simulation by identifying the boundary particles and reconstructing the phase surface in two dimensions (2D) or three dimensions (3D). Breinlinger et al. [8] were able to model wetting effects in two-phase SPH simulations by prescribing the contact angle with respect to a rigid substrate. Thermo-capillary effects based on a temperaturedependent surface tension were studied by Tong and Browne [9] using a two-phase SPH approach. For this purpose, they used a projection of the gradient operator onto the phase interface. Yeganehdoust et al. [10] developed a method to extrapolate the color function into a substrate in order to model both static and dynamic wetting phenomena. Huber et al. [11] introduced a balance equation for the contact line between the two fluid phases and the substrate to accurately model advancing and receding contact angles in dynamic wetting situations. Ordoubadi et al. [12] used a single-phase SPH formulation for the CSF approach in which they identified the surface by a combination of a geometric interface tracking and the calculation of the position divergence field [13]. Hirschler et al. [14] formulated a single-phase model which obtains the surface normal and the surface delta function directly from the SPH kernel gradient. The surface curvature is obtained using a corrected kernel gradient [15]. Krimi et al. [16] used an interface stress tensor formulation to model multi-phase flow situations with large density ratios of the involved phases. Hopp-Hirschler et al. [17] investigated thermo-capillary phenomena in internal flow situations or for two-phase systems, respectively. They stabilized the SPH simulations by a particle shifting technique [18]. Lin et al. [19] proposed the detection of an interface between two phases as well as the calculations of the interface normal and curvature based on a Delaunay triangulation. Fürstenau et al. [20] formulated a single-phase description using the position divergence to identify the surface and SPH kernel gradient to compute the surface normal. They improved the surface curvature calculation by using a curvature tensor. Liu et al. [21] proposed a machine learning approach in order to obtain a more accurate local interface curvature in two-phase simulations. Härdi et al. [22] investigated wetting phenomena of thin films by using an SPH discretization of the shallow water equation and introducing a contact line force. Liu et al. [23] proposed a single-phase surface tension model based on the δ + SPH scheme [24]. They calculated the surface normal from the kernel gradient and used a normalization procedure to improve the accuracy of the surface delta function.
On the other hand, a surface tension formulation can be based on an additional pairwise force which is repulsive at short range, attractive at intermediate range and vanishing at long range. Such models were developed by Nugent and Posch [25] and refined by Tartakovsky and Meakin [26] as well as Akinci et al. [27]. A pairwise force mimics the molecular origin of surface tension [28] and results in a net force towards the bulk for particles at a curved, convex surface [29]. The pairwise force approach originally required calibration by adjusting model parameters which have no correspondence in experiments. Tartakovsky and Panchenko [30] made inroads to overcome this drawback by deriving relations between the pairwise model parameters and physical properties for multi-phase systems. Nair and Pöschel [31] extended this approach to single-phase simulations. Bao et al. [32] improved the pairwise force formulation for dynamic wetting phenomena by introducing an apparent viscosity at the interface to a rigid substrate which mimics the effect of a finite slip length.
In this work a CSF single-phase surface tension scheme is presented which is an extension of an earlier work by the author [33]. The approach for modeling wetting behavior is adapted from Breinlinger et al. [8] and the method to incorporate thermo-capillarity from Tong and Browne [9]. The novelty of the present work is the consistent formulation which is used to compute all surface related properties. The surface identification scheme from Marrone et al. [34], relying on the renormalization tensor, is used as the basis to obtain the surface normal, surface curvature and surface delta function in a consistent and accurate manner. In previous works, at least one of these three surface properties was not based on the renormalization tensor, leading to an inconsistency in the definition of the actual surface. The surface delta function is also used to formulate thermal boundary conditions which take into account heat transfer, heat flux and radiation.
The remaining part of this paper is organized as follows. Section 2 outlines the governing equations and presents all relevant ingredients of the numerical model. Section 3 uses several benchmark cases to assess the predictive quality of the model. The obtained results are discussed in Section 4. Section 5 concludes the main part of the paper. Details of the surface delta function computation are given in Appendix A.

Materials and Methods
The governing equations with respect to mass balance, momentum balance and energy balance are presented in Section 2.1 with focus on the surface related terms. The discretization of the governing equations using SPH is described in detail in Section 2.2. The main contribution of the present work is the SPH surface treatment presented in Section 2.2.2. It forms the basis for modeling surface thermal boundary conditions (Section 2.2.4), surface tension (Section 2.2.10) and surface wetting (Section 2.2.11).

Governing Equations
The evolution of density ρ over time t is given by the mass balance equation, where v is the velocity. The Navier-Stokes momentum equation provides the time evolution of velocity, with pressure P, dynamic viscosity η, rate of strain tensor E = (∇v + (∇v) T )/2, volumetric surface tension force F S and gravitational acceleration g. The time evolution of the position r is simply Dr The surface tension force can be decomposed into components normal and tangential, respectively, to the surface, The CSF model converts a force per area, i.e., pressure or tension, into a force per volume by multiplication with a delta function δ S which marks the location of the surface in space. This approach yields the following expression for the normal component of the volumetric surface tension force, with σ n as the surface tension, κ as the curvature and n as the surface unit normal vector. A corresponding tangential volumetric force caused by a gradient of the surface tension is given by with ∇ S as the Nabla operator along the surface. One reason for a non-constant surface tension is a dependency on temperature T which might be expressed as ∇ S σ n = σ t ∇ S T with σ t = ∂σ n /∂T which is typically negative and causes Marangoni flow along the surface from warm towards cold regions. The delta function is constructed as the magnitude of the gradient of a color function which is 1 on one side of the surface (or interface) and 0 on the other side. The evolution for the thermal energy per unit mass H is given by Here, k is the thermal conductivity, q is the heat flux, f is the heat transfer coefficient, is the emissivity, σ B is the Stefan-Boltzmann constant and T a is the ambient temperature.
The heat capacity c relates the temperature to the specific thermal energy,

SPH Kernel Approximation
A quantity A at the position r is interpolated using a weighted sum over the contributions A(r j ) of all neighboring particles j at positions r j with masses m j and densities ρ j , where W(r, h) is a kernel function. The particle mass is related to the rest density ρ 0 and the initial particle spacing ∆x via m j = ρ 0 (∆x) d where d is the number of spatial dimensions.
is the quintic C2 Wendland kernel [35] with h as smoothing length and ζ = r/h. The normalization factor is ν = 7/(4π) for d = 2 and ν = 21/(16π) for d = 3. Equation (9) can be used to obtain a smoothed approximation of a field quantity at the position of particle i via where A j = A(r j ), W ij = W(|r ij |, h), and r ij = r i − r j . In order to discretize the governing equations in the SPH scheme, spatial gradients of field quantities can analogously be expressed via where the kernel gradient ∇W ij is explicitly given by In comparison to the canonical formulation given by Equation (12), the expressions and prove to be numerically favorable [36] and are therefore used in the following SPH discretizations.

SPH Surface Modeling
A discretization of the governing equations requires the computation of several surface related quantities such as the surface delta function δ S , the surface unit normal n, the surface curvature κ and the gradient operating along the surface ∇ S . The following section describes our approach to obtain these quantities.
The present SPH formulation makes use of several aspects introduced in [14,24,34]. A central feature is given by the method to identify the surface by means of the renormalization tensor L. The inverse renormalization tensor for particle i is given by The minimum eigenvalue λ i of L −1 i defines a scalar field which is 1 inside the volume and approaches 0 at the surface of the particle distribution [34]. This property makes the λ field similar to the color function on which the CSF model is based. This analogy is used as key ingredient of the presented surface tension model. The corrected gradient of the λ field for particle i is given by and the surface unit normal vector follows simply as Following [24] we identify those particles with λ i ≤ 0.2 as free surface particles. For each particle i with λ i > 0.2 and λ i ≤ 0.75 an umbrella-shaped region with cone half angle 45 • and range 2 ∆x in the direction of the normal n i is scanned. If there is no particle j within this umbrella-shaped region, the particle i is also identified as free surface particle. We define the set of particles which are at most a distance of 1.5 ∆x away from the free surface particles as surface region particles, including the free surface particles. Thus, the surface region is typically composed of two neighboring layers of particles. In order to avoid numerical problems caused by isolated particles, we exclude particles with λ i < 0.1 from the surface region. All further surface related calculations are carried out only for the surface region particles defined in this paragraph.
Furthermore, we define the surface delta function by where Ξ λ is a dimensionless normalization factor which depends on the ratio h/∆x and the actual function used as SPH kernel (see Appendix A). An alternative approach to obtain the surface delta function follows [14], where ∇1 i is Equation (12) applied to unity at the position of each particle j, i.e., A j = 1. Differences between the approaches given by Equations (19) and (20) are discussed in Appendix A. The local curvature of the surface is then calculated following [14] using a corrected kernel gradient, The sum in Equation (21) is restricted to the surface region particles. The curvature is limited by the spatial resolution, 2 h |κ i | ≤ 1.

SPH Mass Balance
The mass balance (Equation (1)) is discretized using Equation (14) in the form where v ij = v i − v j is the relative particle velocity.

SPH Thermal Energy Balance
The thermal energy balance (Equation (7)) for each SPH particle i takes the form The first term on the right hand side representing thermal diffusion is calculated using the approach from [37]. The next three terms for heat flux, heat transfer and radiation, respectively, make use of the surface delta function given by Equation (19). The temperature for each particle i is then simply obtained as T i = H i /c.

SPH Temperature Gradient
In order to model thermo-capillary effects, we need to obtain the temperature gradient within our SPH formulation similarly to [9]. A corrected gradient formulation [15] is used in order to get reliable values at the surface,

SPH Momentum Balance
The Navier-Stokes momentum balance (Equation (2)) for each SPH particle i takes the form with the accelerations caused by pressure, f P i , by viscosity, f V i , by artificial viscosity, f A i and by surface tension, f S i .

SPH Pressure
The pressure acceleration [1] is obtained by using the gradient discretization according to Equation (15), The pressure P i of particle i is related to its density ρ i by Tait's equation of state, where s is the numerical speed of sound and γ is the dimensionless isentropic exponent.

SPH Viscosity
Depending on the details of the actual simulation scenario, one of two approaches to model the effect of viscosity on the flow is used. Cleary [37] introduced a formulation based on the relative normal velocity of each pair or particles i and j, An alternative formulation is given by Sigalotti et al. [38] which is structurally similar to the pressure term (Equation (27)), The rate of strain tensor E i = (∇v i + (∇v i ) T )/2 of particle i is based on its velocity gradient tensor,

SPH Artificial Viscosity
An artificial viscosity formulation is in some situations required to stabilize the numerical simulation by damping high frequency noise. It should not influence the dynamics of the simulation to a statistically significant extent. The artificial viscosity term is given by Monaghan [1], where with α and β as dimensionless control parameters of the order of 10 −2 and φ ij as a rescaled normal velocity,

SPH Surface Tension
Using the quantities derived in Section 2.2.2 we are able to formulate the expressions for the normal surface tension acceleration in analogy to Equation (5), and for the tangential surface tension acceleration in analogy to Equation (6), The sum f S i = f n i + f t i enters the SPH momentum balance (Equation (26)).

SPH Wetting
Following [3,8] we use a modification of the surface normal n i of a fluid particle i in the vicinity of a wall in order to enforce an equilibrium contact angle. To do so, we first extrapolate the surface normal n j of the wall particles j on each fluid particle i, A corresponding tangent is then evaluated by Similarly to Equation (37) the wall contact angle θ j is mapped onto the fluid particles, In many situations the wall contact angle will be a constant but Equation (39) allows to model also graded or sharp variations of the wetting behavior. Next, the distance to the wall is calculated for each fluid particle by This wall distance is used to obtain a normalized weight, Finally, the smooth modification is applied, and n * i is used in Equations (35) and (36) instead of n i .

SPH Particle Shifting
A particle shifting technique as introduced by Sun [24] can be used to homogenize the particle distribution. The position shifting for a particle i inside the material, i.e., not in the surface region, within a time step ∆t is defined as For particles within the surface region and with λ i ≥ 0.4, the position shifting is modified into (I − n i ⊗ n i )∆r PST i where I is the identity matrix. In the vicinity of a wall, n i is replaced by n * i . For surface region particles with λ i < 0.4, no position shifting is applied.

SPH Time Integration
The numerical integration time step is based on the fastest physical mechanism in the simulation. The competing mechanisms are inertia, sound propagation, surface tension, viscous diffusion and thermal diffusion. Time step conditions leading to stable simulations involving these mechanisms are given in [1,4,[39][40][41], respectively, where f max is the magnitude of the maximum acceleration of all particles. It is usually suggested to choose the numerical speed of sound s about 10 times faster than the maximum particle velocity in the simulation [2]. In our visco-capillary simulations, we found that also the inequality should be fulfilled in order to maintain stability. Time integration is done using a velocity Verlet scheme [42]. The density is propagated by The specific thermal energy is propagated by The velocity is propagated by Furthermore, the position is propagated including particle shifting by

Results
A series of benchmark cases are simulated for which analytical solutions are available for comparison. In Sections 3.1-3.7 capillary or thermo-capillary benchmarks are assessed. In Sections 3.8-3.11 the free surface thermal boundary conditions are benchmarked. Sections 3.12 and 3.13 contain more complex cases which provide suggestions for further areas of application of the proposed numerical model.
Unless otherwise stated we are using an ethylene glycol parametrization for the present SPH simulations, i.e., a rest density ρ 0 = 1113 kg m −3 , a viscosity η = 16.1 mPa s and a surface tension σ n = 48.4 mN m −1 . Ethylene glycol is chosen here as an example for a real fluid which is used in many technical applications. If gravity is used, then g = (0, 0, −g) with g = 9.81 m s −2 . Furthermore, the speed of sound is set to s = 10 m s −1 unless stated otherwise and the isentropic exponent of the equation of state is γ = 1. We typically use a ratio of smoothing length and initial particle spacing of h/∆x = 1.5 and accordingly Ξ λ = 1.88 for d = 2 and Ξ λ = 1.82 for d = 3 (compare Appendix A). In Sections 3.5 and 3.6 the viscosity model according to Equation (30) is used. For all other fluidic benchmarks, the viscosity formulation according to Equation (29) is used. For the benchmark cases in Sections 3.5, 3.6, 3.12 and 3.13 artifical viscosity with parameters α = 0.01 and β = 0.02 is used. For all other cases artificial viscosity is not activated.
All simulations are carried out using the SimPARTIX software developed by Fraunhofer IWM [43].

Droplet Oscillation
As a first benchmark case we analyze the oscillation of a circular droplet with an initial divergence-free velocity field [4][5][6], with v 0 = 1 m s −1 and r 0 = 1 mm. The lateral and vertical coordinates are x and y, respectively, and r = x 2 + y 2 . The origin is placed at the center of the droplet. The theoretical oscillation frequency for a 2D droplet is given by The initial radius of the droplet is R = 4 mm and the resolution is varied between ∆x = 500 µm and ∆x = 125 µm which corresponds to 8 and 32 particles per radius. Snapshots from the simulation with the finest resolution are shown in Figure 1. The lateral expansion L x of the droplet as a function of time is shown in Figure 2. The regular oscillation with decay in amplitude due to viscous dissipation is clearly perceptible. A Fourier analysis is carried out in order to assess the oscillation behavior in more detail. The Fourier transform of the droplet expansion, F (L x ), is shown in Figure 3 as a function of frequency f sim . For the coarsest spatial resolution, the oscillation is too slow. However, for the two finer resolutions the main frequency peak is in good agreement with the theoretical prediction given by Equation (52).

Droplet Formation
The second benchmark is the formation of a droplet from an initial square arrangement of SPH particles as proposed by Adami et al. [6]. The edge length of the square is L = 6 mm and resolutions between 15 and 60 particles per edge length are used yielding a particle spacing between ∆x = 400 µm and ∆x = 100 µm. Figure 4 shows a sequence of characteristic simulation snapshots using 60 2 particles. The temporal evolution of the two-dimensional kinetic energy of the droplet for the different spatial resolutions is shown in Figure 5. Following a sharp initial increase of kinetic energy, a decay can be observed for all cases. The decay amounts to about three orders of magnitude for the coarsest resolution while it amounts to more than five orders of magnitude for the finest resolution. The pressure profile in the eventually formed droplet is evaluated using the SPH approximation according to Equation (9). The results are shown in Figure 6. The finer the spatial resolution the better is the pressure jump at the droplet surface approximated. The agreement with the Young-Laplace equation, which relates the pressure P in the droplet to its radius R, also gets better for finer spatial resolutions. The convergence of the numerical solution P sim with respect to the theoretical result P theo for the pressure in the droplet is analyzed in Figure 7. It can be observed that the convergence behavior is significantly better than linear and only slightly worse than quadratic. As a final test, the Young-Laplace scaling with R or L, respectively, is tested. In this case, L is increased by a factor of 10 or decreased by a factor of 0.1 but the particle number of 30 2 is kept constant. The droplet radius of the reference case is referred to as R 0 . Figure 8 shows that the Young-Laplace scaling is well reproduced.

Horizontal Substrate Wetting
The third benchmark assesses the wetting behavior of the model similarly to [8]. An initial circular droplet of particles with a radius of R = 3 mm is placed directly above a planar substrate. The contact angle parameter is varied between θ = 20°and θ = 160°. No gravity is applied is order to compare the results with simple analytical solutions. Gravity, of course, would have an effect on a droplet of this size. The resulting droplet shapes on the substrate are shown in Figure 9 for a particle spacing of ∆x = 100 µm. The predictive quality of these simulations is analyzed by measuring the equilibrium height H eq and base diameter D eq of the resulting droplets and comparing them with the analytical results and D eq = 2 sin θ π R 2 θ − sin θ cos θ .
(55) Figure 10 shows the according simulation results in comparison to the theoretical predictions. The non-trivial shape of both analytical functions is well reproduced up to a contact angle of θ = 100°. Noticeable deviations occur for larger contact angles and the error increases with increasing contact angle.

Vertical Plate Wetting
As fourth benchmark, the wetting of a vertically oriented plate immersed into a fluid is studied. The system is 15 mm wide with periodic boundary conditions in lateral direction. The initial fluid filling height is 4 mm. A vertical plate is immersed into the fluid down to a distance of 1.5 mm from the ground. The acceleration of gravity acts in downward direction and the contact angle is varied between θ = 30°and θ = 150°. A particle spacing of ∆x = 25 µm is used. Figure 11 shows simulation snapshots after the terminal height of capillary rise or depression, respectively, is reached. The accuracy of the simulations is evaluated by comparing the contact angle θ which enters the numerical model via Equations (39) and (42) with the effective, measured wetting angle with respect to the plate. Figure 12 shows that up to θ = 90°the agreement is very good. For larger contact angles deviations become more pronounced similarly to the horizontal substrate wetting case in Section 3.3.

Flow between Parallel Plates
In the fifth benchmark the steady-state flow of a fluid column between parallel plates is analyzed [26]. The length of the column is h 0 = 10 cm and distance of the plates is either L = 0.5 mm or L = 1 mm. The spatial resolution is ∆x = 25 µm. Different pressure profiles and steady-state velocities of the fluid column are enforced by variation of the advancing contact angle θ adv and the receding contact angle θ rec shown in Figure 13. The pressure at the advancing end of the fluid column is given by and at the receding end of the fluid column by An overview of the pressure drop along the fluid column for all simulated cases in comparison to the theory is presented in Figure 14. The obtained pressure profiles are generally in agreement with the analytical predictions. Yet, some finite curvature can be observed in the simulation results while the pressure gradient should be constant in theory.
The average steady-state velocity of the fluid in the channel is v = L σ n 6 h 0 η (cos θ adv − cos θ rec ).
This theoretical prediction is compared to the simulation results in Figure 15. Within the error bars of the simulation, which represent one standard deviation, agreement with the theory is obtained. Note that in order to obtain stable simulations for this benchmark case it is necessary to increase the numerical speed of sound to s = 30 m s −1 . For smaller values of s, fragmentation of the particle arrangement in the vicinity of the concave surface is observed. This phenomenon is discussed in Section 4.

Flow into a Nozzle
The sixth benchmark is similar to the previous one with the difference that the parallel plates are attached to a fluid reservoir resulting in a nozzle geometry. The distance of the parallel plates is either L = 0.5 mm or L = 1 mm and the spatial resolution is ∆x = 25 µm. The capillary action of the surface tension with a contact angle θ = 30°causes the fluid from the reservoir to penetrate the space between the parallel plates. The rise velocity of the fluid decreases over time as shown in Figure 16.
The rise height h r of the fluid column between the plates and the average rise velocity v r are obtained by solving the set of differential equations [44] The according solutions are explicitly given by These theoretical predictions are compared with the data obtained from the numerical simulations in Figure 17. The transient behavior of both height and velocity are satisfactorily reproduced in the simulations although the actual values are systematically underestimated to a small extent. Note that in order to obtain stable simulations for this benchmark case it is necessary to increase the numerical speed of sound to s = 50 m s −1 . The reason for this is the avoidance of fragmentation close to the surface similarly to the case in Section 3.5. The even higher speed of sound in this case is required because of the strong acceleration during the initial phase of the flow into the nozzle.

Marangoni Flow
The seventh benchmark focuses on the effect of a surface tension gradient caused by a temperature gradient [9]. A fluid initially rests in a basin of L = 10 mm width. The left wall of the basin has a constant temperature of T min = 300 K and the right wall a constant temperature of T max = T min + ∆T = 400 K. The fluid temperature initially describes a linear profile interpolating between the values of the walls. The basin bottom is adiabatic. The specific heat capacity of the fluid is c = 1000 J kg −1 K −1 , the thermal conductivity is k = 10 W m −1 K −1 and σ t = −1 mN m −1 K −1 . The particle spacing is varied between ∆x = 100 µm and ∆x = 400 µm. The contact angle with respect to both side walls is θ = 90°. Images from the simulation using the spatial resolution L/∆x = 100 are displayed in Figure 18. The fluid strongly rises on the left wall due to the Marangoni effect which drives the fluid at the surface from the warm region towards the cold region. The last snapshot shows the steady state where an eddy has formed below the surface. The model is quantitatively assessed by evaluating the Marangoni tension at the surface at the beginning of the simulation using different spatial resolutions as well as different ratios of smoothing length and particle spacing. For this purpose, the tangential volumetric surface force (Equation (6)) is evaluated using SPH interpolation and then integrated along a line perpendicular to the surface. The result τ should be equal to |σ t ∇ S T| = 10 Pa. Figure 19 shows that the theoretical result is met for all studied resolutions and smoothing length ratios.

Plate with Double-Sided Heat Transfer
The eighth to eleventh benchmark case do not involve any fluid flow but assess the predictive quality of the thermal boundary conditions at the surface as introduced in Section 2.2.4. A convenient parameter for these cases is the thermal diffusivity, a = k/(ρ 0 c), as it allows to normalize the time and a characteristic length scale in terms of a Fourier number. In all of these four cases the thermal conductivity is k = 1 W m −1 K −1 , the specific heat capacity is c = 100 J kg −1 K −1 and the density is ρ 0 = 7800 kg m −3 .
As eighth benchmark an infinite plate with heat transfer on both surfaces is analyzed. The plate thickness is 2L = 2 mm and its initial temperature is T 0 = 300 K. The heat transfer coefficient is f = 1 × 10 4 W m −2 K −1 and the ambient temperature is T a = 400 K.
The analytical solution for this case is given by the series where γ n is the nth root of the implicit equation The solution is calculated up to n max = 500 and compared with the simulation results for different spatial resolutions in Figure 20. While for the lowest spatial resolution the transient temperature is overestimated, differences between theory and simulation vanish for a sufficiently fine resolution.

Plate with Double-Sided Heat Flux
The ninth benchmark is similar to the eighth one with the difference that instead of heat transfer a heat flux q = 1 × 10 4 W m −2 is prescribed at both surfaces of the plate.
The solution for this situation is given by with γ n defined by γ n = n π.
Again, the series is evaluated up to n max = 500 and simulation and theory are compared in Figure 21 for different spatial resolutions. Perfect agreement is found independent of the actual resolution.

Plate with Heat Transfer and Heat Flux
The tenth benchmark is a combination of the two previous ones. An infinite plate of thickness L = 1 mm experiences heat transfer with f = 1 × 10 4 W m −2 K −1 and T a = 300 K at x/L = 0 and heat flux with q = 1 × 10 4 W m −2 at x/L = 1. The initial temperature of the plate is T 0 = T a .
The analytical series solution is given by cos(γ n x/L) + f L k γ n sin(γ n x/L) 2 sin(γ n ) exp(−γ 2 n a t/L 2 ) γ n + sin(γ n ) cos(γ n ) , with the roots γ n given by f L k = γ n tan(γ n ).
A comparison between simulation and theory (n max = 500) for different spatial resolutions is provided in Figure 22. For too coarse resolutions and Fourier numbers ta/L 2 1 the temperature is underestimated by the simulation. Yet, similarly to the benchmark case with double-sided heat transfer (Section 3.8), the analytical solution gets better approximated with finer spatial resolution.

Rod with Heat Flux
In the eleventh benchmark an infinite rod of radius R = 1 mm and initial temperature T 0 = 300 K is heated by the surface heat flux q = 1 × 10 4 W m −2 .
The according analytical solution for the temperature distribution along the radial coordinate r is given by with the nth root γ n given by Due to the rotational symmetry of this benchmark case, the Bessel functions of the first kind J 0 and J 1 appear in the solution. The comparison between the simulation results and the analytical solution in Figure 23 reveals that the temperature is underestimated for too coarse spatial resolutions. However, for finer resolution the analytical result is well predicted. For this case, the accuracy of the surface heat transfer in the simulation is assessed by evaluating the rate of change of thermal energy in the simulation. As heat transfer is the only active heating mechanism, the rate of thermal energy change should be directly proportional to the applied heat flux. An according convergence analysis is plotted in Figure 24. The effective heat flux in the simulation converges approximately linear with the spatial resolution.

Droplet Deformation and Splitting
The twelfth case is intended to study the usage of the presented model in a more complex three-dimensional situation. A droplet of fluid with a diameter of 1 mm is discretized using ∆x = 20 µm. The droplet floats above a substrate initially. Gravity pushes the droplet onto the substrate which consists of hydrophilic regions with a contact angle of θ = 45°and hydrophobic regions with a contact angle of θ = 135°. Two different substrates are used. One is composed of square regions having a certain contact angle. The other substrate consists of striped regions with a width of 400 µm. The numerical speed of sound in the simulation is s = 3 m s −1 . Figure 25 shows the temporal evolution of the fluid droplet on both substrates. For the square regions, the wetting and de-wetting behavior leads eventually to a splitting of the droplet into two. For the striped regions, the droplet does not split up but ends up with a very peculiar shape of its perimeter.

Thermo-Visco-Capillary Droplet Coalescence
The thirteenth and final case contains the most complex simulation setup. Two droplets of fluid with diameters of 1 mm are discretized using ∆x = 20 µm. The droplets touch a substrate below as well as each other initially. Gravity acts normal to the substrate surface. The fluid has an initial temperature of T min = 300 K and the substrate a constant temperature of T max = T min + ∆T = 500 K. The specific heat capacity of the fluid is c = 1000 J kg −1 K −1 and the thermal conductivity of both fluid and substrate is k = 10 W m −1 K −1 . The surface tension of the fluid is withσ n = 48.4 mN m −1 and σ t = −0.2 mN m −1 K −1 . The contact angle with respect to the substrate is θ = 45°. The ambient temperature is T a = T min and the fluid is cooled by surface heat transfer with f = 1 × 10 2 W m −2 K −1 and by radiation with ε = 0.5. The viscosity of the fluid is temperature-dependent in the form of an Arrhenius relation,

Discussion
The presented consistent modeling of thermo-capillary effects and thermal boundary conditions within the SPH method could be verified and validated by the extensive analyses in Section 3. In the following, we would like to address some details which could be highlighted by the analyses.
For free droplets, the visco-capillary oscillation frequency and the Laplace pressure are adequately reproduced by the simulations with sufficient spatial resolution. For the Laplace pressure, a near quadratic convergence with spatial resolution has been found. At the finest resolution used, the relative error is less than 10 −2 , which is typically an acceptable tolerance for SPH simulations [2].
The investigations on static wetting have shown that contact angles up to θ ≈ 100°a re well reproduced in the simulation. Larger contact angles are increasingly poorly reproduced. A comparable observation was made by Breinlinger et al. [8]. This is apparently a systematic weakness of the submodel for imprinting the desired contact angle in Section 2.2.11. The reason for this is that in the approach used to determine the free surface and the quantities derived from it, the fluid and the substrate are initially considered as a common phase. As a result, it is not possible to identify free surfaces that are only separated by small gaps. However, precisely this situation exists in the case of large contact angles. Alternatively, the fluid and the substrate could be considered as distinct phases from the beginning. In that case, however, the reproduction of small contact angles would become increasingly inaccurate, which has been observed in corresponding tests.
In the dynamic wetting scenarios, a satisfactory overall agreement of the simulation with theoretical predictions for the pressure and velocity has been observed. A slight, systematic deviation in the dynamic wetting case in Section 3.6 can be explained by the fact that the momentum of the reservoir is not considered in the respective balance equation (Equation (59)). The strong tensile stresses caused by the surface tension in the dynamic wetting scenarios are seen as the reason why it is necessary to slightly increase the numerical speed of sound in these cases. Without such an increase, fragmentation of the particles near the surface has been observed. In this context, the used approach of a surface delta function based on the minimum eigenvalue λ of the inverse renormalization tensor has proved to be helpful. As shown in the Appendix A, the surface delta function is more evenly distributed among the particle layers on the surface than when using a delta function which is based only on the kernel gradient. This attenuates the tensile forces acting on the individual particles, which helps to stabilize the simulation. The use of a weak artificial viscosity and the viscosity formulation of Sigalotti et al. [38] have also proved to stabilize the dynamic wetting scenarios.
A fundamental strength of the presented model is that only one phase has to be modeled. This saves computational time, since especially for systems consisting of a fluid and a gas, the low density of the gas requires a very small time step, as shown by Colagrossi and Landrini [45]. However, even the presented model may require a relatively small time step due to a high numerical speed of sound, as described in the previous paragraph. For such a situation, the use of an implicit SPH formulation can be useful, such as the ones used by Nair and Pöschel [31] or Fürstenau et al. [20] in the context of capillary phenomena.
The investigations on the heat transfer boundary condition have shown that agreements with analytical results are obtained by simulations with sufficiently fine spatial resolution. In the case of a pure heat flux boundary condition on a flat surface, excellent accuracy of the numerical simulations has been found regardless of the spatial resolution. This can probably be attributed to the fact that the heat input is independent of the temperature computed in the numerical simulation. In the case of a curved surface, a linear convergence of the relative error of the heat flux with the spatial discretization could be determined. This convergence behavior can be considered a measure of how well the surface delta function describes a curved surface.
The two cases in Sections 3.12 and 3.13 are not intended to validate the model but to provide ideas for application in more complex scenarios. Use cases for thermo-capillary SPH simulations offer, for example, welding processes (Hu and Eberhard [46]) or additive manufacturing processes (Russell et al. [47], Fürstenau et al. [48], Meier et al. [49]). The presented method has already been applied by the author and coworkers to derive universal process maps for laser powder bed fusion of polymers [50].
Based on the identified weakness of the presented model, a more accurate description of large contact angles is considered as a future field of development. In addition, the modeling of dynamic wetting phenomena can be improved by considering the dynamic contact angle as in Huber et al. [11]. A fundamental problem for single-phase models based on the CSF approach is the lack of possibility to formulate the surface tension forces in a momentum-preserving way. As a consequence, sooner or later a non-physical drift of free-floating droplets or of droplets on flat substrates occurs in the simulations. A solution to this problem would be highly desirable in terms of a more universal applicability of the method.

Conclusions
The main contributions of the present work are summarized in the following list: • The surface delta function, surface normal and surface curvature are derived in a consistent manner from the renormalization tensor in single-phase SPH simulations. • Thermo-capillarity, wetting and thermal boundary conditions are formulated in a uniform way based on these surface properties. • The quantitative agreement of corresponding simulations with a large number of analytical test cases demonstrates the validity of the presented approach. Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Surface Delta Function
Different approaches exist in order to evaluate a suitable surface delta function. Equations (19) and (20) are two possibilities for this purpose. We consider using Equation (19) to be a consistent choice, because this way the surface normal (Equation (18)) and the surface delta function are derived from the same quantity, namely λ. However, as other authors [14,20,23] use Equation (20) to obtain the surface delta function, we compare these two approaches in the following.
The half space x < 0 is filled with SPH particles at positions x n = (0.5 − n)∆x, n ∈ N. In Figure A1 the SPH interpolations of λ, and of unity, are plotted as functions of a coordinate |r| perpendicular to the surface for different ratios of smoothing length and particle spacing. Both for λ and unity the functions drop from 1 to 0 in the vicinity of the surface. Yet, there are subtle differences. The descent appears to be steeper for the case of unity. Furthermore, the inflection point is located exactly at |r| = 0 for the unity field while it is shifted somewhat towards the half space filled with particle for the λ field. Figure A1. Functions that are approximately 1 inside and 0 outside the particle distribution for different ratios of smoothing length and particle spacing. The 3D Wendland kernel is used as an example. The particle positions are indicated by the circles. (a) Minimum eigenvalue of the inverse renormalization tensor (Equation (A1)). (b) Unity field (Equation (A2)).
The functional forms of the corresponding gradients are obtained by applying Equation (9) to either |∇λ i | defined by Equation (17) or to |∇1 i | defined by Equation (20). The results are plotted in normalized form in Figure A2. The overall shape of the gradients for both fields is similar. While the gradients are symmetric for the unity field, they are a bit skewed for the λ field. In addition, finite values for the gradients are more localized at the surface for the unity field. As to be expected, the smaller the ratio h/∆x the more localized are the gradients of both fields.
The skewedness in case of the λ field is another reason why we prefer it over the unity field. Due to this property the magnitude of the delta function is more evenly distributed over the two layers of surface region particles. The delta function of the unity field is more localized on the outermost particle layer. We found that it is beneficial for the stability of the numerical simulations if the action of the surface tension is more evenly distributed over two particle layers. In order to obtain the normalization factors Ξ λ in Equation (19) and Ξ 1 in Equation (20) the gradient magnitudes |∇λ i | and |∇1 i | are summed over a column of particles i perpendicular to the surface similar to the in [23]. The inverse sums yield the normalization factors shown in Figure A3 for four different kinds of SPH kernel functions. For the case of the λ field, the normalization factor increases with the ratio h/∆x and appears to reach a plateau for large ratios. For the case of the unity field, the normalization factor oscillates around 2 for ratios h/∆x < 1.5 and is close to 2 for larger ratios. As mentioned in Section 2.2.2, we restrict all surface related calculations to those particles which are no more than 1.5 ∆x away from the free surface particles. Consequently, the normalization factors must reflect this restriction. In Figure A4 the corresponding results are presented. The normalization factors are in good approximation linear functions of the ratio h/∆x for the case of the λ field. For the unity field the normalization factors deviate systematically from 2 and increase with h/∆x if this ratio is larger than ∼1.4. The linear function to approximate the normalization factor for the λ field as plotted in Figure A4 is given by The slope and intercept values for different SPH kernel functions in two and three dimensions are listed in Table A1. Table A1. Parameters for the Ξ λ approximation function (Equation (A3)).