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