tsigma-sim1-starter.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tsigma-sim1-starter.py (5430B)
---
1 #!/usr/bin/env python
2 import sphere
3 import numpy
4 import sys
5
6 # launch with:
7 # $ ipython sigma-sim1-starter.py <device> <fluid> <c_phi> <k_c> <sigma_0> <mu> <velfac>
8
9 # start with
10 # ipython sigma-sim1-starter.py 0 1 1.0 2.0e-16 10000.0 2.080e-7 1.0
11
12 device = int(sys.argv[1])
13 wet = int(sys.argv[2])
14 c_phi = float(sys.argv[3])
15 k_c = float(sys.argv[4])
16 sigma0 = float(sys.argv[5])
17 mu = float(sys.argv[6])
18 velfac = float(sys.argv[7])
19
20 if wet == 1:
21 fluid = True
22 else:
23 fluid = False
24
25 #sim = sphere.sim('halfshear-sigma0=' + str(sigma0), fluid=False)
26 sim = sphere.sim('shear-sigma0=' + str(sigma0), fluid=False)
27 sim.readlast()
28 #sim.readbin('../input/shear-sigma0=10000.0-new.bin')
29 #sim.scaleSize(0.01) # from 1 cm to 0.01 cm = 100 micro m (fine sand)
30
31 sim.fluid = fluid
32 if fluid:
33 #sim.id('halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
34 #'-mu=' + str(mu) + '-velfac=' + str(velfac) + '-shear')
35 sim.id('s1-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
36 '-mu=' + str(mu) + '-velfac=' + str(velfac) + '-noflux-shear')
37 else:
38 sim.id('s1-darcy-sigma0=' + str(sigma0) + '-velfac=' + str(velfac) + \
39 '-noflux-shear')
40
41 #sim.checkerboardColors(nx=6,ny=3,nz=6)
42 sim.checkerboardColors(nx=6,ny=6,nz=6)
43 sim.cleanup()
44 sim.adjustUpperWall()
45 sim.zeroKinematics()
46
47 # customized shear function for linear velocity gradient
48 def shearVelocityGradient(sim, shear_strain_rate = 1.0, shear_stress = False):
49 '''
50 Setup shear experiment either by a constant shear rate or a constant
51 shear stress. The shear strain rate is the shear velocity divided by
52 the initial height per second. The shear movement is along the positive
53 x axis. The function zeroes the tangential wall viscosity (gamma_wt) and
54 the wall friction coefficients (mu_ws, mu_wn).
55
56 :param shear_strain_rate: The shear strain rate [-] to use if
57 shear_stress isn't False.
58 :type shear_strain_rate: float
59 :param shear_stress: The shear stress value to use [Pa].
60 :type shear_stress: float or bool
61 '''
62
63 sim.nw = 1
64
65 # Find lowest and heighest point
66 z_min = numpy.min(sim.x[:,2] - sim.radius)
67 z_max = numpy.max(sim.x[:,2] + sim.radius)
68
69 # the grid cell size is equal to the max. particle diameter
70 cellsize = sim.L[0] / sim.num[0]
71
72 # make grid one cell heigher to allow dilation
73 sim.num[2] += 1
74 sim.L[2] = sim.num[2] * cellsize
75
76 # zero kinematics
77 sim.zeroKinematics()
78
79 # Adjust grid and placement of upper wall
80 sim.wmode = numpy.array([1])
81
82 # Fix horizontal velocity to 0.0 of lowermost particles
83 d_max_below = numpy.max(sim.radius[numpy.nonzero(sim.x[:,2] <
84 (z_max-z_min)*0.3)])*2.0
85 I = numpy.nonzero(sim.x[:,2] < (z_min + d_max_below))
86 sim.fixvel[I] = 1
87 sim.angvel[I,0] = 0.0
88 sim.angvel[I,1] = 0.0
89 sim.angvel[I,2] = 0.0
90 sim.vel[I,0] = 0.0 # x-dim
91 sim.vel[I,1] = 0.0 # y-dim
92 sim.color[I] = -1
93
94 # Fix horizontal velocity to specific value of uppermost particles
95 d_max_top = numpy.max(sim.radius[numpy.nonzero(sim.x[:,2] >
96 (z_max-z_min)*0.7)])*2.0
97 I = numpy.nonzero(sim.x[:,2] > (z_max - d_max_top))
98 sim.fixvel[I] = 1
99 sim.angvel[I,0] = 0.0
100 sim.angvel[I,1] = 0.0
101 sim.angvel[I,2] = 0.0
102 if shear_stress == False:
103 prefactor = sim.x[I,1]/sim.L[1]
104 sim.vel[I,0] = prefactor*(z_max-z_min)*shear_strain_rate
105 else:
106 sim.vel[I,0] = 0.0
107 sim.wmode[0] = 3
108 sim.w_tau_x[0] = float(shear_stress)
109 sim.vel[I,1] = 0.0 # y-dim
110 sim.color[I] = -1
111
112 # Set wall tangential viscosity to zero
113 sim.gamma_wt[0] = 0.0
114
115 # Set wall friction coefficients to zero
116 sim.mu_ws[0] = 0.0
117 sim.mu_wd[0] = 0.0
118 return sim
119
120 sim = shearVelocityGradient(sim, 1.0/20.0 * velfac)
121 K_q_real = 36.4e9
122 K_w_real = 2.2e9
123
124
125 #K_q_sim = 1.16e9
126 K_q_sim = 1.16e6
127 sim.setStiffnessNormal(K_q_sim)
128 sim.setStiffnessTangential(K_q_sim)
129 K_w_sim = K_w_real/K_q_real * K_q_sim
130
131
132 if fluid:
133 #sim.num[2] *= 2
134 sim.num[:] /= 2
135 #sim.L[2] *= 2.0
136 #sim.initFluid(mu = 1.787e-6, p = 600.0e3, cfd_solver = 1)
137 sim.initFluid(mu = mu, p = 0.0, cfd_solver = 1)
138 sim.setFluidTopFixedPressure()
139 #sim.setFluidTopFixedFlow()
140 sim.setFluidBottomNoFlow()
141 #sim.setFluidBottomFixedPressure()
142 #sim.setDEMstepsPerCFDstep(10)
143 sim.setMaxIterations(2e5)
144 sim.setPermeabilityPrefactor(k_c)
145 sim.setFluidCompressibility(1.0/K_w_sim)
146
147
148 # frictionless side boundaries
149 sim.periodicBoundariesX()
150
151 # rearrange particle assemblage to accomodate frictionless side boundaries
152 sim.x[:,1] += numpy.abs(numpy.min(sim.x[:,1] - sim.radius[:]))
153 sim.L[1] = numpy.max(sim.x[:,1] + sim.radius[:])
154
155
156 sim.w_sigma0[0] = sigma0
157 sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
158
159 #sim.setStiffnessNormal(36.4e9 * 0.1 / 2.0)
160 #sim.setStiffnessTangential(36.4e9/3.0 * 0.1 / 2.0)
161 sim.setStiffnessNormal(K_q_sim)
162 sim.setStiffnessTangential(K_q_sim)
163 sim.mu_s[0] = 0.5
164 sim.mu_d[0] = 0.5
165 sim.setDampingNormal(0.0)
166 sim.setDampingTangential(0.0)
167 #sim.deleteAllParticles()
168 #sim.fixvel[:] = -1.0
169
170 sim.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
171 #sim.time_dt[0] *= 1.0e-2
172 #sim.initTemporal(total = 1.0e-4, file_dt = 1.0e-5, epsilon=0.07)
173
174 # Fix lowermost particles
175 #dz = sim.L[2]/sim.num[2]
176 #I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
177 #sim.fixvel[I] = 1
178
179 sim.run(dry=True)
180
181 sim.run(device=device)
182 sim.writeVTKall()