Transcript Notes - Department of Civil Engineering

IGERT: SPH and Free Surface
Flows
Robert A. Dalrymple
Civil Engineering
Wave Theories
Flat Bottom Theories
Small Amplitude (Airy Theory, 1845)
Shallow Water (Boussinesq, 1871, KdV, 1895)
Intermediate and Deep Depths (Stokes, 1856, Stokes V)
Wave Modeling
N
u
m
e
r
i
c
a
l
m
o
d
e
l
i
n
g
Stream Function on the Web
Finite Difference and Finite
Elements in the 70’s
• 2-D to 3-D
• Multiple numbers of waves
• Shoaling, refraction, diffraction
Parabolic Modeling (REF/DIF)
Scripps Canyon--NCEX
http://science.whoi.edu/users/elgar/NCEX/wp-refdif.html
Boussinesq Model
Much more computing time
Empirical breaking
Wave-induced currents!
Wave Theory & Modeling
Perturbation
Theory
O(/εn)
• Cross validation
• Interfacial boundary
conditions
• Parameters for
simulations
• Dissipation model
development
• Cross validation
• Shear instability characteristics
• Parameter for simulations
Direct Wavefield
Simulation
O(100)
DNS/LES
O() ~ O(10)
• Kinematic and dynamic
boundary conditions on
interface
• Mud energy dissipation
for wavefield simulation
Laboratory Experiments O(10) / Field Measurements O(100)
HOS/SNOW (MIT, JHU)
Meshfree Lagrangian Numerical Method
for Wave Breaking
Fluid is described by quantities
at discrete nodes
Radius of
Kernel function, W
Water
Particles
r
2h
Boundary
Particles
Approximated by a
summation interpolant;
other options: MLS
Topics
•
•
•
•
•
Meshfree methods
Interpolation methods
SPH modeling
Results
GPU : the future
Mesh-Free Methods
Smoothed particle hydrodynamics (SPH) (1977)
Diffuse element method (DEM) (1992)
Element-free Galerkin method (EFG / EFGM) (1994)
Reproducing kernel particle method (RKPM) (1995)
hp-clouds
Natural element method (NEM)
Material point method (MPM)
Meshless local Petrov Galerkin (MLPG)
Generalized finite difference method (GFDM)
Particle-in-cell (PIC)
Moving particle finite element method (MPFEM)
Finite cloud method (FCM)
Boundary node method (BNM)
Boundary cloud method (BCM)
Method of finite spheres (MFS)
Radial Basis Functions (RBF)
Particle Methods
Discrete Element Method
Molecular Dynamics
SPH
Vortex Methods
Why Interpolation?
• For discrete models of continuous systems, we
need the ability to interpolate values in
between discrete points.
• Half of the SPH technique involves
interpolation of values known at particles (or
nodes).
Interpolation
• To find the value of a function between known
values.
Consider the two pairs of values (x,y):
(0.0, 1.0), (1.0, 2.0)
What is y at x = 0.5? That is, what’s (0.5, y)?
Linear Interpolation
Given two points, (x1,y1), (x2,y2):
Fit a straight line between the points.
y(x) = a x +b
y1
y1= a x1 + b
y2= a x2 + b
x1
y2
x2
a=(y2-y1)/(x2-x1), b= (y1 x2-y2 x1)/(x2-x1),
Use this equation to find y values for any
x1 < x < x 2
Polynomial Interpolants
Given N (=4) data points,
Find the interpolating function that goes through the points:
If there were N+1 data points, the function would be
with N+1 unknown values, ai, of the Nth order polynomial
Polynomial Interpolant
Force the interpolant through the four points to get four equations:
Rewriting:
The solution for a is found by inverting p
Example
Data are: (0,2), (1,0.3975), (2, -0.1126), (3, -0.0986).
Fitting a cubic polynomial through the four points gives:
Polynomial Fit to Example
Exact: red
Polynomial fit: blue
Beware of Extrapolation
Exact: red
An Nth order polynomial has N roots!
Least Squares Interpolant
For N points, we will use a lower order fitting polynomial of order m < (N-1).
The least squares fitting polynomial is similar to the exact fit form:
Now p is N x m matrix. Since we have fewer unknown coefficients as data points, the
interpolant cannot go through each point. Define the error as the amount of “miss”
Sum of the (errors)2:
Least Squares Interpolant
Minimizing the sum with respect to the coefficients a:
Solving,
This can be rewritten in this form,
which introduces a pseudo-inverse.
Reminder:
for cubic fit
Question
Show that the equation above leads to the
following expression for the best fit straight
line:
Cubic Least Squares Example
x: -0.2 .44 1.0 1.34 1.98 2.50 2.95 3.62 4.13 4.64 4.94
y: 2.75 1.80 -1.52 -2.83 -1.62 1.49 2.98 0.81 -2.14 -2.93 -1.81
Data
irregularly
spaced
Least Squares Interpolant
Cubic Least Squares Fit: * is the fitting polynomial
o is the given data
Exact
Piecewise Interpolation
Piecewise polynomials: fit all points
Linear: continuity in y+, y- (fit pairs of points)
Quadratic: +continuity in slope
Cubic splines: +continuity in second derivative
RBF
All of the above, but smoother
Moving Least Squares Interpolant
are monomials in x for 1D (1, x, x2, x3)
x,y in 2D, e.g. (1, x, y, x2, xy, y2 ….)
Note aj are functions of x
Moving Least Squares Interpolant
Define a weighted mean-squared error:
where W(x-xi) is a weighting function that decays
with increasing x-xi.
Same as previous least squares approach, except for W(x-xi)
Weighting Function
q=x/h
Moving Least Squares Interpolant
Minimizing the weighted squared errors for the coefficients:
,
,
,
Moving Least Squares Interpolant
Solving
The final locally valid interpolant is:
MLS Fit to (Same) Irregular Data
h=0.51
Given data: circles; MLS: *; exact: line
1.0
.3
Varying h Values
.5
1.5
Basis for Smoothed Particle Hydrodynamics
SPH
From astrophysics (Lucy,1977; Gingold and Monaghan, 1977))
An integral interpolant
(an approximation to a Dirac delta function):
kernel
The Kernel (or Weighting Function)
h
Compact support: 2D-circle of radius h
3D-sphere of radius h
1D-line of length 2h
Fundamental Equation of SPH
where W(s-x,h) is a kernel, which behaves like Dirac  function.
Delta Function (1 D)
Kernel Requirements (Monaghan)
Monotonically decreasing with distance |s-x|
Symmetric with distance
Numerical Approximations of the Integrals
Partition of unity
The incremental volume: mj/j , where the mass is constant.
which is an interpolant!
Kernels
Gaussian:
Not compactly supported -- extends to infinity.
Kernels
Moving Particle Semi-implicit (MPS):
(discontinuous slope at q=1)
Quadratic:
Spline Kernel
Same kernel as used
in MLS interpolant.
Gradients
Given SPH interpolant:
Find the gradient directly and analytically:
O.K., but when uj is a constant, there is a problem.
Gradients (2)
(A)
To fix problem, recall partition of unity equation:
The gradient of this equation is zero:
So we can multiply by ui = u(x) and subtract from Eq.(A)
Consistency
Taylor series expansion of u(x,t) about point s (1-D):
Consistency conditions are then:
Integrals of all higher moments must be zero.
MLS and Shepard Interpolant
Governing Equations
Kinematics
Conservation of Mass
Conservation of Momentum
Equation of State: Pressure = f()
Kinematic Equation
for i=1, np
Mass Conservation
Conservation of Mass
Integrate both sides by the kernel and integrate over the domain:
Next use Gauss theorem:
Conservation of Mass (2)
0
Put in discrete form:
j
Conservation of Mass (3)
From before (consistency condition):
The derivative of this expression is
Multiplying by ui and adding to previous conservation
equation:
Conservation of Mass (4)
To determine the density at a given point s:
A more accurate approach is:
During computations, the density can be regularly redetermined through
this expression, due to Randles and Liberski.
Conservation of Momentum
Equations of motion:
Multiply by  and the kernel and integrate as before:
Subtract the zero sum:
Conservation of Momentum
Equations of motion:
Multiply by
and the kernel and integrate as before:
Conservation of Momentum (2)
Another formulation (used in JHU SPH):
Equation of State: Weak Compressibility
No need to solve PDE for pressure
which implies a speed of sound, CS:
where g is usually 7 and B is a constant
Closure Submodels
Viscosity generally accounted for by an artificial empirical term (Monaghan 1992):
  cij ij

 ij   ij
 0

v ij .rij  0
v ij .rij  0
ij 
hvij .rij
rij2   2
Sub-Particle Scale (SPS) Turbulence Model
•
Spatial-filter over the governing equations (compressible):
(Favre-averaging)
~
f  f 
= SPS stress tensor with elements:
~ 2~
*
Sij = strain tensor
 ij   2 t Sij  3 Skkij  23 kij
 t  Cs l  S
• Eddy viscosity:
• Smagorinsky constant: Cs  0.12
2
S  2Sij Sij 
See Lo & Shao (2002), Gotoh et al. (2002) for incompressible SPS
1/ 2
Neighbor Lists
Which particles within 2h of particle i?
In principle, each particle interacts with every other particle
in the domain: N x N calculations (“all pair search”)
Linked List:
“Cells”: Examine particles in neighboring cells;
only 9 cells to examine.
Tree Search
Develop a tree. Search geometrically.
Neighbor List
Link list is reconstructed each time step. Cells are 2h x 2h in size.
from SPHysics Users Manual
Time Stepping
Discretize time:
Take an Euler step (first order derivative expression):
Use this to get a value of
Then make a corrector step:
Monaghan; more
accurate than Euler step
Beeman’s Method
Better velocity estimate (but more storage, and needs an+1)
Predictor-corrector version (for velocity):
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. 1994)
Leonard-Jones forces (Monaghan 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)
y
(slip BC)
b
(Can use kernel normalization techniques to reduce interpolation errors at
the boundaries, Bonet and Lok 2001)
SPHYSICS
DEVELOPED JOINTLY BY
RESEARCHERS AT:
Johns Hopkins
University
(USA)
Universidade de
Vigo
(Spain)
University of
Manchester
(U.K.)
Universita di
Roma “La Sapienza”
(Italy)
Model Validation – hydrodynamics model
– Star points are experiment data
Model Validation – hydrodynamics model
– Star points are experiment data
Model Validation – hydrodynamics model
– Star points are experiment data
Bore in a Box
Tsunami represented by
a dam break
Bore in a Box: 2004
3D Weak Plunger
Coherent Structures
, velocity gradient tensor
The rate of strain tensor and the rotation tensor are
Q criterion
Second invariant of the velocity gradient tensor
Positive values of Q imply vorticity greater than strain
Coherent Structures: Q Criterion
Paradigm Shift: GPU Models
CPU vs GPU
Floating point operations per second: GPU vs CPU
Nvidia CUDA Programming Guide, 1.1
CPU/GPU Architectures
Nvidia CUDA Programming Guide, 1.1
Bore in a Box
Real time: CPU vs GPU on a laptop
Bore in a Box
677,360 particles
Nvidia Tesla Computing Card
4 Gb/240 = 16 Mb/ core