GPU-accelerated smoothed particle hydrodynamics modeling of granular flow

Journal Pre-proof GPU-accelerated smoothed particle hydrodynamics modeling of granular flow Jian-Yu Chen, Fue-Sang Lien, Chong Peng, Eugene Yee PII:

S0032-5910(19)30851-4

DOI:

https://doi.org/10.1016/j.powtec.2019.10.017

Reference:

PTEC 14765

To appear in:

Powder Technology

Received Date: 5 June 2019 Revised Date:

1 October 2019

Accepted Date: 4 October 2019

Please cite this article as: J.-Y. Chen, F.-S. Lien, C. Peng, E. Yee, GPU-accelerated smoothed particle hydrodynamics modeling of granular flow, Powder Technology (2019), doi: https://doi.org/10.1016/ j.powtec.2019.10.017. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Head on collision of granular jets

Granular impact with wave structure

GPU-accelerated smoothed particle hydrodynamics modeling of granular flow Jian-Yu Chena , Fue-Sang Lienb , Chong Pengc,∗, Eugene Yeeb a School

of Mechanical Engineering, Nanjing University of Science and Technology, 200 Xiaolingwei Street, Xuanwu District, Nanjing, 210094, China b Department of Mechanical and Mechatronics Engineering, University of Waterloo, 200 University Avenue West, Waterloo, Ontario, N2L 3G1, Canada c Institut f¨ ur Geotechnik, Universit¨at f¨ur Bodenkultur, Feistmantelstrasse 4, 1180 Vienna, Austria

Abstract In this paper, the dynamical behavior of granular media is investigated numerically using smoothed particle hydrodynamics (SPH) with an elastoplastic and a phase-change constitutive model. A GPU-accelerated SPH code for granular flow modeling is implemented and validated using four test cases: (1) the collapse of an axisymmetric sand pile; (2) the impact force of sand on a rigid wall; (3) the head-on collisions of dense granular jets; and, (4) the impact of a granular jet with a wave structure on the granular film. The results of SPH simulations with the two models are compared to available experimental and numerical results. The validation showed that the dynamics of granular media are well predicted using the SPH methodology with an appropriate material model. An SPH simulation of granular flow run on a rather outdated GPU was already about 160 times as efficient as the same simulation run on a mainstream single CPU. Keywords: elastoplastic, CUDA, GPU, granular flow, phase change, SPH

1

1. Introduction

2

Various geotechnical problems associated with large deformations, such as landslides and granular impacts, are

3

difficult to address using conventional grid-based methods. More specifically, the large deformations that characterize

4

this class of problems frequently result in severe mesh distortions that provoke instabilities in the numerical algorithm,

5

which in turn cause the simulation to terminate prematurely. Owing to these difficulties, we will investigate the ap-

6

plication of an alternative mesh-free methodology for the simulation of granular impact. To this purpose, smoothed

7

particle hydrodynamics (SPH) implemented using the compute unified device architecture (CUDA) programming

8

interface for parallel computing on a graphics processing unit (GPU) is used to address the problem of large defor-

9

mations that generally characterize three-dimensional (3D) granular impact problems. The rationale for using SPH

10

that harnesses the power of a GPU for the simulation of granular impact are as follows: (1) the Lagrangian and

11

meshless property of the SPH methodology facilitates the simulation of granular flow and impact problems involving ∗ Corresponding

author Email address: [email protected] (Chong Peng )

Preprint submitted to Powder Technology

October 5, 2019

12

large deformations as compared to the more conventional grid-based methods; and, (2) the implementation of SPH

13

using CUDA allows the parallelizable parts of the computation to off-loaded to a GPU to give a high floating-point

14

computing performance.

15

The SPH methodology was originally proposed by Gingold and Monaghan [1] and Lucy [2] and applied to address

16

various astrophysics problems. Subsequently, this methodology was applied to problems in other areas of scientific

17

endeavor (e.g., explosive detonations, geomechanics problems). Bui et al. [3, 4] report some seminal work concerning

18

the implementation of an elastoplastic constitutive model for the simulation of the elastic and plastic behavior of soil

19

within the SPH framework. Peng et al. [5, 6] introduced the hypoplastic constitutive model into the SPH methodology

20

and applied it to the simulation of granular flow. This effort demonstrated that the hypoplastic constitutive model

21

can provide good predictions in granular flow problems involving large deformations. Fan et al. [7] and Fan and

22

Li [8] simulated landmine detonations using three different soil constitutive models within the SPH methodology in

23

conjunction with the peridynamics method. Chen and Lien [9] and Chen et al. [10] applied the SPH methodology with

24

an elastoplastic and a hypoplastic constitutive model to investigate phenomoneology associated with the detonation

25

of landmines and its effect on structures.

26

Recently, a number of experiments and simulations related to granular impacts have been conducted using a mesh-

27

free methodology. This work concerning granular material impacts is of significant importance for various industrial

28

applications including ink-jet printing and inpinging jet milling. Pacheco-V´azquez and Ruiz-Su´arez [11] carried out

29

impact cratering experiments using granular projectiles with different porosities. Ellowitz [12] investigated the phe-

30

nomenon of two-dimensional (2D) head-on collisions of dense granular jets using discrete particle simulations . Shi

31

et al. [13] studied the dynamic behavior of a dense granular jet impacting on a circular target using the Discrete Ele-

32

ment Method (DEM). To date, most of the modeling effort concerning granular impact modeling has been conducted

33

using the DEM method as described in Cundall and Strack [14]. While accurate, the DEM method involves solving

34

Newton’s equation of motion for each grain individually. In consequence, this method is computationally prohibitive

35

for problems involving a large physical domain [15]. In comparison to DEM, the SPH methodology is advantageous

36

(and hence, preferred) for the simulation of granular flows for two reasons. Firstly, it is straightforward to implement

37

different constitutive models for the representation of the mechanical behavior of granular materials within the SPH

38

framework. Secondly, the material domain in the SPH methodology is discretized using particles (discrete moving

39

Lagrangian elements) that carry the relevant physical variables (properties) of the flow. Unlike DEM, these particles

40

are not physical grains of the material and, as a result, their size can be chosen in accordance to the problem that needs

41

to be addressed. Indeed, for certain problems, the size of the SPH particles can be much larger than the physical grains

42

of the material, which can result in a significant reduction in the computational cost of SPH simulations per number

43

of particles (in comparison to DEM).

44

Although numerous investigations of the simulation of granular flows have been conducted, the application of the

45

SPH methodology to the modeling and simulation of 3D granular flow is still in its infancy. In this paper, the phase-

46

change constitutive model will be incorporated into the SPH framework and used to investigate various phenomena 2

47

associated with granular flow for the first time. This material constitutive model was recently proposed by Dunatunga

48

and Karmrin [15] who used it to investigate 2D granular flow problems within the material point method (MPM)

49

framework. In comparison with the more commonly used rate-independent elastoplastic and Mohr–Coulomb models,

50

the phase-change constitutive model is more suitable for addressing problems involving dense zones and rapid flows.

51

Furthermore, the effects of particle size have also been incorporated into the phase-change constitutive model [15].

52

The search for the neighbors of a particle and the calculation of the rates of change of various physical variables

53

are extremely time-consuming operations in the SPH methodology. As a consequence, the parallelization of the

54

SPH methodology is necessary, especially for simulations of 3D granular dynamic behavior that necessarily involve

55

large numbers of particles. To this purpose, it is noted that CUDA, which is a parallel computing platform and

56

programming interface for GPUs developed by NVIDIA, is especially well-suited for the efficient implementation

57

of the SPH methodology. The implementation of the SPH methodology on GPUs was first described by Harada

58

et al. [16]. Subsequently, Crespo et al. [17] reported GPU-accelerated SPH simulations for the investigation of

59

complex free surface flows. These researchers showed that their SPH simulations on a single GPU were two orders

60

of magnitude faster than that conducted on a single-core CPU. Recently, GPU-accelerated SPH models have been

61

applied to the simulation of particle-fluid flows and to the solution of the shallow water equations [18, 19], and the

62

computational efficiency of these models has been further improved through the incorporation of a adaptive particle

63

splitting and merging methodology [20]. Furthermore, a GPU-accelerated DEM methodology has been implemented

64

to tackle the challenges associated with modelling powder compaction [21, 22].

65

In this paper, the CUDA programming interface has been used for the implementation of a new SPH code for sim-

66

ulation of granular flows on a GPU. The in-house parallelized SPH code has been implemented in C++ and designed

67

specifically for the computationally efficient and accurate simulation of 3D granular impact problems involving large

68

numbers of particles. To the best of our knowledge, this is the first time that the phase-change constitutive model has

69

been incorporated in the SPH methodology and used to investigate the phenomenology associated with 3D granular

70

impact problems. Furthermore, phenomena associated with 3D head-on collisions and the impact of granular jets with

71

a wave structure on a granular film are modeled and simulated within the SPH framework for the first time.

72

The paper is organized into six sections. Section 2 contains a brief introduction of the SPH methodology. Section

73

3 describes the phase-change and elastoplastic constitutive models. Section 4 presents the computational details of

74

our implementation of the SPH methodology on a GPU. The two constitutive models mentioned above are validated

75

using four test cases in Section 5: namely, the collapse of an axisymmetric sand pile, the impact force of sand on

76

a rigid wall, head-on collisions of dense granular jets, and the impact of a granular jet with a wave structure on a

77

granular film. Finally, the key results and novelty of the research are summarized in Section 6.

3

78

2. Fundamentals of the SPH method

79

2.1. Function approximations in SPH

80

The derivation of the SPH methodology can be divided into two steps. The first step involves a kernel approxi-

81

mation, whereby all the physical variables are approximated using an integral transform involving the convolution of

82

the physical variable with a chosen kernel function. To this purpose, the kernel approximation of a function f (x) at

83

position x = (x, y, z) in a physical domain Ω is h f (x)i =

84

85

Z Ω

f (x0 )W(x − x0 , h)dx0 ,

(1)

where W is the kernel function and h is the smoothing length. The kernel approximation of the gradient of the function ∇ f (x) can also be obtained in the same manner to give Z h∇ f (x)i = − f (x0 )∇W(x − x0 , h)dx0 . Ω

(2)

86

The Wendland C2 function is used as the kernel function W for the SPH simulations conducted herein. This kernel

87

function assumes the form

88

89

     q 4     1 − 2 (2q + 1) 0 ≤ q < 2 ; W(q, h) = αd      0 q≥2,

(3)


 where q is the Euclidean distance r normalized by the kernel width h so q = r/h. The coefficient αd = 7 4πh2 in
 a two-dimensional space and αd = 21 16πh3 in a three-dimensional space. In SPH, the kernel representation for a

90

function and its gradient as exhibited above can be approximated by a Riemann summation over the particles used in

91

the simulation, so h f (xi )i =

N X

f (x j )Wi j

j=1 92

and h∇ f (xi )i = −

N X j=1

mj , ρj

f (x j )∇Wi j

(4)

mj , ρj

(5)

93

where ρ j and m j are the density and mass of particle j, respectively, and N is the number of particles. For notational

94

convenience, W(xi − x j , h) is abbreviated as Wi j .

95

The second step of the SPH methodology involves solving the equations governing granular impact; namely, the

96

continuity and momentum equations. Within the SPH framework, these equations can be discretized with the kernel

97

approximation, respectively, as follows:  N X  ∂Wi j   dρi   m j vαi j · = , (a)  dt   ∂xiα   j=1   N  σαβ + σαβ  ∂W X   i j ij dvαi    = m + Πi j δαβ , (b)  j dt  β  ρi ρ j ∂x j=1

(6)

i

98

 where ρi , vi , pi , mi are the density, velocity, pressure, and mass of particle i; d( · ) dt denotes the time derivative of a

99

physical quantity; and, vαi j = vαi −vαj . The term Πi j in the momentum equation is the Monaghan-type artificial viscosity. 4

100

2.2. Artificial viscosity

101

In order to smooth potential unphysical oscillations, to prevent unphysical particle-particle penetrations, and to

102

stabilize the numerical solutions, a Monaghan-type artificial viscosity is incorporated into the momentum equation as

103

follows [23]:   −αci j φi j +βφ2i j      ρi j Πi j =      0

(7)

vij · xij ≥ 0 ,

+ c j ), ρi j =

+ ρ j ), hi j =

+ h j ), vij = vi − vj , and

104

where φi j = hi j vij · xij

105

xij = xi − xj . Here, ci is the speed of sound associated with the particle i; α and β are constant coefficients that have

106

values of about 1.0; and, ϕ = 0.1hi j is applied in order to prevent any unphysical penetration of two particles.

107

2.3. Neighboring particle search





|xij |2 + ϕ2 , ci j =

vij · xij < 0 ;

1 2 (ci

1 2 (ρi

1 2 (hi

108

Owing to the fact that calculation of any physical variable at a given position in the SPH methodology involves the

109

summation over nearest neighbor particles to that position, a critical issue related to the computational performance

110

of SPH involves the formulation of an efficient algorithm for searching for the neighboring particles at a specified lo-

111

cation in a problem domain. A number of different approaches has been proposed for determination of a dynamically

112

evolving list of neighbouring particles required for the computation of sums in SPH. In this paper, the cell-linked list

113

(CLL) is used for the search of these neighborhood sets or lists (see Fig. 1) [24]. In the CLL approach for a neighbor-

114

hood search, the computational domain is divided into cells with a side length of 2h. A particle in the simulation is

115

referenced with respect to the cell that it belongs to. In consequence, this action results in a list of particles that have

116

been reordered with respect to the uniform grid of cells used in the discretization of the computational domain. For the

117

calculation of the rates of change of a physical variable in SPH, each particle i in a particular cell of the computational

118

domain only needs to query the cells adjacent to this particular cell in order to find all the nearest neighbor particles

119

j of particle i. More specifically, if the distance between particle i and particle j is less than 2h (cut-off limit), then

120

particle j is a genuine neighbor of the given particle i.

121

2.4. Boundary conditions and impact force calculation

122

Because SPH is a Lagrangian method, the imposition of a specific boundary condition is problematic. A number

123

of approaches has been proposed for the specification of boundary conditions in SPH [25, 26, 27]. In this paper,

124

a methodology involving the introduction of dummy particles is used for the imposition of the no-slip boundary

125

condition along a solid (wall) surface [28] (see Fig. 2). The velocity of the solid wall boundary vw is assigned to the

126

velocity of the dummy particles vd (viz., vd = vw ). For the calculation of the rate of change of a physical quantity

127

associated with a real particle κ, the components of the stress tensor of the dummy particles are extrapolated using the

5

128

values of the neighboring granular particles, so N X

σidj =

σiκj W(rdκ ) − (gi − aid )

κ=1

N X

i ρκ rdκ W(rdκ )δi j

κ=1 N X

,

(8)

W(rdκ )

κ=1 129

where subscripts d and κ are used to denote the dummy and granular particles, respectively. Moreover, ad is the

130

prescribed acceleration of the dummy particles; g is the gravitational acceleration; rdκ is the displacement vector

131

from the granular particle κ to the dummy particle d; rdκ is the Euclidean distance between the granular particle κ

132

and the dummy particle d, so rdκ = ||rdκ ||; δi j is the Kronecker delta function (equivalent to the identity matrix in

133

three dimensions); and, N denotes the number of granular particles within the support domain (taken to be a spherical

134

region with a radius that is some multiple k of the kernel smoothing length h) of the dummy particle d.

135

When calculating the impact force of a granular material on the boundary, we first use the extrapolated stresses

136

of the dummy particles given above to determine the interaction force on each dummy particle through application of

137

the momentum equation as follows [29]: Fdα = md

X

mi

 σαβ + σαβ i

d

ρd ρi

i

+ Πdi δαβ

 ∂W

di

∂xdβ

.

(9)

138

Following this calculation, the impact force FW exerted on the solid wall boundary can be calculated as the summation

139

of the forces Fd on the dummy particles, so FW =

XX d

md mi nα

 σαβ + σαβ i

d

ρd ρi

i

+ Πdi δαβ

 ∂W

di

∂xdβ

,

(10)

140

where the subscripts i and d denote the granular and dummy particles, respectively; and, n is the normal vector to the

141

solid wall boundary.

142

2.5. Physical variables update

143

The evolution of the various physical variables is obtained by integrating equation (6) with respect to the time t.

144

To this purpose, a second-order predictor-corrector scheme is used to update the physical variables. For the prediction

145

step, the position xi , velocity vi , and density ρi associated with particle i are updated as follows: t+ ∆t2

xi

= xti +

∆t t v , 2 i

(11)

= vit +

∆t t a , 2 i

(12)

146

t+ ∆t2

vi 147

and t+ ∆t2

ρi

= ρti + 6

∆t dρi 2 dt

!t .

(13)

The stress tensor associated with particle i can be updated in a similar manner, so !t ∆t dσi t+ ∆t2 t = σi + . σi 2 dt

(14)

For the correction step, the various physical variables are updated in accordance to

148

t+ ∆t2

= xti + ∆tvi xt+∆t i

,

(15)

,

(16)

149

t+ ∆t2

vit+∆t = vit + ∆tai 150

and

!t+ ∆t dρi 2 , (17) = + ∆t dt where ai is the acceleration of particle i, which is calculated using equation (6(b)). By contrast, in the correction step, ρit+∆t

ρti

the stress tensor is updated using the value of the stress tensor obtained at end to the prediction step (half-time step) in order to simplify the computation and minimize the storage required in the GPU memory. As a consequence, !t+ ∆t ∆t dσi 2 t+ ∆t2 t+∆t . (18) σi = σi + 2 dt 151

Finally, the Courant-Friedrichs-Lewy (CFL) condition is imposed on this explicit time integration scheme in order

152

to guarantee the stability of this second-order predictor-corrector scheme. To this purpose, the time step for the

153

integration is chosen as ! 1/2 h ∆t = 0.2 min , h/c , (19) max ||ai || where c is the speed of sound for the granular material. The speed of sound is determined in accordance to c = 

154

155

(E/ρ)1/2 , where E is Young’s modulus of the granular material.

156

3. Constitutive equations for the modeling of granular impact

157

158

A constitutive equation is required in order to calculate the stress tensor σ and to close the governing equations. Towards this objective, the stress tensor is decomposed into two parts as follows: σαβ = −pδαβ + S αβ ,

(20)

where −p and S are the pressure and the deviatoric stress tensor, respectively. Normally, the pressure −p is calculated using an equation of state (EOS), while the deviatoric stress S is obtained through a constitutive equation. However, in the simulations conducted herein, the stress tensor σ is calculated directly using a constitutive equation and, subsequently, the pressure is obtained from the mean stress without having to use an EOS. More specifically, p=−


σαα 1 = − σ11 + σ22 + σ33 , 3 3

(21)

159

where σ11 , σ22 , and σ33 are the components of the normal stress tensor in x, y and z directions, respectively. Both the

160

phase-change and elastoplastic constitutive models are adapted for use in our SPH simulations in order to update the

161

stress tensor σ. 7

162

3.1. Phase-change constitutive model

163

The phase-change constitutive model proposed by Dunatunga and Kamrin [15] is implemented within the SPH

164

framework for the first time. This constitutive model allows the granular media to change through several com-

165

mon phases (viz., the grains can behave like a solid, liquid or gas). The granular matter is considered to be elasto-

166

viscoplastic when the media is sufficiently dense and is treated as a disconnected material when the density of the

167

material is smaller than a critical value. In order to calculate an objective rate of change for the stress tensor σ, ˚ the

168

Jaumann frame of reference is used: σ˚ = σ˙ + σ ω ˙ − ωσ ˙ ,

169

170

(22)

where σ˙ is the stress rate tensor and ω ˙ is the spin tensor. Now, the components of the trial stress tensor σαβ tr can be calculated as ! ! 1 αβ γγ αγ βγ γβ αγ αβ γγ αβ σαβ = ∆t σ ω ˙ + σ ω ˙ + 2G ε ˙ − δ ε ˙ + K ε ˙ δ + σαβ , tr 3

(23)

where ε˙ is the strain rate tensor. The shear modulus G and the elastic bulk modulus K are determined as follows: K=

E , 3(1 − 2ν)

(24)

G=

E , 2(1 + ν)

(25)

and

171

where E is Young’s modulus and ν is Poisson’s ratio.

172

Using the trial stress tensor determined above, the procedure for updating the stress consists of the following steps.

173

1. If the particle density ρn+1 is less than the critical value ρc , then set the stress tensor of the particle to σ n+1 = 0

174

175

for the next time step.
2. For all other cases, if the trial pressure ptr is negative, then σ n+1 = 0 and ptr = − 31 tr σtr . 1/2  3. Otherwise, the stress tensor is computed as follows. If τtr ≤ S 0 where τtr = 21 (Str : Str ) and S 0 = µ s ptr , the granular material is in the elastic regime, so set σ n+1 = σtr . If τtr > S 0 , the granular material is in the plastic flow regime, and τn+1 is determined from 

τn+1

2

− Bτn+1 + H = 0 .

(26)

This quadratic equation for τn+1 can be solved explicitly to give τn+1 =

2H B+

B2

− 4H


1/2 ,

(27)

where B and H are obtained from B = S 2 + τtr + α and H = S 2 τtr + S 0 α. Furthermore, S 2 and α are determined from S 2 = µ2 ptr and α = ξG∆tptr , where ξ is a parameter of the granular media, µ s is the friction coefficient at zero shear rate, and µ2 is the maximum value of the friction coefficient. The stress tensor at the next time step is given by σ n+1 =

τn+1 Str − ptr I . τtr 8

(28)

176

3.2. Elastoplastic constitutive model

177

The classical elastoplastic constitutive model has also been incorporated into SPH framework to model the gran-

178

ular impact. The results obtained from this constitutive model will be compared with those obtained from the phase-

179

change constitutive model. In the elastoplastic constitutive model, the plastic deformation occurs only if the Drucker-

180

Prager criteria is satisfied. This criteria is given by f (I1 , J2 ) = J21/2 − ϕφ I1 − kc = 0 ,

(29)

181

where I1 = σαα and J2 = 21 S αβ S αβ are the first invariant of the stress tensor and the second invariant of the deviatoric

182

stress tensor, respectively. The constants ϕφ and kc in the Drucker-Prager model are calculated from the cohesion c

183

and the internal friction φ as ϕφ =

184

tan φ 9 + 12 tan2 φ

(30)

and 3c

kc =

(31)
. 9 + 12 tan2 φ 1/2 The Jaumann frame of reference is used also with respect to the elastoplastic constitutive model. The final stress-

185

186


1/2 ,

strain relationship of the granular flow model evolves in accordance to      dσαβ 1 αβ γγ  G  αβ αγ βγ γβ αγ αβ γγ αβ i αβ = σi ω ˙ i + σi ω ˙ i + 2G ε˙ i − δi ε˙ i + K ε˙ i δi − λ˙ i 9K sin ψδ + 1/2 S i  , dt 3 J2

(32)

187

where ψ is the dilatancy angle. The parameters K and G can also be obtained from Eqs (24) and (25), respectively.

188

The plastic multiplier λ˙ i is determined from λ˙ i =

189

1/2
αβ αβ 3ϕφ K ε˙ γγ S i ε˙ i i + G/J2 . 27ϕφ K sin ψ + G

The strain rate tensor and spin tensor are approximated within the SPH framework as follows:   N  1 X m j  α ∂Wi j αβ β ∂Wi j  ε˙ i = v ji β + v ji α  , 2 j=1 ρ j ∂xi ∂x

(33)

(34)

i

190

and ω ˙ αβ i

N

1 X mj = 2 j=1 ρ j

   α ∂Wi j  β ∂Wi j  v ji − v ji α  . β ∂xi ∂xi

(35)

191

The reader is referred to Bui et al. [3] for a more detailed description of the stress-strain relationship in the context of

192

the elastoplastic constitutive model.

193

4. GPU implementation

194

4.1. GPU solution procedure

195

The application of SPH for the simulation of granular flow in three dimensions is computationally intensive owing

196

to the large number of particles required for such simulations. As a consequence, to enable these types of simulations 9

197

to complete in a reasonable execution time, it is imperative to accelerate the computational workload by parallelization

198

of the SPH simulations on a graphics processing unit (GPU). The parallelization of the SPH methodology on a GPU

199

is depicted in Fig. 3. There are a number of steps involved in the implementation of the SPH methodology on a GPU.

200

Firstly, the data describing the initial configuration of the granular impact is loaded into the central processing unit

201

(CPU). Next, this information is transferred from the CPU to the GPU. The GPU is used to facilitate the following

202

operations: building the neighborhood lists and searching these neighborhood sets (which necessarily involve dynam-

203

ically changing lists of neighboring particles), determination of the particle interactions, and update of the physical

204

variables which involve the computation of sums over the dynamically evolving neighborhood sets. At the end of

205

this computational process, the information computed by the GPU is transferred from the GPU to the CPU and the

206

relevant numerical data is saved. All the SPH simulations of granular flow reported herein were conducted on a PC

207

with a single NVIDIA Quadro K620 graphics processor. The performance of the in-house SPH code used for our

208

simulations, such as the accuracy and computational efficiency, is described in detail in [29, 30].

209

4.2. Pseudocode for determination of particle interactions on GPU

210

The computation of the time development of the physical variables is a compute-intensive process that must be

211

parallelized in order to obtain a good computational performance for an SPH simulation. In consequence, the use of

212

the CUDA programming interface in our in-house SPH code is a key aspect in the acceleration of the computational

213

workload on a GPU. A pseudocode for the algorithm used for the computation of the particle interactions on a GPU

214

is summarized in Algorithm 1 below. The implementation utilizes the CUDA platform (viz., the parallel computing

215

platform and application programming interface developed by NVIDIA for their GPUs).

216

5. Numerical cases

217

5.1. Axisymmetric collapse of dry sand

218

We first test the phase-change and elastoplastic constitutive models for SPH simulations of the 3D axisymmetric

219

collapse of dry sand. Lube et al. [31] conducted experiments of the axisymmetric collapse of columns of dry sand

220

that was initially contained within a cylindrical column resting on a flat surface (see Fig. 4). The nature and mode of

221

collapse of these granular columns is determined by the initial column height-to-halfwidth ratio a = h0 /d0 .

222

Four different cases involving different ratios of a were simulated using our SPH code and the results of these

223

simulations were compared to the experimental results in order to validate the two constitutive models. Table 1 sum-

224

marizes the basic information for each simulation. The parameters for the elastoplastic and phase-change constitutive

225

models are listed in Table 2. These parameters were used for all subsequent SPH simulations reported herein, unless

226

noted otherwise.

227

The magnitude of the velocity of the granular collapse at two different ratios a = 0.5 and 4.0, obtained using the

228

phase-change constitutive model, is shown in Fig. 5. It can be seen that in the final phase of the granular collapse, 10

229

the top surface of the shorter and stouter sand column is flatter than that of the taller and more slender sand column.

230

Qualitatively, it is noted that the overall general shape of the sand pile in the final phase (steady state) of collapse

231

predicted by the simulation is similar to that obtained from the experimental measurements conducted by Lube et

232

al. [31].

233

The stress distributions of the sand collapse using the elastoplastic constitutive model are exhibited in Fig. 6.

234

The stress in vertical direction yy is proportional to the depth of the sand. This is consistent with the vertical stress

235

distribution of the sand confined in a cylindrical column and subjected to a gravitational force. Furthermore, the

236

numerical results for the runout distance d∞ and the final height h∞ obtained using the phase-change and elastoplastic

237

constitutive models are summarized in Figs 7 and 8, respectively. Superimposed on these two figures are some

238

experimental measurements of d∞ and h∞ obtained by Lube et al. [31]. A perusal of Figs 7 and 8 shows that the

239

numerical predictions for d∞ and h∞ are generally in good conformance with the experimental results. Finally, it is

240

noted that the quality of the current numerical simulations of 3D granular collapse obtained using the SPH framework

241

with both the phase-change and elastoplastic constitutive models are comparable to the 2D sand collapse simulations

242

obtained by Dunatunga and Kamrin [15] using the MPM framework with the phase-change constitutive model.

243

We close this subsection with an evaluation of the computational efficiency of our GPU-accelerated SPH code.

244

To this purpose, we will measure the compututational efficiency of the implementation of our GPU-accelerated SPH

245

code in terms of frames per second (FPS), which is the number of computational time steps executed in one second

246

of wall-clock time (see Table 3 which summarizes the FPS and total wall-clock execution time T for four test cases

247

involving the simulation of the collapse of a sand pile with an initial aspect ratio of a = 0.5). The total physical time

248

of the simulation for all the test cases summarized in Table 3 is fixed at a value 0.1 s. A careful examination of this

249

table shows that the FPS decreases as the number of particles used in the simulation increases. More specifically, it

250

is seen that for the case that uses nearly 500,000 particles in the simulation, the FPS achieved (≈ 6.8 time steps per

251

second) is still very good. As mentioned previously, the simulations reported here were conducted on a PC with a

252

single NVIDIA Quadro K620 graphics processor. Note that this graphic card was released in 2014, five years ago. It

253

is an outdated card and only has 384 cuda cores. It is expected that the computational performance reported here will

254

be significantly improved on a more powerful NVIDIA graphics processor, such as the NVIDIA GeForce RTX 2080

255

Ti.

256

In order to characterize the acceleration in computational performance provided by the use of a GPU, the SPH

257

simulations summarized in Table 3 were also conducted on a PC with a single CPU (AMD Opteron Processor 6320).

258

The speed-up in computational performance between the GPU and CPU is defined as Sp =

FPS on a single GPU . FPS on a single CPU

(36)

259

Note that the time required for the input of the data files and the output of the simulation results is excluded from the

260

execution time. The speed-up in the computational performance of the SPH code executed on a single GPU versus

261

that on a single CPU is displayed in Fig. 9 for five test cases (four of which have been summarized in Table 3). It 11

262

is seen that a speed-up of up to about a factor of 160 is achieved as the number of particles N used in the simulation

263

increases to about 750,000.

264

We can compare our GPU-accelerated SPH implementation with an open-source SPH code. To this purpose, it

265

is noted that a granular flow cannot be simulated directly using either the open-source code SPHysics [32] or the

266

code developed by Liu and Liu [33] owing to the fact that there is no constitutive model for soil in these codes.

267

In consequence, we have implemented a soil constitutive model in the serial open-source SPH code of Liu and Liu

268

[33] and used it as the basis for comparison with the computational efficiency of our GPU-accelerated SPH code

269

implementation. Towards this objective, we found that the computational efficiency of our GPH-accelerated SPH

270

code is 170 times faster that the serial open-source code of Liu and Liu [33] when 170,000 particles were used

271

in the simulation. This speedup suggests that our GPU-accelerated SPH implementation has achieved a high parallel

272

computational efficiency, in spite of the fact that the simulations reported herein were conducted with a rather outdated

273

graphics card.

274

5.2. Impact force of sand on a rigid wall

275

Having completed the validation of the phase-change and elastoplastic constitutive models in the SPH framework

276

using a case study involving the axisymmetric collapse of dry sand, this numerical framework will now be applied

277

to investigate the phenomenon of granular impacts. Moriguchi et al. [34] conducted experimental measurements of

278

the impact force exerted on a rigid wall resulting from the impingement of dry Toyoura fine sand released from an

279

inclined flume. The sand was initially confined inside a box at the top of the flume. This material was released through

280

a side door, and the sand impacted a rigid wall at the bottom of the flume. The impact force exerted by the sand on

281

the rigid wall was measured as a function of time using a sensor attached to the wall barrier.

282

In this subsection, we will conduct a 3D numerical simulation of this experiment using SPH with the phase-change

283

constitutive model. The initial configuration for the simulation is shown in Fig. 10. The length and height of the sand

284

box are 50 cm and 30 cm, respectively. The total mass of the sand contained initially in the box is 50 kg. The density

285

of the granular material (sand) is ρ0 = 1340 kg m−3 and the critical density is ρc = 500 kg m−3 . The initial particle

286

spacing used in the simulation is 10 mm. A total of 123,012 particles were used in the SPH simulation. Three layers

287

of dummy particles are used to impose the no-slip boundary condition along the solid wall of the incline. The total

288

physical time simulated using the SPH methdology was 3 s. The computational (wall-clock) time required to complete

289

this simulation was 72 minutes.

290

The results for the simulation of granular impact down a 45◦ inclined flume obtained using SPH with the phase-

291

change constitutive model is presented in Fig. 11. The sand moves downward gradually owing to the gravitational

292

force. The magnitude of the velocity (speed) of the sand increases from zero to almost 4.5 m s−1 during this movement.

293

After about 0.70 s after the initial release, the sand impacts the rigid wall and stops almost immediately after the

294

impact. In consequence, the sand at the bottom of the flume accumulates at the rigid wall and this pile of sand grows

295

in size with increasing time. At the end of the process, the pile of sand reaches the top of the wall barrier and flows 12

296

out of the flume. This dynamic behavior of the sand obtained from the numerical simulations is in general agreement

297

with both the experimental measurements reported by Moriguchi et al. [34] and the simulations conducted by Neto

298

and Borja [35]. The dynamic behavior of the sand obtained in the SPH framework using the elastoplastic constitutive

299

model is similar to that obtained using the phase-change constitutive model. However, it is noted that using the

300

elastoplastic constitutive model, the sand is predicted to impact the rigid wall at 0.65 s after the initial release, which

301

occurs earlier than the prediction provided by the phase-change constitutive model.

302

303

304

The impact force exerted by the sand on the rigid wall is calculated using equation (10). This equation required the specification of the normal vector n = (nα ). For the current case, the vector normal to the rigid wall at the bottom √ √
of the flume for an inclination angle of 45◦ is n = 22 , 0, 22 . Fig. 12 exhibits the predicted time variation of the

305

impact force obtained from our SPH simulations using both the phase-change and elastoplastic constitutive models.

306

The predicted peak impact force obtained using the phase-change constitutive model is 192 N. This predicted value

307

for the peak impact force is larger than the measured value of the peak impact force obtained by Moriguchi et al. [34].

308

More specifically, the overprediction of the peak impact force is about 3.2%. It is noted that the prediction of the

309

peak impact force obtained from our SPH simulations with the phase-change constitutive model is significantly better

310

that obtained from either our SPH simulations with the elastoplastic constitutive model or the numerical simulations

311

conducted by Neto and Borja [35]. In the latter two cases, the predicted peak value of the impact force was about

312

212 N which significantly overpredicts the measured value of the peak impact force.

313

After the rapid initial increase to its peak value, the impact force decreases gradually to value of about 160 N.

314

This decrease occurs because some of the sand flows out of the bottom boundary and the velocity of the sand is

315

reduced after its initial impact with the rigid wall. In general, it is noted that SPH simulations with the phase-change

316

constitutive model is capable of capturing the dynamical behavior of the impact of the sand against the rigid wall.

317

Furthermore, these simulations give generally good predictions for the time variation of the impact force. The slight

318

overestimation of peak impact force is probably due to a force contribution from the soil viscosity and the boundary

319

implementation used here. The effects of the soil viscosity, which has the tendency to slow down the granular flow

320

and reduce the impact force, has not been incorporated in our numerical simulations. In addition, it is expected once

321

a more appropriate model for the calculation of the stress tensor of the boundary particles has been developed, its

322

inclusion in our numerical simulation will lead probably to more accurate predictions of the impact force.

323

5.3. Head-on collisions of dense granular jets

324

In this subsection, we investigate the 3D head-on collisions of dense granular jets. A study of the physical mech-

325

anisms involved in the collisions of granular jets will enhance our understanding of the loads applied to structures by

326

the impact of granular matter on them at high velocity. Ellowitz [12] conducted 2D simulations of head-on collisions

327

of dense granular jets having equal speeds using DEM. Here, we extend this effort to consider fully 3D simulations

328

of dense granular jet collisions using the SPH framework. To this purpose, two granular jets having different sizes

329

will be considered in our simulations: namely, the diameter of the first jet is taken to be 2R and that of the second 13

330

jet is taken to be 2kR, where 0 ≤ k ≤ 1. The centers of these two jets are assumed to be aligned, and their initial

331

velocities are both v0 in magnitude. These two jet velocities are in opposite directions, so that the two jets undergo

332

a head-on collision. The ejecta angle ψ0 , which is the maximum angle between the tangent line of the surface of the

333

jet to the horizontal direction, and the velocity u0 of the impact center point at the interface of the two jets are used to

334

characterize the impact process (see Fig. 13). The initial configurations used for three different simulation scenarios

335

for head-on collisions of the two granular jets are summarized in Table 4.

336

Fig. 14 shows the predicted results obtained using SPH with the phase-change constitutive model for the head-on

337

collision of two granular jets with the initial configuration specified by scenario 1 (cf. Table 4). It is noted that the

338

ejecta angle of the granular jet with the larger diameter increases during the impact process, and its configuration after

339

the impact is similar to the 2D geometry described in Ellowitz [12]. The impact center drifts at a constant speed of

340

about u0 = 0.14 m s−1 , which is consistent with the simulations described in [12]. Fig. 15 exhibits the stress and strain

341

distributions obtained using SPH with the phase-change constitutive model at a time of 0.02 s. A persual of this figure

342

shows that the stress in the horizontal direction assumes the largest values within the impact area. The shear strain

343

and stress in the yz direction are antisymmetric with respect to horizontal x direction.

344

The effect of different velocities of the two jets (as characterized by the coefficient k) on the velocity of the impact

345

center has also been studied. To this purpose, the initial configurations for the colliding granular jets are summarized

346

in scenarios 2 and 3 (cf. Table 4). The constant velocity u0 of the impact center is 0 m s−1 for scenario 2 where the

347

radii of the two jets were equal (viz., k = 1) and 3.96 m s−1 for scenario 3 where the radius of one jet is 1/4 that of

348

the other jet (viz., k = 0.25). It seen that a smaller value for k results in a larger value for the velocity u0 of the impact

349

center, a result that is in conformance with the numerical simulations conducted by Ellowitz [12]. The results of our

350

numerical simulations have shown that the application of the SPH methodology with the phase-change constitutive

351

model provides a good alternative scheme to study the impact of two granular jets.

352

In order to demonstrate the ability of the proposed SPH methodology to capture phase changes in the collision

353

of two granular jets, we set the initial density of the disconnected granular jets to have a value of ρ = 1200 kg

354

m−3 . This value for the initial jet density is less than the critical density which has a value of ρc = 1500 kg m−3 .

355

Otherwise, initial conditions for this case are exactly the same as those of scenario 3. The results of the simulation

356

of the impact of the two granular jets for this case obtained using SPH with the phase-change constitutive model

357

are displayed in Fig. 16. An examination of this figure shows that the two disconnected jets undergoes a process of

358

densification after they collide with the each other. This simulation highlights the ability of the SPH methodology

359

to model large inhomogeneous deformations of the granular material. These results suggest that the effects of the

360

granular particle size and fracture on the interaction of high-velocity granular slugs can be investigated using the

361

proposed SPH methodology.

14

362

5.4. The impact of granular jet with a wave structure on a granular film

363

An experimental study of the liquid-like wave structures on a granular film resulting from a granular jet impact

364

has been described by Shi et al. [36]. In this study, the initial configuration of the particles is a sinusoidal wave. In

365

particular, the motion of the center particles at the jet exit is characterized by y = A0 sin (2π f0 t) ,

(37)

366

where A0 and f0 are the amplitude and frequency of the initial granular jet, respectively, and t is the time. In this

367

subsection, we simulate the liquid-like wave structure on a granular film using the SPH methodology with both the

368

phase-change and elastoplastic constitutive models.

369

To this purpose, the wave structure generated by an amplitude of A0 = 0.002 m and a frequency of f0 = 160 Hz

370

for the initial granular jet has been investigated using the SPH methdology. The initial particle spacing is 0.002 m and

371

800,000 particles were used in this simulation. The initial velocity of the granular jet is 4.0 m s−1 . Figure 17 exhibits

372

the experimental wave structure on the granular film measured by Shi et al. [36] and that obtained from the SPH

373

simulations using both the elastoplastic and phase-change constitutive models. It is noted that the wave structures

374

obtained from the two simulations are in very good qualitative agreement with the experimental results.

375

The reason for the wave structures seen in Fig.

17 is that the initial configuration of the granular media is

376

antisymmetric with respect to vertical z-direction and the periodic oscillation of the resulting disturbance results in

377

the generation of wave structure. Furthermore, the effects of the amplitude and the frequency of the disturbance

378

in the granular jet on the wave structure have also been investigated using SPH similations. The results of these

379

investigations are summarized in Fig. 18. In particular, Fig. 18 (a1–a4) exhibit the wave structures for a frequency of

380

f0 = 160 Hz for a series of increasing values of the amplitude A0 : namely, for A0 = 0.0005 m, 0.001 m, 0.002 m,

381

and 0.003 m. Similarly, Fig. 18 (b1–b4) display the wave structure for an amplitude of A0 = 0.002 m for a series

382

of increasing frequencies f0 : namely, for f0 = 80 Hz, 120 Hz, 160 Hz, and 240 Hz. A careful examination of these

383

results shows that the wave structure disappears when the amplitude of the disturbance decreases below a certain

384

value as the disturbance then becomes too small to generate a wave structure (cf. Fig. 18 (a1)). The wavelength of the

385

wave structure on the granular film increases with a decrease in the frequency of the disturbance as is evident from a

386

perusal Fig. 18 (b1–b4). The results reported here demonstrate that the SPH methodology used in conjunction with

387

the phase-change and elastoplatic constitutive models can be used to predict the nature of the wave structures on a

388

granular film.

389

6. Conclusion

390

Several investigations were conducted in the application of a GPU-accelerated SPH methodology for the 3D

391

simulation of granular impact problems involving a large number of particles. To this purpose, the elastoplastic

392

and phase-change constitutive models have been incorporated in the SPH code. The CUDA parallel programming 15

393

interface for an NVIDIA GPU has been used to improve the computational efficiency of the SPH computations for

394

granular flows. To demonstrate the feasibility of the proposed approach, SPH simulations were conducted for four

395

different cases involving granular flow. The results of these simulations have been compared to available experimental

396

data and to numerical data obtained from alternative simulation methdologies.

397

The key conclusions of this study are as follows:

398

• Fully 3D granular impact problems involving large deformations can be addressed successfully using the SPH

399

methodology with either the phase-change or the elastoplastic constitutive models. To the best of our knowl-

400

edge, this is the first time that the phase-change constitutive model has been implemented within SPH frame-

401

work and successfully used for the simulation of 3D granular impact problems.

402

• The SPH methodology can be used to efficiently simulate granular impact problems with a large number of

403

particles on readily available commodity GPUs. The 3D impact of a granular jet with a wave structure on the

404

granular film and the head-on collisions of dense granular jets have been successfully simulated using the SPH

405

methodology for the first time.

406

407

• An SPH simulation is about 160 times as computationally efficient as the same simulation run on a single CPU even with a seriously outdated graphic card.

408

It is expected that a number of alternative constitutive models (e.g., Mohr-Coulomb model) will be implemented

409

in the future in order to investigate the influence and impact of the constitutive models in obtaining high-fidelity SPH

410

simulations of the dynamical behavior of granular media. Finally, our current implementation of the SPH methodology

411

for granular flow can be extended so that it can be executed either on multiple GPUs on a single computer or on

412

multiple GPUs on different computational nodes.

413

Acknowledgements

414

The first author is supported by the China Scholarship Council (No. 201506030072) and the Natural Sciences

415

and Engineering Research Council of Canada (NSERC). This work was made possible by the facilities of the Shared

416

Hierarchical Academic Research Computing Network (SHARCNET) and Compute/Calcul Canada.

417

418

419 420

References [1] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: Theory and application to non-spherical stars, Monthly Notices of the Royal Astronomical Society 181 (3) (1977) 375–389.

421

[2] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.

422

[3] H. H. Bui, R. Fukagawa, K. Sako, S. Ohno, Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geo-

423

material using elastic–plastic soil constitutive model, International Journal for Numerical and Analytical Methods in Geomechanics 32 (12)

424

(2008) 1537–1570.

16

425

[4] H. H. Bui, K. Sako, R. Fukagawa, J. C. Wells, SPH-based numerical simulations for large deformation of geomaterial considering soil-

426

structure interaction, in: Proceeding of the 12th International Conference of International Association for Computer Methods and Advances

427

in Geomechanics, IACMAG, 2008, pp. 570–578.

428

[5] C. Peng, W. Wu, H. S. Yu, Hypoplastic constitutive model in SPH, in: Computer Methods and Recent Advances in Geomechanics: Pro-

429

ceedings of the 14th International Conference of International Association for Computer Methods and Recent Advances in Geomechanics,

430 431 432 433 434 435 436 437 438 439 440 441 442

IACMAG, 2014, pp. 1863–1868. [6] C. Peng, X. G. Guo, W. Wu, Y. Q. Wang, Unified modelling of granular media with smoothed particle hydrodynamics, Acta Geotechnica 11 (6) (2016) 1231–1247. [7] H. F. Fan, G. L. Bergel, S. F. Li, A hybrid peridynamics – SPH simulation of soil fragmentation by blast loads of buried explosive, International Journal of Impact Engineering 87 (2016) 14–27. [8] H. F. Fan, S. F. Li, A peridynamics-SPH modeling and simulation of blast fragmentation of soil under buried explosive loads, Computer Methods in Applied Mechanics and Engineering 318 (2017) 349–381. [9] J. Y. Chen, F. S. Lien, Simulations for soil explosion and its effects on structures using SPH method, International Journal of Impact Engineering 112 (2018) 41–51. [10] J. Y. Chen, C. Peng, F. S. Lien, Simulations for three-dimensional landmine detonation using the SPH method, International Journal of Impact Engineering 126 (2019) 40–49. [11] F. Pacheco-V´azquez, J. Ruiz-Su´arez, Impact craters in granular media: Grains against grains, Physical Review Letters 107 (21) (2011) 218001.

443

[12] J. Ellowitz, Head-on collisions of dense granular jets, Physical Review E 93 (1) (2016) 012907.

444

[13] Z. H. Shi, W. F. Li, Y. Wang, H. F. Liu, F. C. Wang, DEM study of liquid-like granular film from granular jet impact, Powder Technology 336

445

(2018) 199–209.

446

[14] P. A. Cundall, O. D. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65.

447

[15] S. Dunatunga, K. Kamrin, Continuum modelling and simulation of granular flows through their many phases, Journal of Fluid Mechanics

448 449 450 451 452 453 454 455 456 457 458 459 460 461 462

779 (2015) 483–513. [16] T. Harada, S. Koshizuka, Y. Kawaguchi, Smoothed particle hydrodynamics on GPUs, in: Computer Graphics International, Vol. 40, SBC Petropolis, 2007, pp. 63–70. [17] A. C. Crespo, J. M. Dominguez, A. Barreiro, M. G´omez-Gesteira, B. D. Rogers, GPUs, a new tool of acceleration in CFD: Efficiency and reliability on smoothed particle hydrodynamics methods, PloS One 6 (6) (2011) e20685. [18] Y. He, A. E. Bayly, A. Hassanpour, F. Muller, K. Wu, D. Yang, A GPU-based coupled SPH-DEM method for particle-fluid flow with free surfaces, Powder technology 338 (2018) 548–562. [19] X. Xia, Q. Liang, A GPU-accelerated smoothed particle hydrodynamics (SPH) model for the shallow water equations, Environmental modelling and software 75 (2016) 28–43. [20] Q. Xiong, B. Li, J. Xu, GPU-accelerated adaptive particle splitting and merging in SPH, Computer Physics Communications 184 (7) (2013) 1701–1707. [21] N. Govender, D. N. Wilke, S. Kok, Blaze-DEMGPU: Modular high performance DEM framework for the GPU architecture, SoftwareX 5 (2016) 62–66. [22] Y. He, T. J. Evans, A. B. Yu, R. Y. Yang, A GPU-based DEM for modelling large scale powder compaction with wide size distributions, Powder Technology 333 (2018) 219–228.

463

[23] J. J. Monaghan, An introduction to SPH, Computer Physics Communications 48 (1) (1988) 89–96.

464

[24] J. M. Dom´ınguez, A. J. Crespo, M. G´omez-Gesteira, J. C. Marongiu, Neighbour lists in smoothed particle hydrodynamics, International

465

Journal for Numerical Methods in Fluids 67 (12) (2011) 2026–2042.

466

[25] J. J. Monaghan, Simulating free surface flows with SPH, Journal of Computational Physics 110 (2) (1994) 399–406.

467

[26] L. D. Libersky, A. G. Petschek, T. C. Carney, J. R. Hipp, F. A. Allahdadi, High strain Lagrangian hydrodynamics: A three-dimensional SPH

17

468 469 470 471 472 473 474 475 476 477 478 479 480

code for dynamic material response, Journal of Computational Physics 109 (1) (1993) 67–75. [27] S. Adami, X. Y. Hu, N. A. Adams, A generalized wall boundary condition for smoothed particle hydrodynamics, Journal of Computational Physics 231 (21) (2012) 7057–7075. [28] C. Peng, W. Wu, H. S. Yu, C. Wang, A SPH approach for large deformation analysis with hypoplastic constitutive model, Acta Geotechnica 10 (6) (2015) 703–717. [29] L. Zhan, C. Peng, B. Y. Zhang, W. Wu, Three-dimensional modeling of granular flow impact on rigid and deformable structures, Computers and Geotechnics 112 (2019) 257–271. [30] C. Peng, S. Wang, W. Wu, H. S. Yu, C. Wang, J. Y. Chen, LOQUAT: An open-source GPU-accelerated SPH solver for geotechnical modeling, Acta Geotechnica. [31] G. Lube, H. E. Huppert, R. S. J. Sparks, M. A. Hallworth, Axisymmetric collapses of granular columns, Journal of Fluid Mechanics 508 (2004) 175–199. [32] M. G´omez-Gesteira, A. J. Crespo, B. D. Rogers, R. A. Dalrymple, J. M. Dominguez, A. Barreiro, SPHysics–development of a free-surface fluid solver–part 2: Efficiency and test cases, Computers and Geosciences 48 (2012) 300–307.

481

[33] G. R. Liu, M. B. Liu, Smoothed particle hydrodynamics: a meshfree particle method, World scientific, 2003.

482

[34] S. Moriguchi, R. I. Borja, A. Yashima, K. Sawada, Estimating the impact force generated by granular flow on a rigid obstruction, Acta

483 484 485 486

Geotechnica 4 (1) (2009) 57–71. [35] A. H. F. Neto, R. I. Borja, Continuum hydrodynamics of dry granular flows employing multiplicative elastoplasticity, Acta Geotechnica 13 (5) (2018) 1027–1040. [36] Z. Shi, W. Li, H. Liu, F. Wang, Liquid-like wave structure on granular film from granular jet impact, AIChE Journal 63 (8) (2017) 3276–3285.

18

Algorithm 1 Computation of all particle interactions for particle i Output: Update the physical variables for particle i 1:

global void ParticleInteraction cuk(. . . );

2:

index = blockIdx.x ∗ blockDim.x + threadIdx.x;

3:

if (index < particle number) then

4:

Find the neighbour particles;

5:

Calculate the rate of change of the physical quantities;

6:

Update the physical variables;

19

Table 1: Summary of the initial configuration used for four different simulations of the 3D axisymmetric collapse of dry sand.

Case

h0 (mm)

d0 (mm)

a = h0 /d0

Initial particle spacing (mm)

Number of particles

1

0.10

0.20

0.5

4.00

482,430

2

0.20

0.20

1.0

5.00

451,969

3

0.40

0.20

2.0

6.25

530,109

4

0.80

0.20

4.0

8.00

475,329

Table 2: The parameters of the elastoplastic and phase-change constitutive models for the sand collapse problem [15].

Elastoplastic constants

Phase-change parameters

ρ0 (kg m−3 )

E (MPa)

ν

φ (◦ )

ψ (◦ )

c (kPa)

ρc (kg m−3 )

ξ (m1/2 kg−1/2 )

µs

µ2

2450

20.0

0.30

37.0

0.0

0.0

1500

1.1233

0.6419

0.8435

Table 3: Computational efficiency of the GPU-accelerated SPH code for simulation of sand collapse with an aspect ratio of a = 0.5.

Case

Initial particle spacing (mm)

Number of particles

FPS (s−1 )

T (s)

1

6.25

170,221

22.32

127

2

5.00

286,333

12.19

288

3

4.00

482,430

6.78

653

4

3.00

752,187

3.96

1324

Table 4: Initial configuration for three different simulations of head-on collisions of dense granular jets.

Scenario

R (mm)

k

v0 (m s−1 )

Initial particle spacing (mm)

Number of particles

1

0.10

0.50

10.0

4.00

369,000

2

0.10

1.00

10.0

4.00

592,800

3

0.10

0.25

10.0

4.00

314,400

20

487

488

489

List of figures: Figure 1. Sketch of the cell-linked list method for determination of the potential neighbors of a particle of interest (highlighted in blue).

490

Figure 2. Sketch of the dummy particles used to represent a solid wall boundary.

491

Figure 3. Flow diagram depicting the implementation of a 3D SPH methodology on a GPU.

492

Figure 4. Initial setup of the experiment and simulation for the 3D axisymmetric collapse of dry sand.

493

Figure 5. Predictions of time development of 3D sand collapse obtained using the phase-change model for a = 0.5

494

495

496

497

498

499

500

501

502

(a1–a4) and 4.0 (b1–b4). Figure 6. The predicted stress distributions obtained in the SPH framework using the elastoplastic constitutive model for a = 0.5. Figure 7. Predicted and experimental non-dimensionalized runout distance d∞ versus initial aspect ratio a. Experimental results are taken from Lube et al. [31]. Figure 8. Predicted and experimental non-dimensionalized runout distance h∞ versus initial aspect ratio a. Experimental results are taken from Lube et al. [31]. Figure 9. The speed-up on a single GPU compared to that on a single CPU of SPH simulations as a function of the number of particles N used in the simulation.

503

Figure 10. The initial configuration of sand on an inclined flume. All the dimensions shown here are in cm.

504

Figure 11. The magnitude of the velocity of granular impact on 45◦ inclined flume obtained using SPH with the

505

phase-change constitutive model at six different times: namely, at 0.2 s (a1), at 0.6 s (a2), at 1.0 s (a3), at 1.4 s (a4), at

506

1.8 s (a5), and at 2.2 s (a6).

507

Figure 12. Time variation of the impact force for a flume inclination of 45◦ .

508

Figure 13. A sketch of the ejecta angle and the velocity u0 of the impact center.

509

Figure 14. The magnitude of the velocity arising from the head-on collision of two dense granular jets predicted

510

511

512

513

514

using SPH with the phase-change constitutive model for an initial configuration given by scenario 1. Figure 15. The stress and strain distributions arising from head-on granular impact obtained using SPH with the phase-change constitutive model. Figure 16. Numerical simulation of the process of densification of two colliding granular jets obtained using SPH with the phase-change constitutive model.

515

Figure 17. Wave structures on a granular film obtained from (a) an experimental study [13] and from SPH sim-

516

ulations with (b) the phase-change and (c) the elastoplastic constitutive models. The numerical simulations were

517

conducted with initial amplitude of A0 = 0.002 m and a frequency of f0 = 160 Hz.

518

Figure 18. Wave structures on a granular film obtained at an initial frequency of f0 = 160 Hz and for amplitudes

519

of A0 = 0.0005 m (a1), 0.001 m (a2), 0.002 m (a3), and 0.003 m (a4) and at an initial amplitude of A0 = 0.002 m and

520

for frequencies of f0 = 80 Hz (b1), 120 Hz (b2), 160 Hz (b3), and 240 Hz (b4).

21

Highlight • A GPU-accelerated SPH code has been developed to simulate granular flows and impact. • An elastoplastic and a phase-change models were used to model granular materials. • Four cases have been tested to validate the GPU-accelerated SPH. • The numerical results were compared against a variety of published data. • A case run on a outdated GPU was 160 times more efficient than run on a single CPU.

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: