IGERT: SPH and Free Surface Flows

1. IGERT: SPH and Free Surface Flows Robert A. Dalrymple Civil Engineering

2. Wave Theories Flat Bottom Theories Small Amplitude (Airy Theory, 1845) Shallow Water (Boussinesq, 1871, KdV, 1895) Intermediate and Deep Depths (Stokes, 1856, Stokes V)

3. Wave Modeling Numerical modeling began in the 1960’s Numerically assisted analytical methods Finite difference modeling (marker and cell)

4. Stream Function on the Web

5. Finite Difference and FiniteElements in the 70’s • 2-D to 3-D • Multiple numbers of waves • Shoaling, refraction, diffraction

6. Parabolic Modeling (REF/DIF) Scripps Canyon--NCEX http://science.whoi.edu/users/elgar/NCEX/wp-refdif.html

7. Boussinesq Model Much more computing time Empirical breaking Wave-induced currents!

8. 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)

9. 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

10. Topics • Meshfree methods • Interpolation methods • SPH modeling • Results • GPU : the future

11. 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)

12. Particle Methods Discrete Element Method Molecular Dynamics SPH Vortex Methods

13. 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).

14. 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)?

15. 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

16. 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

17. Polynomial Interpolant Force the interpolant through the four points to get four equations: Rewriting: The solution for a is found by inverting p

18. Example Data are: (0,2), (1,0.3975), (2, -0.1126), (3, -0.0986). Fitting a cubic polynomial through the four points gives:

19. Polynomial Fit to Example Exact: red Polynomial fit: blue

20. Beware of Extrapolation Exact: red An Nth order polynomial has N roots!

21. 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:

22. 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

23. Question Show that the equation above leads to the following expression for the best fit straight line:

24. 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

25. Least Squares Interpolant Cubic Least Squares Fit: * is the fitting polynomial o is the given data Exact

26. 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

27. 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

28. 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)

29. Weighting Function q=x/h

30. Moving Least Squares Interpolant Minimizing the weighted squared errors for the coefficients: , , ,

31. Moving Least Squares Interpolant Solving The final locally valid interpolant is:

32. MLS Fit to (Same) Irregular Data h=0.51 Given data: circles; MLS: *; exact: line

33. Varying h Values 1.0 .3 1.5 .5

34. 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

35. The Kernel (or Weighting Function) h Compact support: 2D-circle of radius h 3D-sphere of radius h 1D-line of length 2h

36. Fundamental Equation of SPH where W(s-x,h) is a kernel,which behaves like Dirac  function.

37. Delta Function (1 D)

38. Kernel Requirements (Monaghan) Monotonically decreasing with distance |s-x| Symmetric with distance

39. Numerical Approximations of the Integrals Partition of unity The incremental volume: mj/j , where the mass is constant. which is an interpolant!

40. Kernels Gaussian: Not compactly supported -- extends to infinity.

41. Kernels Moving Particle Semi-implicit (MPS): Quadratic: (discontinuous slope at q=1)

42. Spline Kernel Same kernel as used in MLS interpolant.

43. Gradients Given SPH interpolant: Find the gradient directly and analytically: O.K., but when uj is a constant, there is a problem.

44. 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)

45. 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.

46. MLS and Shepard Interpolant

47. Governing Equations Kinematics Conservation of Mass Conservation of Momentum Equation of State: Pressure = f()

48. Kinematic Equation for i=1, np Mass Conservation

49. Conservation of Mass Integrate both sides by the kernel and integrate over the domain: Next use Gauss theorem:

50. Conservation of Mass (2) 0 Put in discrete form: j