File: /nfs/home/0/users/jenkins/mfix.git/model/check_data/check_solids_dem.f

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     
12     ! Global Variables:
13     !---------------------------------------------------------------------//
14     ! Runtime Flag: Calculate clusters during a DEM simulation
15           use discretelement, only: FACTOR_RLM
16           use discretelement, only: MAX_RADIUS
17     
18     ! Global Parameters:
19     !---------------------------------------------------------------------//
20           use param1, only: UNDEFINED
21     
22     ! Use the error manager for posting error messages.
23     !---------------------------------------------------------------------//
24           use error_manager
25     
26           implicit none
27     
28     ! Initialize the error manager.
29           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM")
30     
31     ! Particle-particle collision parameters.
32           CALL CHECK_SOLIDS_DEM_COLLISION
33     ! DES cohesion model parameters.
34           CALL CHECK_SOLIDS_DEM_COHESION
35     ! Particle-particle conduction model parameters.
36           CALL CHECK_SOLIDS_DEM_ENERGY
37     
38           CALL FINL_ERR_MSG
39     
40           RETURN
41     
42      1000 FORMAT('Error 1000: Required input not specified: ',A,/'Please ',&
43              'correct the mfix.dat file.')
44     
45      1001 FORMAT('Error 1001: Illegal or unknown input: ',A,' = ',A,/   &
46              'Please correct the mfix.dat file.')
47     
48           END SUBROUTINE CHECK_SOLIDS_DEM
49     
50     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
51     !                                                                      !
52     !  SUBROUTINE CHECK_SOLIDS_DEM_ENERGY                                  !
53     !  Author: J.Musser                                   Date: 02-FEB-14  !
54     !                                                                      !
55     !  Purpose:                                                            !
56     !                                                                      !
57     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
58           SUBROUTINE CHECK_SOLIDS_DEM_ENERGY
59     
60           use des_thermo, only: DES_MIN_COND_DIST
61           use run, only: UNITS
62     
63     ! Global Parameters:
64     !---------------------------------------------------------------------//
65           use param1, only: UNDEFINED
66     
67           use error_manager
68     
69           IMPLICIT NONE
70     
71     
72     ! Initialize the error manager.
73           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM_ENERGY")
74     
75     ! Conduction Equations:
76     !'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
77     
78     ! Set the default value for the minimum distance separating particles'
79     ! surfaces.
80           IF(DES_MIN_COND_DIST == UNDEFINED)THEN
81              DES_MIN_COND_DIST = 1.0D-04 ! cm
82              IF (UNITS == 'SI') DES_MIN_COND_DIST = &
83                 DES_MIN_COND_DIST/100.0  ! m
84           ENDIF
85     
86           CALL FINL_ERR_MSG
87     
88           RETURN
89     
90           END SUBROUTINE CHECK_SOLIDS_DEM_ENERGY
91     
92     
93     
94     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
95     !                                                                      !
96     !  Module name: CHECK_SOLIDS_DEM_COHESION                              !
97     !  Author: J.Musser                                   Date: 11-DEC-13  !
98     !                                                                      !
99     !  Purpose: Check/set parameters for DES cohesion models.              !
100     !                                                                      !
101     !  Comments: Original code moved from CHECK_DES_DATA                   !
102     !                                                                      !
103     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
104           SUBROUTINE CHECK_SOLIDS_DEM_COHESION
105     
106     ! Global Variables:
107     !---------------------------------------------------------------------//
108     
109     ! Runtime Flag: Invoke a cohesion model for DES simulation
110           use discretelement, only: USE_COHESION
111     ! Largest discrete particle diameter.
112           use discretelement, only: MAX_RADIUS
113     
114     ! Runtime Flag: Invoke Square Well
115           use discretelement, only: SQUARE_WELL
116     ! Runtime Flag: Invoke Van der Waals model.
117           use discretelement, only: VAN_DER_WAALS
118     ! User specified parameter for increase neighbor search area.
119           use discretelement, only: FACTOR_RLM
120     
121     ! User specified parameters for Van der Waals model.
122           use discretelement, only: VDW_INNER_CUTOFF
123           use discretelement, only: VDW_OUTER_CUTOFF
124           use discretelement, only: HAMAKER_CONSTANT
125           use discretelement, only: WALL_VDW_INNER_CUTOFF
126           use discretelement, only: WALL_VDW_OUTER_CUTOFF
127           use discretelement, only: WALL_HAMAKER_CONSTANT
128           use discretelement, only: ASPERITIES
129           use discretelement, only: SURFACE_ENERGY
130           use discretelement, only: WALL_SURFACE_ENERGY
131     
132     
133     ! Global Parameters:
134     !---------------------------------------------------------------------//
135           use param1, only: UNDEFINED
136           use param1, only: ZERO
137           use constant, only: Pi
138     
139           use error_manager
140     
141           IMPLICIT NONE
142     
143     ! Local Variables:
144     !---------------------------------------------------------------------//
145     ! Neighborhood size for Van der Waals force.
146           DOUBLE PRECISION :: VDW_NEIGHBORHOOD
147     
148     
149     ! Override the following settings if cohesion not used.
150           IF(.NOT.USE_COHESION) THEN
151     !No more square well in the code
152              SQUARE_WELL = .FALSE.
153              VAN_DER_WAALS = .FALSE.
154              WALL_VDW_OUTER_CUTOFF = ZERO
155              RETURN
156           ENDIF
157     
158     
159     ! Initialize the error manager.
160           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM_COHESION")
161     
162     
163     ! Verify that only one cohesion model is specified.
164           IF (SQUARE_WELL .AND. VAN_DER_WAALS) THEN
165              WRITE(ERR_MSG,1002)
166              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
167     
168     ! Verify that at a cohesion model is specified.
169           ELSEIF(.NOT.SQUARE_WELL .AND. .NOT.VAN_DER_WAALS) THEN
170              WRITE(ERR_MSG,1003)
171              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
172           ENDIF
173     
174     
175     ! Van der Waals model checks.
176           IF (VAN_DER_WAALS) THEN
177     
178              IF (VDW_INNER_CUTOFF .EQ. UNDEFINED) THEN
179                 WRITE(ERR_MSG,1201) 'VDW_INNER_CUTOFF'
180                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
181              ENDIF
182     
183              IF(VDW_OUTER_CUTOFF .EQ. UNDEFINED) THEN
184                 WRITE(ERR_MSG,1201) 'VDW_OUTER_CUTOFF'
185                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
186              ENDIF
187     
188              IF(HAMAKER_CONSTANT .EQ. UNDEFINED) THEN
189                 WRITE(ERR_MSG,1201) 'HAMAKER_CONSTANT'
190                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
191              ENDIF
192     
193              IF (WALL_VDW_INNER_CUTOFF .EQ. UNDEFINED)THEN
194                 WRITE(ERR_MSG,1201) 'WALL_VDW_INNER_CUTOFF'
195                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
196              ENDIF
197     
198              IF (WALL_VDW_OUTER_CUTOFF .EQ. UNDEFINED)THEN
199                 WRITE(ERR_MSG,1201) 'WALL_VDW_OUTER_CUTOFF'
200                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
201              ENDIF
202     
203              IF(WALL_HAMAKER_CONSTANT .EQ. UNDEFINED) THEN
204                 WRITE(ERR_MSG,1201) 'WALL_HAMAKER_CONSTANT'
205                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
206              ENDIF
207     
208              VDW_NEIGHBORHOOD = 1.0d0 + (VDW_OUTER_CUTOFF/(2.d0*MAX_RADIUS))
209              IF (FACTOR_RLM < VDW_NEIGHBORHOOD) THEN
210                 WRITE(ERR_MSG,1202)
211                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
212              ENDIF
213     
214              IF (ASPERITIES < ZERO) THEN
215                 WRITE(ERR_MSG,1001) 'ASPERITIES', trim(iVal(ASPERITIES))
216                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
217              ENDIF
218     
219              SURFACE_ENERGY=HAMAKER_CONSTANT/&
220                 (24.d0*Pi*VDW_INNER_CUTOFF**2)
221     
222              WALL_SURFACE_ENERGY=WALL_HAMAKER_CONSTANT/&
223                 (24.d0*Pi*WALL_VDW_INNER_CUTOFF**2)
224     
225           ENDIF
226     
227     
228           CALL FINL_ERR_MSG
229     
230           RETURN
231     
232      1001 FORMAT('Error 1001: Illegal or unknown input: ',A, ' = ',A,/     &
233                  'Please correct the mfix.dat file.')
234     
235      1002 FORMAT('Error 1000: Cannot use SQUARE_WELL and VAN_DER_WAALS ',  &
236              'cohesion',/'models simultaneously.')
237     
238      1003 FORMAT('Error 1001: A cohesion model was not selected. Specify ',&
239              'one of the available models in the mfix.dat file.')
240     
241     
242     !<------------------- Van der Waals model messages. ----------------->!
243     
244      1201 FORMAT('Error 1201: Missing input data for Van der Waals ',      &
245              'cohesion model.',/'Input parameter ',A,' is UNDEFINED.')
246     
247      1202 FORMAT('Error 1202: VDW_OUTER_CUTOFF outside of the neighbor ',  &
248              'search distance.',/'Increase FACTOR_RLM to increase the ',   &
249              'search distance.')
250     
251           END SUBROUTINE CHECK_SOLIDS_DEM_COHESION
252     
253     
254     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
255     !                                                                      !
256     !  Subroutine: CHECK_SOLIDS_DEM_COLLISION                              !
257     !  Author: J.Musser                                   Date: 11-Dec-13  !
258     !                                                                      !
259     !  Purpose: Check user input data for DES collision calculations.      !
260     !                                                                      !
261     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
262           SUBROUTINE CHECK_SOLIDS_DEM_COLLISION
263     
264     
265     ! Global Variables:
266     !---------------------------------------------------------------------//
267     
268     ! User specified collision model
269           USE discretelement, only: DES_COLL_MODEL
270           USE discretelement, only: DES_COLL_MODEL_ENUM
271           USE discretelement, only: LSD
272           USE discretelement, only: HERTZIAN
273     ! Number of discrete solids phases
274           USE discretelement, only: DES_MMAX, DES_D_P0, DES_RO_s
275     ! Particle and wall friction coeff.
276           USE discretelement, only: MEW, MEW_W
277     ! Parameter constatns.
278           USE param1, only: ONE, ZERO, UNDEFINED
279     
280     !      USE mpi_utility
281           use error_manager
282     
283           IMPLICIT NONE
284     
285     ! Initialize the error manager.
286           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM_COLLISION")
287     
288     ! Check coefficient friction
289           IF(MEW == UNDEFINED) THEN
290              WRITE(ERR_MSG,1000) 'MEW'
291              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
292           ELSEIF (MEW < ZERO .OR. MEW_W > ONE) THEN
293              WRITE(ERR_MSG,1001) 'MEW', trim(iVal(MEW))
294              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
295           ENDIF
296     
297           IF(MEW_W == UNDEFINED) THEN
298              WRITE(ERR_MSG,1000) 'MEW_W'
299              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
300           ELSEIF(MEW_w < ZERO .OR. MEW_W > ONE) THEN
301              WRITE(ERR_MSG,1001) 'MEW_W', trim(iVal(MEW_W))
302              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
303           ENDIF
304     
305     ! Check collision model specific parameters.
306           SELECT CASE (trim(DES_COLL_MODEL))
307     ! Linear spring-dashpot model.
308           CASE('LSD')
309             DES_COLL_MODEL_ENUM = LSD
310             CALL CHECK_SOLIDS_DEM_COLL_LSD
311     ! Hertzian collision model.
312           CASE('HERTZIAN')
313              DES_COLL_MODEL_ENUM = HERTZIAN
314              CALL CHECK_SOLIDS_DEM_COLL_HERTZ
315     ! Unknown collision model.
316           CASE DEFAULT
317              WRITE(ERR_MSG,2000) TRIM(DES_COLL_MODEL)
318              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
319           END SELECT
320     
321      2000 FORMAT('Error 2000: Invalid particle-particle collision model:',&
322              A,/'Please correct the mfix.dat file.')
323     
324     
325           CALL FINL_ERR_MSG
326     
327           RETURN
328     
329      1000 FORMAT('Error 1000: Required input not specified: ',A,/'Please ',&
330              'correct the mfix.dat file.')
331     
332      1001 FORMAT('Error 1001: Illegal or unknown input: ',A,' = ',A,/      &
333              'Please correct the mfix.dat file.')
334     
335      1002 FORMAT('Error 1002: Required input not specified: ',A,/          &
336              'Description:',A,/'Please correct the mfix.dat file.')
337     
338           END SUBROUTINE CHECK_SOLIDS_DEM_COLLISION
339     
340     
341     
342     
343     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
344     !                                                                      !
345     !  Subroutine: CHECK_SOLIDS_DEM_COLL_LSD                               !
346     !  Author: J.Musser                                   Date: 11-Dec-13  !
347     !                                                                      !
348     !  Purpose: Check user input data for DES collision calculations.      !
349     !                                                                      !
350     !  References:                                                         !
351     !   - Schafer et al., J. Phys. I France, 1996, 6, 5-20 (see page 7&13) !
352     !   -  Van der Hoef et al., Advances in Chemical Engineering, 2006, 31,!
353     !      65-149 (pages 94-95)                                            !
354     !   - Silbert et al., Physical Review E, 2001, 64, 051302 1-14 (page 5)!
355     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
356           SUBROUTINE CHECK_SOLIDS_DEM_COLL_LSD
357     
358     
359     ! Global Variables:
360     !---------------------------------------------------------------------//
361     ! Number of discrete solids phases
362           USE discretelement, only: DES_MMAX, DES_D_P0, DES_RO_s
363     ! Particle and wall friction coeff.
364           USE discretelement, only: MEW, MEW_W
365     ! Particle and wall normal and tangential spring constants
366           USE discretelement, only: KN, KN_W
367           USE discretelement, only: KT, KT_W
368     ! Particle and wall tangential spring factor := KN/KT
369           USE discretelement, only: KT_FAC, KT_W_FAC
370     
371           use discretelement, only: DES_ETAN, DES_ETAN_WALL
372           use discretelement, only: DES_ETAT, DES_ETAT_WALL
373     
374     ! Coefficients of restitution: Normal and Tangential
375           USE discretelement, only: DES_EN_INPUT, DES_EN_WALL_INPUT
376           USE discretelement, only: DES_ET_INPUT, DES_ET_WALL_INPUT
377     
378     ! Tangential damping factors := ET/EN
379           USE discretelement, only: DES_ETAT_FAC, DES_ETAT_W_FAC
380     ! Particle and wall Young's modulus
381           USE discretelement, only: E_YOUNG, Ew_YOUNG
382     ! Particle and wall Poisson ratio
383           USE discretelement, only: V_POISSON, Vw_POISSON
384     ! Shear modulus
385           USE discretelement, only: G_MOD
386     
387           use constant, only: PI
388     
389           use discretelement, only: DTSOLID
390     
391     ! Parameter constatns.
392           USE param1, only: ZERO, HALF, ONE, UNDEFINED
393     
394     !      USE mpi_utility
395           use error_manager
396     
397           IMPLICIT NONE
398     
399     
400     ! Local Variables:
401     !---------------------------------------------------------------------//
402     ! Loop index.
403           INTEGER :: M, L, LC
404     ! Flag to warn user.
405           LOGICAL :: FLAG_WARN
406     ! Collision length scale.
407           DOUBLE PRECISION :: TCOLL, TCOLL_TMP
408     ! Collision length scale.
409           DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
410     ! Alias for coefficient restitution
411           DOUBLE PRECISION :: EN
412     
413     ! Initialize the error manager.
414           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM_COLL_LSD")
415     
416     ! Initialize.
417           TCOLL = UNDEFINED
418     
419     ! Check for particle-particle normal spring constants.
420           IF(KN == UNDEFINED) THEN
421              WRITE(ERR_MSG, 1000) 'KN'
422              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
423           ENDIF
424     
425     ! Check for particle-particle tangential spring constant factors.
426           IF(KT_FAC == UNDEFINED) THEN
427              WRITE (ERR_MSG, 2100) 'KT_FAC'
428              CALL FLUSH_ERR_MSG()
429              KT_FAC = 2.0d0/7.0d0
430           ELSEIF(KT_FAC > ONE .OR. KT_FAC < ZERO) THEN
431              WRITE(ERR_MSG,1001) 'KT_FAC', trim(iVal(KT_FAC))
432              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
433           ENDIF
434     ! Calculate the particle-particle tangential spring factor.
435           KT = KT_FAC*KN
436     
437     ! Check for particle-wall normal spring constants.
438           IF(KN_W == UNDEFINED) THEN
439              WRITE(ERR_MSG, 1000) 'KN_W'
440              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
441           ENDIF
442     
443     ! Check for particle-wall tangential spring constant factors.
444           IF(KT_W_FAC == UNDEFINED) THEN
445              WRITE (ERR_MSG, 2100) 'KT_W_FAC'
446              CALL FLUSH_ERR_MSG()
447              KT_W_FAC = 2.0d0/7.0d0
448           ELSEIF(KT_W_FAC > ONE .OR. KT_W_FAC < ZERO) THEN
449              WRITE(ERR_MSG,1001) 'KT_W_FAC', trim(iVal(KT_W_FAC))
450              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
451           ENDIF
452     ! Calculate the particle-wall tangential spring factor.
453           KT_W = KT_W_FAC*KN_W
454     
455      2100 FORMAT('Warning 2100: Tangential spring factor ',A,' not ',      &
456              'specified in mfix.dat.',/'Setting to default: (2/7).')
457     
458     ! Check for particle-particle tangential damping coefficients
459           IF(DES_ETAT_FAC == UNDEFINED) THEN
460              WRITE (ERR_MSG, 2101) 'DES_ETAT_FAC'
461              CALL FLUSH_ERR_MSG
462              DES_ETAT_FAC = HALF
463           ELSEIF(DES_ETAT_FAC > ONE .OR. DES_ETAT_FAC < ZERO) THEN
464              WRITE(ERR_MSG,1001) 'DES_ETAT_FAC', iVal(DES_ETAT_FAC)
465              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
466           ENDIF
467     
468     ! Check for particle-wall tangential damping coefficients
469           IF(DES_ETAT_W_FAC == UNDEFINED) THEN
470              WRITE (ERR_MSG, 2101) 'DES_ETAT_W_FAC'
471              CALL FLUSH_ERR_MSG
472              DES_ETAT_W_FAC = HALF
473           ELSEIF(DES_ETAT_W_FAC > ONE .OR. DES_ETAT_W_FAC < ZERO) THEN
474              WRITE(ERR_MSG,1001) 'DES_ETAT_W_FAC', iVal(DES_ETAT_W_FAC)
475              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
476           ENDIF
477     
478      2101 FORMAT('Warning 2101: Tangential damping factor ',A,' not ',     &
479              'specified',/'in mfix.dat. Setting to default: (1/2).')
480     
481     
482           LC = 0
483           DO M = 1, DES_MMAX
484     
485     ! Calculate the mass of a phase M particle.
486              MASS_M = (PI/6.d0)*(DES_D_P0(M)**3)*DES_RO_S(M)
487     
488     ! Particle-Particle Collision Parameters ------------------------------>
489              DO L = M, DES_MMAX
490                 LC = LC+1
491     
492     ! Check particle-particle normal restitution coefficient
493                 IF(DES_EN_INPUT(LC) == UNDEFINED) THEN
494                    WRITE(ERR_MSG,1000) trim(iVar('DES_EN_INPUT',LC))
495                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
496                 ELSEIF(DES_EN_INPUT(LC) > ONE .OR.                         &
497                    DES_EN_INPUT(LC) < ZERO) THEN
498                    WRITE(ERR_MSG,1001) trim(iVar('DES_EN_INPUT',LC)),      &
499                       trim(iVal(DES_EN_INPUT(LC)))
500                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
501                 ENDIF
502                 EN = DES_EN_INPUT(LC)
503     
504     ! Calculate masses used for collision calculations.
505                 MASS_L = (PI/6.d0)*(DES_D_P0(L)**3)*DES_RO_S(L)
506                 MASS_EFF = MASS_M*MASS_L/(MASS_M+MASS_L)
507     
508     ! Calculate the M-L normal and tangential damping cofficients.
509                 IF(EN .NE. ZERO) THEN
510                    DES_ETAN(M,L) = 2.0D0*SQRT(KN*MASS_EFF) * ABS(LOG(EN))
511                    DES_ETAN(M,L) = DES_ETAN(M,L)/SQRT(PI*PI + (LOG(EN)**2))
512                 ELSE
513                    DES_ETAN(M,L) = 2.0D0*SQRT(KN*MASS_EFF)
514                 ENDIF
515                 DES_ETAT(M,L) = DES_ETAT_FAC*DES_ETAN(M,L)
516     
517     ! Store the entries in the symmetric matrix.
518                 DES_ETAN(L,M) = DES_ETAN(M,L)
519                 DES_ETAT(L,M) = DES_ETAT(M,L)
520     
521     ! Calculate the collision time scale.
522                 TCOLL_TMP = PI/SQRT(KN/MASS_EFF -                          &
523                    ((DES_ETAN(M,L)/MASS_EFF)**2)/4.d0)
524                 TCOLL = MIN(TCOLL_TMP, TCOLL)
525              ENDDO
526     
527     
528     ! Particle-Wall Collision Parameters ---------------------------------->
529     ! Check particle-wall normal restitution coefficient.
530              IF(DES_EN_WALL_INPUT(M) == UNDEFINED) THEN
531                 WRITE(ERR_MSG,1000) trim(iVar('DES_EN_WALL_INPUT',M))
532                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
533              ELSEIF(DES_EN_WALL_INPUT(M) > ONE .OR.                        &
534                 DES_EN_WALL_INPUT(M) < ZERO) THEN
535                 WRITE(ERR_MSG,1001) trim(iVar('DES_EN_WALL_INPUT',M)),     &
536                    trim(iVal(DES_EN_WALL_INPUT(M)))
537                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
538              ENDIF
539              EN = DES_EN_WALL_INPUT(M)
540     
541     ! Calculate masses used for collision calculations.
542              MASS_EFF = MASS_M
543     
544     ! Calculate the M-Wall normal and tangential damping coefficients.
545              IF(EN .NE. ZERO) THEN
546                 DES_ETAN_WALL(M) = 2.d0*SQRT(KN_W*MASS_EFF)*ABS(LOG(EN))
547                 DES_ETAN_WALL(M) = DES_ETAN_WALL(M)/SQRT(PI*PI+(LOG(EN))**2)
548              ELSE
549                 DES_ETAN_WALL(M) = 2.D0*SQRT(KN_W*MASS_EFF)
550              ENDIF
551              DES_ETAT_WALL(M) = DES_ETAT_W_FAC*DES_ETAN_WALL(M)
552     
553     ! Calculate the collision time scale.
554              TCOLL_TMP = PI/SQRT(KN_W/MASS_EFF -                           &
555                 ((DES_ETAN_WALL(M)/MASS_EFF)**2.d0)/4.d0)
556     !         TCOLL = MIN(TCOLL_TMP, TCOLL)
557           ENDDO
558     
559     ! if following are assigned warn user they are discarded
560           FLAG_WARN = .FALSE.
561           DO M = 1, DES_MMAX+DES_MMAX*(DES_MMAX-1)/2
562              IF(DES_ET_INPUT(M) .NE. UNDEFINED) FLAG_WARN = .TRUE.
563           ENDDO
564           IF (FLAG_WARN) THEN
565              WRITE(ERR_MSG,2102) 'DES_ET_INPUT'
566              CALL FLUSH_ERR_MSG
567           ENDIF
568     
569           FLAG_WARN = .FALSE.
570           DO M = 1, DES_MMAX
571              IF(DES_ET_WALL_INPUT(M) .NE. UNDEFINED) FLAG_WARN = .TRUE.
572           ENDDO
573           IF (FLAG_WARN)THEN
574              WRITE(ERR_MSG,2102) 'DES_ET_WALL_INPUT'
575              CALL FLUSH_ERR_MSG
576           ENDIF
577     
578      2102 FORMAT('Warning 2102: ',A,' values are not used ',/' with the',  &
579              ' linear spring-dashpot collision model.')
580     
581     ! Store the smalled calculated collision time scale. This value is used
582     ! in time-marching the DEM solids.
583           DTSOLID = TCOLL/50.d0
584     
585     
586           CALL FINL_ERR_MSG
587     
588           RETURN
589     
590      1000 FORMAT('Error 1000: Required input not specified: ',A,/'Please ',&
591              'correct the mfix.dat file.')
592     
593      1001 FORMAT('Error 1001: Illegal or unknown input: ',A,' = ',A,/      &
594              'Please correct the mfix.dat file.')
595     
596      1002 FORMAT('Error 1002: Required input not specified: ',A,/          &
597              'Description:',A,/'Please correct the mfix.dat file.')
598     
599           END SUBROUTINE CHECK_SOLIDS_DEM_COLL_LSD
600     
601     
602     
603     
604     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
605     !                                                                      !
606     !  Subroutine: CHECK_SOLIDS_DEM_COLL_HERTZ                             !
607     !  Author: J.Musser                                   Date: 11-Dec-13  !
608     !                                                                      !
609     !  Purpose: Check user input data for Hertzian collisions.             !
610     !                                                                      !
611     !  References:                                                         !
612     !   - Schafer et al., J. Phys. I France, 1996, 6, 5-20 (see page 7&13) !
613     !   -  Van der Hoef et al., Advances in Chemical Engineering, 2006, 31,!
614     !      65-149 (pages 94-95)                                            !
615     !   - Silbert et al., Physical Review E, 2001, 64, 051302 1-14 (page 5)!
616     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
617           SUBROUTINE CHECK_SOLIDS_DEM_COLL_HERTZ
618     
619     
620     ! Global Variables:
621     !---------------------------------------------------------------------//
622     ! Number of discrete solids phases, diameters and densities
623           USE discretelement, only: DES_MMAX, DES_D_P0, DES_RO_s
624     ! User defined coefficients of restitution: Normal and Tangential
625           USE discretelement, only: DES_EN_INPUT, DES_EN_WALL_INPUT
626           USE discretelement, only: DES_ET_INPUT, DES_ET_WALL_INPUT
627     ! Damping coefficients: Normal and Tangential
628           use discretelement, only: DES_ETAN, DES_ETAN_WALL
629           use discretelement, only: DES_ETAT, DES_ETAT_WALL
630     ! Hertzian spring constants: Normal and Tangential
631           use discretelement, only: HERT_KN, HERT_Kwn
632           use discretelement, only: HERT_KT, HERT_Kwt
633     ! Tangential damping factors := ET/EN
634           USE discretelement, only: DES_ETAT_FAC, DES_ETAT_W_FAC
635     ! Particle and wall Young's modulus and Shear modulus
636           USE discretelement, only: E_YOUNG, Ew_YOUNG, G_MOD
637     ! Particle and wall Poisson ratio
638           USE discretelement, only: V_POISSON, Vw_POISSON
639     ! Solids time step-size.
640           use discretelement, only: DTSOLID
641     ! Particle and wall normal spring constants
642           USE discretelement, only: KN, KN_W
643     ! Particle and wall tangential spring factor := KN/KT
644           USE discretelement, only: KT_FAC, KT_W_FAC
645     ! The constant PI
646           use constant, only: PI
647     
648     ! Parameter constatns.
649           USE param1, only: ZERO, HALF, ONE, UNDEFINED
650     
651     !      USE mpi_utility
652           use error_manager
653     
654           IMPLICIT NONE
655     
656     
657     ! Local Variables:
658     !---------------------------------------------------------------------//
659     ! Loop index.
660           INTEGER :: M, L, LC
661     ! Message for formatted output.
662           CHARACTER(len=64) :: MSG
663     ! Collision length scale.
664           DOUBLE PRECISION :: TCOLL, TCOLL_TMP
665     ! Particle and effective mass.
666           DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
667     ! Effective physical quantities. Radius, Youngs, Shear
668           DOUBLE PRECISION :: R_EFF, E_EFF, G_MOD_EFF, RED_MASS_EFF
669     ! Alias for coefficient restitution
670           DOUBLE PRECISION :: EN, ET
671     ! Shear modules for wall
672           DOUBLE PRECISION :: G_MOD_WALL
673     
674     ! Initialize the error manager.
675           CALL INIT_ERR_MSG("CHECK_SOLIDS_DEM_COLL_HERTZ")
676     
677     ! Initialize.
678           TCOLL = UNDEFINED
679     
680     ! check young's modulus and poisson ratio
681           IF(Ew_YOUNG == UNDEFINED ) THEN
682              MSG='Wall value for Youngs modulus'
683              WRITE(ERR_MSG,1002) 'Ew_YOUNG', MSG
684              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
685           ENDIF
686     
687           IF(Vw_POISSON == UNDEFINED) THEN
688              MSG='Wall value for Poissons ratio'
689              WRITE(ERR_MSG,1002) 'Vw_POISSON', MSG
690              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
691           ELSEIF (Vw_POISSON > 0.5d0 .OR. Vw_POISSON <= -ONE) THEN
692              WRITE(ERR_MSG,1001) 'Vw_POISSON',iVal(Vw_POISSON)
693              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
694           ENDIF
695     
696           G_MOD_WALL = 0.5d0*Ew_YOUNG/(1.d0+Vw_POISSON)
697     
698           LC = 0
699           DO M=1,DES_MMAX
700     
701              IF(E_YOUNG(M) == UNDEFINED) THEN
702                 MSG=''; WRITE(MSG,"('Phase ',I2,' Youngs modulus')") M
703                 WRITE(ERR_MSG,1002) 'E_YOUNG', MSG
704                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
705              ENDIF
706              IF(V_POISSON(M) == UNDEFINED) THEN
707                 MSG=''; WRITE(MSG,"('Phase ',I2,' Poissons ratio')") M
708                 WRITE(ERR_MSG,1002) 'V_POISSON', MSG
709                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
710              ELSEIF(V_POISSON(M) > 0.5d0 .OR.                              &
711                 V_POISSON(M) <= -ONE) THEN
712                 WRITE(ERR_MSG,1001) trim(iVar('V_POISSON',M)),             &
713                    iVal(V_POISSON(M))
714                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
715              ENDIF
716     ! Calculate the shear moduleus for phase M.
717              G_MOD(M) = 0.5d0*E_YOUNG(M)/(1.d0+V_POISSON(M))
718     
719     ! Calculate the mass of a phase M particle.
720              MASS_M = (PI/6.d0)*(DES_D_P0(M)**3)*DES_RO_S(M)
721     
722     ! Particle-Particle Collision Parameters ------------------------------>
723              DO L=M,DES_MMAX
724                 LC = LC+1
725     
726     ! Check particle-particle normal restitution coefficient
727                 IF(DES_EN_INPUT(LC) == UNDEFINED) THEN
728                    WRITE(ERR_MSG,1000) trim(iVar('DES_EN_INPUT',LC))
729                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
730     
731                 ELSEIF(DES_EN_INPUT(LC) > ONE .OR.                         &
732                    DES_EN_INPUT(LC) < ZERO) THEN
733                    WRITE(ERR_MSG,1001) trim(iVar('DES_EN_INPUT',LC)),      &
734                       trim(iVal(DES_EN_INPUT(LC)))
735                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
736                 ENDIF
737                 EN = DES_EN_INPUT(LC)
738     
739     ! Check particle-particle tangential restitution coefficient
740                 IF(DES_ET_INPUT(M) == UNDEFINED) THEN
741                    WRITE(ERR_MSG,1000) trim(iVar('DES_ET_INPUT',M))
742                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
743                 ELSEIF(DES_ET_INPUT(M) > ONE .OR.                          &
744                    DES_ET_INPUT(M) < ZERO) THEN
745                    WRITE(ERR_MSG,1001) trim(iVar('DES_ET_INPUT',M)),       &
746                       iVal(DES_ET_INPUT(M))
747                    CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
748                 ENDIF
749                 ET = DES_EN_INPUT(LC)
750     
751     
752     ! Calculate masses used for collision calculations.
753                 MASS_L = (PI/6.d0)*(DES_D_P0(L)**3)*DES_RO_S(L)
754                 MASS_EFF = (MASS_M*MASS_L)/(MASS_M+MASS_L)
755                 RED_MASS_EFF = (2.d0/7.d0)*MASS_EFF
756     ! Calculate the effective radius, Youngs modulus, and shear modulus.
757                 R_EFF = 0.5d0*(DES_D_P0(M)*DES_D_P0(L)/                    &
758                    (DES_D_P0(M) + DES_D_P0(L)))
759                 E_EFF = E_YOUNG(M)*E_YOUNG(L) /                            &
760                    (E_YOUNG(M)*(1.d0 - V_POISSON(L)**2) +                  &
761                     E_YOUNG(L)*(1.d0 - V_POISSON(M)**2))
762                 G_MOD_EFF = G_MOD(M)*G_MOD(L)/                             &
763                    (G_MOD(M)*(2.d0 - V_POISSON(L)) +                       &
764                     G_MOD(L)*(2.d0 - V_POISSON(M)))
765     
766     ! Calculate the spring properties and store in symmetric matrix format.
767                 HERT_KN(M,L)=( 4.d0/3.d0)*SQRT(R_EFF)*E_EFF
768                 HERT_KT(M,L)=(16.d0/3.d0)*SQRT(R_EFF)*G_MOD_EFF
769                 HERT_KN(L,M) = HERT_KN(M,L)
770                 HERT_KT(L,M) = HERT_KT(M,L)
771     
772     ! Calculate the normal coefficient.
773                 IF(EN .NE. ZERO) THEN
774                    DES_ETAN(M,L) = 2.d0*SQRT(HERT_KN(M,L)*MASS_EFF)*       &
775                       ABS(LOG(EN))
776                    DES_ETAN(M,L) = DES_ETAN(M,L)/                          &
777                       SQRT(PI*PI + (LOG(EN))**2)
778                 ELSE
779                    DES_ETAN(M,L) = 2.d0*SQRT(HERT_KN(M,L)*MASS_EFF)
780                 ENDIF
781                 DES_ETAN(L,M) = DES_ETAN(M,L)
782     
783     ! Calculate the tangential coefficients.
784                 IF(ET .NE. ZERO) THEN
785                    DES_ETAT(M,L) = 2.d0*SQRT(HERT_KT(M,L)*RED_MASS_EFF)*   &
786                       ABS(LOG(ET))
787                    DES_ETAT(M,L) = DES_ETAT(M,L)/ SQRT(PI*PI + (LOG(ET))**2)
788                 ELSE
789                    DES_ETAT(M,L) = 2.d0*SQRT(HERT_KT(M,L)*RED_MASS_EFF)
790                 ENDIF
791                 DES_ETAT(L,M) = DES_ETAT(M,L)
792     
793                 TCOLL_TMP = PI/SQRT(HERT_KN(M,L)/MASS_EFF -                &
794                    ((DES_ETAN(M,L)/MASS_EFF)**2)/4.d0)
795                 TCOLL = MIN(TCOLL_TMP, TCOLL)
796              ENDDO
797     
798     
799     
800     ! Particle-Wall Collision Parameters ---------------------------------->
801     ! Check particle-wall normal restitution coefficient.
802              IF(DES_EN_WALL_INPUT(M) == UNDEFINED) THEN
803                 WRITE(ERR_MSG,1000) trim(iVar('DES_EN_WALL_INPUT',M))
804                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
805              ELSEIF(DES_EN_WALL_INPUT(M) > ONE .OR.                        &
806                 DES_EN_WALL_INPUT(M) < ZERO) THEN
807                 WRITE(ERR_MSG,1001) trim(iVar('DES_EN_WALL_INPUT',M)),     &
808                    trim(iVal(DES_EN_WALL_INPUT(M)))
809                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
810              ENDIF
811              EN = DES_EN_WALL_INPUT(M)
812     
813     ! Check particle-wall tangential restitution coefficient
814              IF(DES_ET_WALL_INPUT(M) == UNDEFINED) THEN
815                 WRITE(ERR_MSG,1000) trim(iVar('DES_ET_WALL_INPUT',M))
816                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
817              ELSEIF(DES_ET_WALL_INPUT(M) > ONE .OR.                        &
818                 DES_ET_WALL_INPUT(M) < ZERO) THEN
819                 WRITE(ERR_MSG,1001) trim(iVar('DES_ET_WALL_INPUT',M)),     &
820                    trim(iVal(DES_ET_WALL_INPUT(M)))
821                 CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
822              ENDIF
823              ET = DES_ET_WALL_INPUT(M)
824     
825     
826     ! Calculate masses used for collision calculations.
827              MASS_EFF = MASS_M
828              RED_MASS_EFF = (2.d0/7.d0)*MASS_EFF
829     ! Calculate the effective radius, Youngs modulus, and shear modulus.
830              R_EFF = 0.5d0*DES_D_P0(M)
831              E_EFF = E_YOUNG(M)*Ew_YOUNG /                                 &
832                 (E_YOUNG(M)*(1.d0-Vw_POISSON**2) +                         &
833                  Ew_YOUNG  *(1.d0-V_POISSON(M)**2))
834              G_MOD_EFF = G_MOD(M)*G_MOD_WALL /                             &
835                 (G_MOD(M)*(2.d0 - Vw_POISSON) +                            &
836                  G_MOD_WALL*(2.d0 - V_POISSON(M)))
837     
838     ! Calculate the spring properties.
839              HERT_Kwn(M) = ( 4.d0/3.d0)*SQRT(R_EFF)*E_EFF
840              HERT_Kwt(M) = (16.d0/3.d0)*SQRT(R_EFF)*G_MOD_EFF
841     
842     ! Calculate the tangential coefficients.
843              DES_ETAN_WALL(M) = 2.d0*SQRT(HERT_Kwn(M)*MASS_EFF)*ABS(LOG(EN))
844              DES_ETAN_WALL(M) = DES_ETAN_WALL(M)/SQRT(PI*PI + (LOG(EN))**2)
845              DES_ETAT_WALL(M) = 2.d0*SQRT(HERT_Kwt(M)*RED_MASS_EFF)*       &
846                  ABS(LOG(ET))
847              DES_ETAT_WALL(M) = DES_ETAT_WALL(M)/SQRT(PI*PI + (LOG(ET))**2)
848     ! Calculate the collision time scale.
849              TCOLL_TMP = PI/SQRT(HERT_Kwn(M)/MASS_EFF -                    &
850                 ((DES_ETAN_WALL(M)/MASS_EFF)**2)/4.d0)
851           ENDDO
852     
853     
854     ! If following are assigned warn user they are discarded.
855            IF(KN .NE. UNDEFINED) THEN
856               WRITE(ERR_MSG, 2200) 'KN'
857               CALL FLUSH_ERR_MSG
858            ENDIF
859            IF(KN_W .NE. UNDEFINED) THEN
860               WRITE(ERR_MSG, 2200) 'KN_W'
861               CALL FLUSH_ERR_MSG
862            ENDIF
863            IF(KT_FAC .NE. UNDEFINED) THEN
864               WRITE(ERR_MSG, 2200) 'KT_FAC'
865               CALL FLUSH_ERR_MSG
866            ENDIF
867            IF(KT_W_FAC .NE. UNDEFINED) THEN
868               WRITE(ERR_MSG, 2200) 'KT_W_FAC'
869               CALL FLUSH_ERR_MSG
870            ENDIF
871            IF(DES_ETAT_FAC .NE. UNDEFINED) THEN
872               WRITE(ERR_MSG, 2200) 'DES_ETAT_FAC'
873               CALL FLUSH_ERR_MSG
874            ENDIF
875            IF(DES_ETAT_W_FAC .NE. UNDEFINED) THEN
876               WRITE(ERR_MSG, 2200) 'DES_ETAT_W_FAC'
877               CALL FLUSH_ERR_MSG
878            ENDIF
879     
880      2200 FORMAT('Warning 2200: ',A,' values are not used ',/' with the',  &
881              ' linear spring-dashpot collision model.')
882     
883     
884     ! Store the smalled calculated collision time scale. This value is used
885     ! in time-marching the DEM solids.
886           DTSOLID = TCOLL/50.d0
887     
888     
889           CALL FINL_ERR_MSG
890     
891           RETURN
892     
893      1000 FORMAT('Error 1000: Required input not specified: ',A,/'Please ',&
894              'correct the mfix.dat file.')
895     
896      1001 FORMAT('Error 1001: Illegal or unknown input: ',A,' = ',A,/      &
897              'Please correct the mfix.dat file.')
898     
899      1002 FORMAT('Error 1002: Required input not specified: ',A,/          &
900              'Description:',A,/'Please correct the mfix.dat file.')
901     
902           END SUBROUTINE CHECK_SOLIDS_DEM_COLL_HERTZ
903     
904     
905     
906     
907     
908     
909     
910     
911     
912     
913