fluid.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
---
fluid.py (29940B)
---
1 import numpy
2 from .common import V_sphere, plt, py_mpl
3
4
5 class SimFluid:
6 'Fluid phase setup and CFD solver parameters for sim objects.'
7
8 def largestFluidTimeStep(self, safety=0.5, v_max=-1.0):
9 '''
10 Finds and returns the largest time step in the fluid phase by von
11 Neumann and Courant-Friedrichs-Lewy analysis given the current
12 velocities. This ensures stability in the diffusive and advective parts
13 of the momentum equation.
14
15 The value of the time step decreases with increasing fluid viscosity
16 (`self.mu`), and increases with fluid cell size (`self.L/self.num`)
17 and fluid velocities (`self.v_f`).
18
19 NOTE: The fluid time step with the Darcy solver is an arbitrarily
20 large value. In practice, this is not a problem since the short
21 DEM time step is stable for fluid computations.
22
23 :param safety: Safety factor which is multiplied to the largest time
24 step.
25 :type safety: float
26 :param v_max: The largest anticipated absolute fluid velocity [m/s]
27 :type v_max: float
28
29 :returns: The largest timestep stable for the current fluid state.
30 :return type: float
31 '''
32
33 if self.fluid:
34
35 # Normalized velocities
36 v_norm = numpy.empty(self.num[0]*self.num[1]*self.num[2])
37 idx = 0
38 for x in numpy.arange(self.num[0]):
39 for y in numpy.arange(self.num[1]):
40 for z in numpy.arange(self.num[2]):
41 v_norm[idx] = numpy.sqrt(self.v_f[x, y, z, :]\
42 .dot(self.v_f[x, y, z, :]))
43 idx += 1
44
45 v_max_obs = numpy.amax(v_norm)
46 if v_max_obs == 0:
47 v_max_obs = 1.0e-7
48 if v_max < 0.0:
49 v_max = v_max_obs
50
51 dx_min = numpy.min(self.L/self.num)
52 dt_min_cfl = dx_min/v_max
53
54 # Navier-Stokes
55 if self.cfd_solver[0] == 0:
56 dt_min_von_neumann = 0.5*dx_min**2/(self.mu[0] + 1.0e-16)
57
58 return numpy.min([dt_min_von_neumann, dt_min_cfl])*safety
59
60 # Darcy
61 elif self.cfd_solver[0] == 1:
62
63 return dt_min_cfl
64
65 '''
66 # Determine on the base of the diffusivity coefficient
67 # components
68 #self.hydraulicPermeability()
69 #alpha_max = numpy.max(self.k/(self.beta_f*0.9*self.mu))
70 k_max = 2.7e-10 # hardcoded in darcy.cuh
71 phi_min = 0.1 # hardcoded in darcy.cuh
72 alpha_max = k_max/(self.beta_f*phi_min*self.mu)
73 print(alpha_max)
74 return safety * 1.0/(2.0*alpha_max)*1.0/(
75 1.0/(self.dx[0]**2) + \
76 1.0/(self.dx[1]**2) + \
77 1.0/(self.dx[2]**2))
78 '''
79
80 '''
81 # Determine value on the base of the hydraulic conductivity
82 g = numpy.max(numpy.abs(self.g))
83
84 # Bulk modulus of fluid
85 K = 1.0/self.beta_f[0]
86
87 self.hydraulicDiffusivity()
88
89 return safety * 1.0/(2.0*self.D)*1.0/( \
90 1.0/(self.dx[0]**2) + \
91 1.0/(self.dx[1]**2) + \
92 1.0/(self.dx[2]**2))
93 '''
94
95 def hydraulicConductivity(self, phi=0.35):
96 '''
97 Determine the hydraulic conductivity (K) [m/s] from the permeability
98 prefactor and a chosen porosity. This value is stored in `self.K_c`.
99 This function only works for the Darcy solver (`self.cfd_solver == 1`)
100
101 :param phi: The porosity to use in the Kozeny-Carman relationship
102 :type phi: float
103 :returns: The hydraulic conductivity [m/s]
104 :return type: float
105 '''
106 if self.cfd_solver[0] == 1:
107 k = self.k_c * phi**3/(1.0 - phi**2)
108 self.K_c = k*self.rho_f*numpy.abs(self.g[2])/self.mu
109 return self.K_c[0]
110 else:
111 raise Exception('This function only works for the Darcy solver')
112
113 def hydraulicPermeability(self):
114 '''
115 Determine the hydraulic permeability (k) [m*m] from the Kozeny-Carman
116 relationship, using the permeability prefactor (`self.k_c`), and the
117 range of valid porosities set in `src/darcy.cuh`, by default in the
118 range 0.1 to 0.9.
119
120 This function is only valid for the Darcy solver (`self.cfd_solver ==
121 1`).
122 '''
123 if self.cfd_solver[0] == 1:
124 self.findPermeabilities()
125 else:
126 raise Exception('This function only works for the Darcy solver')
127
128 def hydraulicDiffusivity(self):
129 '''
130 Determine the hydraulic diffusivity (D) [m*m/s]. The result is stored in
131 `self.D`. This function only works for the Darcy solver
132 (`self.cfd_solver[0] == 1`)
133 '''
134 if self.cfd_solver[0] == 1:
135 self.hydraulicConductivity()
136 phi_bar = numpy.mean(self.phi)
137 self.D = self.K_c/(self.rho_f*self.g[2]
138 *(self.k_n[0] + phi_bar*self.K))
139 else:
140 raise Exception('This function only works for the Darcy solver')
141
142 def dry(self):
143 '''
144 Set the simulation to be dry (no fluids).
145
146 See also :func:`wet()`
147 '''
148 self.fluid = False
149
150 def wet(self):
151 '''
152 Set the simulation to be wet (total fluid saturation).
153
154 See also :func:`dry()`
155 '''
156 self.fluid = True
157 self.initFluid()
158
159 def initFluid(self, mu=8.9e-4, rho=1.0e3, p=0.0, hydrostatic=False,
160 cfd_solver=0):
161 '''
162 Initialize the fluid arrays and the fluid viscosity. The default value
163 of ``mu`` equals the dynamic viscosity of water at 25 degrees Celcius.
164 The value for water at 0 degrees Celcius is 17.87e-4 kg/(m*s).
165
166 :param mu: The fluid dynamic viscosity [kg/(m*s)]
167 :type mu: float
168 :param rho: The fluid density [kg/(m^3)]
169 :type rho: float
170 :param p: The hydraulic pressure to initialize the cells to. If the
171 parameter `hydrostatic` is set to `True`, this value will apply to
172 the fluid cells at the top
173 :param hydrostatic: Initialize the fluid pressures to the hydrostatic
174 pressure distribution. A pressure gradient with depth is only
175 created if a gravitational acceleration along :math:`z` previously
176 has been specified
177 :type hydrostatic: bool
178 :param cfd_solver: Solver to use for the computational fluid dynamics.
179 Accepted values: 0 (Navier Stokes, default) and 1 (Darcy).
180 :type cfd_solver: int
181 '''
182 self.fluid = True
183 self.mu = numpy.ones(1, dtype=numpy.float64) * mu
184 self.rho_f = numpy.ones(1, dtype=numpy.float64) * rho
185
186 self.p_f = numpy.ones((self.num[0], self.num[1], self.num[2]),
187 dtype=numpy.float64) * p
188
189 if hydrostatic:
190
191 dz = self.L[2]/self.num[2]
192 # Zero pressure gradient from grid top to top wall, linear pressure
193 # distribution from top wall to grid bottom
194 if self.nw == 1:
195 wall0_iz = int(self.w_x[0]/(self.L[2]/self.num[2]))
196 self.p_f[:, :, wall0_iz:] = p
197
198 for iz in numpy.arange(wall0_iz - 1):
199 z = dz*iz + 0.5*dz
200 depth = self.w_x[0] - z
201 self.p_f[:, :, iz] = p + (depth-dz) * rho * -self.g[2]
202
203 # Linear pressure distribution from grid top to grid bottom
204 else:
205 for iz in numpy.arange(self.num[2] - 1):
206 z = dz*iz + 0.5*dz
207 depth = self.L[2] - z
208 self.p_f[:, :, iz] = p + (depth-dz) * rho * -self.g[2]
209
210
211 self.v_f = numpy.zeros((self.num[0], self.num[1], self.num[2], self.nd),
212 dtype=numpy.float64)
213 self.phi = numpy.ones((self.num[0], self.num[1], self.num[2]),
214 dtype=numpy.float64)
215 self.dphi = numpy.zeros((self.num[0], self.num[1], self.num[2]),
216 dtype=numpy.float64)
217
218 self.p_mod_A = numpy.zeros(1, dtype=numpy.float64) # Amplitude [Pa]
219 self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
220 self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
221
222 self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
223 self.bc_top = numpy.zeros(1, dtype=numpy.int32)
224 self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
225 self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
226 self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
227 self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
228
229 self.p_f_constant = numpy.zeros((self.num[0], self.num[1], self.num[2]),
230 dtype=numpy.int32)
231
232 # Fluid solver type
233 # 0: Navier Stokes (fluid with inertia)
234 # 1: Stokes-Darcy (fluid without inertia)
235 self.cfd_solver = numpy.ones(1)*cfd_solver
236
237 if self.cfd_solver[0] == 0:
238 self.gamma = numpy.array(0.0)
239 self.theta = numpy.array(1.0)
240 self.beta = numpy.array(0.0)
241 self.tolerance = numpy.array(1.0e-3)
242 self.maxiter = numpy.array(1e4)
243 self.ndem = numpy.array(1)
244
245 self.c_phi = numpy.ones(1, dtype=numpy.float64)
246 self.c_v = numpy.ones(1, dtype=numpy.float64)
247 self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
248
249 self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
250 self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
251 self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
252 self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
253
254 elif self.cfd_solver[0] == 1:
255 self.tolerance = numpy.array(1.0e-3)
256 self.maxiter = numpy.array(1e4)
257 self.ndem = numpy.array(1)
258 self.c_phi = numpy.ones(1, dtype=numpy.float64)
259 self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
260 self.beta_f = numpy.ones(1, dtype=numpy.float64)*4.5e-10
261 self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
262 self.k_c = numpy.ones(1, dtype=numpy.float64)*4.6e-10
263
264 self.bc_xn = numpy.ones(1, dtype=numpy.int32)*2
265 self.bc_xp = numpy.ones(1, dtype=numpy.int32)*2
266 self.bc_yn = numpy.ones(1, dtype=numpy.int32)*2
267 self.bc_yp = numpy.ones(1, dtype=numpy.int32)*2
268
269 else:
270 raise Exception('Value of cfd_solver not understood (' + \
271 str(self.cfd_solver[0]) + ')')
272
273 def setFluidBottomNoFlow(self):
274 '''
275 Set the lower boundary of the fluid domain to follow the no-flow
276 (Neumann) boundary condition with free slip parallel to the boundary.
277
278 The default behavior for the boundary is fixed value (Dirichlet), see
279 :func:`setFluidBottomFixedPressure()`.
280 '''
281 self.bc_bot[0] = 1
282
283 def setFluidBottomNoFlowNoSlip(self):
284 '''
285 Set the lower boundary of the fluid domain to follow the no-flow
286 (Neumann) boundary condition with no slip parallel to the boundary.
287
288 The default behavior for the boundary is fixed value (Dirichlet), see
289 :func:`setFluidBottomFixedPressure()`.
290 '''
291 self.bc_bot[0] = 2
292
293 def setFluidBottomFixedPressure(self):
294 '''
295 Set the lower boundary of the fluid domain to follow the fixed pressure
296 value (Dirichlet) boundary condition.
297
298 This is the default behavior for the boundary. See also
299 :func:`setFluidBottomNoFlow()`
300 '''
301 self.bc_bot[0] = 0
302
303 def setFluidBottomFixedFlux(self, specific_flux):
304 '''
305 Define a constant fluid flux normal to the boundary.
306
307 The default behavior for the boundary is fixed value (Dirichlet), see
308 :func:`setFluidBottomFixedPressure()`.
309
310 :param specific_flux: Specific flux values across boundary (positive
311 values upwards), [m/s]
312 '''
313 self.bc_bot[0] = 4
314 self.bc_bot_flux[0] = specific_flux
315
316 def setFluidTopNoFlow(self):
317 '''
318 Set the upper boundary of the fluid domain to follow the no-flow
319 (Neumann) boundary condition with free slip parallel to the boundary.
320
321 The default behavior for the boundary is fixed value (Dirichlet), see
322 :func:`setFluidTopFixedPressure()`.
323 '''
324 self.bc_top[0] = 1
325
326 def setFluidTopNoFlowNoSlip(self):
327 '''
328 Set the upper boundary of the fluid domain to follow the no-flow
329 (Neumann) boundary condition with no slip parallel to the boundary.
330
331 The default behavior for the boundary is fixed value (Dirichlet), see
332 :func:`setFluidTopFixedPressure()`.
333 '''
334 self.bc_top[0] = 2
335
336 def setFluidTopFixedPressure(self):
337 '''
338 Set the upper boundary of the fluid domain to follow the fixed pressure
339 value (Dirichlet) boundary condition.
340
341 This is the default behavior for the boundary. See also
342 :func:`setFluidTopNoFlow()`
343 '''
344 self.bc_top[0] = 0
345
346 def setFluidTopFixedFlux(self, specific_flux):
347 '''
348 Define a constant fluid flux normal to the boundary.
349
350 The default behavior for the boundary is fixed value (Dirichlet), see
351 :func:`setFluidBottomFixedPressure()`.
352
353 :param specific_flux: Specific flux values across boundary (positive
354 values upwards), [m/s]
355 '''
356 self.bc_top[0] = 4
357 self.bc_top_flux[0] = specific_flux
358
359 def setFluidXFixedPressure(self):
360 '''
361 Set the X boundaries of the fluid domain to follow the fixed pressure
362 value (Dirichlet) boundary condition.
363
364 This is not the default behavior for the boundary. See also
365 :func:`setFluidXFixedPressure()`,
366 :func:`setFluidXNoFlow()`, and
367 :func:`setFluidXPeriodic()` (default)
368 '''
369 self.bc_xn[0] = 0
370 self.bc_xp[0] = 0
371
372 def setFluidXNoFlow(self):
373 '''
374 Set the X boundaries of the fluid domain to follow the no-flow
375 (Neumann) boundary condition.
376
377 This is not the default behavior for the boundary. See also
378 :func:`setFluidXFixedPressure()`,
379 :func:`setFluidXNoFlow()`, and
380 :func:`setFluidXPeriodic()` (default)
381 '''
382 self.bc_xn[0] = 1
383 self.bc_xp[0] = 1
384
385 def setFluidXPeriodic(self):
386 '''
387 Set the X boundaries of the fluid domain to follow the periodic
388 (cyclic) boundary condition.
389
390 This is the default behavior for the boundary. See also
391 :func:`setFluidXFixedPressure()` and
392 :func:`setFluidXNoFlow()`
393 '''
394 self.bc_xn[0] = 2
395 self.bc_xp[0] = 2
396
397 def setFluidYFixedPressure(self):
398 '''
399 Set the Y boundaries of the fluid domain to follow the fixed pressure
400 value (Dirichlet) boundary condition.
401
402 This is not the default behavior for the boundary. See also
403 :func:`setFluidYNoFlow()` and
404 :func:`setFluidYPeriodic()` (default)
405 '''
406 self.bc_yn[0] = 0
407 self.bc_yp[0] = 0
408
409 def setFluidYNoFlow(self):
410 '''
411 Set the Y boundaries of the fluid domain to follow the no-flow
412 (Neumann) boundary condition.
413
414 This is not the default behavior for the boundary. See also
415 :func:`setFluidYFixedPressure()` and
416 :func:`setFluidYPeriodic()` (default)
417 '''
418 self.bc_yn[0] = 1
419 self.bc_yp[0] = 1
420
421 def setFluidYPeriodic(self):
422 '''
423 Set the Y boundaries of the fluid domain to follow the periodic
424 (cyclic) boundary condition.
425
426 This is the default behavior for the boundary. See also
427 :func:`setFluidYFixedPressure()` and
428 :func:`setFluidYNoFlow()`
429 '''
430 self.bc_yn[0] = 2
431 self.bc_yp[0] = 2
432
433 def setPermeabilityGrainSize(self, verbose=True):
434 '''
435 Set the permeability prefactor based on the mean grain size (Damsgaard
436 et al., 2015, eq. 10).
437
438 :param verbose: Print information about the realistic permeabilities
439 hydraulic conductivities to expect with the chosen permeability
440 prefactor.
441 :type verbose: bool
442 '''
443 self.setPermeabilityPrefactor(k_c=numpy.mean(self.radius*2.0)**2.0/180.0,
444 verbose=verbose)
445
446 def setPermeabilityPrefactor(self, k_c, verbose=True):
447 '''
448 Set the permeability prefactor from Goren et al 2011, eq. 24. The
449 function will print the limits of permeabilities to be simulated. This
450 parameter is only used in the Darcy solver.
451
452 :param k_c: Permeability prefactor value [m*m]
453 :type k_c: float
454 :param verbose: Print information about the realistic permeabilities and
455 hydraulic conductivities to expect with the chosen permeability
456 prefactor.
457 :type verbose: bool
458 '''
459 if self.cfd_solver[0] == 1:
460 self.k_c[0] = k_c
461 if verbose:
462 phi = numpy.array([0.1, 0.35, 0.9])
463 k = self.k_c * phi**3/(1.0 - phi**2)
464 K = k * self.rho_f*numpy.abs(self.g[2])/self.mu
465 print('Hydraulic permeability limits for porosity phi=' + \
466 str(phi) + ':')
467 print('\tk=' + str(k) + ' m*m')
468 print('Hydraulic conductivity limits for porosity phi=' + \
469 str(phi) + ':')
470 print('\tK=' + str(K) + ' m/s')
471 else:
472 raise Exception('setPermeabilityPrefactor() only relevant for the '
473 'Darcy solver (cfd_solver=1)')
474
475 def findPermeabilities(self):
476 '''
477 Calculates the hydrological permeabilities from the Kozeny-Carman
478 relationship. These values are only relevant when the Darcy solver is
479 used (`self.cfd_solver=1`). The permeability pre-factor `self.k_c`
480 and the assemblage porosities must be set beforehand. The former values
481 are set if a file from the `output/` folder is read using
482 `self.readbin`.
483 '''
484 if self.cfd_solver[0] == 1:
485 phi = numpy.clip(self.phi, 0.1, 0.9)
486 self.k = self.k_c * phi**3/(1.0 - phi**2)
487 else:
488 raise Exception('findPermeabilities() only relevant for the '
489 'Darcy solver (cfd_solver=1)')
490
491 def findHydraulicConductivities(self):
492 '''
493 Calculates the hydrological conductivities from the Kozeny-Carman
494 relationship. These values are only relevant when the Darcy solver is
495 used (`self.cfd_solver=1`). The permeability pre-factor `self.k_c`
496 and the assemblage porosities must be set beforehand. The former values
497 are set if a file from the `output/` folder is read using
498 `self.readbin`.
499 '''
500 if self.cfd_solver[0] == 1:
501 self.findPermeabilities()
502 self.K = self.k*self.rho_f*numpy.abs(self.g[2])/self.mu
503 else:
504 raise Exception('findPermeabilities() only relevant for the '
505 'Darcy solver (cfd_solver=1)')
506
507 def setFluidCompressibility(self, beta_f):
508 '''
509 Set the fluid adiabatic compressibility [1/Pa]. This value is equal to
510 `1/K` where `K` is the bulk modulus [Pa]. The value for water is 5.1e-10
511 for water at 0 degrees Celcius. This parameter is used for the Darcy
512 solver exclusively.
513
514 :param beta_f: The fluid compressibility [1/Pa]
515 :type beta_f: float
516
517 See also: :func:`setFluidDensity()` and :func:`setFluidViscosity()`
518 '''
519 if self.cfd_solver[0] == 1:
520 self.beta_f[0] = beta_f
521 else:
522 raise Exception('setFluidCompressibility() only relevant for the '
523 'Darcy solver (cfd_solver=1)')
524
525 def setFluidViscosity(self, mu):
526 '''
527 Set the fluid dynamic viscosity [Pa*s]. The value for water is
528 1.797e-3 at 0 degrees Celcius. This parameter is used for both the Darcy
529 and Navier-Stokes fluid solver.
530
531 :param mu: The fluid dynamic viscosity [Pa*s]
532 :type mu: float
533
534 See also: :func:`setFluidDensity()` and
535 :func:`setFluidCompressibility()`
536 '''
537 self.mu[0] = mu
538
539 def setFluidDensity(self, rho_f):
540 '''
541 Set the fluid density [kg/(m*m*m)]. The value for water is 1000. This
542 parameter is used for the Navier-Stokes fluid solver exclusively.
543
544 :param rho_f: The fluid density [kg/(m*m*m)]
545 :type rho_f: float
546
547 See also: :func:`setFluidViscosity()` and
548 :func:`setFluidCompressibility()`
549 '''
550 self.rho_f[0] = rho_f
551
552 def setTopWallNormalStressModulation(self, A, f, plot=False):
553 '''
554 Set the parameters for the sine wave modulating the normal stress
555 at the top wall. Note that a cos-wave is obtained with phi=pi/2.
556
557 :param A: Fluctuation amplitude [Pa]
558 :type A: float
559 :param f: Fluctuation frequency [Hz]
560 :type f: float
561 :param plot: Show a plot of the resulting modulation
562 :type plot: bool
563
564 See also: :func:`setFluidPressureModulation()` and
565 :func:`disableTopWallNormalStressModulation()`
566 '''
567 self.w_sigma0_A[0] = A
568 self.w_sigma0_f[0] = f
569
570 if plot and py_mpl:
571 self.plotSinFunction(self.w_sigma0[0], A, f, phi=0.0,
572 xlabel='$t$ [s]', ylabel='$\\sigma_0$ [Pa]')
573
574 def disableTopWallNormalStressModulation(self):
575 '''
576 Set the parameters for the sine wave modulating the normal stress
577 at the top dynamic wall to zero.
578
579 See also: :func:`setTopWallNormalStressModulation()`
580 '''
581 self.setTopWallNormalStressModulation(A=0.0, f=0.0)
582
583 def setFluidPressureModulation(self, A, f, phi=0.0, plot=False):
584 '''
585 Set the parameters for the sine wave modulating the fluid pressures
586 at the top boundary. Note that a cos-wave is obtained with phi=pi/2.
587
588 :param A: Fluctuation amplitude [Pa]
589 :type A: float
590 :param f: Fluctuation frequency [Hz]
591 :type f: float
592 :param phi: Fluctuation phase shift (default=0.0) [rad]
593 :type phi: float
594 :param plot: Show a plot of the resulting modulation
595 :type plot: bool
596
597 See also: :func:`setTopWallNormalStressModulation()` and
598 :func:`disableFluidPressureModulation()`
599 '''
600 self.p_mod_A[0] = A
601 self.p_mod_f[0] = f
602 self.p_mod_phi[0] = phi
603
604 if plot:
605 self.plotSinFunction(self.p_f[0, 0, -1], A, f, phi=0.0,
606 xlabel='$t$ [s]', ylabel='$p_f$ [kPa]')
607
608 def disableFluidPressureModulation(self):
609 '''
610 Set the parameters for the sine wave modulating the fluid pressures
611 at the top boundary to zero.
612
613 See also: :func:`setFluidPressureModulation()`
614 '''
615 self.setFluidPressureModulation(A=0.0, f=0.0)
616
617 def plotPrescribedFluidPressures(self, graphics_format='png',
618 verbose=True):
619 '''
620 Plot the prescribed fluid pressures through time that may be
621 modulated through the class parameters p_mod_A, p_mod_f, and p_mod_phi.
622 The plot is saved in the output folder with the file name
623 '<simulation id>-pres.<graphics_format>'.
624 '''
625 if not py_mpl:
626 print('Error: matplotlib module not found ' +
627 '(plotPrescribedFluidPressures).')
628 return
629
630 fig = plt.figure()
631
632 plt.title('Prescribed fluid pressures at the top in "' + self.sid + '"')
633 plt.xlabel('Time [s]')
634 plt.ylabel('Pressure [Pa]')
635 t = numpy.linspace(0, self.time_total, self.time_total/self.time_file_dt)
636 p = self.p_f[0, 0, -1] + self.p_mod_A * \
637 numpy.sin(2.0*numpy.pi*self.p_mod_f*t + self.p_mod_phi)
638 plt.plot(t, p, '.-')
639 plt.grid()
640 filename = '../output/' + self.sid + '-pres.' + graphics_format
641 plt.savefig(filename)
642 if verbose:
643 print('saved to ' + filename)
644 plt.clf()
645 plt.close(fig)
646
647 def acceleration(self, idx=-1):
648 '''
649 Returns the acceleration of one or more particles, selected by their
650 index. If the index is equal to -1 (default value), all accelerations
651 are returned.
652
653 :param idx: Index or index range of particles
654 :type idx: int, list or numpy.array
655 :returns: n-by-3 matrix of acceleration(s)
656 :return type: numpy.array
657 '''
658 if idx == -1:
659 idx = range(self.np)
660 return self.force[idx, :]/(V_sphere(self.radius[idx])*self.rho[0]) + \
661 self.g
662
663 def setGamma(self, gamma):
664 '''
665 Gamma is a fluid solver parameter, used for smoothing the pressure
666 values. The epsilon (pressure) values are smoothed by including the
667 average epsilon value of the six closest (face) neighbor cells. This
668 parameter should be in the range [0.0;1.0[. The higher the value, the
669 more averaging is introduced. A value of 0.0 disables all averaging.
670
671 The default and recommended value is 0.0.
672
673 :param theta: The smoothing parameter value
674 :type theta: float
675
676 Other solver parameter setting functions: :func:`setTheta()`,
677 :func:`setBeta()`, :func:`setTolerance()`,
678 :func:`setDEMstepsPerCFDstep()` and :func:`setMaxIterations()`
679 '''
680 self.gamma = numpy.asarray(gamma)
681
682 def setTheta(self, theta):
683 '''
684 Theta is a fluid solver under-relaxation parameter, used in solution of
685 Poisson equation. The value should be within the range ]0.0;1.0]. At a
686 value of 1.0, the new estimate of epsilon values is used exclusively. At
687 lower values, a linear interpolation between new and old values is used.
688 The solution typically converges faster with a value of 1.0, but
689 instabilities may be avoided with lower values.
690
691 The default and recommended value is 1.0.
692
693 :param theta: The under-relaxation parameter value
694 :type theta: float
695
696 Other solver parameter setting functions: :func:`setGamma()`,
697 :func:`setBeta()`, :func:`setTolerance()`,
698 :func:`setDEMstepsPerCFDstep()` and :func:`setMaxIterations()`
699 '''
700 self.theta = numpy.asarray(theta)
701
702 def setBeta(self, beta):
703 '''
704 Beta is a fluid solver parameter, used in velocity prediction and
705 pressure iteration 1.0: Use old pressures for fluid velocity prediction
706 (see Langtangen et al. 2002) 0.0: Do not use old pressures for fluid
707 velocity prediction (Chorin's original projection method, see Chorin
708 (1968) and "Projection method (fluid dynamics)" page on Wikipedia. The
709 best results precision and performance-wise are obtained by using a beta
710 of 0 and a low tolerance criteria value.
711
712 The default and recommended value is 0.0.
713
714 Other solver parameter setting functions: :func:`setGamma()`,
715 :func:`setTheta()`, :func:`setTolerance()`,
716 :func:`setDEMstepsPerCFDstep()` and
717 :func:`setMaxIterations()`
718 '''
719 self.beta = numpy.asarray(beta)
720
721 def setTolerance(self, tolerance):
722 '''
723 A fluid solver parameter, the value of the tolerance parameter denotes
724 the required value of the maximum normalized residual for the fluid
725 solver.
726
727 The default and recommended value is 1.0e-3.
728
729 :param tolerance: The tolerance criteria for the maximal normalized
730 residual
731 :type tolerance: float
732
733 Other solver parameter setting functions: :func:`setGamma()`,
734 :func:`setTheta()`, :func:`setBeta()`, :func:`setDEMstepsPerCFDstep()` and
735 :func:`setMaxIterations()`
736 '''
737 self.tolerance = numpy.asarray(tolerance)
738
739 def setMaxIterations(self, maxiter):
740 '''
741 A fluid solver parameter, the value of the maxiter parameter denotes the
742 maximal allowed number of fluid solver iterations before ending the
743 fluid solver loop prematurely. The residual values are at that point not
744 fulfilling the tolerance criteria. The parameter is included to avoid
745 infinite hangs.
746
747 The default and recommended value is 1e4.
748
749 :param maxiter: The maximum number of Jacobi iterations in the fluid
750 solver
751 :type maxiter: int
752
753 Other solver parameter setting functions: :func:`setGamma()`,
754 :func:`setTheta()`, :func:`setBeta()`, :func:`setDEMstepsPerCFDstep()`
755 and :func:`setTolerance()`
756 '''
757 self.maxiter = numpy.asarray(maxiter)
758
759 def setDEMstepsPerCFDstep(self, ndem):
760 '''
761 A fluid solver parameter, the value of the maxiter parameter denotes the
762 number of DEM time steps to be performed per CFD time step.
763
764 The default value is 1.
765
766 :param ndem: The DEM/CFD time step ratio
767 :type ndem: int
768
769 Other solver parameter setting functions: :func:`setGamma()`,
770 :func:`setTheta()`, :func:`setBeta()`, :func:`setTolerance()` and
771 :func:`setMaxIterations()`.
772 '''
773 self.ndem = numpy.asarray(ndem)