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