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 Numerical modeling began in the 1960’s Numerically assisted analytical methods Finite difference modeling (marker and cell)
Finite Difference and FiniteElements 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!
HOS/SNOW (MIT, JHU) Wave Theory & Modeling • Cross validation • Shear instability characteristics • Parameter for simulations Perturbation Theory O(/εn) DNS/LES O() ~ O(10) • Kinematic and dynamic boundary conditions on interface • Mud energy dissipation for wavefield simulation • Cross validation • Interfacial boundary conditions • Parameters for simulations • Dissipation model development Direct Wavefield Simulation O(100) Laboratory Experiments O(10) / Field Measurements O(100)
Meshfree Lagrangian Numerical Method for Wave Breaking Fluid is described by quantities at discrete nodes Approximated by a summation interpolant; other options: MLS Radius of Kernel function, W Water Particles r 2h Boundary Particles
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= a x1 + b y2= a x2 + b a=(y2-y1)/(x2-x1), b= (y1 x2-y2 x1)/(x2-x1), Use this equation to find y values for any x1 < x < x2 y2 y1 x1x2
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 Data irregularly spaced 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
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 ofx
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
Varying h Values 1.0 .3 1.5 .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.
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): Quadratic: (discontinuous slope at q=1)
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.
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