Scaling laws for implicit viscosities in smoothed particle hydrodynamics

Smoothed particle hydrodynamics (SPH) is a particle-based method which solves continuum equations such as the Navier-Stokes equations. A periodic fluidic system under homogeneous shear is studied using SPH in the present work. The total pressure of the system and the shear stress contributions from the SPH interaction terms for pressure and viscosity as well as the contribution caused by velocity fluctuations are analyzed. It is found that the pressure and the shear stress contributions obey certain scaling laws depending on physical properties of the system such as compressibility, viscosity and shear rate as well as the spatial resolution. Some of the identified relations resemble scaling laws for the rheology of dense granular flows. These findings render an assessment of the convergence with respect to the spatial resolution of SPH simulations possible. Furthermore, the similarities between numerical SPH particles and physical grains in dense flow provide a deeper understanding of the nature of the SPH method.


Introduction
Since its conception in 1977 [1] the SPH method has been used in a wide range of scientific and industrial applications due to its simplicity to implement and its flexibility in modelling (multi-phase) fluids with free surfaces [2,3]. Despite the successful application of the method, it is still lacking a theoretical proof of its convergence with respect to the spatial discretization, i.e., the SPH particle size.
Several studies focused on the accuracy of SPH shear flow results. Posch et al. [4] described implicit SPH viscosities, i.e., artificial shear stress contributions which have their origin in the particle-based nature of the method. Ellero et al. [5] investigated this phenomenon further and quantified a kinetic contribution which acts as an average Reynolds stress and a potential contribution without continuum analog. Both contributions can become larger than the specified shear viscosity of the system for low viscous flows. Previous work of the present authors [6] showed that the magnitudes of the implicit viscosities depend on the used SPH formulation of the Navier-Stokes viscosity term. Imaeda et al. [7] reported large density errors which were attributed to the local shear deformation of the neighborhood of each SPH particle.
In the present work we analyze the resulting density, pressure and shear stress of SPH shear flow simulations from a different perspective, i.e., by interpreting SPH particles as grains in dense granular flow. This approach enables the mapping of several relations which are valid for granular rheology to SPH simulations and, thus, a detailed quantification of the implicit viscosities. e-mail: claas.bierwisch@iwm.fraunhofer.de

SPH method
The motion of an incompressible Newtonian fluid in the absence of body forces is described in the moving frame of reference by the Navier-Stokes momentum equation with D/Dt as the substantial derivative, v as the velocity, ρ as the density, p as the pressure and η 0 as the dynamic viscosity. Using the SPH summation formalism [8], the value of the quantity A at the position r is given by the summation over the contribution A j from each particle j at position r j with fixed mass m and variable density ρ j within range of the kernel function in three dimensions, where h is the smoothing length. We are using the common cubic spline kernel function, Let r i j = r i − r j , r i j = |r i j |, W i j = W(r i j , h) and ∇W i j = r i j (d f (q)/dq)/(r i j h 4 ). By applying the SPH summation formalism, Eq. (2), the accelaration of particle i due to the pressure term of the momentum equation can be written as yielding a conservative (potential) force f pot i j between the particles. A weakly-compressible approach is used for closure by calculating the particle pressure from the equation of state, with ρ 0 as the reference density, γ as the isentropic exponent and B as a compressibility which limits density fluctuations. The particle mass and the reference density are related by the initial particle spacing s via m = ρ 0 s 3 . The SPH summation formalism yields the density of an SPH particle i via ρ i = m j W i j . For the viscous term of the momentum equation the form proposed by Cleary et al. [9] is used which is based on dissipation due to normal SPH particle collisions, ξ is a dimensionless coefficient in the order of 20.
= 0.01 is used to prevent zero division. Equation (6) yields a dissipative force f diss i j between the particles. The explicit time integration with time step Δt of the total force on particle i, , is based on the Velocity Verlet scheme. To ensure numerical stability a time step condition is enforced, Δt ≤ min(0.25 h/c s , 0.05 h 2 ρ 0 /η 0 ), with the speed of sound c s = √ (γB/ρ 0 ).

Simulation setup and analysis
The simulation setup consists of a cubic three-dimensional cell of volume V = L 3 . The SPH particles are initially positioned on a simple cubic lattice with lattice constant s filling the volume. Then, an equilibration scheme is used to obtain vanishing kernel derivatives before running the actual simulations [10]. Lees-Edwards boundary conditions are used to impose a simple shear flow with rateγ on the system [11]. The average stress tensor is evaluated using the Irving-Kirkwood relation [12], The first term in Eq. (7) is the kinetic contribution with being the deviation of the velocity of particle i from the homogeneous shear flow velocity field. The second term in Eq. (7) is the potential contribution caused by the Navier-Stokes pressure term and the third term in Eq. (7) is the dissipative contribution caused by the Navier-Stokes viscosity term. The shear stress contributions in the simulation are denoted as with kind ∈ {kin, pot, diss}. The total pressure in the system is given by The density of the system, ρ, is evaluated as the ensemble average of the SPH particle densities ρ i . An overview of the ranges of simulation parameters is given in Table 1. The ratio of the simulation volume edge length and the initial particle spacing is L/s = 10. The ratio of the smoothing length and the particle spacing is h/s = 1. The reference density is ρ 0 = 10 3 kg/m 3 and the isentropic exponent is γ = 7. The viscosity model parameter ξ is set to 23.9 according to numerical calibration results. It should be noted that ξ was found to be a function of h/s but not of η 0 . A total number of 45 simulations were carried out using only parameter combinations which yield a system Reynolds number ρ 0γ L 2 /η 0 ≤ 10 in order to limit the study to the non-turbulent regime. All simulations are run until a shear strain, i.e., the product of shear rate and simulation time, of 50 is reached. For the following analysis all dynamic quantities are time averaged over the stationary second half of each simulation. The software SimPARTIX is used for the SPH simulations.

Results
Let us define three stress or energy density scales, respectively, based on the simulation input parameters, Equation (10) is the kinetic energy density of the relative shear movement of the SPH particles, Eq. (11) is the viscous shear stress and Eq. (12) is a measure of the inverse compressibility of the model fluid. Based on these stress scales and the total pressure in the system we are able to define a set of dimensionless numbers, The dimensionless numbers are inspired by the particlebased nature of the SPH simulation method. Equation (  defines a particle Reynolds number in the shear flow while Eq. (14) defines a particle Mach number by relating the relative particle velocity and the speed of sound in the system. Finally, Eq. (15) defines a particle inertial number which resembles the inertial number of dense granular flows [13,14]. Note that Re and Ma are only based on input parameters while I is also based on the simulation result via P tot . Figure 1 displays the linear correlation of the density deviation and the normalized total pressure in the system. This result is not surprising as a Taylor expansion of Eq. (5) around ρ 0 yields P i ≈ γB (ρ i /ρ 0 − 1). Thus, the following relation apparently holds, with P as the time and ensemble averaged pressure form the equation of state.
Furthermore, Fig. 1 shows that the deviation of the dissipative viscosity contribution from the input viscosity is below 0.5 % for a normalized total pressure below ∼ 10 −3 while the dissipative viscosity deviation increases with increasing normalized pressure above ∼ 10 −3 . Ideally, ρ/ρ 0 = S diss /Σ shear = 1.
A power-law relation between the relative kinetic shear stress contribution and the ratio of particle Reynolds and particle inertial number was found as presented in Fig. 2. This result is a first indicator for the connection of SPH simulations and dense granular flow. GDR Midi [13] showed that the mean velocity fluctuations in granular plane shear flow scale as δv 2 ∝γ 2 d 2 /I g , with d as the grain diameter and I g as the granular inertial number, I g =γd √ (ρ 0 /P). Recalling Eqs. (7), (8), (13) and substituting I g → I and d → s, we can write   A power-law scaling of the relative kinetic viscosity with Re/I is discernible in Fig. 2, however, with an exponent of ∼ 2 and not 1 as predicted by Eq. (17). Figure 3 displays a scaling relation for the ratio of the potential shear stress and the total pressure as a function of the normalized total pressure. Similar regimes as found for the scaling bahavior of the dissipative shear stress can be identified (compare Fig. 1). Below a normalized total pressure of about 10 −3 the potential shear stress ratio is constant. Above this value the shear stress ratio increases with the normalized total pressure. Jop et al. [14] used the following expression to model the dependency of the ratio of shear stress and pressure on the granular inertial number I g for dense granular flow, S /P = μ s + (μ 2 − μ s )/(I 0 /I g + 1), with μ s and μ 2 as the limiting values at low and high I g , respectively, and I 0 as a model parameter. A substitution using S → S pot , P → P tot , I g → (Ma/I) 2 and setting μ 2 = 1 yields an expression which describes the data in  Figure 4. Scaling of total pressure. See text for details of the analytic expression fitting the data. Fig. 3 well. The substitution of I g means that the ratio of pressure and stiffness is relevant for the transition in the present system while the ratio of kinetic energy density and pressure is relevent for the granular system. So far, all of the discovered scaling relations include I or P tot , respectively, as a factor on the abscissa. Thus, it would be desirable to obtain a dependency of I only on input parameters of the simulation. Such a relation exists as shown in Fig. 4. The analytic expression which accurately describes the data collapse was used in a very similar form for modelling the rheology of soft, frictional, noncohesive spheres [15,16]. Omitting factors related to the solid volume fraction which has no correspondence in the present system, Gu et al. [16] related the normalized total pressure P to the shear rateγ via with the grain diameter d, the speed of sound c g ,γ =γd/c g and m int ≈ 2/3. By interpreting the SPH particles again as grains of diameter s and substituting P → P tot , c g → c s andγ → Ma/Re 1/2 , an analytic expression is found which describes the simulation data in Fig. 4 accurately. Equation (19) is valid for sheared granular systems below a critical solid volume fraction and describes a transition from an inertial flow regime with Bagnold scaling, i.e., P ∝γ 2 , to an intermediate flow regime where the finite grain stiffness plays a role. Such a transition in the stress scaling behavior with respect to the shear rate was previously observed by Campbell [17] in simulations of soft granular media. The substitution ofγ indicates that the role of inertial shear motion in the granular system is replaced by viscous shear motion in the SPH simulations with a scaling of P tot ∝γ. Interestingly, the transitions in the dissipative shear stress scaling (Fig. 1) and the potential shear stress scaling (Fig. 3) at (Ma/I) 2 ≈ 10 −3 correspond to the transition between flow regimes in Fig. 4 at Ma/Re 1/2 ≈ 3 · 10 −3 .

Conclusion
It is demonstrated that deviations of pressure and shear stress in SPH shear flow simulations from their theoretical values can be expressed as functions of dimensionless system parameters, namely a Reynolds number (Re), a Mach number (Ma) and an inertial number (I).
Interestingly, the SPH simulation results contain signatures observed in dense granular flow. The SPH kinetic shear stress scales as power-law of Re/I which corresponds to experimental measurements of the scaling of granular velocity fluctuations albeit with a different powerlaw exponent [13]. The ratio of SPH potential shear stress and total pressure shows a transition from a constant ratio to a dependency on the total pressure which is also a fundamental feature in the rheology of dense granular flows [14]. The functional dependency of the density or, equivalently, the total pressure on SPH simulation parameters resembles a similar relation found in discrete element simulations of granular flow [16].
The identified scaling relations render quantitative predictions of the accuracy of the stress in SPH simulations possible. Moreover, the influence of the spatial SPH discretization, i.e., the particle size, on the accuracy can be quantified.