analysis.py - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
analysis.py (23402B)
---
1 import math
2 import subprocess
3 import numpy
4 from .common import V_sphere
5
6
7 class SimAnalysis:
8 'Analysis of energies, stresses, contacts and porosity.'
9
10 def currentNormalStress(self, type='defined'):
11 '''
12 Calculates the current magnitude of the defined or effective top wall
13 normal stress.
14
15 :param type: Find the 'defined' (default) or 'effective' normal stress
16 :type type: str
17
18 :returns: The current top wall normal stress in Pascal
19 :return type: float
20 '''
21 if type == 'defined':
22 return self.w_sigma0[0] \
23 + self.w_sigma0_A[0] \
24 *numpy.sin(2.0*numpy.pi*self.w_sigma0_f[0]\
25 *self.time_current[0])
26 elif type == 'effective':
27 return self.w_force[0]/(self.L[0]*self.L[1])
28 else:
29 raise Exception('Normal stress type ' + type + ' not understood')
30
31 def surfaceArea(self, idx):
32 '''
33 Returns the surface area of a particle.
34
35 :param idx: Particle index
36 :type idx: int
37 :returns: The surface area of the particle [m^2]
38 :return type: float
39 '''
40 return 4.0*numpy.pi*self.radius[idx]**2
41
42 def volume(self, idx):
43 '''
44 Returns the volume of a particle.
45
46 :param idx: Particle index
47 :type idx: int
48 :returns: The volume of the particle [m^3]
49 :return type: float
50 '''
51 return V_sphere(self.radius[idx])
52
53 def mass(self, idx):
54 '''
55 Returns the mass of a particle.
56
57 :param idx: Particle index
58 :type idx: int
59 :returns: The mass of the particle [kg]
60 :return type: float
61 '''
62 return self.rho[0]*self.volume(idx)
63
64 def totalMass(self):
65 '''
66 Returns the total mass of all particles.
67
68 :returns: The total mass in [kg]
69 '''
70 m = 0.0
71 for i in range(self.np):
72 m += self.mass(i)
73 return m
74
75 def smallestMass(self):
76 '''
77 Returns the mass of the leightest particle.
78
79 :param idx: Particle index
80 :type idx: int
81 :returns: The mass of the particle [kg]
82 :return type: float
83 '''
84 return V_sphere(numpy.min(self.radius))
85
86 def largestMass(self):
87 '''
88 Returns the mass of the heaviest particle.
89
90 :param idx: Particle index
91 :type idx: int
92 :returns: The mass of the particle [kg]
93 :return type: float
94 '''
95 return V_sphere(numpy.max(self.radius))
96
97 def momentOfInertia(self, idx):
98 '''
99 Returns the moment of inertia of a particle.
100
101 :param idx: Particle index
102 :type idx: int
103 :returns: The moment of inertia [kg*m^2]
104 :return type: float
105 '''
106 return 2.0/5.0*self.mass(idx)*self.radius[idx]**2
107
108 def kineticEnergy(self, idx):
109 '''
110 Returns the (linear) kinetic energy for a particle.
111
112 :param idx: Particle index
113 :type idx: int
114 :returns: The kinetic energy of the particle [J]
115 :return type: float
116 '''
117 return 0.5*self.mass(idx) \
118 *numpy.sqrt(numpy.dot(self.vel[idx, :], self.vel[idx, :]))**2
119
120 def totalKineticEnergy(self):
121 '''
122 Returns the total linear kinetic energy for all particles.
123
124 :returns: The kinetic energy of all particles [J]
125 '''
126 esum = 0.0
127 for i in range(self.np):
128 esum += self.kineticEnergy(i)
129 return esum
130
131 def rotationalEnergy(self, idx):
132 '''
133 Returns the rotational energy for a particle.
134
135 :param idx: Particle index
136 :type idx: int
137 :returns: The rotational kinetic energy of the particle [J]
138 :return type: float
139 '''
140 return 0.5*self.momentOfInertia(idx) \
141 *numpy.sqrt(numpy.dot(self.angvel[idx, :], self.angvel[idx, :]))**2
142
143 def totalRotationalEnergy(self):
144 '''
145 Returns the total rotational kinetic energy for all particles.
146
147 :returns: The rotational energy of all particles [J]
148 '''
149 esum = 0.0
150 for i in range(self.np):
151 esum += self.rotationalEnergy(i)
152 return esum
153
154 def viscousEnergy(self, idx):
155 '''
156 Returns the viscous dissipated energy for a particle.
157
158 :param idx: Particle index
159 :type idx: int
160 :returns: The energy lost by the particle by viscous dissipation [J]
161 :return type: float
162 '''
163 return self.ev[idx]
164
165 def totalViscousEnergy(self):
166 '''
167 Returns the total viscous dissipated energy for all particles.
168
169 :returns: The normal viscous energy lost by all particles [J]
170 :return type: float
171 '''
172 esum = 0.0
173 for i in range(self.np):
174 esum += self.viscousEnergy(i)
175 return esum
176
177 def frictionalEnergy(self, idx):
178 '''
179 Returns the frictional dissipated energy for a particle.
180
181 :param idx: Particle index
182 :type idx: int
183 :returns: The frictional energy lost of the particle [J]
184 :return type: float
185 '''
186 return self.es[idx]
187
188 def totalFrictionalEnergy(self):
189 '''
190 Returns the total frictional dissipated energy for all particles.
191
192 :returns: The total frictional energy lost of all particles [J]
193 :return type: float
194 '''
195 esum = 0.0
196 for i in range(self.np):
197 esum += self.frictionalEnergy(i)
198 return esum
199
200 def energy(self, method):
201 '''
202 Calculates the sum of the energy components of all particles.
203
204 :param method: The type of energy to return. Possible values are 'pot'
205 for potential energy [J], 'kin' for kinetic energy [J], 'rot' for
206 rotational energy [J], 'shear' for energy lost by friction,
207 'shearrate' for the rate of frictional energy loss [W], 'visc_n' for
208 viscous losses normal to the contact [J], 'visc_n_rate' for the rate
209 of viscous losses normal to the contact [W], and finally 'bondpot'
210 for the potential energy stored in bonds [J]
211 :type method: str
212 :returns: The value of the selected energy type
213 :return type: float
214 '''
215
216 if method == 'pot':
217 m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
218 return numpy.sum(m*math.sqrt(numpy.dot(self.g, self.g))*self.x[:, 2])
219
220 elif method == 'kin':
221 m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
222 esum = 0.0
223 for i in range(self.np):
224 esum += 0.5*m[i]*math.sqrt(\
225 numpy.dot(self.vel[i, :], self.vel[i, :]))**2
226 return esum
227
228 elif method == 'rot':
229 m = numpy.ones(self.np)*4.0/3.0*math.pi*self.radius**3*self.rho
230 esum = 0.0
231 for i in range(self.np):
232 esum += 0.5*2.0/5.0*m[i]*self.radius[i]**2 \
233 *math.sqrt(\
234 numpy.dot(self.angvel[i, :], self.angvel[i, :]))**2
235 return esum
236
237 elif method == 'shear':
238 return numpy.sum(self.es)
239
240 elif method == 'shearrate':
241 return numpy.sum(self.es_dot)
242
243 elif method == 'visc_n':
244 return numpy.sum(self.ev)
245
246 elif method == 'visc_n_rate':
247 return numpy.sum(self.ev_dot)
248
249 elif method == 'bondpot':
250 if self.nb0 > 0:
251 R_bar = self.lambda_bar*numpy.minimum(\
252 self.radius[self.bonds[:, 0]],\
253 self.radius[self.bonds[:, 1]])
254 A = numpy.pi*R_bar**2
255 I = 0.25*numpy.pi*R_bar**4
256 J = I*2.0
257 bondpot_fn = numpy.sum(\
258 0.5*A*self.k_n*numpy.abs(self.bonds_delta_n)**2)
259 bondpot_ft = numpy.sum(\
260 0.5*A*self.k_t*numpy.linalg.norm(self.bonds_delta_t)**2)
261 bondpot_tn = numpy.sum(\
262 0.5*J*self.k_t*numpy.abs(self.bonds_omega_n)**2)
263 bondpot_tt = numpy.sum(\
264 0.5*I*self.k_n*numpy.linalg.norm(self.bonds_omega_t)**2)
265 return bondpot_fn + bondpot_ft + bondpot_tn + bondpot_tt
266 else:
267 return 0.0
268 else:
269 raise Exception('Unknownw energy() method "' + method + '"')
270
271 def voidRatio(self):
272 '''
273 Calculates the current void ratio
274
275 :returns: The void ratio (pore volume relative to solid volume)
276 :return type: float
277 '''
278
279 # Find the bulk volume
280 V_t = (self.L[0] - self.origo[0]) \
281 *(self.L[1] - self.origo[1]) \
282 *(self.w_x[0] - self.origo[2])
283
284 # Find the volume of solids
285 V_s = numpy.sum(4.0/3.0 * math.pi * self.radius**3)
286
287 # Return the void ratio
288 e = (V_t - V_s)/V_s
289 return e
290
291 def bulkPorosity(self, trim=True):
292 '''
293 Calculates the bulk porosity of the particle assemblage.
294
295 :param trim: Trim the total volume to the smallest axis-parallel cube
296 containing all particles.
297 :type trim: bool
298
299 :returns: The bulk porosity, in [0:1]
300 :return type: float
301 '''
302
303 V_total = 0.0
304 if trim:
305 min_x = numpy.min(self.x[:, 0] - self.radius)
306 min_y = numpy.min(self.x[:, 1] - self.radius)
307 min_z = numpy.min(self.x[:, 2] - self.radius)
308 max_x = numpy.max(self.x[:, 0] + self.radius)
309 max_y = numpy.max(self.x[:, 1] + self.radius)
310 max_z = numpy.max(self.x[:, 2] + self.radius)
311 V_total = (max_x - min_x)*(max_y - min_y)*(max_z - min_z)
312
313 else:
314 if self.nw == 0:
315 V_total = self.L[0] * self.L[1] * self.L[2]
316 elif self.nw == 1:
317 V_total = self.L[0] * self.L[1] * self.w_x[0]
318 if V_total <= 0.0:
319 raise Exception("Could not determine total volume")
320
321 # Find the volume of solids
322 V_solid = numpy.sum(V_sphere(self.radius))
323 return (V_total - V_solid) / V_total
324
325 def porosity(self, slices=10, verbose=False):
326 '''
327 Calculates the porosity as a function of depth, by averaging values in
328 horizontal slabs. Returns porosity values and their corresponding depth.
329 The values are calculated using the external ``porosity`` program.
330
331 :param slices: The number of vertical slabs to find porosities in.
332 :type slices: int
333 :param verbose: Show the file name of the temporary file written to
334 disk
335 :type verbose: bool
336 :returns: A 2d array of depths and their averaged porosities
337 :return type: numpy.array
338 '''
339
340 # Write data as binary
341 self.writebin(verbose=False)
342
343 # Run porosity program on binary
344 pipe = subprocess.Popen(["../porosity",\
345 "-s", "{}".format(slices),
346 "../input/" + self.sid + ".bin"],
347 stdout=subprocess.PIPE)
348 output, err = pipe.communicate()
349
350 if err:
351 print(err)
352 raise Exception("Could not run external 'porosity' program")
353
354 # read one line of output at a time
355 s2 = output.split(b'\n')
356 depth = []
357 porosity = []
358 for row in s2:
359 if row != '\n' or row != '' or row != ' ': # skip blank lines
360 s3 = row.split(b'\t')
361 if s3 != '' and len(s3) == 2: # make sure line has two vals
362 depth.append(float(s3[0]))
363 porosity.append(float(s3[1]))
364
365 return numpy.array(porosity), numpy.array(depth)
366
367 def shearDisplacement(self):
368 '''
369 Calculates and returns the current shear displacement. The displacement
370 is found by determining the total x-axis displacement of the upper,
371 fixed particles.
372
373 :returns: The total shear displacement [m]
374 :return type: float
375
376 See also: :func:`shearStrain()` and :func:`shearVelocity()`
377 '''
378
379 # Displacement of the upper, fixed particles in the shear direction
380 #xdisp = self.time_current[0] * self.shearVel()
381 fixvel = numpy.nonzero(self.fixvel > 0.0)
382 return numpy.max(self.xyzsum[fixvel, 0])
383
384 def shearVelocity(self):
385 '''
386 Calculates and returns the current shear velocity. The displacement
387 is found by determining the total x-axis velocity of the upper,
388 fixed particles.
389
390 :returns: The shear velocity [m/s]
391 :return type: float
392
393 See also: :func:`shearStrainRate()` and :func:`shearDisplacement()`
394 '''
395 # Displacement of the upper, fixed particles in the shear direction
396 #xdisp = self.time_current[0] * self.shearVel()
397 fixvel = numpy.nonzero(self.fixvel > 0.0)
398 return numpy.max(self.vel[fixvel, 0])
399
400 def shearVel(self):
401 '''
402 Alias of :func:`shearVelocity()`
403 '''
404 return self.shearVelocity()
405
406 def shearStrain(self):
407 '''
408 Calculates and returns the current shear strain (gamma) value of the
409 experiment. The shear strain is found by determining the total x-axis
410 displacement of the upper, fixed particles.
411
412 :returns: The total shear strain [-]
413 :return type: float
414
415 See also: :func:`shearStrainRate()` and :func:`shearVel()`
416 '''
417
418 # Current height
419 w_x0 = self.w_x[0]
420
421 # Displacement of the upper, fixed particles in the shear direction
422 xdisp = self.shearDisplacement()
423
424 # Return shear strain
425 return xdisp/w_x0
426
427 def shearStrainRate(self):
428 '''
429 Calculates the shear strain rate (dot(gamma)) value of the experiment.
430
431 :returns: The value of dot(gamma)
432 :return type: float
433
434 See also: :func:`shearStrain()` and :func:`shearVel()`
435 '''
436 #return self.shearStrain()/self.time_current[1]
437
438 # Current height
439 w_x0 = self.w_x[0]
440 v = self.shearVelocity()
441
442 # Return shear strain rate
443 return v/w_x0
444
445 def inertiaParameterPlanarShear(self):
446 '''
447 Returns the value of the inertia parameter $I$ during planar shear
448 proposed by GDR-MiDi 2004.
449
450 :returns: The value of $I$
451 :return type: float
452
453 See also: :func:`shearStrainRate()` and :func:`shearVel()`
454 '''
455 return self.shearStrainRate() * numpy.mean(self.radius) \
456 * numpy.sqrt(self.rho[0]/self.currentNormalStress())
457
458 def findOverlaps(self):
459 '''
460 Find all particle-particle overlaps by a n^2 contact search, which is
461 done in C++. The particle pair indexes and the distance of the overlaps
462 is saved in the object itself as the ``.pairs`` and ``.overlaps``
463 members.
464
465 See also: :func:`findNormalForces()`
466 '''
467 self.writebin(verbose=False)
468 subprocess.call('cd .. && ./sphere --contacts input/' + self.sid
469 + '.bin > output/' + self.sid + '.contacts.txt',
470 shell=True)
471 contactdata = numpy.loadtxt('../output/' + self.sid + '.contacts.txt')
472 self.pairs = numpy.array((contactdata[:, 0], contactdata[:, 1]),
473 dtype=numpy.int32)
474 self.overlaps = numpy.array(contactdata[:, 2])
475
476 def findCoordinationNumber(self):
477 '''
478 Finds the coordination number (the average number of contacts per
479 particle). Requires a previous call to :func:`findOverlaps()`. Values
480 are stored in ``self.coordinationnumber``.
481 '''
482 self.coordinationnumber = numpy.zeros(self.np, dtype=int)
483 for i in numpy.arange(self.overlaps.size):
484 self.coordinationnumber[self.pairs[0, i]] += 1
485 self.coordinationnumber[self.pairs[1, i]] += 1
486
487 def findMeanCoordinationNumber(self):
488 '''
489 Returns the coordination number (the average number of contacts per
490 particle). Requires a previous call to :func:`findOverlaps()`
491
492 :returns: The mean particle coordination number
493 :return type: float
494 '''
495 return numpy.mean(self.coordinationnumber)
496
497 def findNormalForces(self):
498 '''
499 Finds all particle-particle overlaps (by first calling
500 :func:`findOverlaps()`) and calculating the normal magnitude by
501 multiplying the overlaps with the elastic stiffness ``self.k_n``.
502
503 The result is saved in ``self.f_n_magn``.
504
505 See also: :func:`findOverlaps()` and :func:`findContactStresses()`
506 '''
507 self.findOverlaps()
508 self.f_n_magn = self.k_n * numpy.abs(self.overlaps)
509
510 def contactSurfaceArea(self, i, j, overlap):
511 '''
512 Finds the contact surface area of an inter-particle contact.
513
514 :param i: Index of first particle
515 :type i: int or array of ints
516 :param j: Index of second particle
517 :type j: int or array of ints
518 :param d: Overlap distance
519 :type d: float or array of floats
520 :returns: Contact area [m*m]
521 :return type: float or array of floats
522 '''
523 r_i = self.radius[i]
524 r_j = self.radius[j]
525 d = r_i + r_j + overlap
526 contact_radius = 1./(2.*d)*((-d + r_i - r_j)*(-d - r_i + r_j)*
527 (-d + r_i + r_j)*(d + r_i + r_j)
528 )**0.5
529 return numpy.pi*contact_radius**2.
530
531 def contactParticleArea(self, i, j):
532 '''
533 Finds the average area of an two particles in an inter-particle contact.
534
535 :param i: Index of first particle
536 :type i: int or array of ints
537 :param j: Index of second particle
538 :type j: int or array of ints
539 :param d: Overlap distance
540 :type d: float or array of floats
541 :returns: Contact area [m*m]
542 :return type: float or array of floats
543 '''
544 r_bar = (self.radius[i] + self.radius[j])*0.5
545 return numpy.pi*r_bar**2.
546
547 def findAllContactSurfaceAreas(self):
548 '''
549 Finds the contact surface area of an inter-particle contact. This
550 function requires a prior call to :func:`findOverlaps()` as it reads
551 from the ``self.pairs`` and ``self.overlaps`` arrays.
552
553 :returns: Array of contact surface areas
554 :return type: array of floats
555 '''
556 return self.contactSurfaceArea(self.pairs[0, :], self.pairs[1, :],
557 self.overlaps)
558
559 def findAllAverageParticlePairAreas(self):
560 '''
561 Finds the average area of an inter-particle contact. This
562 function requires a prior call to :func:`findOverlaps()` as it reads
563 from the ``self.pairs`` and ``self.overlaps`` arrays.
564
565 :returns: Array of contact surface areas
566 :return type: array of floats
567 '''
568 return self.contactParticleArea(self.pairs[0, :], self.pairs[1, :])
569
570 def findContactStresses(self, area='average'):
571 '''
572 Finds all particle-particle uniaxial normal stresses (by first calling
573 :func:`findNormalForces()`) and calculating the stress magnitudes by
574 dividing the normal force magnitude with the average particle area
575 ('average') or by the contact surface area ('contact').
576
577 The result is saved in ``self.sigma_contacts``.
578
579 :param area: Area to use: 'average' (default) or 'contact'
580 :type area: str
581
582 See also: :func:`findNormalForces()` and :func:`findOverlaps()`
583 '''
584 self.findNormalForces()
585 if area == 'average':
586 areas = self.findAllAverageParticlePairAreas()
587 elif area == 'contact':
588 areas = self.findAllContactSurfaceAreas()
589 else:
590 raise Exception('Contact area type "' + area + '" not understood')
591
592 self.sigma_contacts = self.f_n_magn/areas
593
594 def findLoadedContacts(self, threshold):
595 '''
596 Finds the indices of contact pairs where the contact stress magnitude
597 exceeds or is equal to a specified threshold value. This function calls
598 :func:`findContactStresses()`.
599
600 :param threshold: Threshold contact stress [Pa]
601 :type threshold: float
602 :returns: Array of contact indices
603 :return type: array of ints
604 '''
605 self.findContactStresses()
606 return numpy.nonzero(self.sigma_contacts >= threshold)
607
608 def momentum(self, idx):
609 '''
610 Returns the momentum (m*v) of a particle.
611
612 :param idx: The particle index
613 :type idx: int
614 :returns: The particle momentum [N*s]
615 :return type: numpy.array
616 '''
617 return self.rho*V_sphere(self.radius[idx])*self.vel[idx, :]
618
619 def totalMomentum(self):
620 '''
621 Returns the sum of particle momentums.
622
623 :returns: The sum of particle momentums (m*v) [N*s]
624 :return type: numpy.array
625 '''
626 m_sum = numpy.zeros(3)
627 for i in range(self.np):
628 m_sum += self.momentum(i)
629 return m_sum
630
631 def ReynoldsNumber(self):
632 '''
633 Estimate the per-cell Reynolds number by: Re=rho * ||v_f|| * dx/mu.
634 This value is returned and also stored in `self.Re`.
635
636 :returns: Reynolds number
637 :return type: Numpy array with dimensions like the fluid grid
638 '''
639
640 # find magnitude of fluid velocity vectors
641 self.v_f_magn = numpy.empty_like(self.p_f)
642 for z in numpy.arange(self.num[2]):
643 for y in numpy.arange(self.num[1]):
644 for x in numpy.arange(self.num[0]):
645 self.v_f_magn[x, y, z] = \
646 self.v_f[x, y, z, :].dot(self.v_f[x, y, z, :])
647
648 Re = self.rho_f*self.v_f_magn*self.L[0]/self.num[0]/(self.mu + \
649 1.0e-16)
650 return Re
651
652 def convergence(self):
653 '''
654 Read the convergence evolution in the CFD solver. The values are stored
655 in `self.conv` with iteration number in the first column and iteration
656 count in the second column.
657
658 See also: :func:`plotConvergence()`
659 '''
660 return numpy.loadtxt('../output/' + self.sid + '-conv.log', dtype=numpy.int32)
661
662 def shearStress(self, type='effective'):
663 '''
664 Calculates the sum of shear stress values measured on any moving
665 particles with a finite and fixed velocity.
666
667 :param type: Find the 'defined' or 'effective' (default) shear stress
668 :type type: str
669
670 :returns: The shear stress in Pa
671 :return type: numpy.array
672 '''
673
674 if type == 'defined':
675 return self.w_tau_x[0]
676
677 elif type == 'effective':
678
679 fixvel = numpy.nonzero(self.fixvel > 0.0)
680 force = numpy.zeros(3)
681
682 # Summation of shear stress contributions
683 for i in fixvel[0]:
684 if self.vel[i, 0] > 0.0:
685 force += -self.force[i, :]
686
687 return force[0]/(self.L[0]*self.L[1])
688
689 else:
690 raise Exception('Shear stress type ' + type + ' not understood')