Modelling Tsunami Waves using Smoothed Particle Hydrodynamics (SPH)

R.A. DALRYMPLE and B.D. ROGERS Department of Civil Engineering, Johns Hopkins University

Introduction

• Motivation  multiply-connected free-surface flows • Mathematical formulation of Smooth Particle Hydrodynamics (SPH) • Inherent Drawbacks of SPH • Modifications - Slip Boundary Conditions - Sub-Particle-Scale (SPS) Model

Numerical Basis of SPH

• SPH describes a fluid by replacing its continuum properties with locally (smoothed) quantities at discrete Lagrangian locations  meshless • SPH is based on integral interpolants (Lucy 1977, Gingold & Monaghan 1977, Liu 2003)

A

 

  

A

  

r



r

' ,

h



d

r

' (

W

is the smoothing kernel) • These can be approximated discretely by a summation interpolant

A

  

j N

  1

A

  

j

r



r

j

,

h



m j



j

The Kernel (or Weighting Function)

Radius of influence

W

(

r

-

r

’,

h

) Water Particles

r

2

h

Compact support of kernel • Quadratic Kernel

W

 3 2 

h

2 1 4

q

2 

q

 1 

q



r h

,

r

 |

r

a



r

b

|

SPH Gradients

• Spatial gradients are approximated using a summation containing the gradient of the chosen kernel function 

i



A i

  .

u



i

  

j



j m j m j



j



u

i A j



i W ij



u

j

 . 

i W ij

• Advantages are: – spatial gradients of the data are calculated analytically – the characteristics of the method can be changed by using a different kernel

Equations of Motion

• Navier-Stokes equations: d      .

v

d

t

d

v

d

t

  1  

p

 

o

 2

u



F

i

• Recast in particle form as d

m i

d

t

 0 d

r

i

d

t

d 

i

d

t



v

i

  

j m j

   

v



ij ji

   

W ij

 

j m j



v

i



v

j

  

i W ij

d

v

i

d

t

  

j m j

(XSPH) 

p i i

2 

p j



j

2  

ij

 

i W ij



F

i

Closure Submodels

• Equation of state (Batchelor 1974):

p



B

     

o

    1   accounts for incompressible flows by setting B such that speed of sound is

c

 d d

p

  10

v

max Compressibility O(M 2 ) • Viscosity generally accounted for by an artificial empirical term (Monaghan 1992): 

ij

   

c ij

 0

ij



ij

v v

ij ij

.

r

ij

.

r

ij

 0  0 

ij



h

v

ij

.

r

ij r ij

2   2

Dissipation and the need for a

•

Sub-Particle-Scale (SPS) Model

• Description of shear and vorticity in conventional SPH is empirical

Π ij

  

c ij



ij



h



u

i



r ij

2

u



j

 .

r

i

 0 .

01

h

2

r

j

 is needed for stability for free-surface flows, but is too dissipative, e.g. vorticity behind foil

Sub-Particle Scale (SPS) Turbulence Model

• Spatial-filter over the governing equations: (Favre-averaging)

f

~  

f

 D      .

~ D

t D Dt

  1  

P



g

 

o

 2  1   .

τ

S ij

= SPS stress tensor with elements: = strain tensor 

ij

   2 

t

~

S ij

   2 3 ~

S kk



ij

  2 3 

k



ij

 2

S ij S ij

 1 / 2 • Smagorinsky constant:

C s

 0.12 (not dynamic!)

Boundary Conditions

Boundary conditions are problematic in SPH due to: – – the boundary is not well defined kernel sum deficiencies at boundaries, e.g. density • Ghost (or virtual) particles (Takeda

et al.

• Leonard-Jones forces (Monaghan 1994) 1994) • Boundary particles with repulsive forces (Monaghan 1999) • Rows of fixed particles that masquerade as interior flow particles (Dalrymple & Knio 2001) a

f

=

n

R(

y

) P(

x

) (slip BC)

y

b (Can use kernel normalisation techniques to reduce interpolation errors at the boundaries, Bonet and Lok 2001)

Determination of the free-surface

Free-surface defined by   1 2 

water

x

2

h

where   

j



j W



x



x

j



V j

 

j



j W



x



x

j



m j



j

Far from perfect!!

Free-surface

g Caveats

: • SPH is inherently a multiply-connected • Each particle represents an interpolation location of the governing equations

JHU-SPH - Test Case 3

R.A. DALRYMPLE and B.D. ROGERS Department of Civil Engineering, Johns Hopkins University

SPH: Test 3 - case A -



= 0.01

• Geometry aspect-ratio proved to be very heavy computationally to the point where meaningful resolution could not be obtained without

high-performance computing ===>

real disadvantage of SPH • hence, work at JHU is focusing on coupling a depth averaged model with SPH e.g. Boussinesq FUNWAVE scheme • Have not investigated using 

z

<< 

x

, 

y

for particles

SPH: Test 3 - case B



= 0.1

• Modelled the landslide by moving the SPH bed particles (similar to a wavemaker)

t 1 t 2

• Involves run-time calculation of boundary normal vectors and velocities, etc.

• Water particles are initially arranged in a grid-pattern …

Test 3 - case B



= 0.1

• • SPH settings: 

x

= 0.196m, 

t

= 0.0001s, Cs = 0.12

• 34465 particles • Machine Info: – Machine: 2.5GHz

– RAM: 512 MB – Compiler: g77 – cpu time: 71750s ~ 20 hrs

Test 3B



= 0.1 animation

t ND = 0.5

Test 3 comparisons with analytical solution

2 1.5

1 0.5

0 -0.5

0 -1 2 1.5

1 0.5

0 -0.5

0 -1 20 40 60 80 100 120

x (m) t ND = 1.0

SPH Analytical SPH Analytical 20 40 60 80 100 120

x (m)

t ND = 2.5

Test 3 comparisons with analytical solution

2 1.5

1 0.5

0 -0.5

0 -1 Free-surface fairly constant with different resolutions 2 1.5

1 0.5

0 -0.5

0 -1 SPH Analytical 20 40 60 80 100 120

x (m) t ND = 4.5

SPH Analytical 20 40 60

x (m)

80 100 120

Points to note:

• Separation of the bottom particles from the bed near the shoreline • Magnitude of SPH shoreline from SWL depended on resolution • Influence of scheme’s viscosity

JHU-SPH - Test Case 4

R.A. DALRYMPLE and B.D. ROGERS Department of Civil Engineering, Johns Hopkins University

JHU-SPH: Test 4

• Modelled the landslide by moving a wedge of rigid particles over a fixed slope according to the prescribed motion of the wedge • Downstream wall in the simulations • 2-D: SPS with repulsive force Monaghan BC • 3-D: artificial viscosity Double layer Particle BC • did not do a comparison with run-up data

2-D, run 30, coarse animation

8600 particles, 

y

= 0.12m, cpu time ~ 3hrs

2-D, run 30, wave gage 1 data

0.2

0.1

0 -0.1

0 -0.2

0.5

1 1.5

2 2.5

-0.3

-0.4

experimental data SPH -0.5

time (s)

• Huge drawdown • little change with higher resolution  lack of 3-D effects 3

2-D, run 32, coarse animation

10691 particles, 

y

= 0.08m, cpu time ~ 4hrs breaking is reduced at higher resolution

2-D, run 32, wave gage 1 data

0.2

0.1

0 -0.1

0 -0.2

-0.3

-0.4

-0.5

0.5

1 1.5

experimental data SPH

time (s)

2 2.5

3 • Huge drawdown & phase difference • Magnitude of max free-surface displacements is reduced • lack of 3-D effects

3-D, run 30, animation 38175 Ps,



x = 0.1m (desktop) cpu time ~ 20hrs

Conclusions and Further Work

• Many of these benchmark problems are inappropriate for the application of SPH as the scales are too large • Described some inherent problems & limitations of SPH • Develop hybrid Boussinesq-SPH code, so that SPH is used solely where detailed flow is needed