1# SPH example
2# Andreas Aigner (CD Lab Particulate FLow Modelling, JKU)
3# andreas.aigner@jku.at
4
5#variables
6variable      mass universe 0.086 #0.001
7variable      h universe 0.012 #0.0012
8
9variable      lat equal $h/1.2 # = 0.01
10variable      lathalf equal ${lat}*0.5
11variable      wallpos equal 0.1+${lathalf}
12variable      skin equal $h*0.25
13variable      eta equal 0.01*$h*$h
14
15variable      ins01 equal 0.1-${lat}
16variable      ins05 equal 0.5-${lat}
17
18
19atom_style    sph
20atom_modify   map array sort 0 0
21communicate   single vel yes
22
23boundary      f f p   #periodic in z; fixed in y-direction
24newton        off
25
26dimension	    2
27
28units         si
29lattice       sq ${lat}
30
31region        reg block -0.021 0.521     -0.021 0.5    -0.001 0.001  units box
32create_box    2 reg
33
34region        regWallOut block 0.01 0.49   0.01 0.5   -0.001 0.001 side out units box
35region        regWall intersect 2 reg regWallOut
36region        insreg block ${lat} 0.1  ${lat} 0.1  -0.001 0.001 units box
37region        insreg2 block 0.10001 ${ins05}  ${lat} 0.05  -0.001 0.001 units box
38create_atoms  1 region insreg
39create_atoms  1 region insreg2
40create_atoms  2 region regWall
41mass          * ${mass}
42
43group         fluid type 1
44group         wall  type 2
45
46neighbor      ${skin} bin
47
48fix           m1 all property/global speedOfSound peratomtype 20. 20.
49fix           m2 all property/global sl peratomtype $h $h
50fix           m3 all property/global artViscAlpha peratomtype 4.1666e-3 4.1666e-3
51fix           m4 all property/global artViscBeta peratomtype 0. 0.
52fix           m5 all property/global artViscEta  scalar ${eta}
53fix           m6 all property/global tensCorrEpsilon scalar 0.2
54fix           m7 all property/global tensCorrDeltaP peratomtype ${lat} ${lat}
55
56#sph pair style
57# mhu = 1e-3 (water) ~ alpha * cAB * h
58# for c=20 (10*vMax) and h=0.012 --> alpha=0.004166667
59# pair_style    sph kernel_style h [artVisc alpha beta cAB eta] [tensCorr epsilon]
60# artifical viscosity [Monaghan and Gingold (1983)] and tensile correction [Monaghan (2000)]
61
62pair_style 	  sph/artVisc/tensCorr cubicspline2d $h artVisc tensCorr
63#pair_coeff	  1 1 sph/artVisc/tensCorr
64#pair_coeff    1 2 sph/artVisc/tensCorr
65#pair_coeff    2 2 none
66pair_coeff    * *
67
68#sph fixes
69# density
70fix           density all sph/density/continuity
71fix           corr all sph/density/corr shepard every 30
72set           group all meso_rho 1000
73
74# pressure
75# fix         id group style type [if Tait: B rho0 gamma] (according to Monaghan 1994)
76# B = c^2*rho0/gamma
77# for c=20 --> B~60000
78fix		        pressure all sph/pressure Tait 60000. 1000. 7.
79
80#wall
81#region        boxw block 0. 0.5  0. 0.1  0. 0.5 units box
82# fix         id group wall/region/sph region_id r0 D (repulsive force similar to Lennard-Jones potential p1=4, p2=2)
83#fix           boxwall_reg all wall/region/sph boxw ${lat} 5.0
84#fix         id group wall/sph args (e.g: xplane lo hi) r0 D
85fix           wall all wall/sph xplane  ${wallpos} NULL ${lathalf} 1.0
86
87# time integration
88# for c=20 --> dt~1e-5
89timestep      1e-5
90fix           integr fluid nve/sph
91fix           integrWall wall nve/sph/stationary
92
93#gravity
94fix           gravi all gravity 9.81 vector 0.0 -1.0 0.0
95
96#2d simulation
97fix           only2d all enforce2d
98
99# set sph kernels for all fix after all fixes are defined
100set        group all sphkernel cubicspline2d
101
102#output settings, include total thermal energy
103thermo_style	custom step atoms ke vol cpu
104thermo          1000
105thermo_modify	lost ignore norm no
106
107run           1
108dump	      dmp all custom/vtk 1000 post/sph_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz vx vy vz p p rho
109run           5000
110
111unfix         wall
112
113run           200000
114
115