66 CHARACTER(len=64) :: MSG
67 INTEGER :: lent, lend, lenc, LC
69 LOGICAL :: ANY_CONDUCTION = .false.
71 DOUBLE PRECISION :: MASS_L, MASS_M, MASS_EFF
72 DOUBLE PRECISION :: R_EFF, E_EFF
82 IF(des_min_cond_dist == undefined)
THEN 83 des_min_cond_dist = 1.0d-04
84 IF (units ==
'SI') des_min_cond_dist = &
85 des_min_cond_dist/100.0
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
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
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
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))
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
132 IF(vw_poisson_actual == undefined)
THEN 133 msg=
'';
WRITE(msg,
"(' Actual poisson ratio (wall)')")
134 WRITE(
err_msg,2002)
'VW_POISSON_ACTUAL', msg
155 mass_eff = (mass_m*mass_l)/(mass_m+mass_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))
165 hert_kn_actual(m,l)=(4.d0/3.d0)*sqrt(r_eff)*e_eff
168 tau_c_base_actual(m,l)=3.21d0*(mass_eff/hert_kn_actual(m,l))
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
176 tau_c_base_sim(m,l)=
pi/sqrt(kn/mass_eff - &
177 ((des_etan(m,l)/mass_eff)**2)/4.d0)
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))
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))*
191 IF (des_coll_model_enum .EQ. hertzian)
THEN 192 tauw_c_base_sim(m)=3.21d0*(mass_eff/hert_kwn(m))**0.4
194 tauw_c_base_sim(m)=
pi/sqrt(kn_w/mass_eff - &
195 ((des_etan_wall(m)/mass_eff)**2)/4.d0)
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.')
229 use discretelement
, only: use_cohesion
231 use discretelement
, only: max_radius
233 use discretelement
, only: square_well
235 use discretelement
, only: van_der_waals
237 use discretelement
, only: factor_rlm
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
260 DOUBLE PRECISION :: VDW_NEIGHBORHOOD
265 IF(.NOT.use_cohesion)
THEN 267 square_well = .false.
268 van_der_waals = .false.
269 wall_vdw_outer_cutoff =
zero 279 IF (square_well .AND. van_der_waals)
THEN 284 ELSEIF(.NOT.square_well .AND. .NOT.van_der_waals)
THEN 291 IF (van_der_waals)
THEN 293 IF (vdw_inner_cutoff .EQ. undefined)
THEN 294 WRITE(
err_msg,1201)
'VDW_INNER_CUTOFF' 298 IF(vdw_outer_cutoff .EQ. undefined)
THEN 299 WRITE(
err_msg,1201)
'VDW_OUTER_CUTOFF' 303 IF(hamaker_constant .EQ. undefined)
THEN 304 WRITE(
err_msg,1201)
'HAMAKER_CONSTANT' 308 IF (wall_vdw_inner_cutoff .EQ. undefined)
THEN 309 WRITE(
err_msg,1201)
'WALL_VDW_INNER_CUTOFF' 313 IF (wall_vdw_outer_cutoff .EQ. undefined)
THEN 314 WRITE(
err_msg,1201)
'WALL_VDW_OUTER_CUTOFF' 318 IF(wall_hamaker_constant .EQ. undefined)
THEN 319 WRITE(
err_msg,1201)
'WALL_HAMAKER_CONSTANT' 323 vdw_neighborhood = 1.0d0 + (vdw_outer_cutoff/(2.d0*max_radius))
324 IF (factor_rlm < vdw_neighborhood)
THEN 329 IF (asperities <
zero)
THEN 330 WRITE(
err_msg,1001)
'ASPERITIES', trim(
ival(asperities))
334 surface_energy=hamaker_constant/&
335 (24.d0*
pi*vdw_inner_cutoff**2)
337 wall_surface_energy=wall_hamaker_constant/&
338 (24.d0*
pi*wall_vdw_inner_cutoff**2)
347 1001
FORMAT(
'Error 1001: Illegal or unknown input: ',a,
' = ',a,/ &
348 'Please correct the mfix.dat file.')
350 1002
FORMAT(
'Error 1000: Cannot use SQUARE_WELL and VAN_DER_WAALS ', &
351 'cohesion',/
'models simultaneously.')
353 1003
FORMAT(
'Error 1001: A cohesion model was not selected. Specify ',&
354 'one of the available models in the mfix.dat file.')
359 1201
FORMAT(
'Error 1201: Missing input data for Van der Waals ', &
360 'cohesion model.',/
'Input parameter ',a,
' is UNDEFINED.')
362 1202
FORMAT(
'Error 1202: VDW_OUTER_CUTOFF outside of the neighbor ', &
363 'search distance.',/
'Increase FACTOR_RLM to increase the ', &
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
389 USE discretelement
, only: mew, mew_w
405 ELSEIF (mew <
zero .OR. mew_w >
one)
THEN 413 ELSEIF(mew_w <
zero .OR. mew_w >
one)
THEN 419 SELECT CASE (trim(des_coll_model))
422 des_coll_model_enum = lsd
426 des_coll_model_enum = hertzian
430 WRITE(
err_msg,2000) trim(des_coll_model)
434 2000
FORMAT(
'Error 2000: Invalid particle-particle collision model:',&
435 a,/
'Please correct the mfix.dat file.')
442 1000
FORMAT(
'Error 1000: Required input not specified: ',a,/
'Please ',&
443 'correct the mfix.dat file.')
445 1001
FORMAT(
'Error 1001: Illegal or unknown input: ',a,
' = ',a,/ &
446 'Please correct the mfix.dat file.')
469 USE discretelement
, only: des_mmax
471 USE discretelement
, only: kn, kn_w
472 USE discretelement
, only: kt, kt_w
474 USE discretelement
, only: kt_fac, kt_w_fac
476 use discretelement
, only: des_etan, des_etan_wall
477 use discretelement
, only: des_etat, des_etat_wall
480 USE discretelement
, only: des_en_input, des_en_wall_input
481 USE discretelement
, only: des_et_input, des_et_wall_input
484 USE discretelement
, only: des_etat_fac, des_etat_w_fac
486 use discretelement
, only: dtsolid
500 INTEGER :: M, L, LC, MMAX_TOT
503 INTEGER :: lent, lend, lenc
507 DOUBLE PRECISION :: TCOLL, TCOLL_TMP
509 DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
511 DOUBLE PRECISION :: EN
531 ELSEIF(kt_fac >
one .OR. kt_fac <
zero)
THEN 546 WRITE (
err_msg, 2100)
'KT_W_FAC' 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))
556 2100
FORMAT(
'Warning 2100: Tangential spring factor ',a,
' not ', &
557 'specified in mfix.dat.',/
'Setting to default: (2/7).')
561 WRITE (
err_msg, 2101)
'DES_ETAT_FAC' 564 ELSEIF(des_etat_fac >
one .OR. des_etat_fac <
zero)
THEN 565 WRITE(
err_msg,1001)
'DES_ETAT_FAC',
ival(des_etat_fac)
571 WRITE (
err_msg, 2101)
'DES_ETAT_W_FAC' 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)
579 2101
FORMAT(
'Warning 2101: Tangential damping factor ',a,
' not ', &
580 'specified',/
'in mfix.dat. Setting to default: (1/2).')
586 mmax_tot = des_mmax+
mmax 587 des_en_wall_input((
mmax+1):mmax_tot) = des_en_wall_input(1:des_mmax
596 DO m =
mmax+1, mmax_tot
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)))
617 en = des_en_input(lc)
621 mass_eff = mass_m*mass_l/(mass_m+mass_l)
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))
628 des_etan(m,l) = 2.0d0*sqrt(kn*mass_eff)
630 des_etat(m,l) = des_etat_fac*des_etan(m,l)
633 des_etan(l,m) = des_etan(m,l)
634 des_etat(l,m) = des_etat(m,l)
637 tcoll_tmp =
pi/sqrt(kn/mass_eff - &
638 ((des_etan(m,l)/mass_eff)**2)/4.d0)
639 tcoll = min(tcoll_tmp, tcoll)
645 IF(des_en_wall_input(m) ==
undefined)
THEN 646 WRITE(
err_msg,1000) trim(
ivar(
'DES_EN_WALL_INPUT',m))
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)))
654 en = des_en_wall_input(m)
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)
664 des_etan_wall(m) = 2.d0*sqrt(kn_w*mass_eff)
666 des_etat_wall(m) = des_etat_w_fac*des_etan_wall(m)
669 tcoll_tmp =
pi/sqrt(kn_w/mass_eff - &
670 ((des_etan_wall(m)/mass_eff)**2.d0)/4.d0)
677 DO m =
mmax+1, mmax_tot
681 IF(des_et_input(m) .NE.
undefined) flag_warn = .true.
685 WRITE(
err_msg,2102)
'DES_ET_INPUT' 690 DO m =
mmax+1, mmax_tot
692 IF(des_et_wall_input(m) .NE.
undefined) flag_warn = .true.
695 WRITE(
err_msg,2102)
'DES_ET_WALL_INPUT' 699 2102
FORMAT(
'Warning 2102: ',a,
' values are not used ',/
' with the', &
700 ' linear spring-dashpot collision model.')
704 dtsolid = tcoll/50.d0
711 1000
FORMAT(
'Error 1000: Required input not specified: ',a,/
'Please ',&
712 'correct the mfix.dat file.')
714 1001
FORMAT(
'Error 1001: Illegal or unknown input: ',a,
' = ',a,/ &
715 'Please correct the mfix.dat file.')
739 USE discretelement
, only: des_mmax
741 USE discretelement
, only: des_en_input, des_en_wall_input
742 USE discretelement
, only: des_et_input, des_et_wall_input
744 use discretelement
, only: des_etan, des_etan_wall
745 use discretelement
, only: des_etat, des_etat_wall
747 use discretelement
, only: hert_kn, hert_kwn
748 use discretelement
, only: hert_kt, hert_kwt
750 USE discretelement
, only: des_etat_fac, des_etat_w_fac
752 USE discretelement
, only: e_young, ew_young
754 USE discretelement
, only: v_poisson, vw_poisson
756 use discretelement
, only: dtsolid
758 USE discretelement
, only: kn, kn_w
760 USE discretelement
, only: kt_fac, kt_w_fac
776 INTEGER :: M, L, LC, MMAX_TOT
779 INTEGER :: lent, lend, lenc
781 CHARACTER(len=64) :: MSG
783 DOUBLE PRECISION :: TCOLL, TCOLL_TMP
785 DOUBLE PRECISION :: MASS_M, MASS_L, MASS_EFF
787 DOUBLE PRECISION :: R_EFF, E_EFF, G_MOD_EFF, RED_MASS_EFF
789 DOUBLE PRECISION :: EN, ET
791 DOUBLE PRECISION :: G_MOD(
dim_m)
793 DOUBLE PRECISION :: G_MOD_WALL
804 msg=
'Wall value for Youngs modulus' 805 WRITE(
err_msg,1002)
'Ew_YOUNG', msg
810 msg=
'Wall value for Poissons ratio' 811 WRITE(
err_msg,1002)
'Vw_POISSON', msg
813 ELSEIF (vw_poisson > 0.5d0 .OR. vw_poisson <= -
one)
THEN 818 g_mod_wall = 0.5d0*ew_young/(1.d0+vw_poisson)
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
839 msg=
'';
WRITE(msg,
"('Phase ',I2,' Youngs modulus')") m
840 WRITE(
err_msg,1002)
'E_YOUNG', msg
844 msg=
'';
WRITE(msg,
"('Phase ',I2,' Poissons ratio')") m
845 WRITE(
err_msg,1002)
'V_POISSON', msg
847 ELSEIF(v_poisson(m) > 0.5d0 .OR. &
848 v_poisson(m) <= -
one)
THEN 854 g_mod(m) = 0.5d0*e_young(m)/(1.d0+v_poisson(m))
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)))
880 en = des_en_input(lc)
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))
892 et = des_et_input(lc)
897 mass_eff = (mass_m*mass_l)/(mass_m+mass_l)
898 red_mass_eff = (2.d0/7.d0)*mass_eff
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)))
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
912 hert_kn(l,m) = hert_kn(m,l)
913 hert_kt(l,m) = hert_kt(m,l)
916 IF(en .NE.
zero)
THEN 917 des_etan(m,l) = 2.d0*sqrt(hert_kn(m,l)*mass_eff)* &
919 des_etan(m,l) = des_etan(m,l)/ &
920 sqrt(
pi*
pi + (log(en))**2)
922 des_etan(m,l) = 2.d0*sqrt(hert_kn(m,l)*mass_eff)
924 des_etan(l,m) = des_etan(m,l)
927 IF(et .NE.
zero)
THEN 928 des_etat(m,l) = 2.d0*sqrt(hert_kt(m,l)*red_mass_eff)* &
930 des_etat(m,l) = des_etat(m,l)/ sqrt(
pi*
pi+(log(et))**2)
932 des_etat(m,l) = 2.d0*sqrt(hert_kt(m,l)*red_mass_eff)
934 des_etat(l,m) = des_etat(m,l)
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)
943 IF(des_en_wall_input(m) ==
undefined)
THEN 944 WRITE(
err_msg,1000) trim(
ivar(
'DES_EN_WALL_INPUT',m))
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)))
952 en = des_en_wall_input(m)
955 IF(des_et_wall_input(m) ==
undefined)
THEN 956 WRITE(
err_msg,1000) trim(
ivar(
'DES_ET_WALL_INPUT',m))
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)))
964 et = des_et_wall_input(m)
968 red_mass_eff = (2.d0/7.d0)*mass_eff
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)))
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
984 des_etan_wall(m) = 2.d0*sqrt(hert_kwn(m)*mass_eff)*&
986 des_etan_wall(m) = des_etan_wall(m)/&
987 sqrt(
pi*
pi + (log(en))**2)
989 des_etan_wall(m) = 2.d0*sqrt(hert_kwn(m)*mass_eff)
993 des_etat_wall(m) = 2.d0*sqrt(hert_kwt(m)*red_mass_eff)* &
995 des_etat_wall(m) = des_etat_wall(m)/sqrt(
pi*
pi+(log(et))**2)
997 des_etat_wall(m) = 2.d0*sqrt(hert_kwt(m)*red_mass_eff)
1001 tcoll_tmp =
pi/sqrt(hert_kwn(m)/mass_eff - &
1002 ((des_etan_wall(m)/mass_eff)**2)/4.d0)
1020 WRITE(
err_msg, 2200)
'KT_W_FAC' 1024 WRITE(
err_msg, 2200)
'DES_ETAT_FAC' 1028 WRITE(
err_msg, 2200)
'DES_ETAT_W_FAC' 1032 2200
FORMAT(
'Warning 2200: ',a,
' values are not used ',/
' with the', &
1033 ' linear spring-dashpot collision model.')
1038 dtsolid = tcoll/50.d0
1045 1000
FORMAT(
'Error 1000: Required input not specified: ',a,/
'Please ',&
1046 'correct the mfix.dat file.')
1048 1001
FORMAT(
'Error 1001: Illegal or unknown input: ',a,
' = ',a,/ &
1049 'Please correct the mfix.dat file.')
1051 1002
FORMAT(
'Error 1002: Required input not specified: ',a,/ &
1052 'Description:',a,/
'Please correct the mfix.dat file.')
subroutine check_solids_dem
double precision, dimension(dim_m) d_p0
character(len=32) function ivar(VAR, i1, i2, i3)
double precision, parameter one
logical, dimension(dim_m) calc_cond_des
double precision des_min_cond_dist
character(len=3), dimension(dim_m) solids_model
double precision, parameter undefined
subroutine init_err_msg(CALLER)
logical do_area_correction
subroutine check_solids_dem_collision
subroutine check_solids_dem_energy
double precision, parameter half
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
double precision, parameter pi
double precision, parameter zero
subroutine flush_err_msg(DEBUG, HEADER, FOOTER, ABORT, LOG, CALL_TREE)