MFIX  2016-1
drag_gs.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Subroutine: DRAG_gs C
4 ! Purpose: Calculate the gas-solids drag coefficient C
5 ! C
6 ! Author: M. Syamlal Date: 20-MAY-96 C
7 ! Reviewer: Date: C
8 ! C
9 ! Literature/Document References: C
10 ! C
11 ! Comments: C
12 ! If changes are made to this file, especially in regard to the C
13 ! structure of any drag subroutine call, then such changes must C
14 ! also be made to the analgous routine des_drag_gp in the file C
15 ! des/drag_fgs.f for consistency! C
16 ! C
17 ! C
18 ! C
19 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20 
21  SUBROUTINE drag_gs(M, IER)
22 
23 !-----------------------------------------------
24 ! Modules
25 !-----------------------------------------------
26  USE compar
27  USE constant
28  USE discretelement
29  USE drag
30  USE fldvar
31  USE fun_avg
32  USE functions
33  USE funits
34  USE geometry
35  USE indices
36  USE machine, only: start_log, end_log
37  USE mms
38  USE parallel
39  USE param
40  USE param1
41  USE physprop
42  USE run, only: syam_obrien, gidaspow, gidaspow_pcf, gidaspow_blend, gidaspow_blend_pcf, wen_yu, wen_yu_pcf, koch_hill, koch_hill_pcf, bvk, user_drag, hys, drag_type, drag_type_enum, ghd_2007, igci, kt_type_enum, milioli, model_b, subgrid_type_enum, undefined_subgrid_type
43  USE sendrecv
44  USE ur_facs
45  IMPLICIT NONE
46 !-----------------------------------------------
47 ! Dummy arguments
48 !-----------------------------------------------
49 ! Solids phase index
50  INTEGER, INTENT(IN) :: M
51 ! Error index
52  INTEGER, INTENT(INOUT) :: IER
53 !-----------------------------------------------
54 ! Local parameters
55 !-----------------------------------------------
56 
57 !-----------------------------------------------
58 ! Local variables
59 !-----------------------------------------------
60 ! indices
61  INTEGER :: I, IJK, IMJK, IJMK, IJKM
62 ! cell center value of U_sm, V_sm, W_sm
63  DOUBLE PRECISION :: USCM, VSCM, WSCM
64 ! cell center value of x, y and z-particle velocity
65  DOUBLE PRECISION :: USCM_HYS, VSCM_HYS, WSCM_HYS
66 ! cell center value of U_g, V_g, W_g
67  DOUBLE PRECISION :: UGC, VGC, WGC
68 ! magnitude of gas-solids relative velocity
69  DOUBLE PRECISION :: VREL
70 ! gas laminar viscosity redefined here to set viscosity at pressure
71 ! boundaries
72  DOUBLE PRECISION :: Mu
73 ! drag coefficient
74  DOUBLE PRECISION :: DgA
75 ! current value of F_gs (i.e., without underrelaxation)
76  DOUBLE PRECISION F_gstmp
77 ! indices of solids phases (continuous, discrete)
78  INTEGER :: CM, DM, L
79 ! temporary shift of total number of solids phases to account for both
80 ! discrete and continuous solids phases used for the hybrid mdoel
81  INTEGER :: MAXM
82 ! tmp local variable for the particle diameter of solids
83 ! phase M (continuous or discrete)
84  DOUBLE PRECISION :: DP_loc(dim_m)
85 ! tmp local variable for the solids volume fraction of solids
86 ! phase M (continuous or discrete)
87  DOUBLE PRECISION :: EPs_loc(dim_m)
88 ! tmp local variable for the particle density of solids
89 ! phase M (continuous or discrete)
90  DOUBLE PRECISION :: ROs_loc(dim_m)
91 ! correction factors for implementing polydisperse drag model
92 ! proposed by van der Hoef et al. (2005)
93  DOUBLE PRECISION :: F_cor, tmp_sum, tmp_fac
94 ! average particle diameter in polydisperse systems
95  DOUBLE PRECISION :: DPA
96 ! diameter ratio in polydisperse systems
97  DOUBLE PRECISION :: Y_i
98 ! total solids volume fraction
99  DOUBLE PRECISION :: phis
100 ! temporary local variables to use for dummy arguments in subroutines
101 ! void fraction, gas density, gas bulk density, solids volume fraction
102 ! particle diameter, particle density
103  DOUBLE PRECISION :: EPG, ROg, ROPg, EP_SM, DPM, ROs
104 !-----------------------------------------------
105 
106 !!$omp parallel do default(shared) &
107 !!$omp private( I, IJK, IMJK, IJMK, IJKM, DM, MAXM, CM, L, &
108 !!$omp UGC, VGC, WGC, USCM, VSCM, WSCM, VREL, USCM_HYS, &
109 !!$omp VSCM_HYS, WSCM_HYS, tmp_sum, tmp_fac, Y_i, F_cor, &
110 !!$omp EP_SM, EPs_loc, ROs_loc, DP_loc, DPA, DPM, ROs, &
111 !!$omp phis, EPg, ROg, ROPg, Mu, DgA, F_gstmp)
112 
113 
114  DO ijk = ijkstart3, ijkend3
115 
116  IF (fluidorp_flow_at(ijk)) THEN
117 
118  i = i_of(ijk)
119  imjk = im_of(ijk)
120  ijmk = jm_of(ijk)
121  ijkm = km_of(ijk)
122 
123  DO l = 1,des_mmax+mmax
124  IF(kt_type_enum == ghd_2007 .AND. m==mmax) THEN
125 ! set this to avoid issues later in routine
126  dp_loc(l) = one
127  eps_loc(l) = zero
128  ros_loc(l) = zero
129  ENDIF
130  dp_loc(l) = d_p(ijk,l)
131  eps_loc(l) = ep_s(ijk,l)
132  ros_loc(l) = ro_s(ijk,l)
133  ENDDO
134 
135 ! Calculate velocity components at i, j, k
136  ugc = avg_x_e(u_g(imjk),u_g(ijk),i)
137  vgc = avg_y_n(v_g(ijmk),v_g(ijk))
138  wgc = avg_z_t(w_g(ijkm),w_g(ijk))
139 
140  uscm = avg_x_e(u_s(imjk,m),u_s(ijk,m),i)
141  vscm = avg_y_n(v_s(ijmk,m),v_s(ijk,m))
142  wscm = avg_z_t(w_s(ijkm,m),w_s(ijk,m))
143 
144 ! magnitude of gas-solids relative velocity
145  vrel = sqrt((ugc-uscm)**2 + (vgc-vscm)**2 + (wgc-wscm)**2)
146 
147 ! Laminar viscosity at a pressure boundary is given the value of the
148 ! fluid cell next to it. This applies just to the calculation of the
149 ! drag, in other routines the value of viscosity at a pressure boundary
150 ! always has a zero value.
151  IF (p_flow_at(ijk)) THEN
152  IF( fluid_at(east_of(ijk)) ) THEN
153  mu = mu_g(east_of(ijk))
154  ELSEIF ( fluid_at(west_of(ijk)) ) THEN
155  mu = mu_g(west_of(ijk))
156  ELSEIF ( fluid_at(north_of(ijk)) ) THEN
157  mu = mu_g(north_of(ijk))
158  ELSEIF ( fluid_at(south_of(ijk)) ) THEN
159  mu = mu_g(south_of(ijk))
160  ELSEIF ( fluid_at(top_of(ijk)) ) THEN
161  mu = mu_g(top_of(ijk))
162  ELSEIF ( fluid_at(bottom_of(ijk)) ) THEN
163  mu = mu_g(bottom_of(ijk))
164  ENDIF
165  ELSE
166  mu = mu_g(ijk)
167  ENDIF
168 
169 ! calculate the total solids volume fraction
170  phis = zero
171  DO l = 1, des_mmax+mmax
172 ! this is slightly /= one-ep_g due to round-off
173  phis = phis + eps_loc(l)
174  ENDDO
175 
176 ! calculate the average particle diameter and particle ratio
177  dpa = zero
178  tmp_sum = zero
179  tmp_fac = zero
180  DO l = 1, des_mmax+mmax
181  IF (phis .GT. zero) THEN
182  tmp_fac = eps_loc(l)/phis
183  tmp_sum = tmp_sum + tmp_fac/dp_loc(l)
184  ELSE
185  tmp_sum = tmp_sum + one/dp_loc(l) ! not important, but will avoid NaN's in empty cells
186  ENDIF
187  ENDDO
188  dpa = one / tmp_sum
189  y_i = dp_loc(m)/dpa
190 
191 ! assign aliases for easy reference
192  epg = ep_g(ijk)
193  rog = ro_g(ijk)
194  ropg = rop_g(ijk)
195  ep_sm = eps_loc(m)
196  dpm = dp_loc(m)
197  ros = ros_loc(m)
198 
199 
200 ! determine the drag coefficient
201  IF (ep_sm <= zero) THEN
202  dga = zero
203  ELSEIF (epg == zero) THEN
204 ! this case will already be caught in most drag subroutines whenever
205 ! RE==0 (for correlations in which RE includes EPg). however, this will
206 ! prevent potential divisions by zero in some models by setting it now.
207  dga = zero
208 ! Force a ZERO drag coefficient for MMS cases.
209  ELSEIF (use_mms) THEN
210  dga = zero
211  ELSE
212  SELECT CASE(drag_type_enum)
213 
214  CASE (syam_obrien)
215  CALL drag_syam_obrien(dga,epg,mu,rog,vrel,dpm)
216 
217  CASE (gidaspow)
218  CALL drag_gidaspow(dga,epg,mu,rog,ropg,vrel,dpm)
219 
220  CASE (gidaspow_pcf)
221  CALL drag_gidaspow(dga,epg,mu,rog,ropg,vrel,dpa)
222 
223  CASE (gidaspow_blend)
224  CALL drag_gidaspow_blend(dga,epg,mu,rog,ropg,vrel,dpm)
225 
226  CASE (gidaspow_blend_pcf)
227  CALL drag_gidaspow_blend(dga,epg,mu,rog,ropg,vrel,dpa)
228 
229  CASE (wen_yu)
230  CALL drag_wen_yu(dga,epg,mu,ropg,vrel,dpm)
231 
232  CASE (wen_yu_pcf)
233  CALL drag_wen_yu(dga,epg,mu,ropg,vrel,dpa)
234 
235  CASE (koch_hill)
236  CALL drag_koch_hill(dga,epg,mu,ropg,vrel,dpm,dpm,phis)
237 
238  CASE (koch_hill_pcf)
239  CALL drag_koch_hill(dga,epg,mu,ropg,vrel,dpm,dpa,phis)
240 
241  CASE (bvk)
242  CALL drag_bvk(dga,epg,mu,ropg,vrel,dpm,dpa,phis)
243 
244  CASE (user_drag)
245  CALL drag_usr(ijk, m, dga, epg, mu, rog, vrel, dpm, &
246  ros, ugc, vgc, wgc)
247 
248  CASE (hys)
249 ! only over the continuous two fluid phases
250  maxm = smax
251 ! calculate velocity components of each solids phase
252  uscm_hys = zero
253  vscm_hys = zero
254  wscm_hys = zero
255  IF(phis > zero) THEN
256  DO l = 1, maxm
257  uscm_hys = uscm_hys + eps_loc(l)*(ugc - &
258  avg_x_e(u_s(imjk,l),u_s(ijk,l),i))
259  vscm_hys = vscm_hys + eps_loc(l)*(vgc - &
260  avg_y_n(v_s(ijmk,l),v_s(ijk,l)))
261  wscm_hys = wscm_hys + eps_loc(l)*(wgc - &
262  avg_z_t(w_s(ijkm,l),w_s(ijk,l)))
263  ENDDO
264  uscm_hys = uscm_hys/phis
265  vscm_hys = vscm_hys/phis
266  wscm_hys = wscm_hys/phis
267  ENDIF
268 ! magnitude of gas-solids relative velocity
269  vrel = sqrt(uscm_hys**2 +vscm_hys**2 +wscm_hys**2)
270 
271  CALL drag_hys(dga,epg,mu,ropg,vrel,&
272  dp_loc(:),dpa,y_i,eps_loc(:),phis,m,maxm,ijk)
273 
274  CASE DEFAULT
275  CALL start_log
276  IF(.NOT.dmp_log) call open_pe_log(ier)
277  IF(dmp_log) WRITE (*, '(A,A)') &
278  'Unknown DRAG_TYPE: ', drag_type
279  WRITE (unit_log, '(A,A)') 'Unknown DRAG_TYPE: ', drag_type
280  CALL end_log
281  CALL mfix_exit(mype)
282  END SELECT ! end selection of drag_type
283  ENDIF ! end if/elseif/else (ep_sm <= zero, ep_g==0)
284 
285 ! modify the drag coefficient to account for subgrid domain effects
286  IF (subgrid_type_enum .ne. undefined_subgrid_type) THEN
287  IF (subgrid_type_enum .EQ. igci) THEN
288  CALL subgrid_drag_igci(dga,epg,mu,rog,dpm,ros,ijk)
289  ELSEIF (subgrid_type_enum .EQ. milioli) THEN
290  CALL subgrid_drag_milioli(dga,epg,mu,rog,vrel,&
291  dpm,ros,ijk)
292  ENDIF
293  ENDIF
294 
295 
296 ! Modify drag coefficient to account for possible corrections and
297 ! for differences between Model B and Model A
298  IF(drag_type_enum == hys) THEN
299 ! this drag model is handled differently than the others
300  IF(model_b)THEN
301  f_gstmp = dga/epg
302  ELSE
303  f_gstmp = dga
304  ENDIF
305  ELSE
306  IF(drag_type_enum == gidaspow_pcf .OR. &
307  drag_type_enum == gidaspow_blend_pcf .OR. &
308  drag_type_enum == wen_yu_pcf .OR. &
309  drag_type_enum == koch_hill_pcf .OR. &
310  drag_type_enum == bvk) THEN
311 ! see erratum by Beetstra et al. (2007) : the correction factor differs
312 ! for model A versus model B.
313 ! application of the correction factor for model A is found from
314 ! the correction factor for model B and neglects the Y_i**3 term
315  IF(model_b) THEN
316  IF (m == 1) THEN
317  f_cor = (epg*y_i + phis*y_i**2)
318  ELSE
319  f_cor = (epg*y_i + phis*y_i**2 + &
320  0.064d0*epg*y_i**3)
321  ENDIF
322  ELSE
323  f_cor = y_i
324  ENDIF
325  dga = one/(y_i*y_i) * dga * f_cor
326  ENDIF
327 
328 ! Calculate the drag coefficient (Model B coeff = Model A coeff/EP_g); all other models, eg., Wen_Yu
329  IF(model_b)THEN
330  f_gstmp = dga * ep_sm/epg
331  ELSE
332  f_gstmp = dga * ep_sm !3D buoyancy model
333  ENDIF
334  ENDIF !end if/else trim(drag_type=='hys')
335 
336 ! Determine drag force coefficient accounting for any under relaxation
337  f_gs(ijk,m) = (one - ur_f_gs)*f_gs(ijk, m) + &
338  ur_f_gs*f_gstmp
339 
340  IF(kt_type_enum == ghd_2007) THEN
341  IF(m==1) THEN
342  f_gs(ijk,mmax) = f_gs(ijk,m)
343  ELSE
344  f_gs(ijk,mmax) = f_gs(ijk,mmax) + f_gs(ijk,m)
345  ENDIF
346  ENDIF
347 
348  ELSE ! .not.(fluidorp_flow_at(ijk)) branch
349 
350  f_gs(ijk,m) = zero
351  IF(kt_type_enum == ghd_2007) f_gs(ijk, mmax) = zero
352 
353  ENDIF ! end if (fluidorp_flow_at(ijk))
354 
355  ENDDO ! end do (ijk=ijkstart3,ijkend3)
356 
357  RETURN
358  END SUBROUTINE drag_gs
359 
360 !-----------------------------------------------------------------<<<
361 
362 ! Turton and Levenspiel (1986)
363 !----------------------------------------------------------------->>>
364  DOUBLE PRECISION FUNCTION c_dsxre_tl(RE)
365  USE param
366  USE param1
367  IMPLICIT NONE
368  DOUBLE PRECISION, INTENT(IN) :: RE ! Reynolds number
369 
370  c_dsxre_tl = 24.d0*(1.d0 + 0.173d0*re**0.657d0) + &
371  0.413d0*re**2.09d0/(re**1.09d0 + 16300.d0)
372  RETURN
373  END FUNCTION c_dsxre_tl
374 !-----------------------------------------------------------------<<<
375 
376 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
377 ! C
378 ! Subroutine: DRAG_SYAM_OBRIEN C
379 ! Purpose: Calculate the gas-solids drag coefficient C
380 ! C
381 ! Literature/Document References: C
382 ! Syamlal M, O'Brien TJ (1988). International Journal of C
383 ! Multiphase Flow 14: 473-481. C
384 ! C
385 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
386 
387  SUBROUTINE drag_syam_obrien(lDgA,EPg,Mug,ROg,VREL,&
388  dpm)
389 
390 !-----------------------------------------------
391 ! Modules
392 !-----------------------------------------------
393  USE constant, only : drag_c1, drag_d1
394  USE drag
395  USE param
396  USE param1
397  IMPLICIT NONE
398 !-----------------------------------------------
399 ! Dummy arguments
400 !-----------------------------------------------
401 ! drag coefficient
402  DOUBLE PRECISION, INTENT(OUT) :: ldGA
403 ! gas volume fraction
404  DOUBLE PRECISION, INTENT(IN) :: EPg
405 ! gas laminar viscosity
406  DOUBLE PRECISION, INTENT(IN) :: Mug
407 ! gas density
408  DOUBLE PRECISION, INTENT(IN) :: ROg
409 ! Magnitude of gas-solids relative velocity
410  DOUBLE PRECISION, INTENT(IN) :: VREL
411 ! particle diameter of solids phase M
412  DOUBLE PRECISION, INTENT(IN) :: DPM
413 !-----------------------------------------------
414 ! Local parameters
415 !-----------------------------------------------
416 ! Parameters in the Cluster-effect model
417 ! PARAMETER (a1 = 250.) !for G_s = 98 kg/m^2.s
418 ! PARAMETER (a1 = 1500.) !for G_s = 147 kg/m^2.s
419 ! a1 depends upon solids flux. It has been represented by C(1)
420 ! defined in the data file.
421 ! DOUBLE PRECISION, PARAMETER :: A2 = 0.005D0
422 ! DOUBLE PRECISION, PARAMETER :: A3 = 90.0D0
423 ! DOUBLE PRECISION, PARAMETER :: RE_C = 5.D0
424 ! DOUBLE PRECISION, PARAMETER :: EP_C = 0.92D0
425 !-----------------------------------------------
426 ! Local variables
427 !-----------------------------------------------
428 ! Variables which are function of EP_g
429  DOUBLE PRECISION :: A, BB
430 ! Ratio of settling velocity of a multiparticle system to
431 ! that of a single particle
432  DOUBLE PRECISION :: V_rm
433 ! Reynolds number
434  DOUBLE PRECISION :: RE
435 !-----------------------------------------------
436 
437  IF(mug > zero) THEN
438  re = dpm*vrel*rog/mug
439  ELSE
440  re = large_number
441  ENDIF
442 
443 ! Calculate V_rm
444  a = epg**4.14d0
445  IF (epg <= 0.85d0) THEN
446  bb = drag_c1*epg**1.28d0
447  ELSE
448  bb = epg**drag_d1
449  ENDIF
450 
451  v_rm=half*(a-0.06d0*re+&
452  sqrt((3.6d-3)*re*re+0.12d0*re*(2.d0*bb-a)+a*a) )
453 
454 !------------------Begin cluster correction --------------------------
455 ! uncomment the following four lines ...
456 ! V_RM = V_RM * (ONE + C(1)*&
457 ! EXP(-A2*(RE-RE_C)**2 - A3*(EPg-EP_C)**2)* &
458 ! RE*(1. - EPg))
459 !------------------End cluster correction ----------------------------
460 
461  ldga = 0.75d0*mug*epg*c_dsxre_dv(re/v_rm)/(v_rm*dpm*dpm)
462 
463  IF (re == zero) ldga = zero
464 
465  RETURN
466 
467  END SUBROUTINE drag_syam_obrien
468 
469 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
470 ! C
471 ! Subroutine: DRAG_GIDASPOW C
472 ! Purpose: Calculate the gas-solids drag coefficient C
473 ! C
474 ! Literature/Document References: C
475 ! Ding J, Gidaspow D (1990). AIChE Journal 36: 523-538. C
476 ! C
477 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
478 
479  SUBROUTINE drag_gidaspow(lDgA,EPg,Mug,ROg,ROPg,VREL, DPM)
480 
481 !-----------------------------------------------
482 ! Modules
483 !-----------------------------------------------
484  USE drag
485  USE param
486  USE param1
487  IMPLICIT NONE
488 !-----------------------------------------------
489 ! Dummy arguments
490 !-----------------------------------------------
491 ! drag coefficient
492  DOUBLE PRECISION, INTENT(OUT) :: lDgA
493 ! gas volume fraction
494  DOUBLE PRECISION, INTENT(IN) :: EPg
495 ! gas laminar viscosity
496  DOUBLE PRECISION, INTENT(IN) :: Mug
497 ! gas density
498  DOUBLE PRECISION, INTENT(IN) :: ROg
499 ! gas density*EP_g
500  DOUBLE PRECISION, INTENT(IN) :: ROPg
501 ! Magnitude of gas-solids relative velocity
502  DOUBLE PRECISION, INTENT(IN) :: VREL
503 ! particle diameter of solids phase M or
504 ! average particle diameter if PCF
505  DOUBLE PRECISION, INTENT(IN) :: DPM
506 
507 !-----------------------------------------------
508 ! Local variables
509 !-----------------------------------------------
510 ! Reynolds number
511  DOUBLE PRECISION :: RE
512 ! Single sphere drag coefficient
513  DOUBLE PRECISION :: C_d
514 !-----------------------------------------------
515 
516 ! Note the presence of gas volume fraction in ROPG
517  re = merge(dpm*vrel*ropg/mug, large_number, mug > zero)
518 
519 ! Dense phase
520  IF(epg <= 0.8d0) THEN
521  ldga = 150d0*(one-epg)*mug / (epg*dpm**2) + &
522  1.75d0*rog*vrel/dpm
523  ELSE
524 ! Dilute phase - EP_g >= 0.8
525  IF(re <= 1000d0)THEN
526 ! this could be replaced with the function C_DS_SN
527  c_d = c_ds_sn(re)
528  ELSE
529  c_d = 0.44d0
530  ENDIF
531  ldga = 0.75d0*c_d*vrel*ropg*epg**(-2.65d0) / dpm
532  ENDIF
533 
534  IF (re == zero) ldga = zero
535 
536  RETURN
537 
538  END SUBROUTINE drag_gidaspow
539 
540 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
541 ! C
542 ! Subroutine: DRAG_GIDASPOW_BLEND C
543 ! Purpose: Calculate the gas-solids drag coefficient C
544 ! C
545 ! Author: Charles E.A. Finney Date: 23-Mar-06 C
546 ! Reviewer: Sreekanth Pannala Date: 24-Mar-06 C
547 ! C
548 ! Literature/Document References: C
549 ! original source unknown: C
550 ! Lathouwers D, Bellan J (2000). Proceedings of the 2000 U.S. DOE C
551 ! Hydrogen Program Review NREL/CP-570-28890. Available from C
552 ! http://www.eere.energy.gov/hydrogenandfuelcells/pdfs/28890k.pdf. C
553 ! C
554 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
555 
556  SUBROUTINE drag_gidaspow_blend(lDgA,EPg,Mug,ROg,ROPg,VREL,&
557  dpm)
558 
559 !-----------------------------------------------
560 ! Modules
561 !-----------------------------------------------
562  USE param
563  USE param1
564  USE constant, only : pi
565  IMPLICIT NONE
566 !-----------------------------------------------
567 ! Dummy arguments
568 !-----------------------------------------------
569 ! drag coefficient
570  DOUBLE PRECISION, INTENT(OUT) :: lDgA
571 ! gas volume fraction
572  DOUBLE PRECISION, INTENT(IN) :: EPg
573 ! gas laminar viscosity
574  DOUBLE PRECISION, INTENT(IN) :: Mug
575 ! gas density
576  DOUBLE PRECISION, INTENT(IN) :: ROg
577 ! gas density*EP_g
578  DOUBLE PRECISION, INTENT(IN) :: ROPg
579 ! Magnitude of gas-solids relative velocity
580  DOUBLE PRECISION, INTENT(IN) :: VREL
581 ! particle diameter of solids phase M or
582 ! average particle diameter if PCF
583  DOUBLE PRECISION, INTENT(IN) :: DPM
584 !-----------------------------------------------
585 ! Local variables
586 !-----------------------------------------------
587 ! Reynolds number
588  DOUBLE PRECISION :: RE
589 ! Single sphere drag coefficient
590  DOUBLE PRECISION :: C_d
591 ! Gidaspow switch function variables
592  DOUBLE PRECISION :: Ergun
593  DOUBLE PRECISION :: WenYu
594  DOUBLE PRECISION :: PHI_gs
595 !-----------------------------------------------
596 
597  IF(mug > zero) THEN
598 ! Note the presence of gas volume fraction in ROPG
599  re = dpm*vrel*ropg/mug
600  ELSE
601  re = large_number
602  ENDIF
603 
604 ! Dense phase - EP_g <= 0.8
605  ergun = 150d0*(one-epg)*mug / (epg*dpm**2) + &
606  1.75d0*rog*vrel/dpm
607 ! Dilute phase - EP_g >= 0.8
608  IF(re <= 1000d0)THEN
609  c_d = (24.d0/(re+small_number)) * &
610  (one + 0.15d0 * re**0.687d0)
611  ELSE
612  c_d = 0.44d0
613  ENDIF
614  wenyu = 0.75d0*c_d*vrel*ropg*epg**(-2.65d0) / dpm
615 
616 ! Switch function
617  phi_gs = atan(150.d0*1.75d0*(epg - 0.8d0))/pi + 0.5d0
618 
619 ! Blend the models
620  ldga = (1.d0-phi_gs)*ergun + phi_gs*wenyu
621  IF (re == zero) ldga = zero
622 
623  RETURN
624  END SUBROUTINE drag_gidaspow_blend
625 
626 
627 
628 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
629 ! C
630 ! Subroutine: DRAG_WEN_YU C
631 ! Purpose: Calculate the gas-solids drag coefficient C
632 ! C
633 ! Literature/Document References: C
634 ! Wen CY, Yu YH (1966). Chemical Engineering Progress Symposium C
635 ! Series 62: 100-111. C
636 ! C
637 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
638 
639  SUBROUTINE drag_wen_yu(lDgA,EPg,Mug,ROPg,VREL,DPM)
640 
641 !-----------------------------------------------
642 ! Modules
643 !-----------------------------------------------
644  USE param
645  USE param1
646  IMPLICIT NONE
647 !-----------------------------------------------
648 ! Dummy arguments
649 !-----------------------------------------------
650 ! drag coefficient
651  DOUBLE PRECISION, INTENT(OUT) :: lDgA
652 ! gas volume fraction
653  DOUBLE PRECISION, INTENT(IN) :: EPg
654 ! gas laminar viscosity
655  DOUBLE PRECISION, INTENT(IN) :: Mug
656 ! gas density*EP_g
657  DOUBLE PRECISION, INTENT(IN) :: ROPg
658 ! Magnitude of gas-solids relative velocity
659  DOUBLE PRECISION, INTENT(IN) :: VREL
660 ! particle diameter of solids phase M or
661 ! average particle diameter if PCF
662  DOUBLE PRECISION, INTENT(IN) :: DPM
663 !-----------------------------------------------
664 ! Local variables
665 !-----------------------------------------------
666 ! Reynolds number
667  DOUBLE PRECISION :: RE
668 ! Single sphere drag coefficient
669  DOUBLE PRECISION :: C_d
670 !-----------------------------------------------
671 
672  IF(mug > zero) THEN
673 ! Note the presence of gas volume fraction in ROPG
674  re = dpm*vrel*ropg/mug
675  ELSE
676  re = large_number
677  ENDIF
678 
679  IF(re <= 1000.0d0)THEN
680  c_d = (24.d0/(re+small_number)) * (one + 0.15d0*re**0.687d0)
681  ELSE
682  c_d = 0.44d0
683  ENDIF
684 
685  ldga = 0.75d0 * c_d * vrel * ropg * epg**(-2.65d0) / dpm
686  IF (re == zero) ldga = zero
687 
688  RETURN
689  END SUBROUTINE drag_wen_yu
690 
691 
692 
693 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
694 ! C
695 ! Subroutine: SUBGRID_DRAG_IGCI C
696 ! Purpose: Calculate subgrid correction to the gas-solids drag C
697 ! coefficient developed by Wen-Yu C
698 ! C
699 ! Author: Sebastien Dartevelle, LANL, May 2013 C
700 ! C
701 ! Revision: 1 C
702 ! Purpose: Minor changes, e.g., fix inconsistenty with analogous C
703 ! calls in des_drag_gp and with new variable density feature C
704 ! Author: Janine Carney, June 2013 C
705 ! C
706 ! Literature/Document References: C
707 ! Igci, Y., Pannala, S., Benyahia, S., & Sundaresan S., C
708 ! Validation studies on filtered model equations for gas- C
709 ! particle flows in risers, Industrial & Engineering Chemistry C
710 ! Research, 2012, 51(4), 2094-2103 C
711 ! C
712 ! C
713 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
714 
715  SUBROUTINE subgrid_drag_igci(lDgA,EPg,Mug,ROg,DPM,ROs,IJK)
716 
717 !-----------------------------------------------
718 ! Modules
719 !-----------------------------------------------
720  USE param
721  USE param1
722  USE run, only : filter_size_ratio, subgrid_wall
723  USE constant, only : gravity
724  USE geometry, only : vol,axy,do_k
725  IMPLICIT NONE
726 !-----------------------------------------------
727 ! Dummy arguments
728 !-----------------------------------------------
729 ! drag coefficient
730  DOUBLE PRECISION, INTENT(INOUT) :: lDgA
731 ! gas volume fraction
732  DOUBLE PRECISION, INTENT(IN) :: EPg
733 ! gas laminar viscosity
734  DOUBLE PRECISION, INTENT(IN) :: Mug
735 ! gas density
736  DOUBLE PRECISION, INTENT(IN) :: ROg
737 ! particle diameter of solids phase M
738  DOUBLE PRECISION, INTENT(IN) :: DPM
739 ! particle density of solids phase M
740  DOUBLE PRECISION, INTENT(IN) :: ROs
741 ! current cell index
742  INTEGER, INTENT(IN) :: IJK
743 !-----------------------------------------------
744 ! Local variables
745 !-----------------------------------------------
746 ! factor to correct the drag for subgrid domain effects
747  DOUBLE PRECISION :: F_Subgrid
748 ! factor to correct the drag for subgrid domain effects arising from
749 ! wall
750  DOUBLE PRECISION :: F_SubGridWall
751 ! particle terminal settling velocity from stokes' formulation
752  DOUBLE PRECISION :: vt
753 ! filter size which is a function of each grid cell volume
754  DOUBLE PRECISION :: filtersize
755 ! inverse Froude number, or dimensionless filter size
756  DOUBLE PRECISION :: Inv_Froude
757 ! total solids volume fraction
758  DOUBLE PRECISION :: EPs
759 ! Variables for Igci model
760  DOUBLE PRECISION :: GG_phip, h_phip, h_phip2, c_function,&
761  f_filter
762 !-----------------------------------------------
763 
764 ! initialize
765  f_subgrid = one
766  f_subgridwall = one
767 
768 ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
769  vt = gravity*dpm*dpm*(ros - rog) / (18.0d0*mug)
770 ! filter size calculation for each specific gridcell volume
771  IF(do_k) THEN
772  filtersize = filter_size_ratio * (vol(ijk)**(one/3.0d0))
773  ELSE
774  filtersize = filter_size_ratio * dsqrt(axy(ijk))
775  ENDIF
776 
777 ! dimensionless inverse of Froude number
778  IF(abs(vt) > small_number) THEN
779  inv_froude = filtersize * gravity / vt**2
780  ELSE
781  inv_froude = large_number
782  ENDIF
783 
784 ! total solids volume fraction
785  eps = one - epg
786 
787  IF (eps .LT. 0.0012d0) THEN
788  h_phip = 2.7d0*(eps**0.234)
789  ELSEIF (eps .LT. 0.014d0) THEN
790  h_phip = -0.019d0/(eps**0.455) + 0.963d0
791  ELSEIF (eps .LT. 0.25d0) THEN
792  h_phip = 0.868d0*exp((-0.38*eps)) - &
793  0.176d0*exp((-119.2*eps))
794  ELSEIF (eps .LT. 0.455d0) THEN
795  h_phip = -4.59d-5*exp((19.75*eps)) + &
796  0.852d0*exp((-0.268*eps))
797  ELSEIF (eps .LE. 0.59d0) THEN
798  h_phip = (eps - 0.59d0) * (-1501.d0*(eps**3) + &
799  2203.d0*(eps**2) - 1054.d0*eps + 162.d0)
800  ELSE
801  h_phip=zero
802  ENDIF
803 
804  IF (eps .LT. 0.18d0) THEN
805  gg_phip = (eps**0.24)*(1.48d0 + exp(-18.0*eps))
806  ELSE
807  gg_phip = one
808  ENDIF
809 
810 ! a filter function needed in Igci Filtered/subgrid Model [dimensionless]
811  f_filter = (inv_froude**1.6) / ((inv_froude**1.6)+0.4d0)
812  h_phip2=h_phip*gg_phip
813  c_function=-h_phip2*f_filter
814  f_subgrid =(one + c_function)
815 
816  IF (subgrid_wall) THEN
817  CALL subgrid_drag_wall(f_subgridwall,vt,ijk)
818  ENDIF
819 
820  ldga = f_subgridwall*f_subgrid * ldga
821 
822  RETURN
823  END SUBROUTINE subgrid_drag_igci
824 
825 
826 
827 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
828 ! C
829 ! Subroutine: SUBGRID_DRAG_MILIOLI C
830 ! Purpose: Calculate subgrid correction to the gas-solids drag C
831 ! coefficient developed by Wen-Yu C
832 ! C
833 ! Author: Sebastien Dartevelle, LANL, May 2013 C
834 ! C
835 ! Revision: 1 C
836 ! Purpose: Minor changes, e.g., fix inconsistenty with analogous C
837 ! calls in des_drag_gp and with new variable density feature C
838 ! Author: Janine Carney, June 2013 C
839 ! C
840 ! Literature/Document References: C
841 ! Milioli, C. C., et al., Filtered two-fluid models of fluidized C
842 ! gas-particle flows: new constitutive relations, AICHE J, C
843 ! doi: 10.1002/aic.14130 C
844 ! C
845 ! Comments: C
846 ! Still needs to be reviewed for accuracy with source material C
847 ! C
848 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
849 
850  SUBROUTINE subgrid_drag_milioli(lDgA,EPg,Mug,ROg,VREL,DPM,ROs,&
851  ijk)
852 
853 !-----------------------------------------------
854 ! Modules
855 !-----------------------------------------------
856  USE param
857  USE param1
858  USE run, only : filter_size_ratio, subgrid_wall
859  USE constant, only : gravity
860  USE geometry, only : vol,axy,do_k
861  IMPLICIT NONE
862 !-----------------------------------------------
863 ! Dummy arguments
864 !-----------------------------------------------
865 ! drag coefficient
866  DOUBLE PRECISION, INTENT(INOUT) :: lDgA
867 ! gas volume fraction
868  DOUBLE PRECISION, INTENT(IN) :: EPg
869 ! gas laminar viscosity
870  DOUBLE PRECISION, INTENT(IN) :: Mug
871 ! gas density
872  DOUBLE PRECISION, INTENT(IN) :: ROg
873 ! Magnitude of gas-solids relative velocity
874  DOUBLE PRECISION, INTENT(IN) :: VREL
875 ! particle diameter of solids phase M
876  DOUBLE PRECISION, INTENT(IN) :: DPM
877 ! particle density of solids phase M
878  DOUBLE PRECISION, INTENT(IN) :: ROs
879 ! current cell index
880  INTEGER, INTENT(IN) :: IJK
881 !-----------------------------------------------
882 ! Local variables
883 !-----------------------------------------------
884 ! factor to correct the drag for subgrid domain effects
885  DOUBLE PRECISION :: F_Subgrid
886 ! factor to correct the drag for subgrid domain effects arising from
887 ! wall
888  DOUBLE PRECISION :: F_SubGridWall
889 ! particle terminal settling velocity from stokes' formulation
890  DOUBLE PRECISION :: vt
891 ! filter size which is a function of each grid cell volume
892  DOUBLE PRECISION :: filtersize
893 ! inverse Froude number, or dimensionless filter size
894  DOUBLE PRECISION :: Inv_Froude
895 ! dimensionless slip velocity = VREL/vt
896  DOUBLE PRECISION :: vslip
897 ! total solids volume fraction
898  DOUBLE PRECISION :: EPs
899 ! Variables for Milioli model
900  DOUBLE PRECISION :: h1, henv, hlin
901 !-----------------------------------------------
902 
903 ! initialize
904  f_subgrid = one
905  f_subgridwall = one
906 
907 ! particle terminal settling velocity: vt = g*d^2*(Rho_s - Rho_g) / 18 * Mu_g
908  vt = gravity*dpm*dpm*(ros - rog) / (18.0d0*mug)
909 ! filter size calculation for each specific gridcell volume
910  IF(do_k) THEN
911  filtersize = filter_size_ratio * (vol(ijk)**(one/3.0d0))
912  ELSE
913  filtersize = filter_size_ratio * dsqrt(axy(ijk))
914  ENDIF
915 ! dimensionless inverse of Froude number
916  IF(abs(vt) > small_number) THEN
917  inv_froude = filtersize * gravity / vt**2
918  ELSE
919  inv_froude = large_number
920  ENDIF
921 ! total solids volume fractionn
922  eps = one - epg
923 ! dimensionless slip velocity between gas and solids phase M
924  vslip = vrel / vt
925 
926  IF (inv_froude .LE. 1.028d0) THEN
927  h1 = (1.076d0 + 0.12d0*vslip - (0.02d0/(vslip+0.01d0)))*eps + &
928  (0.084d0 + 0.09d0*vslip - (0.01d0/(0.1d0*vslip+0.01d0)))
929  IF (eps .LE. 0.53d0) THEN
930  henv = (6.8d0*(one+eps)*(eps**0.3)) / &
931  (10.d0*(eps**1.5) + 5.d0)
932  ELSEIF (eps .GT. 0.53d0 .AND. eps .LE. 0.65d0) THEN
933  henv = (2.23d0*((0.65d0-eps)**(0.45))) / &
934  ((one/eps)-one)
935  ELSEIF (eps .GT. 0.65d0) THEN
936  henv=zero
937  ENDIF
938 
939  ELSEIF (inv_froude .GT. 1.028d0 .AND. &
940  inv_froude .LE. 2.056d0) THEN
941  h1 = (1.268d0 - (0.2d0*vslip) + (0.14d0/(vslip+0.01d0)))*eps + &
942  (0.385d0 + 0.09d0*vslip - (0.05d0/(0.2d0*vslip+0.01d0)))
943  IF (eps .LE. 0.53d0) THEN
944  henv = (8.6d0*(one+eps)*(eps**0.2)) / (10.d0*eps + 6.3d0)
945  ELSEIF (eps .GT. 0.53d0 .AND. eps .LE. 0.65d0) THEN
946  henv = (0.423d0*((0.65d0-eps)**0.3)) / (one-(eps**0.4))
947  ELSEIF (eps .GT. 0.65d0) THEN
948  henv=zero
949  ENDIF
950 
951  ELSEIF (inv_froude .GT. 2.056d0 .AND. &
952  inv_froude .LE. 4.112d0) THEN
953  h1 = ((0.018d0*vslip + 0.1d0)/(0.14d0*vslip + 0.01d0))*eps + &
954  (0.9454d0 - (0.09d0/(0.2d0*vslip + 0.01d0)))
955  IF (eps .LE. 0.5d0) THEN
956  henv = (7.9d0*(one+eps)*(eps**0.2)) / &
957  (10.d0*(eps**0.9) + 5.d0)
958  ELSEIF (eps .GT. 0.5d0 .AND. eps .LE. 0.63d0) THEN
959  henv = (0.705d0*((0.63d0-eps)**0.3)) / (one-(eps**0.7))
960  ELSEIF (eps .GT. 0.63d0) THEN
961  henv=zero
962  ENDIF
963 
964  ELSEIF (inv_froude .GT. 4.112d0 .AND. &
965  inv_froude .LE. 8.224d0) THEN
966  h1 = ((0.05d0*vslip+0.3d0)/(0.4d0*vslip+0.06d0))*eps + &
967  (0.9466d0 - (0.05d0/(0.11d0*vslip+0.01d0)))
968  IF (eps .LE. 0.45d0) THEN
969  henv = (7.9d0*(one+eps)*(eps**0.2)) / &
970  ((10.d0*(eps**0.6)) + 3.6d0)
971  ELSEIF (eps .GT. 0.45d0 .AND. eps .LE. 0.57d0) THEN
972  henv = (0.78d0*((0.57d0-eps)**0.2)) / (one-(eps**0.9))
973  ELSEIF (eps .GT. 0.57d0) THEN
974  henv=zero
975  ENDIF
976 
977  ELSEIF (inv_froude .GT. 8.224d0 .AND. &
978  inv_froude .LE. 12.336d0) THEN
979  h1 = ((1.3d0*vslip+2.2d0)/(5.2d0*vslip+0.07d0))*eps + &
980  (0.9363d0-(0.11d0/(0.3d0*vslip+0.01d0)))
981  IF (eps .LE. 0.35d0) THEN
982  henv = (7.6d0*(one+eps)*(eps**0.2)) / &
983  ((10.d0*(eps**0.6)) + 3.3d0)
984  ELSEIF (eps .GT. 0.35d0 .AND. eps .LE. 0.55d0) THEN
985  henv = (0.81d0*((0.55d0-eps)**0.3)) / (one-(eps**0.7))
986  ELSEIF (eps .GT. 0.55d0) THEN
987  henv=zero
988  ENDIF
989 
990  ELSEIF (inv_froude .GT. 12.336d0 .AND. &
991  inv_froude .LE. 16.448d0) THEN
992  h1 = ((2.6d0*vslip+4.d0)/(10.d0*vslip+0.08d0))*eps + &
993  (0.926d0-(0.17d0/(0.5d0*vslip+0.01d0)))
994  IF (eps .LE. 0.25d0) THEN
995  henv = (8.4d0*(one+eps)*(eps**0.2)) / &
996  ((10.d0*(eps**0.5)) + 3.3d0)
997  ELSEIF (eps .GT. 0.25d0 .AND. eps .LE. 0.52d0) THEN
998  henv = (1.01d0*((0.52d0-eps)**0.03))/(one-(eps**0.9))
999  ELSEIF (eps .GT. 0.52d0) THEN
1000  henv=zero
1001  ENDIF
1002 
1003  ELSEIF (inv_froude .GT. 16.448d0 .AND. &
1004  inv_froude .LE. 20.56d0) THEN
1005  h1 = ((2.5d0*vslip+4.d0)/(10.d0*vslip+0.08d0))*eps + &
1006  (0.9261d0-(0.17d0/(0.5d0*vslip+0.01d0)))
1007  IF (eps .LE. 0.25d0) THEN
1008  henv = (8.4d0*(one+eps)*(eps**0.2)) / &
1009  ((10.d0*(eps**0.5)) + 3.3d0)
1010  ELSEIF (eps .GT. 0.25d0 .AND. eps .LE. 0.52d0) THEN
1011  henv = (1.065d0*((0.52d0-eps)**0.3))/(one-eps)
1012  ELSEIF (eps .GT.0.52d0) THEN
1013  henv=zero
1014  ENDIF
1015 
1016  ELSEIF (inv_froude .GT. 20.56d0) THEN
1017  h1 = ((1.6d0*vslip+4.d0)/(7.9d0*vslip+0.08d0))*eps + &
1018  (0.9394d0 - (0.22d0/(0.6d0*vslip+0.01d0)))
1019  IF (eps .LE. 0.25d0) THEN
1020  henv = (9.d0*(one+eps)*(eps**0.15)) / &
1021  (10.d0*(eps**0.45) + 4.2d0)
1022  ELSEIF (eps .GT. 0.25d0 .AND. eps .LE. 0.52d0) THEN
1023  henv = (0.91d0*((0.52d0-eps)**0.4))/(one-(eps**0.6))
1024  ELSEIF (eps .GT. 0.52d0) THEN
1025  henv=zero
1026  ENDIF
1027  ENDIF
1028 
1029  IF (h1 .GT. zero) THEN
1030  hlin=h1
1031  ELSE
1032  hlin=zero
1033  ENDIF
1034 
1035  IF (inv_froude .LT. 1.028d0) THEN
1036 ! for very small filtered size, the drag wont be changed:
1037 ! F_Subgrid = 1.0 - H where H = 0.0
1038  f_subgrid = one
1039  ELSE
1040 ! MIN(henv,hlin) is H in Milioli paper, 2013
1041  f_subgrid = one - min(henv,hlin)
1042  ENDIF
1043 
1044 ! Filtered drag = (1 - H)*Microscopic_drag; it is strange Milioli takes EPs
1045 ! F_Subgrid = EPs*(ONE-hmili)
1046 
1047  IF (subgrid_wall) THEN
1048  CALL subgrid_drag_wall(f_subgridwall,vt,ijk)
1049  ENDIF
1050 
1051  ldga = f_subgridwall*f_subgrid * ldga
1052 
1053 
1054  RETURN
1055  END SUBROUTINE subgrid_drag_milioli
1056 
1057 
1058 
1059 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1060 ! C
1061 ! Subroutine: SUBGRID_DRAG_WALL C
1062 ! Purpose: Calculate subgrid correction arising from wall to the C
1063 ! gas-solids drag coefficient C
1064 ! C
1065 ! Author: Sebastien Dartevelle, LANL, May 2013 C
1066 ! C
1067 ! Revision: 1 C
1068 ! Author: Janine Carney, June 2013 C
1069 ! C
1070 ! Literature/Document References: C
1071 ! Igci, Y., and Sundaresan, S., Verification of filtered two- C
1072 ! fluid models for gas-particle flows in risers, AICHE J., C
1073 ! 2011, 57 (10), 2691-2707. C
1074 ! C
1075 ! Comments: Currently only valid for free-slip wall but no checks C
1076 ! are made to ensure user has selected free-slip wall when this C
1077 ! option is invoked C
1078 ! C
1079 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1080 
1081  SUBROUTINE subgrid_drag_wall(lSubgridWall,vt,IJK)
1082 
1083 !-----------------------------------------------
1084 ! Modules
1085 !-----------------------------------------------
1086  USE param
1087  USE param1
1088  USE constant, only : gravity
1089  USE cutcell, only : dwall
1090  IMPLICIT NONE
1091 !-----------------------------------------------
1092 ! Dummy arguments
1093 !-----------------------------------------------
1094 ! factor to correct the drag for subgrid domain effects arising from
1095 ! wall
1096  DOUBLE PRECISION, INTENT(OUT) :: lSubGridWall
1097 ! particle terminal settling velocity from stokes' formulation
1098  DOUBLE PRECISION, INTENT(IN) :: vt
1099 ! current cell index
1100  INTEGER, INTENT(IN) :: IJK
1101 !-----------------------------------------------
1102 ! Local parameters
1103 !-----------------------------------------------
1104 ! values are only correct for FREE-Slip walls
1105  DOUBLE PRECISION, PARAMETER :: a22=6.0d0, b22=0.295d0
1106 !-----------------------------------------------
1107 ! Local variables
1108 !-----------------------------------------------
1109 ! dimensionless distance to wall
1110  DOUBLE PRECISION :: x_d
1111 !-----------------------------------------------
1112 ! initialize
1113  lsubgridwall = one
1114 
1115 ! dimensionless distance to the Wall
1116  x_d = dwall(ijk) * gravity / vt**2
1117 
1118 ! decrease exponentionally away from the wall
1119 ! more complex model could be implemented with JJ wall model
1120  lsubgridwall = one / ( one + a22 * (exp(-b22*x_d)) )
1121 
1122  RETURN
1123  END SUBROUTINE subgrid_drag_wall
1124 
1125 
1126 
1127 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1128 ! C
1129 ! Subroutine: DRAG_KOCH_HILL C
1130 ! Purpose: Calculate the gas-solids drag coefficient C
1131 ! C
1132 ! Author: Clay Sutton (Lehigh University) Date: 14-Jul-04 C
1133 ! C
1134 ! Revision: 1 C
1135 ! Author: Sofiane Benyahia Date: 21-Jan-05 C
1136 ! C
1137 ! Literature/Document References: C
1138 ! Benyahia S, Syamlal M, O'Brien TJ (2006). Powder Technology C
1139 ! 162: 166-174. C
1140 ! Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics C
1141 ! 448: 213-241. C
1142 ! Hill RJ, Koch DL, Ladd JC (2001). Journal of Fluid Mechanics C
1143 ! 448: 243-278. C
1144 ! C
1145 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1146 
1147  SUBROUTINE drag_koch_hill(lDgA,EPg,Mug,ROPg,VREL,&
1148  dpm,dpa,phis)
1149 
1150 !-----------------------------------------------
1151 ! Modules
1152 !-----------------------------------------------
1153  USE param
1154  USE param1
1155  IMPLICIT NONE
1156 !-----------------------------------------------
1157 ! Dummy arguments
1158 !-----------------------------------------------
1159 ! drag coefficient
1160  DOUBLE PRECISION, INTENT(OUT) :: lDgA
1161 ! gas volume fraction
1162  DOUBLE PRECISION, INTENT(IN) :: EPg
1163 ! gas laminar viscosity
1164  DOUBLE PRECISION, INTENT(IN) :: Mug
1165 ! gas density*EP_g
1166  DOUBLE PRECISION, INTENT(IN) :: ROPg
1167 ! Magnitude of gas-solids relative velocity
1168  DOUBLE PRECISION, INTENT(IN) :: VREL
1169 ! particle diameter of solids phase M
1170  DOUBLE PRECISION, INTENT(IN) :: DPM
1171 ! average particle diameter if pcf otherwise DPM again
1172  DOUBLE PRECISION, INTENT(IN) :: DPA
1173 ! total solids volume fraction of solids phases
1174  DOUBLE PRECISION, INTENT(IN) :: PHIS
1175 !-----------------------------------------------
1176 ! Local variables
1177 !-----------------------------------------------
1178 ! Reynolds number
1179  DOUBLE PRECISION :: RE
1180 ! transition Reynolds numbers
1181  DOUBLE PRECISION :: Re_Trans_1, Re_Trans_2
1182 ! Stokes Drag Force
1183  DOUBLE PRECISION :: F_STOKES
1184 ! zero Re function for low Reynolds number
1185  DOUBLE PRECISION :: F_0
1186 ! inertial function for low Reynolds number
1187  DOUBLE PRECISION :: F_1
1188 ! zero Re function for high Reynolds number
1189  DOUBLE PRECISION :: F_2
1190 ! inertial function for high Reynolds number
1191  DOUBLE PRECISION :: F_3
1192 ! dimensionless drag force F
1193  DOUBLE PRECISION :: F
1194 ! weighting factor to compute F_0 and F_2
1195  DOUBLE PRECISION :: ww
1196 !-----------------------------------------------
1197 
1198 
1199  IF(mug > zero) THEN
1200 ! Note the presence of gas volume fraction in ROPG and factor of 1/2
1201  re = 0.5d0*dpa*vrel*ropg/mug ! if pcf DPA otherwise DPM
1202  ELSE
1203  re = large_number
1204  ENDIF
1205 
1206  f_stokes = 18.d0*mug*epg*epg/dpm**2 ! use DPM
1207  ww = exp(-10.0d0*(0.4d0-phis)/phis)
1208 
1209  IF(phis > 0.01d0 .AND. phis < 0.4d0) THEN
1210  f_0 = (1.0d0-ww) * (1.0d0 + 3.0d0*dsqrt(phis/2.0d0) + &
1211  135.0d0/64.0d0*phis*log(phis) + 17.14d0*phis) / &
1212  (1.0d0 + 0.681d0*phis - 8.48d0*phis*phis + &
1213  8.16d0*phis**3) + ww*10.0d0*phis/(1.0d0-phis)**3
1214  ELSEIF(phis >= 0.4d0) THEN
1215  f_0 = 10.0d0*phis/(1.0d0-phis)**3
1216  ENDIF
1217 
1218  IF(phis > 0.01d0 .AND. phis <= 0.1d0) THEN
1219  f_1 = dsqrt(2.0d0/phis) / 40.0d0
1220  ELSE IF(phis > 0.1d0) THEN
1221  f_1 = 0.11d0 + 5.1d-04 * exp(11.6d0*phis)
1222  ENDIF
1223 
1224  IF(phis < 0.4d0) THEN
1225  f_2 = (1.0d0-ww) * (1.0d0 + 3.0d0*dsqrt(phis/2.0d0) + &
1226  135.0d0/64.0d0*phis*log(phis) + 17.89d0*phis) / &
1227  (1.0d0 + 0.681d0*phis - 11.03d0*phis*phis + &
1228  15.41d0*phis**3)+ ww*10.0d0*phis/(1.0d0-phis)**3
1229  ELSE
1230  f_2 = 10.0d0*phis/(1.0d0-phis)**3
1231  ENDIF
1232 
1233  IF(phis < 0.0953d0) THEN
1234  f_3 = 0.9351d0*phis + 0.03667d0
1235  ELSE
1236  f_3 = 0.0673d0 + 0.212d0*phis +0.0232d0/(1.0-phis)**5
1237  ENDIF
1238 
1239  re_trans_1 = (f_2 - 1.0d0)/(3.0d0/8.0d0 - f_3)
1240  re_trans_2 = (f_3 + dsqrt(f_3*f_3 - 4.0d0*f_1 &
1241  *(f_0-f_2))) / (2.0d0*f_1)
1242 
1243  IF(phis <= 0.01d0 .AND. re <= re_trans_1) THEN
1244  f = 1.0d0 + 3.0d0/8.0d0*re
1245  ELSEIF(phis > 0.01d0 .AND. re <= re_trans_2) THEN
1246  f = f_0 + f_1*re*re
1247  ELSEIF(phis <= 0.01d0 .AND. re > re_trans_1 .OR. &
1248  phis > 0.01d0 .AND. re > re_trans_2) THEN
1249  f = f_2 + f_3*re
1250  ELSE
1251  f = zero
1252  ENDIF
1253 
1254  ldga = f * f_stokes
1255  IF (re == zero) ldga = zero
1256 
1257  RETURN
1258  END SUBROUTINE drag_koch_hill
1259 
1260 
1261 
1262 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1263 ! C
1264 ! Subroutine: DRAG_BVK C
1265 ! Purpose: Calculate the gas-solids drag coefficient C
1266 ! C
1267 ! Literature/Document References: C
1268 ! Beetstra, van der Hoef, Kuipers, Chem. Eng. Science 62 C
1269 ! (Jan 2007) C
1270 ! C
1271 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1272 
1273  SUBROUTINE drag_bvk(lDgA,EPg,Mug,ROPg,VREL,&
1274  dpm,dpa,phis)
1275 
1276 !-----------------------------------------------
1277 ! Modules
1278 !-----------------------------------------------
1279  USE param
1280  USE param1
1281  IMPLICIT NONE
1282 !-----------------------------------------------
1283 ! Dummy arguments
1284 !-----------------------------------------------
1285 ! drag coefficient
1286  DOUBLE PRECISION, INTENT(OUT) :: lDgA
1287 ! gas volume fraction
1288  DOUBLE PRECISION, INTENT(IN) :: EPg
1289 ! gas laminar viscosity
1290  DOUBLE PRECISION, INTENT(IN) :: Mug
1291 ! gas density*EP_g
1292  DOUBLE PRECISION, INTENT(IN) :: ROPg
1293 ! magnitude of gas-solids relative velocity
1294  DOUBLE PRECISION, INTENT(IN) :: VREL
1295 ! particle diameter of solids phase M or
1296  DOUBLE PRECISION, INTENT(IN) :: DPM
1297 ! average particle diameter
1298  DOUBLE PRECISION, INTENT(IN) :: DPA
1299 ! total solids volume fraction of solids phases
1300  DOUBLE PRECISION, INTENT(IN) :: PHIS
1301 !-----------------------------------------------
1302 ! Local variables
1303 !-----------------------------------------------
1304 ! Reynolds number
1305  DOUBLE PRECISION :: RE
1306 ! Stokes Drag Force
1307  DOUBLE PRECISION :: F_STOKES
1308 ! dimensionless drag force F
1309  DOUBLE PRECISION :: F
1310 !-----------------------------------------------
1311 
1312  IF(mug > zero) THEN
1313 ! Note the presence of gas volume fraction in ROPG
1314  re = dpa*vrel*ropg/mug ! use DPA
1315  ELSE
1316  re = large_number
1317  ENDIF
1318 
1319 ! eq(9) BVK J. fluid. Mech. 528, 2005
1320 ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g)
1321  f_stokes = 18d0*mug*epg/dpm**2 ! use DPM
1322 
1323  f = 10d0*phis/epg**2 + epg**2*(one+1.5d0*dsqrt(phis))
1324  f = f + 0.413d0*re/(24.d0*epg**2) * &
1325  (one/epg + 3d0*epg*phis + 8.4d0/re**0.343) / &
1326  (one+10.d0**(3d0*phis)/re**(0.5+2.d0*phis))
1327 
1328  ldga = f*f_stokes
1329  IF (re == zero) ldga = zero
1330 
1331  RETURN
1332  END SUBROUTINE drag_bvk
1333 
1334 
1335 
1336 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1337 ! C
1338 ! Subroutine: DRAG_HYS C
1339 ! Purpose: Calculate the gas-solids drag coefficient C
1340 ! C
1341 ! Literature/Document References: C
1342 ! Yin, X, Sundaresan, S. (2009). AIChE Journal 55: no 6, 1352- C
1343 ! 1368 C
1344 ! C
1345 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1346 
1347  SUBROUTINE drag_hys(lDgA,EPg,Mug,ROPg,VREL,&
1348  dpm,dpa,y_i,ep_sm,phis,m,maxm,ijk)
1349 
1350 !-----------------------------------------------
1351 ! Modules
1352 !-----------------------------------------------
1353  USE param
1354  USE param1
1355  USE drag, only : beta_ij
1356  USE run, only : lam_hys
1357  IMPLICIT NONE
1358 !-----------------------------------------------
1359 ! Dummy arguments
1360 !-----------------------------------------------
1361 ! drag coefficient
1362  DOUBLE PRECISION, INTENT(OUT) :: lDgA
1363 ! gas volume fraction
1364  DOUBLE PRECISION, INTENT(IN) :: EPg
1365 ! gas laminar viscosity
1366  DOUBLE PRECISION, INTENT(IN) :: Mug
1367 ! gas density*EP_g
1368  DOUBLE PRECISION, INTENT(IN) :: ROPg
1369 ! magnitude of gas-solids relative velocity
1370  DOUBLE PRECISION, INTENT(IN) :: VREL
1371 ! local variable for the particle diameter
1372  DOUBLE PRECISION :: DPM(dim_m)
1373 ! average particle diameter
1374  DOUBLE PRECISION, INTENT(IN) :: DPA
1375 ! diameter ratio in polydisperse systems
1376  DOUBLE PRECISION, INTENT(IN) :: Y_i
1377 ! local variable for the solids volume fraction
1378  DOUBLE PRECISION :: EP_sM(dim_m)
1379 ! total solids volume fraction of solids phases
1380  DOUBLE PRECISION, INTENT(IN) :: PHIS
1381 ! current solids phase index and fluid cell index
1382  INTEGER, INTENT(IN) :: M, IJK
1383 ! maximum number of solids phases
1384  INTEGER, INTENT(IN) :: MAXM
1385 !-----------------------------------------------
1386 ! Local variables
1387 !-----------------------------------------------
1388 ! minimum particle diameter in mixture
1389  DOUBLE PRECISION :: DPmin
1390 ! Index for particles of other species
1391  INTEGER :: L
1392 ! Reynolds number
1393  DOUBLE PRECISION :: RE
1394 ! Stokes Drag Force
1395  DOUBLE PRECISION :: F_STOKES
1396 ! dimensionless drag force F
1397  DOUBLE PRECISION :: F
1398 ! Polydisperse correction factor for YS drag relation
1399  DOUBLE PRECISION :: a_YS
1400 ! Lubrication interaction prefactor in YS drag relation
1401  DOUBLE PRECISION :: alpha_YS
1402 ! Friction coefficient for a particle of type i (HYS drag relation)
1403  DOUBLE PRECISION :: beta_i_HYS
1404 ! Friction coefficient for a particle of type j (HYS drag relation)
1405  DOUBLE PRECISION :: beta_j_HYS
1406 ! Stokes drag of a particle of type j
1407  DOUBLE PRECISION :: FSTOKES_j
1408 ! Diameter ratio for particle of type j
1409  DOUBLE PRECISION :: Y_i_J
1410 ! Variable for Beetstra et. al. drag relation
1411  DOUBLE PRECISION :: F_D_BVK
1412 ! Variable for YS drag relation
1413  DOUBLE PRECISION :: F_YS
1414 !-----------------------------------------------
1415 
1416  IF (mug > zero) THEN
1417 ! Note the presence of gas volume fraction in ROPG
1418  re = dpa*vrel*ropg/mug ! use DPA
1419  ELSE
1420  re = large_number
1421  ENDIF
1422 
1423 ! (this F_Stokes is /= of Koch_Hill by a factor of ep_g/ep_sm)
1424 ! (this F_Stokes is /= of BVK by a factor of 1/ep_sm)
1425  f_stokes = 18d0*mug*epg*ep_sm(m)/dpm(m)**2 ! use DPM
1426 
1427 ! Find smallest diameter if number of particle types is greater than 1
1428  dpmin= dpm(1)
1429  IF (maxm > 1) THEN
1430  DO l=2,maxm
1431  dpmin = min(dpmin,dpm(l))
1432  ENDDO
1433  ENDIF
1434 
1435  a_ys = 1d0 - 2.66d0*phis + 9.096d0*phis**2 - 11.338d0*phis**3
1436 
1437 ! Calculate the prefactor of the off-diagonal friction coefficient
1438 ! Use default value of lamdba if there are no particle asparities
1439  alpha_ys = 1.313d0*log10(dpmin/lam_hys) - 1.249d0
1440 
1441 
1442 ! Beetstra correction for monodisperse drag
1443  f_d_bvk = zero
1444  f = 10.d0*phis/epg**2 + epg**2*(one+1.5d0*dsqrt(phis))
1445 
1446  IF(re > zero) f_d_bvk = f + 0.413d0*re/(24.d0*epg**2)*&
1447  (one/epg + 3.d0*epg*phis + 8.4d0/re**0.343d0) / &
1448  (one+10**(3.d0*phis)/re**(0.5d0+2.d0*phis))
1449 
1450 ! YS correction for polydisperse drag
1451  f_ys = 1d0/epg + (f_d_bvk - 1d0/epg)*&
1452  (a_ys*y_i+(1d0-a_ys)*y_i**2)
1453  f_ys = f_ys*f_stokes
1454  beta_i_hys = f_ys
1455 
1456  DO l= 1,maxm
1457  IF (l /= m) THEN
1458  y_i_j = dpm(l)/dpa
1459  beta_j_hys = 1.d0/epg + (f_d_bvk - 1.d0/epg) * &
1460  (a_ys*y_i_j + (1d0-a_ys)*y_i_j**2)
1461  fstokes_j = 18.d0*mug*ep_sm(l)*epg/&
1462  dpm(l)**2
1463 
1464  beta_j_hys = beta_j_hys*fstokes_j
1465 
1466 ! Calculate off-diagonal friction coefficient
1467  beta_ij(ijk,m,l) = zero
1468 
1469 ! This if statement prevents NaN values from appearing for beta_ij
1470  IF (ep_sm(m) > zero .AND. ep_sm(l) > zero) &
1471  beta_ij(ijk,m,l) = (2.d0*alpha_ys*ep_sm(m)*ep_sm(l))/ &
1472  (ep_sm(m)/beta_i_hys + ep_sm(l)/beta_j_hys)
1473  f_ys = f_ys + beta_ij(ijk,m,l)
1474 
1475  ENDIF ! end if (J/=M)
1476  ENDDO ! end do (J=1,MAXM)
1477 
1478 
1479  ldga = f_ys
1480 
1481  RETURN
1482  END SUBROUTINE drag_hys
drag::c_ds_sn
double precision function c_ds_sn(RE)
Definition: drag_mod.f:49
fldvar::v_s
double precision, dimension(:,:), allocatable v_s
Definition: fldvar_mod.f:105
funits::dmp_log
logical dmp_log
Definition: funits_mod.f:6
ur_facs::ur_f_gs
double precision ur_f_gs
Definition: ur_facs_mod.f:17
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
subgrid_drag_igci
subroutine subgrid_drag_igci(lDgA, EPg, Mug, ROg, DPM, ROs, IJK)
Definition: drag_gs.f:716
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param1::one
double precision, parameter one
Definition: param1_mod.f:29
drag_usr
subroutine drag_usr(IJK, M_NP, lDgA, EPg, Mug, ROg, VREL, DPM, ROs, lUg, lVg, lWg)
Definition: usr_drag.f:27
run::lam_hys
double precision lam_hys
Definition: run_mod.f:140
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
fldvar::w_s
double precision, dimension(:,:), allocatable w_s
Definition: fldvar_mod.f:117
run::subgrid_wall
logical subgrid_wall
Definition: run_mod.f:131
drag
Definition: drag_mod.f:11
param::dim_m
integer, parameter dim_m
Definition: param_mod.f:67
drag_bvk
subroutine drag_bvk(lDgA, EPg, Mug, ROPg, VREL, DPM, DPA, PHIS)
Definition: drag_gs.f:1275
fldvar::u_s
double precision, dimension(:,:), allocatable u_s
Definition: fldvar_mod.f:93
fun_avg
Definition: fun_avg_mod.f:1
drag_gidaspow
subroutine drag_gidaspow(lDgA, EPg, Mug, ROg, ROPg, VREL, DPM)
Definition: drag_gs.f:480
fldvar::d_p
double precision, dimension(:,:), allocatable d_p
Definition: fldvar_mod.f:57
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
ur_facs
Definition: ur_facs_mod.f:1
mms
Definition: mms_mod.f:12
subgrid_drag_wall
subroutine subgrid_drag_wall(lSubgridWall, vt, IJK)
Definition: drag_gs.f:1082
drag::c_dsxre_dv
double precision function c_dsxre_dv(RE)
Definition: drag_mod.f:37
run::filter_size_ratio
double precision filter_size_ratio
Definition: run_mod.f:133
fldvar::v_g
double precision, dimension(:), allocatable v_g
Definition: fldvar_mod.f:99
fldvar::w_g
double precision, dimension(:), allocatable w_g
Definition: fldvar_mod.f:111
param1::half
double precision, parameter half
Definition: param1_mod.f:28
funits::unit_log
integer, parameter unit_log
Definition: funits_mod.f:21
run
Definition: run_mod.f:13
param1::large_number
double precision, parameter large_number
Definition: param1_mod.f:23
param
Definition: param_mod.f:2
fldvar::ro_s
double precision, dimension(:,:), allocatable ro_s
Definition: fldvar_mod.f:45
subgrid_drag_milioli
subroutine subgrid_drag_milioli(lDgA, EPg, Mug, ROg, VREL, DPM, ROs, IJK)
Definition: drag_gs.f:852
drag_gidaspow_blend
subroutine drag_gidaspow_blend(lDgA, EPg, Mug, ROg, ROPg, VREL, DPM)
Definition: drag_gs.f:558
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
drag_koch_hill
subroutine drag_koch_hill(lDgA, EPg, Mug, ROPg, VREL, DPM, DPA, PHIS)
Definition: drag_gs.f:1149
compar::mype
integer mype
Definition: compar_mod.f:24
physprop::mu_g
double precision, dimension(:), allocatable mu_g
Definition: physprop_mod.f:68
constant::gravity
double precision gravity
Definition: constant_mod.f:149
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
constant::drag_d1
double precision drag_d1
Definition: constant_mod.f:139
mms::use_mms
logical use_mms
Definition: mms_mod.f:15
drag_hys
subroutine drag_hys(lDgA, EPg, Mug, ROPg, VREL, DPM, DPA, Y_I, EP_sM, PHIS, M, MAXM, IJK)
Definition: drag_gs.f:1349
fldvar::u_g
double precision, dimension(:), allocatable u_g
Definition: fldvar_mod.f:87
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
machine::start_log
subroutine start_log
Definition: machine_mod.f:182
drag::f_gs
double precision, dimension(:,:), allocatable f_gs
Definition: drag_mod.f:14
run::model_b
logical model_b
Definition: run_mod.f:88
physprop::smax
integer smax
Definition: physprop_mod.f:22
cutcell::dwall
double precision, dimension(:), allocatable dwall
Definition: cutcell_mod.f:488
c_dsxre_tl
double precision function c_dsxre_tl(RE)
Definition: drag_gs.f:365
drag_wen_yu
subroutine drag_wen_yu(lDgA, EPg, Mug, ROPg, VREL, DPM)
Definition: drag_gs.f:640
drag::beta_ij
double precision, dimension(:,:,:), allocatable beta_ij
Definition: drag_mod.f:20
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
drag_syam_obrien
subroutine drag_syam_obrien(lDgA, EPg, Mug, ROg, VREL, DPM)
Definition: drag_gs.f:389
constant::drag_c1
double precision drag_c1
Definition: constant_mod.f:139
geometry::vol
double precision, dimension(:), allocatable vol
Definition: geometry_mod.f:212
drag_gs
subroutine drag_gs(M, IER)
Definition: drag_gs.f:22
open_pe_log
subroutine open_pe_log(IER)
Definition: open_files.f:270
fldvar::ro_g
double precision, dimension(:), allocatable ro_g
Definition: fldvar_mod.f:32
fldvar::rop_g
double precision, dimension(:), allocatable rop_g
Definition: fldvar_mod.f:38
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
machine::end_log
subroutine end_log
Definition: machine_mod.f:208