File: N:\mfix\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     ! 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.
21           CALL CHECK_SOLIDS_DEM_COLLISION
22     ! DES cohesion model parameters.
23           CALL CHECK_SOLIDS_DEM_COHESION
24     ! Particle-particle conduction model parameters.
25           CALL CHECK_SOLIDS_DEM_ENERGY
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
54           use des_thermo_cond, only: DO_AREA_CORRECTION
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
225     
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     !<------------------- Van der Waals model messages. ----------------->!
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
378     
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
423             CALL CHECK_SOLIDS_DEM_COLL_LSD
424     ! Hertzian collision model.
425           CASE('HERTZIAN')
426              DES_COLL_MODEL_ENUM = HERTZIAN
427              CALL CHECK_SOLIDS_DEM_COLL_HERTZ
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
464     
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
733     
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     
1065