SPH-WBL

							Simulating complex surface flow by
Smoothed Particle Hydrodynamics
  & Moving Particle Semi-implicit
             methods


       Benlong Wang    Kai Gong     Hua Liu

            benlongwang@sjtu.edu.cn

           Shanghai Jiaotong University
                   Contents
• Introduction
• SPH & MPS methods
• Parallel strategy and approaches
  – SPH:
  – MPS:
• Numerical results
  –   2D dam breaking
  –   2D wedge entry
  –   3D cavity flow
  –   3D dam breaking
   Modeling free surface flows
• Multiphase flows:
        MAC, VOF, LevelSet etc.



• ALE

• Meshless methods
  & particle methods

     SPH & MPS        LBM
                        Kernel function
                                                                (r )
   f (r )   f (r ')W (r  r '; h)dV             h
             

             f (r ')W (r ' r; h)dV                                   dx
                                                      W
            f a   f bWabVb
                   b
                                        0
• Properties:
   – Narrow support                     -3   -2       -1   0       1     2   3

   –    W (r  r '; h)dV  1
       
   –    W (r; h) decreases monotonously as r increase
   – h->0, Dirac delta function
                 expression of derivatives

                         W                   f '(r )   f '(r ')W (r  r '; h)dV
            W’                   h                    
                                                                                         0

0                                                     f (r ')W (r  r '; h) nds   f (r ')W '(r  r '; h)dV
                                                                                   



                                                      f (r ')W '(r  r '; h)dV   f (r ')W '(r ' r; h)dV
                                                                                    




                                                               f a '   f bWab 'Vb
-3     -2    -1     0        1       2   3
                                                                      b



         f ( x)dV   f ( x)V
                    b
                         b       b               h
                                                                        1.3 ~ 1.5
                                                 dx
        Integral   Summation                                              3.0       2h 130+ (2D)
     Trapeze like quadrature formula
              Correction and Consistance
                              ——advanced topic …

                                                             f a '   f bWab 'Vb
                                                                      b


 0
                                                             f  const              fa '  0

                                                             f a '   ( f b  f a )Wab 'Vb
                                                                      b



                                                            f  ax  by  c
 -3      -2      -1       0        1      2       3




 a f a  (  a f a )  f a a   b fbWabVb  f a  bWabVb   b ( fb  f a )WabVb
                               b                      b              b
                        Lists of kernel function
    0.5
                                                                 Cubic spline                     2h
    0.4


    0.3                                                         Quartic spline                  2.5h
w




    0.2
                                                            Fifth order B-spline                  3h

    0.1
                                                           Truncated Gaussian                     
     0
      0     0.5     1     1.5     2        2.5     3
                           h                                              3 2 3 3
                    s2          2                                   1  2 s  4 s    0  s 1
               exp   2   exp   2                                  
             1      h           h                              10  1
    W ( s)                                        3h   W (s)              (2  s) 2   1 s  2
              2 h2   2        2                             7 h 2  4
               h           exp   2 
                       2         h                                            0         2s
                                                                         
                                                                         
      Hydrodynamics governing equations
  dva      p                              dxa
       g  a  0 2 va                        va
   dt      a                               dt
                mb                              2 0 (va  vb )(ra  rb )  aWab
       g         ( pa  pb ) aWab   mb
           b  a b                           a  b          ra  rb
                                                                       2
                                        b


SPH:      weakly compressible method: State Equation
                                                           1
da                                                   
      a va   mb (va  vb )  aWab pa      1
                                                    a
                                                                                Ma < 0.1
 dt              b                               0 
                                                         
                                                          
MPS: projection method: Pressure Poisson Equation
                                                   1     va
                                                  pa  
    va  0                                        a     t
                                 2 ( pa  pb )(ra  rb )  aWab       1   mb
                         mb
                         b      a b         ra  rb
                                                        2
                                                                         (va  vb )aWab
                                                                      t b b
   Link-List neighbour search


   L




                             SPH: the most time consuming part
back ground mesh (L X L)           ~90%

L=2h, 3h, support distance   MPS: generally less than PPE solver
            Boundary Condition
• Sym:        ghost particles,


              p         p
                 0          g
              x         y

             v '  v   vn '   vn

• Free surface, p0
  Identify the surface particle: 95% const. density

                                                                mb
                                                 a   bWab
                                                      b         b
Large Scale Computation
(a few millions particles)
share memory architecture
(NEC SX8: 8 nodes, 128G RAM)
(Dell T5400: 2 Quad cores Xeon 16G RAM)



• SPH
   – Particle lists partition, NOT domain partition


• MPS
   – parallel ICCG method
Black-box   Parallel Sparse Matrix Solver

  Why not Domain decomposition ?

                 SPH Method
                 Lagrangian Method
                 Large deformation
                 Continue changing domain
                 Complex domain structure
                 So, Black-box solver
                 give me a matrix, I will solve it in parallel…
     PPE solver : ICCG method

Ax  b                Direct solver or Iterative solver


    Sparse symmetric positive definite matrix

•   Precondition ILU(0)
•   Forward and backward substitutions
•   Inner products
•   Matrix-vector products          Parallel
•   Vector updates
                Coloring
• COLOR: Unit of independent sets.
• Any two adjacent nodes have different colors. Elements
  grouped in the same ―color‖ are independent from each
  other, thus parallel/vector operation is possible.
• Many colors provide faster convergence, but shorter
  vector length.
  Main Idea of the Coloring
Algebraic Multi-Color Ordering
 The number of the colors has a lower boundary
      the max bandwidth of the sparse matrix
  Any two adjacent nodes have different colors



           2h
                             T. Iwashita & M. Shimasaki
                             2002 IEEE Trans. Magn.

                The connection info could be obtained from the
                distribution of non-zeros in the sparse matrix
    bcsstk14
n=1806,nnz=63454
MC=50




MC=180
       Parallelized ICCG with AMC
   Ax  b                                       A          LDL      T



 C1,1 C1,2            C1,nc        x1        D1                        
C                     C2,nc       x         L          D2              
 2,1 C2,2                          2         2,1                       
                                                                       
                                                                       
                                                                       
Cnc ,1 Cnc ,2         Cnc ,nc      xnc       Lnc ,1   Lnc ,2       Dnc 
                                                                       

 Forward and backward substitutions: parallelized in each color

                                                    ic 1
                                                                
             Ly  r                yic  Dic  ric   Lic,k rk 
                                           1

                                                    k 1       
SPH Parallel Strategy: OpenMP


     Almost linear speedup



MPS Parallel Strategy: OpenMP
Numerical Results

•   2D dam breaking
•   2D wedge water entry
•   3D cavity flow
•   3D dam breaking
Dambreaking Test




                   Surge front location
Water entry of a wedge




    4.5M particles Speed up around 7
       Dell T5400 2 Xeon Quadcores
        8000
                                                                          6.5




        6000                                                               6
Fy(N)




                                                                 v(m/s)
        4000                                                              5.5




        2000                                                               5
                                    SPH Results
                                    Analysis
                                    Experiment                                              SPH Results
                                    Oger' s results                                         Experiment
           0
               0   0.005   0.01          0.015    0.02   0.025            4.5
                                                                                0   0.005   0.01          0.015   0.02   0.025
                                  t(s)                                                             t(s)
         3D Cavity Flow: Re=400

                               45 X 45 X 45 nodes




Yang Jaw-Yen et al. 1998
J. Comput. Phys. 146:464-487
                                   h/dx=1.5
3D Dambreaking Tests




        Kleefsman, K.M.T. et al 2005
        J. Comput. Phys. 206:363-393
                  0.6
                                                               H4
                                                               H3
                                         MARIN Exp. Results    H2
                  0.5
                                                               H1
Water Level (m)
                  0.4

                  0.3

                           SPH Results
                  0.2

                  0.1

                  0.0
                     0.0          0.5       1.0          1.5    2.0
                                          Time (s)
             Conclusions
• 2D code is developed for both SPH and
  MPS methods
• 3D code is developed for complex free
  surface flows
• Computation costs of SPH is generally
  cheaper than MPS method
• Good agreements are obtained, a
  promising method for complex free surface
  flows.