MFIX  2016-1
check_solids_dem.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2 ! !
3 ! SUBROUTINE: CHECK_DES_SOLIDS !
4 ! Author: J.Musser Date: 02-FEB-14 !
5 ! !
6 ! Purpose: !
7 ! !
8 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
9  SUBROUTINE check_solids_dem
10 
11 ! Use the error manager for posting error messages.
12 !---------------------------------------------------------------------//
13  use error_manager
14 
15  implicit none
16 
17 ! Initialize the error manager.
18  CALL init_err_msg("CHECK_SOLIDS_DEM")
19 
20 ! Particle-particle collision parameters.
22 ! DES cohesion model parameters.
24 ! Particle-particle conduction model parameters.
26 
27  CALL finl_err_msg
28 
29  RETURN
30 
31  END SUBROUTINE check_solids_dem
32 
33 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
34 ! !
35 ! SUBROUTINE CHECK_SOLIDS_DEM_ENERGY !
36 ! Author: J.Musser Date: 02-FEB-14 !
37 ! !
38 ! Purpose: !
39 ! !
40 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
41  SUBROUTINE check_solids_dem_energy
42 
43 ! Global Variables
44 !---------------------------------------------------------------------//
45 
46  use des_thermo, only: des_min_cond_dist
47  use run, only: units
48 
49 ! Global Parameters:
50 !---------------------------------------------------------------------//
51  use param1, only: undefined
52  use error_manager
53  use des_thermo, only : calc_cond_des
55  use discretelement
56  use param1, only : one
57  use physprop, only : mmax, ro_s0, d_p0
58  use run, only : solids_model
59  use constant, only : pi
60 
61  IMPLICIT NONE
62 !......................................................................!
63  INTEGER :: M, L ! Loop indices for DEM solids
64 ! Calculate phase index offset for certain inputs until it can be
65 ! addressed in other ways. should not matter unlesss hybrid
66  CHARACTER(len=64) :: MSG
67  INTEGER :: lent, lend, lenc, LC
68  INTEGER :: MMAX_TOT
69  LOGICAL :: ANY_CONDUCTION = .false.
70 
71  DOUBLE PRECISION :: MASS_L, MASS_M, MASS_EFF
72  DOUBLE PRECISION :: R_EFF, E_EFF
73 
74 ! Initialize the error manager.
75  CALL init_err_msg("CHECK_SOLIDS_DEM_ENERGY")
76 
77 ! Conduction Equations:
78 !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
79 
80 ! Set the default value for the minimum distance separating particles'
81 ! surfaces.
82  IF(des_min_cond_dist == undefined)THEN
83  des_min_cond_dist = 1.0d-04 ! cm
84  IF (units == 'SI') des_min_cond_dist = &
85  des_min_cond_dist/100.0 ! m
86  ENDIF
87 
88 ! Setup code for conduction correction terms for artificial softening
89  ! Calculate masses used for collision calculations.
90 
91 ! Shift the phase index for certain inputs to match the global phase
92 ! index until this matter can be addressed otherwise (i.e., require
93 ! the user specify correct indexing in mfix.dat). This should have no
94 ! impact if not running a hybrid case
95  mmax_tot = des_mmax+mmax
96  e_young_actual((mmax+1):mmax_tot) = e_young_actual(1:des_mmax)
97  v_poisson_actual((mmax+1):mmax_tot) = v_poisson_actual(1:des_mmax)
98  lent = mmax_tot+mmax_tot*(mmax_tot-1)/2
99  lend = des_mmax+des_mmax*(des_mmax-1)/2
100  lenc = lent-lend
101  do_area_correction = .true.
102  DO m=mmax+1,mmax_tot
103  IF (solids_model(m) /= 'DEM') cycle
104  IF (.NOT.calc_cond_des(m))cycle
105  any_conduction = .true.
106  IF(e_young_actual(m) == undefined) THEN
107  msg=''; WRITE(msg,"('Phase ',I2,' Actual EYoungs')") m
108  WRITE(err_msg,2002) 'E_YOUNG_ACTUAL', msg
109  do_area_correction = .false.
110  CALL flush_err_msg
111  ENDIF
112  IF(v_poisson_actual(m) == undefined) THEN
113  msg=''; WRITE(msg,"('Phase ',I2,' Actual poissons ratio')") m
114  WRITE(err_msg,2002) 'V_POISSON_ACTUAL', msg
115  do_area_correction = .false.
116  CALL flush_err_msg
117  ELSEIF(v_poisson_actual(m) > 0.5d0 .OR. &
118  v_poisson_actual(m) <= -one) THEN
119  WRITE(err_msg,1001) trim(ivar('V_POISSON_ACTUAL',m)), &
120  ival(v_poisson_actual(m))
121  CALL flush_err_msg(abort=.true.)
122  ENDIF
123  ENDDO
124 
125  IF(any_conduction)THEN
126  IF(ew_young_actual == undefined) THEN
127  msg=''; WRITE(msg,"('Actual EYoungs (wall)')")
128  WRITE(err_msg,2002) 'EW_YOUNG_ACTUAL', msg
129  do_area_correction = .false.
130  CALL flush_err_msg
131  ENDIF
132  IF(vw_poisson_actual == undefined) THEN
133  msg=''; WRITE(msg,"(' Actual poisson ratio (wall)')")
134  WRITE(err_msg,2002) 'VW_POISSON_ACTUAL', msg
135  do_area_correction = .false.
136  CALL flush_err_msg
137  ENDIF
138  ENDIF
139 
140 ! see above index shift
141  lc = lenc
142  DO m=mmax+1,mmax_tot
143  IF(solids_model(m) /='DEM') cycle
144  IF(.NOT.do_area_correction)cycle
145  IF(.NOT.calc_cond_des(m)) cycle
146 ! Calculate the mass of a phase M particle.
147  mass_m = (pi/6.d0)*(d_p0(m)**3)*ro_s0(m)
148 
149 ! Particle-Particle Collision Parameters ------------------------------>
150  DO l=m,mmax_tot
151  lc = lc + 1
152  IF(.NOT.calc_cond_des(l)) cycle
153 
154  mass_l = (pi/6.d0)*(d_p0(l)**3)*ro_s0(l)
155  mass_eff = (mass_m*mass_l)/(mass_m+mass_l)
156 
157 ! Calculate the effective radius, Youngs modulus, and shear modulus.
158  r_eff = 0.5d0*(d_p0(m)*d_p0(l)/(d_p0(m) + d_p0(l)))
159  e_eff = e_young_actual(m)*e_young_actual(l) / &
160  & (e_young_actual(m)*(1.d0 - v_poisson_actual(l)**2) + &
161  & e_young_actual(l)*(1.d0 - v_poisson_actual(m)**2))
162 
163 
164 ! Calculate the spring properties and store in symmetric matrix format.
165  hert_kn_actual(m,l)=(4.d0/3.d0)*sqrt(r_eff)*e_eff
166 
167 ! Compute baseline for Hertzian collision time
168  tau_c_base_actual(m,l)=3.21d0*(mass_eff/hert_kn_actual(m,l))**0.4
169  ! Can compute actual collision time via:
170  ! TAU_C_ACTUAL = TAU_C_BASE_ACTUAL * (1/ImpactVel)^0.2
171 
172 ! Compute base for simulated collision time. If Hertzian, only include base (without impact vel component).
173  IF (des_coll_model_enum .EQ. hertzian)THEN
174  tau_c_base_sim(m,l)=3.21d0*(mass_eff/hert_kn(m,l))**0.4
175  ELSE
176  tau_c_base_sim(m,l)=pi/sqrt(kn/mass_eff - &
177  ((des_etan(m,l)/mass_eff)**2)/4.d0)
178  ENDIF
179  ENDDO
180 
181 ! Do particle-wall calculations
182  mass_eff = mass_m
183  r_eff = 0.5d0*d_p0(m)
184  e_eff = e_young_actual(m)*ew_young_actual / &
185  & (e_young_actual(m)*(1.d0 - vw_poisson_actual**2) + &
186  & ew_young_actual*(1.d0 - v_poisson_actual(m)**2))
187 
188  hert_kwn_actual(m) = (4.d0/3.d0)*sqrt(r_eff)*e_eff
189  tauw_c_base_actual(m) = 3.21d0 * (mass_eff/hert_kwn_actual(m))**0.4
190 
191  IF (des_coll_model_enum .EQ. hertzian)THEN
192  tauw_c_base_sim(m)=3.21d0*(mass_eff/hert_kwn(m))**0.4
193  ELSE
194  tauw_c_base_sim(m)=pi/sqrt(kn_w/mass_eff - &
195  ((des_etan_wall(m)/mass_eff)**2)/4.d0)
196  ENDIF
197 
198  ENDDO
199 
200  CALL finl_err_msg
201 
202  RETURN
203 
204  1001 FORMAT('Warning 2001: Illegal or unknown input: ',a,' = ',a,/ &
205  'Please correct the mfix.dat file.')
206  2002 FORMAT('Warning 2002: Recommended input not specified: ',a,/ &
207  'Description:',a,/'Not correcting contact area.')
208 
209 
210  END SUBROUTINE check_solids_dem_energy
211 
212 
213 
214 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
215 ! !
216 ! Module name: CHECK_SOLIDS_DEM_COHESION !
217 ! Author: J.Musser Date: 11-DEC-13 !
218 ! !
219 ! Purpose: Check/set parameters for DES cohesion models. !
220 ! !
221 ! Comments: Original code moved from CHECK_DES_DATA !
222 ! !
223 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
224  SUBROUTINE check_solids_dem_cohesion
226 ! Global Variables:
227 !---------------------------------------------------------------------//
228 ! Runtime Flag: Invoke a cohesion model for DES simulation
229  use discretelement, only: use_cohesion
230 ! Largest discrete particle diameter.
231  use discretelement, only: max_radius
232 ! Runtime Flag: Invoke Square Well
233  use discretelement, only: square_well
234 ! Runtime Flag: Invoke Van der Waals model.
235  use discretelement, only: van_der_waals
236 ! User specified parameter for increase neighbor search area.
237  use discretelement, only: factor_rlm
238 ! User specified parameters for Van der Waals model.
239  use discretelement, only: vdw_inner_cutoff
240  use discretelement, only: vdw_outer_cutoff
241  use discretelement, only: hamaker_constant
242  use discretelement, only: wall_vdw_inner_cutoff
243  use discretelement, only: wall_vdw_outer_cutoff
244  use discretelement, only: wall_hamaker_constant
245  use discretelement, only: asperities
246  use discretelement, only: surface_energy
247  use discretelement, only: wall_surface_energy
248  use error_manager
249 
250 ! Global Parameters:
251 !---------------------------------------------------------------------//
252  use param1, only: undefined
253  use param1, only: zero
254  use constant, only: pi
255  IMPLICIT NONE
256 
257 ! Local Variables:
258 !---------------------------------------------------------------------//
259 ! Neighborhood size for Van der Waals force.
260  DOUBLE PRECISION :: VDW_NEIGHBORHOOD
261 !......................................................................!
262 
263 
264 ! Override the following settings if cohesion not used.
265  IF(.NOT.use_cohesion) THEN
266 !No more square well in the code
267  square_well = .false.
268  van_der_waals = .false.
269  wall_vdw_outer_cutoff = zero
270  RETURN
271  ENDIF
272 
273 
274 ! Initialize the error manager.
275  CALL init_err_msg("CHECK_SOLIDS_DEM_COHESION")
276 
277 
278 ! Verify that only one cohesion model is specified.
279  IF (square_well .AND. van_der_waals) THEN
280  WRITE(err_msg,1002)
281  CALL flush_err_msg(abort=.true.)
282 
283 ! Verify that at a cohesion model is specified.
284  ELSEIF(.NOT.square_well .AND. .NOT.van_der_waals) THEN
285  WRITE(err_msg,1003)
286  CALL flush_err_msg(abort=.true.)
287  ENDIF
288 
289 
290 ! Van der Waals model checks.
291  IF (van_der_waals) THEN
292 
293  IF (vdw_inner_cutoff .EQ. undefined) THEN
294  WRITE(err_msg,1201) 'VDW_INNER_CUTOFF'
295  CALL flush_err_msg(abort=.true.)
296  ENDIF
297 
298  IF(vdw_outer_cutoff .EQ. undefined) THEN
299  WRITE(err_msg,1201) 'VDW_OUTER_CUTOFF'
300  CALL flush_err_msg(abort=.true.)
301  ENDIF
302 
303  IF(hamaker_constant .EQ. undefined) THEN
304  WRITE(err_msg,1201) 'HAMAKER_CONSTANT'
305  CALL flush_err_msg(abort=.true.)
306  ENDIF
307 
308  IF (wall_vdw_inner_cutoff .EQ. undefined)THEN
309  WRITE(err_msg,1201) 'WALL_VDW_INNER_CUTOFF'
310  CALL flush_err_msg(abort=.true.)
311  ENDIF
312 
313  IF (wall_vdw_outer_cutoff .EQ. undefined)THEN
314  WRITE(err_msg,1201) 'WALL_VDW_OUTER_CUTOFF'
315  CALL flush_err_msg(abort=.true.)
316  ENDIF
317 
318  IF(wall_hamaker_constant .EQ. undefined) THEN
319  WRITE(err_msg,1201) 'WALL_HAMAKER_CONSTANT'
320  CALL flush_err_msg(abort=.true.)
321  ENDIF
322 
323  vdw_neighborhood = 1.0d0 + (vdw_outer_cutoff/(2.d0*max_radius))
324  IF (factor_rlm < vdw_neighborhood) THEN
325  WRITE(err_msg,1202)
326  CALL flush_err_msg(abort=.true.)
327  ENDIF
328 
329  IF (asperities < zero) THEN
330  WRITE(err_msg,1001) 'ASPERITIES', trim(ival(asperities))
331  CALL flush_err_msg(abort=.true.)
332  ENDIF
333 
334  surface_energy=hamaker_constant/&
335  (24.d0*pi*vdw_inner_cutoff**2)
336 
337  wall_surface_energy=wall_hamaker_constant/&
338  (24.d0*pi*wall_vdw_inner_cutoff**2)
339 
340  ENDIF
341 
342 
343  CALL finl_err_msg
344 
345  RETURN
346 
347  1001 FORMAT('Error 1001: Illegal or unknown input: ',a, ' = ',a,/ &
348  'Please correct the mfix.dat file.')
349 
350  1002 FORMAT('Error 1000: Cannot use SQUARE_WELL and VAN_DER_WAALS ', &
351  'cohesion',/'models simultaneously.')
352 
353  1003 FORMAT('Error 1001: A cohesion model was not selected. Specify ',&
354  'one of the available models in the mfix.dat file.')
355 
356 
357 
358 
359  1201 FORMAT('Error 1201: Missing input data for Van der Waals ', &
360  'cohesion model.',/'Input parameter ',a,' is UNDEFINED.')
361 
362  1202 FORMAT('Error 1202: VDW_OUTER_CUTOFF outside of the neighbor ', &
363  'search distance.',/'Increase FACTOR_RLM to increase the ', &
364  'search distance.')
365 
366  END SUBROUTINE check_solids_dem_cohesion
367 
368 
369 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
370 ! !
371 ! Subroutine: CHECK_SOLIDS_DEM_COLLISION !
372 ! Author: J.Musser Date: 11-Dec-13 !
373 ! !
374 ! Purpose: Check user input data for DES collision calculations. !
375 ! !
376 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
377  SUBROUTINE check_solids_dem_collision
379 
380 ! Global Variables:
381 !---------------------------------------------------------------------//
382 
383 ! User specified collision model
384  USE discretelement, only: des_coll_model
385  USE discretelement, only: des_coll_model_enum
386  USE discretelement, only: lsd
387  USE discretelement, only: hertzian
388 ! Particle and wall friction coeff.
389  USE discretelement, only: mew, mew_w
390 ! Parameter constatns.
391  USE param1, only: one, zero, undefined
392 
393 ! USE mpi_utility
394  use error_manager
395 
396  IMPLICIT NONE
397 
398 ! Initialize the error manager.
399  CALL init_err_msg("CHECK_SOLIDS_DEM_COLLISION")
400 
401 ! Check coefficient friction
402  IF(mew == undefined) THEN
403  WRITE(err_msg,1000) 'MEW'
404  CALL flush_err_msg(abort=.true.)
405  ELSEIF (mew < zero .OR. mew_w > one) THEN
406  WRITE(err_msg,1001) 'MEW', trim(ival(mew))
407  CALL flush_err_msg(abort=.true.)
408  ENDIF
409 
410  IF(mew_w == undefined) THEN
411  WRITE(err_msg,1000) 'MEW_W'
412  CALL flush_err_msg(abort=.true.)
413  ELSEIF(mew_w < zero .OR. mew_w > one) THEN
414  WRITE(err_msg,1001) 'MEW_W', trim(ival(mew_w))
415  CALL flush_err_msg(abort=.true.)
416  ENDIF
417 
418 ! Check collision model specific parameters.
419  SELECT CASE (trim(des_coll_model))
420 ! Linear spring-dashpot model.
421  CASE('LSD')
422  des_coll_model_enum = lsd
424 ! Hertzian collision model.
425  CASE('HERTZIAN')
426  des_coll_model_enum = hertzian
428 ! Unknown collision model.
429  CASE DEFAULT
430  WRITE(err_msg,2000) trim(des_coll_model)
431  CALL flush_err_msg(abort=.true.)
432  END SELECT
433 
434  2000 FORMAT('Error 2000: Invalid particle-particle collision model:',&
435  a,/'Please correct the mfix.dat file.')
436 
437 
438  CALL finl_err_msg
439 
440  RETURN
441 
442  1000 FORMAT('Error 1000: Required input not specified: ',a,/'Please ',&
443  'correct the mfix.dat file.')
444 
445  1001 FORMAT('Error 1001: Illegal or unknown input: ',a,' = ',a,/ &
446  'Please correct the mfix.dat file.')
447 
448  END SUBROUTINE check_solids_dem_collision
449 
450 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
451 ! !
452 ! Subroutine: CHECK_SOLIDS_DEM_COLL_LSD !
453 ! Author: J.Musser Date: 11-Dec-13 !
454 ! !
455 ! Purpose: Check user input data for DES collision calculations. !
456 ! !
457 ! References: !
458 ! - Schafer et al., J. Phys. I France, 1996, 6, 5-20 (see page 7&13) !
459 ! - Van der Hoef et al., Advances in Chemical Engineering, 2006, 31,!
460 ! 65-149 (pages 94-95) !
461 ! - Silbert et al., Physical Review E, 2001, 64, 051302 1-14 (page 5)!
462 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
463  SUBROUTINE check_solids_dem_coll_lsd
465 ! Modules
466 !---------------------------------------------------------------------//
467  use constant, only: pi
468 ! Number of discrete solids phases
469  USE discretelement, only: des_mmax
470 ! Particle and wall normal and tangential spring constants
471  USE discretelement, only: kn, kn_w
472  USE discretelement, only: kt, kt_w
473 ! Particle and wall tangential spring factor := KN/KT
474  USE discretelement, only: kt_fac, kt_w_fac
475 
476  use discretelement, only: des_etan, des_etan_wall
477  use discretelement, only: des_etat, des_etat_wall
478 
479 ! Coefficients of restitution: Normal and Tangential
480  USE discretelement, only: des_en_input, des_en_wall_input
481  USE discretelement, only: des_et_input, des_et_wall_input
482 
483 ! Tangential damping factors := ET/EN
484  USE discretelement, only: des_etat_fac, des_etat_w_fac
485 
486  use discretelement, only: dtsolid
487 
488 ! Parameter constants.
489  USE param1, only: zero, half, one, undefined
490 !
491  use physprop, only: mmax, d_p0, ro_s0
492  use run, only: solids_model
493 ! USE mpi_utility
494  use error_manager
495  IMPLICIT NONE
496 
497 ! Local Variables:
498 !---------------------------------------------------------------------//
499 ! Loop index.
500  INTEGER :: M, L, LC, MMAX_TOT
501 ! Calculate phase index offset for certain inputs until it can be
502 ! addressed in other ways. should not matter unlesss hybrid
503  INTEGER :: lent, lend, lenc
504 ! Flag to warn user.
505  LOGICAL :: FLAG_WARN
506 ! Collision length scale.
507  DOUBLE PRECISION :: TCOLL, TCOLL_TMP
508 ! Collision length scale.
509  DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
510 ! Alias for coefficient restitution
511  DOUBLE PRECISION :: EN
512 !......................................................................!
513 
514 ! Initialize the error manager.
515  CALL init_err_msg("CHECK_SOLIDS_DEM_COLL_LSD")
516 
517 ! Initialize.
518  tcoll = undefined
519 
520 ! Check for particle-particle normal spring constants.
521  IF(kn == undefined) THEN
522  WRITE(err_msg, 1000) 'KN'
523  CALL flush_err_msg(abort=.true.)
524  ENDIF
525 
526 ! Check for particle-particle tangential spring constant factors.
527  IF(kt_fac == undefined) THEN
528  WRITE (err_msg, 2100) 'KT_FAC'
529  CALL flush_err_msg()
530  kt_fac = 2.0d0/7.0d0
531  ELSEIF(kt_fac > one .OR. kt_fac < zero) THEN
532  WRITE(err_msg,1001) 'KT_FAC', trim(ival(kt_fac))
533  CALL flush_err_msg(abort=.true.)
534  ENDIF
535 ! Calculate the particle-particle tangential spring factor.
536  kt = kt_fac*kn
537 
538 ! Check for particle-wall normal spring constants.
539  IF(kn_w == undefined) THEN
540  WRITE(err_msg, 1000) 'KN_W'
541  CALL flush_err_msg(abort=.true.)
542  ENDIF
543 
544 ! Check for particle-wall tangential spring constant factors.
545  IF(kt_w_fac == undefined) THEN
546  WRITE (err_msg, 2100) 'KT_W_FAC'
547  CALL flush_err_msg()
548  kt_w_fac = 2.0d0/7.0d0
549  ELSEIF(kt_w_fac > one .OR. kt_w_fac < zero) THEN
550  WRITE(err_msg,1001) 'KT_W_FAC', trim(ival(kt_w_fac))
551  CALL flush_err_msg(abort=.true.)
552  ENDIF
553 ! Calculate the particle-wall tangential spring factor.
554  kt_w = kt_w_fac*kn_w
555 
556  2100 FORMAT('Warning 2100: Tangential spring factor ',a,' not ', &
557  'specified in mfix.dat.',/'Setting to default: (2/7).')
558 
559 ! Check for particle-particle tangential damping coefficients
560  IF(des_etat_fac == undefined) THEN
561  WRITE (err_msg, 2101) 'DES_ETAT_FAC'
562  CALL flush_err_msg
563  des_etat_fac = half
564  ELSEIF(des_etat_fac > one .OR. des_etat_fac < zero) THEN
565  WRITE(err_msg,1001) 'DES_ETAT_FAC', ival(des_etat_fac)
566  CALL flush_err_msg(abort=.true.)
567  ENDIF
568 
569 ! Check for particle-wall tangential damping coefficients
570  IF(des_etat_w_fac == undefined) THEN
571  WRITE (err_msg, 2101) 'DES_ETAT_W_FAC'
572  CALL flush_err_msg
573  des_etat_w_fac = half
574  ELSEIF(des_etat_w_fac > one .OR. des_etat_w_fac < zero) THEN
575  WRITE(err_msg,1001) 'DES_ETAT_W_FAC', ival(des_etat_w_fac)
576  CALL flush_err_msg(abort=.true.)
577  ENDIF
578 
579  2101 FORMAT('Warning 2101: Tangential damping factor ',a,' not ', &
580  'specified',/'in mfix.dat. Setting to default: (1/2).')
581 
582 ! Shift the phase index for certain inputs to match the global phase
583 ! index until this matter can be addressed otherwise (i.e., require
584 ! the user specify correct indexing in mfix.dat). This should have no
585 ! impact if not running a hybrid case
586  mmax_tot = des_mmax+mmax
587  des_en_wall_input((mmax+1):mmax_tot) = des_en_wall_input(1:des_mmax)
588  des_et_wall_input((mmax+1):mmax_tot) = des_et_wall_input(1:des_mmax)
589  lent = mmax_tot+mmax_tot*(mmax_tot-1)/2
590  lend = des_mmax+des_mmax*(des_mmax-1)/2
591  lenc = lent-lend
592  des_en_input((lenc+1):lent) = des_en_input(1:lend)
593  des_et_input((lenc+1):lent) = des_et_input(1:lend)
594  lc = lenc
595 
596  DO m = mmax+1, mmax_tot
597  IF (solids_model(m) /= 'DEM') cycle
598 
599 ! Calculate the mass of a phase M particle.
600  mass_m = (pi/6.d0)*(d_p0(m)**3)*ro_s0(m)
601 
602 ! Particle-Particle Collision Parameters ------------------------------>
603  DO l = m, mmax_tot
604 ! not necessary to cycle here given requirements on ordering of solids input
605  lc = lc+1
606 
607 ! Check particle-particle normal restitution coefficient
608  IF(des_en_input(lc) == undefined) THEN
609  WRITE(err_msg,1000) trim(ivar('DES_EN_INPUT',lc))
610  CALL flush_err_msg(abort=.true.)
611  ELSEIF(des_en_input(lc) > one .OR. &
612  des_en_input(lc) < zero) THEN
613  WRITE(err_msg,1001) trim(ivar('DES_EN_INPUT',lc)), &
614  trim(ival(des_en_input(lc)))
615  CALL flush_err_msg(abort=.true.)
616  ENDIF
617  en = des_en_input(lc)
618 
619 ! Calculate masses used for collision calculations.
620  mass_l = (pi/6.d0)*(d_p0(l)**3)*ro_s0(l)
621  mass_eff = mass_m*mass_l/(mass_m+mass_l)
622 
623 ! Calculate the M-L normal and tangential damping coefficients.
624  IF(en .NE. zero) THEN
625  des_etan(m,l) = 2.0d0*sqrt(kn*mass_eff) * abs(log(en))
626  des_etan(m,l) = des_etan(m,l)/sqrt(pi*pi + (log(en)**2))
627  ELSE
628  des_etan(m,l) = 2.0d0*sqrt(kn*mass_eff)
629  ENDIF
630  des_etat(m,l) = des_etat_fac*des_etan(m,l)
631 
632 ! Store the entries in the symmetric matrix.
633  des_etan(l,m) = des_etan(m,l)
634  des_etat(l,m) = des_etat(m,l)
635 
636 ! Calculate the collision time scale.
637  tcoll_tmp = pi/sqrt(kn/mass_eff - &
638  ((des_etan(m,l)/mass_eff)**2)/4.d0)
639  tcoll = min(tcoll_tmp, tcoll)
640  ENDDO
641 
642 
643 ! Particle-Wall Collision Parameters ---------------------------------->
644 ! Check particle-wall normal restitution coefficient.
645  IF(des_en_wall_input(m) == undefined) THEN
646  WRITE(err_msg,1000) trim(ivar('DES_EN_WALL_INPUT',m))
647  CALL flush_err_msg(abort=.true.)
648  ELSEIF(des_en_wall_input(m) > one .OR. &
649  des_en_wall_input(m) < zero) THEN
650  WRITE(err_msg,1001) trim(ivar('DES_EN_WALL_INPUT',m)), &
651  trim(ival(des_en_wall_input(m)))
652  CALL flush_err_msg(abort=.true.)
653  ENDIF
654  en = des_en_wall_input(m)
655 
656 ! Calculate masses used for collision calculations.
657  mass_eff = mass_m
658 
659 ! Calculate the M-Wall normal and tangential damping coefficients.
660  IF(en .NE. zero) THEN
661  des_etan_wall(m) = 2.d0*sqrt(kn_w*mass_eff)*abs(log(en))
662  des_etan_wall(m) = des_etan_wall(m)/sqrt(pi*pi+(log(en))**2)
663  ELSE
664  des_etan_wall(m) = 2.d0*sqrt(kn_w*mass_eff)
665  ENDIF
666  des_etat_wall(m) = des_etat_w_fac*des_etan_wall(m)
667 
668 ! Calculate the collision time scale.
669  tcoll_tmp = pi/sqrt(kn_w/mass_eff - &
670  ((des_etan_wall(m)/mass_eff)**2.d0)/4.d0)
671 ! TCOLL = MIN(TCOLL_TMP, TCOLL)
672  ENDDO
673 
674 ! if following are assigned warn user they are discarded
675  flag_warn = .false.
676  lc = lenc
677  DO m = mmax+1, mmax_tot
678  IF(solids_model(m) /= 'DEM') cycle
679  DO l = m, mmax_tot
680  lc = lc + 1
681  IF(des_et_input(m) .NE. undefined) flag_warn = .true.
682  ENDDO
683  ENDDO
684  IF (flag_warn) THEN
685  WRITE(err_msg,2102) 'DES_ET_INPUT'
686  CALL flush_err_msg
687  ENDIF
688 
689  flag_warn = .false.
690  DO m = mmax+1, mmax_tot
691  IF (solids_model(m) /= 'DEM') cycle
692  IF(des_et_wall_input(m) .NE. undefined) flag_warn = .true.
693  ENDDO
694  IF (flag_warn)THEN
695  WRITE(err_msg,2102) 'DES_ET_WALL_INPUT'
696  CALL flush_err_msg
697  ENDIF
698 
699  2102 FORMAT('Warning 2102: ',a,' values are not used ',/' with the', &
700  ' linear spring-dashpot collision model.')
701 
702 ! Store the smalled calculated collision time scale. This value is used
703 ! in time-marching the DEM solids.
704  dtsolid = tcoll/50.d0
705 
706 
707  CALL finl_err_msg
708 
709  RETURN
710 
711  1000 FORMAT('Error 1000: Required input not specified: ',a,/'Please ',&
712  'correct the mfix.dat file.')
713 
714  1001 FORMAT('Error 1001: Illegal or unknown input: ',a,' = ',a,/ &
715  'Please correct the mfix.dat file.')
716 
717  END SUBROUTINE check_solids_dem_coll_lsd
718 
719 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
720 ! !
721 ! Subroutine: CHECK_SOLIDS_DEM_COLL_HERTZ !
722 ! Author: J.Musser Date: 11-Dec-13 !
723 ! !
724 ! Purpose: Check user input data for Hertzian collisions. !
725 ! !
726 ! References: !
727 ! - Schafer et al., J. Phys. I France, 1996, 6, 5-20 (see page 7&13) !
728 ! - Van der Hoef et al., Advances in Chemical Engineering, 2006, 31,!
729 ! 65-149 (pages 94-95) !
730 ! - Silbert et al., Physical Review E, 2001, 64, 051302 1-14 (page 5)!
731 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
732  SUBROUTINE check_solids_dem_coll_hertz
734 ! Modules
735 !---------------------------------------------------------------------//
736 ! The constant PI
737  use constant, only: pi
738 ! Number of discrete solids phases, diameters and densities
739  USE discretelement, only: des_mmax
740 ! User defined coefficients of restitution: Normal and Tangential
741  USE discretelement, only: des_en_input, des_en_wall_input
742  USE discretelement, only: des_et_input, des_et_wall_input
743 ! Damping coefficients: Normal and Tangential
744  use discretelement, only: des_etan, des_etan_wall
745  use discretelement, only: des_etat, des_etat_wall
746 ! Hertzian spring constants: Normal and Tangential
747  use discretelement, only: hert_kn, hert_kwn
748  use discretelement, only: hert_kt, hert_kwt
749 ! Tangential damping factors := ET/EN
750  USE discretelement, only: des_etat_fac, des_etat_w_fac
751 ! Particle and wall Young's modulus and Shear modulus
752  USE discretelement, only: e_young, ew_young
753 ! Particle and wall Poisson ratio
754  USE discretelement, only: v_poisson, vw_poisson
755 ! Solids time step-size.
756  use discretelement, only: dtsolid
757 ! Particle and wall normal spring constants
758  USE discretelement, only: kn, kn_w
759 ! Particle and wall tangential spring factor := KN/KT
760  USE discretelement, only: kt_fac, kt_w_fac
761 
762  use param, only: dim_m
763 ! Parameter constatns.
764  USE param1, only: zero, one, undefined
765 
766  USE physprop, only: mmax, ro_s0, d_p0
767  USE run, only: solids_model
768 ! USE mpi_utility
769  use error_manager
770 
771  IMPLICIT NONE
772 
773 ! Local Variables:
774 !---------------------------------------------------------------------//
775 ! Loop index.
776  INTEGER :: M, L, LC, MMAX_TOT
777 ! Calculate phase index offset for certain inputs until it can be
778 ! addressed in other ways. should not matter unlesss hybrid
779  INTEGER :: lent, lend, lenc
780 ! Message for formatted output.
781  CHARACTER(len=64) :: MSG
782 ! Collision length scale.
783  DOUBLE PRECISION :: TCOLL, TCOLL_TMP
784 ! Particle and effective mass.
785  DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
786 ! Effective physical quantities. Radius, Youngs, Shear
787  DOUBLE PRECISION :: R_EFF, E_EFF, G_MOD_EFF, RED_MASS_EFF
788 ! Alias for coefficient restitution
789  DOUBLE PRECISION :: EN, ET
790 ! shear modulus
791  DOUBLE PRECISION :: G_MOD(dim_m)
792 ! Shear modules for wall
793  DOUBLE PRECISION :: G_MOD_WALL
794 !......................................................................!
795 
796 ! Initialize the error manager.
797  CALL init_err_msg("CHECK_SOLIDS_DEM_COLL_HERTZ")
798 
799 ! Initialize.
800  tcoll = undefined
801 
802 ! check young's modulus and poisson ratio
803  IF(ew_young == undefined ) THEN
804  msg='Wall value for Youngs modulus'
805  WRITE(err_msg,1002) 'Ew_YOUNG', msg
806  CALL flush_err_msg(abort=.true.)
807  ENDIF
808 
809  IF(vw_poisson == undefined) THEN
810  msg='Wall value for Poissons ratio'
811  WRITE(err_msg,1002) 'Vw_POISSON', msg
812  CALL flush_err_msg(abort=.true.)
813  ELSEIF (vw_poisson > 0.5d0 .OR. vw_poisson <= -one) THEN
814  WRITE(err_msg,1001) 'Vw_POISSON',ival(vw_poisson)
815  CALL flush_err_msg(abort=.true.)
816  ENDIF
817 
818  g_mod_wall = 0.5d0*ew_young/(1.d0+vw_poisson)
819 
820 ! Shift the phase index for certain inputs to match the global phase
821 ! index until this matter can be addressed otherwise (i.e., require
822 ! the user specify correct indexing in mfix.dat). This should have no
823 ! impact if not running a hybrid case
824  mmax_tot = des_mmax+mmax
825  e_young((mmax+1):mmax_tot) = e_young(1:des_mmax)
826  v_poisson((mmax+1):mmax_tot) = v_poisson(1:des_mmax)
827  des_en_wall_input((mmax+1):mmax_tot) = des_en_wall_input(1:des_mmax)
828  des_et_wall_input((mmax+1):mmax_tot) = des_et_wall_input(1:des_mmax)
829  lent = mmax_tot+mmax_tot*(mmax_tot-1)/2
830  lend = des_mmax+des_mmax*(des_mmax-1)/2
831  lenc = lent-lend
832  des_en_input((lenc+1):lent) = des_en_input(1:lend)
833  des_et_input((lenc+1):lent) = des_et_input(1:lend)
834 
835  DO m=mmax+1,mmax_tot
836  IF (solids_model(m) /= 'DEM') cycle
837 
838  IF(e_young(m) == undefined) THEN
839  msg=''; WRITE(msg,"('Phase ',I2,' Youngs modulus')") m
840  WRITE(err_msg,1002) 'E_YOUNG', msg
841  CALL flush_err_msg(abort=.true.)
842  ENDIF
843  IF(v_poisson(m) == undefined) THEN
844  msg=''; WRITE(msg,"('Phase ',I2,' Poissons ratio')") m
845  WRITE(err_msg,1002) 'V_POISSON', msg
846  CALL flush_err_msg(abort=.true.)
847  ELSEIF(v_poisson(m) > 0.5d0 .OR. &
848  v_poisson(m) <= -one) THEN
849  WRITE(err_msg,1001) trim(ivar('V_POISSON',m)), &
850  ival(v_poisson(m))
851  CALL flush_err_msg(abort=.true.)
852  ENDIF
853 ! Calculate the shear modulus for phase M.
854  g_mod(m) = 0.5d0*e_young(m)/(1.d0+v_poisson(m))
855  ENDDO
856 
857 ! see above index shift
858  lc = lenc
859  DO m=mmax+1,mmax_tot
860  IF(solids_model(m) /='DEM') cycle
861 
862 ! Calculate the mass of a phase M particle.
863  mass_m = (pi/6.d0)*(d_p0(m)**3)*ro_s0(m)
864 
865 ! Particle-Particle Collision Parameters ------------------------------>
866  DO l=m,mmax_tot
867 ! not necessary to cycle here given requirements on ordering of solids input
868  lc = lc+1
869 
870 ! Check particle-particle normal restitution coefficient
871  IF(des_en_input(lc) == undefined) THEN
872  WRITE(err_msg,1000) trim(ivar('DES_EN_INPUT',lc))
873  CALL flush_err_msg(abort=.true.)
874  ELSEIF(des_en_input(lc) > one .OR. &
875  des_en_input(lc) < zero) THEN
876  WRITE(err_msg,1001) trim(ivar('DES_EN_INPUT',lc)), &
877  trim(ival(des_en_input(lc)))
878  CALL flush_err_msg(abort=.true.)
879  ENDIF
880  en = des_en_input(lc)
881 
882 ! Check particle-particle tangential restitution coefficient
883  IF(des_et_input(lc) == undefined) THEN
884  WRITE(err_msg,1000) trim(ivar('DES_ET_INPUT',lc))
885  CALL flush_err_msg(abort=.true.)
886  ELSEIF(des_et_input(lc) > one .OR. &
887  des_et_input(lc) < zero) THEN
888  WRITE(err_msg,1001) trim(ivar('DES_ET_INPUT',lc)), &
889  ival(des_et_input(lc))
890  CALL flush_err_msg(abort=.true.)
891  ENDIF
892  et = des_et_input(lc)
893 
894 
895 ! Calculate masses used for collision calculations.
896  mass_l = (pi/6.d0)*(d_p0(l)**3)*ro_s0(l)
897  mass_eff = (mass_m*mass_l)/(mass_m+mass_l)
898  red_mass_eff = (2.d0/7.d0)*mass_eff
899 ! Calculate the effective radius, Youngs modulus, and shear modulus.
900  r_eff = 0.5d0*(d_p0(m)*d_p0(l)/(d_p0(m) + d_p0(l)))
901  e_eff = e_young(m)*e_young(l) / &
902  (e_young(m)*(1.d0 - v_poisson(l)**2) + &
903  e_young(l)*(1.d0 - v_poisson(m)**2))
904  g_mod_eff = g_mod(m)*g_mod(l)/ &
905  (g_mod(m)*(2.d0 - v_poisson(l)) + &
906  g_mod(l)*(2.d0 - v_poisson(m)))
907 
908 ! Calculate the spring properties and store in symmetric matrix format.
909  hert_kn(m,l)=(4.d0/3.d0)*sqrt(r_eff)*e_eff
910  hert_kt(m,l)= 8.d0*sqrt(r_eff)*g_mod_eff
911 
912  hert_kn(l,m) = hert_kn(m,l)
913  hert_kt(l,m) = hert_kt(m,l)
914 
915 ! Calculate the normal damping coefficient.
916  IF(en .NE. zero) THEN
917  des_etan(m,l) = 2.d0*sqrt(hert_kn(m,l)*mass_eff)* &
918  abs(log(en))
919  des_etan(m,l) = des_etan(m,l)/ &
920  sqrt(pi*pi + (log(en))**2)
921  ELSE
922  des_etan(m,l) = 2.d0*sqrt(hert_kn(m,l)*mass_eff)
923  ENDIF
924  des_etan(l,m) = des_etan(m,l)
925 
926 ! Calculate the tangential coefficients.
927  IF(et .NE. zero) THEN
928  des_etat(m,l) = 2.d0*sqrt(hert_kt(m,l)*red_mass_eff)* &
929  abs(log(et))
930  des_etat(m,l) = des_etat(m,l)/ sqrt(pi*pi+(log(et))**2)
931  ELSE
932  des_etat(m,l) = 2.d0*sqrt(hert_kt(m,l)*red_mass_eff)
933  ENDIF
934  des_etat(l,m) = des_etat(m,l)
935 
936  tcoll_tmp = pi/sqrt(hert_kn(m,l)/mass_eff - &
937  ((des_etan(m,l)/mass_eff)**2)/4.d0)
938  tcoll = min(tcoll_tmp, tcoll)
939  ENDDO
940 
941 ! Particle-Wall Collision Parameters ---------------------------------->
942 ! Check particle-wall normal restitution coefficient.
943  IF(des_en_wall_input(m) == undefined) THEN
944  WRITE(err_msg,1000) trim(ivar('DES_EN_WALL_INPUT',m))
945  CALL flush_err_msg(abort=.true.)
946  ELSEIF(des_en_wall_input(m) > one .OR. &
947  des_en_wall_input(m) < zero) THEN
948  WRITE(err_msg,1001) trim(ivar('DES_EN_WALL_INPUT',m)), &
949  trim(ival(des_en_wall_input(m)))
950  CALL flush_err_msg(abort=.true.)
951  ENDIF
952  en = des_en_wall_input(m)
953 
954 ! Check particle-wall tangential restitution coefficient
955  IF(des_et_wall_input(m) == undefined) THEN
956  WRITE(err_msg,1000) trim(ivar('DES_ET_WALL_INPUT',m))
957  CALL flush_err_msg(abort=.true.)
958  ELSEIF(des_et_wall_input(m) > one .OR. &
959  des_et_wall_input(m) < zero) THEN
960  WRITE(err_msg,1001) trim(ivar('DES_ET_WALL_INPUT',m)), &
961  trim(ival(des_et_wall_input(m)))
962  CALL flush_err_msg(abort=.true.)
963  ENDIF
964  et = des_et_wall_input(m)
965 
966 ! Calculate masses used for collision calculations.
967  mass_eff = mass_m
968  red_mass_eff = (2.d0/7.d0)*mass_eff
969 ! Calculate the effective radius, Youngs modulus, and shear modulus.
970  r_eff = 0.5d0*d_p0(m)
971  e_eff = e_young(m)*ew_young / &
972  (e_young(m)*(1.d0-vw_poisson**2) + &
973  ew_young *(1.d0-v_poisson(m)**2))
974  g_mod_eff = g_mod(m)*g_mod_wall / &
975  (g_mod(m)*(2.d0 - vw_poisson) + &
976  g_mod_wall*(2.d0 - v_poisson(m)))
977 
978 ! Calculate the spring properties.
979  hert_kwn(m) = (4.d0/3.d0)*sqrt(r_eff)*e_eff
980  hert_kwt(m) = 8.0*sqrt(r_eff)*g_mod_eff
981 
982 ! Calculate the tangential coefficients.
983  IF(en /= zero) THEN
984  des_etan_wall(m) = 2.d0*sqrt(hert_kwn(m)*mass_eff)*&
985  abs(log(en))
986  des_etan_wall(m) = des_etan_wall(m)/&
987  sqrt(pi*pi + (log(en))**2)
988  ELSE
989  des_etan_wall(m) = 2.d0*sqrt(hert_kwn(m)*mass_eff)
990  ENDIF
991 
992  IF(et /= zero) THEN
993  des_etat_wall(m) = 2.d0*sqrt(hert_kwt(m)*red_mass_eff)* &
994  abs(log(et))
995  des_etat_wall(m) = des_etat_wall(m)/sqrt(pi*pi+(log(et))**2)
996  ELSE
997  des_etat_wall(m) = 2.d0*sqrt(hert_kwt(m)*red_mass_eff)
998  ENDIF
999 
1000 ! Calculate the collision time scale.
1001  tcoll_tmp = pi/sqrt(hert_kwn(m)/mass_eff - &
1002  ((des_etan_wall(m)/mass_eff)**2)/4.d0)
1003  ENDDO
1004 
1005 
1006 ! If following are assigned warn user they are discarded.
1007  IF(kn .NE. undefined) THEN
1008  WRITE(err_msg, 2200) 'KN'
1009  CALL flush_err_msg
1010  ENDIF
1011  IF(kn_w .NE. undefined) THEN
1012  WRITE(err_msg, 2200) 'KN_W'
1013  CALL flush_err_msg
1014  ENDIF
1015  IF(kt_fac .NE. undefined) THEN
1016  WRITE(err_msg, 2200) 'KT_FAC'
1017  CALL flush_err_msg
1018  ENDIF
1019  IF(kt_w_fac .NE. undefined) THEN
1020  WRITE(err_msg, 2200) 'KT_W_FAC'
1021  CALL flush_err_msg
1022  ENDIF
1023  IF(des_etat_fac .NE. undefined) THEN
1024  WRITE(err_msg, 2200) 'DES_ETAT_FAC'
1025  CALL flush_err_msg
1026  ENDIF
1027  IF(des_etat_w_fac .NE. undefined) THEN
1028  WRITE(err_msg, 2200) 'DES_ETAT_W_FAC'
1029  CALL flush_err_msg
1030  ENDIF
1031 
1032  2200 FORMAT('Warning 2200: ',a,' values are not used ',/' with the', &
1033  ' linear spring-dashpot collision model.')
1034 
1035 
1036 ! Store the smalled calculated collision time scale. This value is used
1037 ! in time-marching the DEM solids.
1038  dtsolid = tcoll/50.d0
1039 
1040 
1041  CALL finl_err_msg
1042 
1043  RETURN
1044 
1045  1000 FORMAT('Error 1000: Required input not specified: ',a,/'Please ',&
1046  'correct the mfix.dat file.')
1047 
1048  1001 FORMAT('Error 1001: Illegal or unknown input: ',a,' = ',a,/ &
1049  'Please correct the mfix.dat file.')
1050 
1051  1002 FORMAT('Error 1002: Required input not specified: ',a,/ &
1052  'Description:',a,/'Please correct the mfix.dat file.')
1053 
1054  END SUBROUTINE check_solids_dem_coll_hertz
1055 
1056 
1057 
1058 
1059 
1060 
1061 
1062 
1063 
1064 
subroutine check_solids_dem
double precision, dimension(dim_m) d_p0
Definition: physprop_mod.f:25
character(len=32) function ivar(VAR, i1, i2, i3)
subroutine finl_err_msg
double precision, parameter one
Definition: param1_mod.f:29
logical, dimension(dim_m) calc_cond_des
integer, parameter dim_m
Definition: param_mod.f:67
double precision des_min_cond_dist
character(len=3), dimension(dim_m) solids_model
Definition: run_mod.f:253
double precision, parameter undefined
Definition: param1_mod.f:18
subroutine init_err_msg(CALLER)
integer mmax
Definition: physprop_mod.f:19
subroutine check_solids_dem_collision
subroutine check_solids_dem_energy
double precision, parameter half
Definition: param1_mod.f:28
Definition: run_mod.f:13
Definition: param_mod.f:2
character(len=16) units
Definition: run_mod.f:30
subroutine check_solids_dem_cohesion
subroutine check_solids_dem_coll_hertz
character(len=line_length), dimension(line_count) err_msg
subroutine check_solids_dem_coll_lsd
double precision, dimension(dim_m) ro_s0
Definition: physprop_mod.f:28
double precision, parameter pi
Definition: constant_mod.f:158
double precision, parameter zero
Definition: param1_mod.f:27
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)