Simulating complex fluids with Smoothed Particle Hydrodynamics

PH669 ABSTRACT Complex fluid dynamics encompasses a large variety of flows, such as fluids with non-Newtonian rheology, multiphase and multi-fluid flows (suspensions, lather, solid/fluid interaction with floating objects, etc.), violent flows breaking waves, dam-breaks, etc.), fluids with thermal dependencies and phase transition or free-surface flows. Correctly modeling the behavior of such flows can be quite challenging, and has led to significant advances in the field of Computational Fluid Dynamics (CFD). Recently, the Smoothed Particle Hydrodynamics (SPH) method has emerged as a powerful alternative to more classic CFD methods (such as finite volumes or finite elements) in many fields, including oceanography, volcanology, structural engineering, nuclear physics and medicine. With SPH, the fluid is discretized by means of particles and thanks to the meshless, Lagrangian nature of the model, it easily allows the modeling and simulation of both simple and complex fluids, simplifying the treatment of aspects that can be challenging with more traditional methods: dynamic free surfaces, large deformations, phase transition, fluid/solid interaction and complex geometries. In addition, the most common SPH formulations are fully parallelizable, which favors implementation on high-performance parallel computing hardware, such as modern Graphics Processing Units (GPUs). We present here how GPUSPH, an implementation of the SPH method that runs on GPUs, can model a variety of complex fluids, highlighting the computational challenges that arise in its applications to problem of great interest in volcanology.


Introduction
The ability to accurately model the behavior of fluids under the most diverse conditions is extremely important in both scientific research and applied industry, and has led to the development of Computational Fluid Dynamics (CFD).
Traditional approaches to CFD have been based on an Eulerian description of the fluid, with fixed reference grids (finite differences, finite volumes).These approaches are often computationally very efficient, but may have issues with complex geometries, moving boundaries or the tracking of the free surface and of internal interfaces (such as the boundaries between fluids in multi-fluid simulations, or the front in phase transitions).Lagrangian, mesh-based approaches such as finite elements can overcome some of these issues, but their application to highly dynamic flows such as dambreaks is limited due to the possible fragmentation of the flow and the computationally expensive need for remeshing when the fluid undergoes large deformations.
These methods have been used to model and study geophysical flows under specific conditions, such as lava tubes [Filippucci and Dragoni 2013], channel flows [Costa and Macedonio 2005], or other simple geometries [Fernandez Nieto and Narbona-Reina 2016].In some cases, it is also possible to find analytical solutions to the equations governing the flow [Tallarico and Dragoni 1999], which can then be used in models still relying on mesh-based approaches [Cappello et al. 2016].However, many geophysical flows present characteristics, such as multiple fluids or phases, non-Newtonian or temperature-dependent rheology, violent or otherwise highly dynamic behavior, which are challenging to model with the classical approaches.
An alternative to mesh-based method has recently emerged in the form of Lagrangian, mesh-free approaches.Among these, Smoothed Particle Hydrodynamics (SPH) has found a large number of applications in many fields [Monaghan 2005] since its inception in 1977.
In SPH the problem is discretized by means of particles which act both as interpolation nodes and to carry physical information (mass, density, temperature, velocity, etc.) about a small volume of the fluid.The evolution of the properties of each particle is therefore determined by a discretized form of the physical equations (Navier-Stokes equations for the motion, thermal equation for the temperature, etc.), so that the properties of the fluid at any point in space can be computed by a smoothed average of the properties of the neighboring particles.
For the simulation of complex fluids, SPH presents a number of advantages over traditional mesh based methods, such as the implicit tracking of the free surface and of any internal interface (such as solidification fronts) and the lack of restrictions in the geometry of the fluid and of the containing geometries.The main downsides are lower accuracy and, for the classical formulation, the use of a weakly-compressible rather than a fully incompressible fluid model, such that the speed of sound can significantly limit the time-step size then used with explicit time-stepping schemes.On the other hand, the standard weakly-compressible SPH formulation has the benefit of being completely parallelizable, since the evolution of each particle at each time-step can be computed independently from that of the others, which makes SPH quite suitable for implementation on massively parallel hardware, such as modern Graphics Processing Units [Hérault et al. 2010].
In this paper we present some applications of GPUSPH [Hérault et al. 2011, Bilotta et al. 2016], a GPU implementation of the SPH method specialized for the modeling of complex fluids in three dimensions.We introduce first the equations governing the fluid motion and heat exchange, the approach used to discretize them using SPH, and finally some applications including validation against standard analytical problems and qualitative results in more sophisticated examples taken from applications of GPUSPH to both industrial topics and geophysics, with a particular attention to volcanology.In the application section we will highlight the computational challenges faced by the model, and discuss the possible approaches to resolve them.

Notation
Throughout this manuscript we use the following conventions: -Vectors will be written in bold; -The reference coordinate system is Cartesian orthogonal with axes x 1 , x 2 and x 3 ; -Latin indices denote coordinates; -Greek indices denote particles; -Physical quantities and mathematical expressions follow the common SPH notation, summarized in Table 1.

Governing equations of fluid flow and heat transfer
We model the dynamics of fluids according to the Navier-Stokes equations for continuity of mass and forces balance: (1) (where ρ is the density, u the velocity and D/Dt the total derivative) and (2) (where P is the pressure, μ the dynamic viscosity, and F represents external forces, such as gravity).The fluid is assumed to be weakly-compressible, so that the pressure can be obtained directly from the density, using Cole's equation of state [Cole 1948, Batchelor 1974]: (3) where ρ 0 is the at-rest density, γ is the polytropic constant and c 0 the speed of sound.A weakly-compressible regime is achieved if c 0 is at least an order of magnitude higher than the maximum velocity experienced during the flow.The thermal evolution is described by the heat equation: where T is the temperature, c p the specific heat at constant pressure, and k the thermal conductivity.
COMPLEx FLUID WITH SPH

3.1.Viscosity models
Equation 2 only applies to Newtonian fluids, for which the shear stress τ is proportional to the strain rate ε ⋅ .The dynamic viscosity is the coefficient of proportionality, and it may depend on local physical and chemical properties of the fluid (such as the temperature, density, composition, etc.), but not on the shear stress or strain rate themselves.
A fluid where the relationship between stress and strain rate is not linear is known as non-Newtonian fluid, and many fluids of great interest in industrial applications, biology and geophysics are indeed non-Newtonian.Among non-Newtonian fluids, there is a large subset of fluids (called generalized Newtonian fluids), for which it is still possible to use the same form (Equation 2) of the Navier-Stokes equations, provided that the dynamic viscosity μ is replaced by an effective viscosity function The effective viscosity may still depend on other physical properties of the fluid.For our model we consider the Herschel-Bulkley rheology, a generalized Newtonian rheology characterized by a yield strength τ 0 , a power law exponent n and a consistency index k such that and (5) otherwise.As special cases of the Herschel-Bulkley rheology, we obtain: Newtonian fluids for τ 0 = 0, n=1 (in which case k is the viscosity); power-law fluids for τ 0 = 0, n≠1; Bingham fluids for τ 0 ≠ 0, n=1.
For an Herschel-Bulkley fluid, the effective viscosity can be written [Zhu et al. 2005] as.

Mathematical basis of Smoothed Particles Hydrodynamics
To solve numerically the Navier-Stokes and thermal equations seen in section 3, we discretize their right-hand sides using the Smoothed Particle Hydrodynamics method.In our presentation of the method and its application to fluid dynamics we follow Liu and Liu [2003] and Monaghan [2005].

SPH discretization of fields
Let us consider a field f defined on a domain Ω.By definition of Dirac's delta distribution δ, the value assumed by f at any location x _ ∈ Ω can be expressed, with typical abuse of notation, as (8) Let us now approximate Dirac's delta with a family of functions W( .,h), parametrized by the length h, and satisfying ( 9) and ( 10) where the limit is to be taken in the sense of distributions.
These W functions are the smoothing kernels and the parameter h is the smoothing length.By substituting the smoothing kernel into (8) we get an initial approximation of f as ( 11) We can further approximate the integral with a summation of a finite set of particles x 1 ,…,x α ,…,x N with volume V α , obtaining: where the summation is extended to all particles.However, if the W are chosen with compact support, the summation can be extended only to the particles in a small neighborhood of x _ .From equations ( 11) and ( 12), we observe that the SPH approximation has two main sources of error: 1. approximation of the Dirac's delta with the smoothing kernels; 2. discretization of the domain by means of a finite set of particles.The first approximation vanishes as h→0, the second error is controlled by the average inter-particle spacing Δp and vanishes for Δp→0.Additionally, Δp must tend to zero faster than h to ensure consistency for the method.Most applications adopt a fixed ratio h/Δp usually in the range [1.3,1.5].

SPH discretization of gradients
The field discretization with SPH can be used to reconstruct a field at any point given its values on a set of particles, but the same principle can also be used to discretize spatial gradients.To this end, we will assume further that the smoothing kernels W have compact support and radial symmetry, so that they only depend on |xx _ |.Let us consider the vector field constituted by ∇f(x); recalling equation ( 11) we can write (13) By applying Green's theorem, we rewrite the righthand side obtaining (14) where Σ=∂Ω is the boundary of the domain and n its normal.Given the compact support for W, the first integral in ( 14) is zero if x _ is far from the boundary; additionally, by symmetry of W, we have ), so that the discretized version can be written as (15) This expression allows us to compute (an approximation of ) the gradient of the field f without having its analytical expression, relying instead on the gradient of the smoothing kernel.When applied to a particle β, the equation takes the form: (16) where W αβ =W(x αx β , h).This can be symmetrized (thus helping preserve conservation properties of the analytical equations) by subtracting the gradient of the function identically equal to f(x β ), obtaining: (17) where Finally, we observe that due to the symmetry of W, its gradient can be written as: (18) where It is therefore convenient to choose a kernel such that ( 19) has an analytical expression, and given

SPH discretization of the fluid equations
We can apply the SPH discretization to the continuity Equation (1), obtaining (20) For the momentum (Equation 2) and thermal (Equation 4) equations, we additionally need a discretization for the Laplacian, for which we follow Brookshaw [1985], Morris et al. [1997] and Cleary and Monaghan [1999].The momentum equation then takes the form: where μ ¯αβ is the harmonic mean of μ α and μ β , and g = F/ρ, while the thermal equation becomes: (22 For non-Newtonian fluids, the viscosity μ β of each particle is computed before the actual forces, using the specific rheological law, from the strain rate: (23) with the (∇u) β tensor computed in SPH form as: (24) The evolution of temperature, density and velocity for each particle is finally obtained by numerical integration of (20), ( 21) and ( 22) over time.GPUSPH uses a predictor-corrector time integration scheme, which can be described as follows, for a given time-step ∇t: 1. (in the non-Newtonian case) compute viscosities 2. compute accelerations and temperature derivatives

(in the non-Newtonian case) compute intermediate viscosities
5. compute corrected accelerations a (n*) = a(x (n*) , u (n*) , ρ (n*) , T (n*) , μ (n*) ), density derivatives and temperature derivatives *) , u (n*) , ρ (n*) , T (n*) ), 6. compute new positions, velocities, densities and temperatures: Adaptive time-stepping can be used, in which case the maximum allowed time-step is computed based on constraints, similar to the Courant-Friedrichs-Lewy conditions, derived from the sound speed, the forces magnitude and the viscous and thermal diffusivity: Some consideration can be done about applications involving highly viscous fluids, where the viscous stability condition may become predominant.The dependence on the square of the smoothing length causes in fact a stronger reduction of the time-step as the resolution increases.This constitutes a relevant problem in simulating geophysical flows, like lava, where large values of viscosity are commonly involved, causing simulation times to be very long.

GPUSPH applications and tests
We now show some applications of the GPUSPH platform, illustrating the wide range of possible applications and the benefits of the method.We will also discuss the computational performance of the model in the different contexts.All simulations were run with a smoothing length h=1.3∆p, which gives about 80 neighbors per particle.The kernel employed is the 5th order smoothing kernel proposed by Wendland defined W(r, h) = W (r, h), general space the form (26) constant in three dimensions is (27) The Wendland kernel has an analytical expression for F = 1/r (∂W/∂r), given by F (r, h) = F ~ (r, h), which in d dimensions is (28) the normalization constant is given by ( 29) In what follows, particular emphasis will be given on discussing simulation performances in reference to the number of particles, the allowed time-step length, and the resulting ratio between simulation run-time and simulated time.The simulations were all run on one or more NVIDIA TITAN x Maxwell GPUs, as discussed in the details of each problem.

High dynamic flows: Dam Break
The dam-break problem reproduces the breaking of a dam, a large body of water that suddenly gets in motion under gravitational effect and impacts the obstacles it encounters (Figure 1).This example shows the ability of GPUSPH to model highly dynamic flows, and makes it possible to appreciate the detail in the reproduction of the free surface, including splashes and droplets.The fluid density is ρ=1000 kg/m 3 and the dynamic viscosity is μ = 10 -4 Pa s.The simulation was run using two GPUs, showing the capability to distribute the computation across multiple domains, even with domain partitions not parallel to the computational axes [Rustico et al. 2012[Rustico et al. , 2014]].With a resolution of 170 particles per meter, about 500,000 particles are involved in the simulation, resulting in a time-step of 10 -4 s, and a time ratio of 1:100 (one second of simulated time takes about 100 seconds of simulation).

Non-Newtonian Fluids
As discussed in subsection 3.1, GPUSPH is able to reproduce any rheology based on the Herschel-Bulkley model.We present here a validation of the implementation reproducing the Plane Poiseuille flow for a Bingham fluid.The distance between the top and bottom plane is 1m, and the domain is periodic in the flow and transverse directions (Figure 2).The employed fluid has consistency index k=0.1 Pa s, density ρ = 1 kg/m 3 and yield strength τ 0 =0.01 Pa, and is driven by an external force conferring an acceleration a = 0.05 m/s 2 result- ing in a Reynolds number of 0.225.The convergence of the numerical method (as discussed in section 4) is illustrated in Figure 3, showing the convergence of the velocity profile (Figure 3a) and the decreasing error (Figure 3b) with increasing resolution.The performance at different resolutions is shown in Table 2.We remark that the time-step in this case is limited by the viscosity, and thus decreases with the square of the linear resolution.This leads to a significant worsening of the ratio of simulated time to runtime, as discussed at the end of section 5, which in the worst case reaches here as low as 50 minutes to simulate 1 second.

Thermodynamics: natural convection
We tested the interaction of the thermal and mechanical models by simulating a thermal convection.We consider a box containing the fluid, with adiabatic walls, a heated bottom plate and a cooled top plate.As effect of the heating, the fluid gets in motion causing the inception of a Rayleigh-Bénard convective cell.
Additionally, the thermal model can be paired with the momentum equation by introducing temperaturedependent rheological parameters, as in Kaddiri et al. [2012], which alters the behavior of the convective cell as shown in Figure 4.The bottom plate is at temperature T b = 100 K and the top plate at T t = 0 K. Consider-ing a Pearson constant m = 10, the kinematic viscosity at temperature T = T b is ν b = 20m 2 /s while at temperature T = T t is ν t ≈ 10 -3 m 2 /s.The simulation is stable despite the large ratio (five orders of magnitude) between the minimum and maximum viscosity.
Again, the time-step is limited by the maximum viscosity, so that with a resolution of 33 particles per meter, resulting in about 55,000 particles and timestep in the order of 10 -5 s, the simulated to run-time is nearly 1:400.

Multi-Fluid and multi-phase flows
GPUSPH allows the simultaneous simulation of several fluids, each one with its own rheological and thermal properties.SPH is quite appropriate for these kind of problems, since interfaces between two fluids are implicitly tracked.

Multiple liquids
Figure 5a illustrates a multi-fluid problem based on a lava-lamp: two fluids with opposite conditions in terms of density, rheology and thermal expansion coefficient are contained within an adiabatic box with heated bottom and cooled top.The red fluid has density ρ 1 =2 kg/m 3 , Bingham rheology with τ 0 = 0.001 Pa, consistency index depending on the temperature as k(T) = exp(-2T) and thermal expansion coefficient α 1 = 0.6 K -1 , while the blue one has ρ 2 = 1 kg/m 3 , Newtonian rheology with kinematic viscosity μ 2 =1 m 2 /s and α 2 = 0.1 K -1 .At a resolution of 33 particles per meter, about 30,800 particles are involved with a timestep in the order of 10 -4 s, and using a single GPU the obtained time-ratio is 1:10.

Liquid/Gas interaction
For multi-phase flows with liquid and gas phases, the entire system can be modeled using SPH.An example with an air bubble rising in water (1:1000 density ratio) is illustrated in Figure 5b, showing the fracture of the bubble surface towards the end of the ascension.As explained in detail in Bilotta [2014], for liquid/gas simulations where the compressibility does not play a major role, it is necessary to give the less dense fluid a higher sound-speed, so that both fluids are at the same quasi-compressible regimen.Due to this phenomenon, such multi-fluid simulations suffer from a time-step limitation dictated by the soundspeed of the less dense fluid.In the raising bubble example, a simulation with 16 particles in the bubble diameter requires over 500,000 particles, and the resulting time-step is in the order of 10 -6 s.The timeratio over 1:7400.It should be noted that in this case the time scale of the entire simulation is short, so that the impact of the high time-ratio is less significant than in larger-scale simulations.

Liquid/Solid interaction
For the fluid/solid interaction, GPUSPH computes the forces and torque acting on each moving object from the interaction of the boundary particles modeling the object with the neighboring fluid particles.The actual object dynamics, including object/object and objects/terrain interaction, is then handled by the ODE library (Smith, 2007).An example application is shown in Figure 5c, with floating tree trunks transported on the spillway of the Goulours Dam, on data provided by EDF [Bilotta et al. 2014].Water density is ρw=1000 kg/m 3 , whereas ρ t =0.7ρ w .Like in the standard dam-break case, this is a highly dynamic problem, but we now introduce a strong interaction with moving objects (the trunks) and complex geometries (the irregular floor of the spillway).The entire simulation consists of over 2,500,000 particles (inter-particle spacing of 0.2 m), and the time-step (controlled by the sound speed of water) is in the order of 10 -4 s.With a single GPU, we achieve a time ratio of about 1:500.This is due partly to the large number of particles, and partly to the use of the semi-analytical boundary model [Mayrhofer et al. 2015] for the fluid/terrain interaction, which requires a larger influence radius and thus a higher number of neighbors per particle.isons between numerical solutions and experimental data [Dietterich et al. 2015] allow both improvements to the numerical model and in an inversion process to better determine which kind of rheology better reproduces the experiment.The dynamical viscosity of the lava depends on the temperature following the Vogel-Fulcher-Tammann equation [Cordonnier et al. 2015], with coefficients determined experimentally [H.Dietterich, private communication]: ( where T is the temperature expressed in Kelvin and the density is ρ = 2350 kg/m 3 .Thermal radiation is modeled as described in Bilotta et al. [2016], using heat capacity c p =1600 J/( kg•K) and thermal emissivity ϵ=0.96.
A further step in our ongoing research is shown in Figure 7, illustrating a simple case of lava-water interaction.In this case, the sharp decrease in temperature caused by the contact with water results in a different emplacement compared to the lava flow behavior in absence of water.In our preliminary tests, simulations have been conducted with a viscosity that is 10 times lower than the experimental viscosity, due to computational constraints.Indeed, the main limitation in performing lava simulations is constituted by very long simulation times due to the high viscosity of lava, especially in the cooling phase.The high values of viscosity can determine time-steps in the order of 10 -7 s or lower.
Considering a resolution of 16 particles per meter, each of the simulations above involved about 300,000 particles, with an effective time ratio of over 1:150,000 (about two days for one simulated second) running on a cluster with four GPUs.Using the experimental law would have led to nearly one month of computational time to simulate a single second.

Conclusions and future work
We have shown how GPUSPH constitutes a very versatile fully three-dimensional platform for CFD, enabling the simulation of natural flows of Newtonian and non-Newtonian fluid with constant as well as temperature-dependent rheology.
The fully explicit nature of the SPH numerical method allows efficient parallelization on current high-performance parallel computing hardware such as GPUs, and can even scale to multiple GPUs.On the other hand, the explicit integration scheme can become a limiting factor in the simulation of highly viscous fluids, for which the time-step decreases with the square of the spatial resolution, quickly leading to very long simulation times.
In the application of GPUSPH to lava flows this can be of crucial importance, particularly in the cooling phase, where the viscosity of lava can grow by more than two orders of magnitude.
A solution to this issue, which we are currently working on, is the use of a semi-implicit integration scheme, in which the inviscid part of the Navier-Stokes equation is integrated with the explicit scheme currently in use, whereas the viscous contribution is computed with an implicit scheme.While this incurs the need to solve a large, sparse implicit system, it provides significant benefits in computational times, while also improving the stability of the numerical method.
c β is the sound speed at density ρ β and the C 1 , C 2 , C 3 , C 4 stability constants in GPUSPH take the values C 1 = C 2 = 0.3, C 3 = 0.125 and C 4 = 0.1.The effective time-step is then computed as ∆t=min β ∆t β .

Figure 1 .
Figure 1.Dam break with obstacle.Particles are colored by device.

Figure 5 .
Figure 5. Multi-phase flows.a) Lava lamp with particles colored by fluid type.b) Rising bubble as example of gas-fluid interaction.c) Spillway of the Goulours dam as example of interaction.

Figure 7 .
Figure 7. Lava flows simulations with GPUSPH.a) Lava-water interaction.b) Comparison between lava-water interaction and dry emplacement.

Table 2 .
Performance of plane Poiseuille flow for a Bingham fluid.