MFIX  2016-1
calc_d_ghd.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: CALC_D_ghd_e(A_m, VxF_gs, d_e, IER) C
4 ! Purpose: calculte coefficients linking velocity correction to C
5 ! pressure correction -- East C
6 ! C
7 ! Note: MFIX convention: center coeff is negative, hence: C
8 ! (-A_M(IJK,0,M)) > or = 0 C
9 ! The EAST face area is AYZ C
10 ! Modifed for GHD theory C
11 ! C
12 ! Author: M. Syamlal Date: 21-JUN-96 C
13 ! Reviewer: Date: C
14 ! C
15 ! Revision Number: 2 C
16 ! Purpose: allow multiparticle in D_E calculation in accounting for C
17 ! the averaged Solid-Solid drag and multi-solid-particle dragC
18 ! C
19 ! Author: S. Dartevelle, LANL Date: 28-FEB-04 C
20 ! Reviewer: Date: dd-mmm-yy C
21 ! C
22 ! Revision Number: 3 C
23 ! Purpose: more accurate formulas in the D_s(M) for Model_A C
24 ! C
25 ! Author: S. Dartevelle, LANL Date: 17-MAR-04 C
26 ! Reviewer: Date: dd-mmm-yy C
27 ! C
28 ! C
29 ! Literature/Document References: C
30 ! C
31 ! Variables referenced: C
32 ! Variables modified: C
33 ! C
34 ! Local variables: C
35 ! C
36 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
37 !
38  SUBROUTINE calc_d_ghd_e(A_M, VXF_GS, D_E)
39 !
40 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98
41 !...Switches: -xf
42 !
43 !-----------------------------------------------
44 ! M o d u l e s
45 !-----------------------------------------------
46  USE param
47  USE param1
48  USE parallel
49  USE fldvar
50  USE geometry
51  USE indices
52  USE physprop
53  USE run
54  USE scales
55  USE compar
56  USE sendrecv
57  USE fun_avg
58  USE functions
59  IMPLICIT NONE
60 !-----------------------------------------------
61 ! G l o b a l P a r a m e t e r s
62 !-----------------------------------------------
63 !-----------------------------------------------
64 ! D u m m y A r g u m e n t s
65 !-----------------------------------------------
66 !
67 ! Septadiagonal matrix A_m
68  DOUBLE PRECISION A_m(dimension_3, -3:3, 0:dimension_m)
69 !
70 ! Volume x average at momentum cell centers
71  DOUBLE PRECISION VxF_gs(dimension_3, dimension_m)
72 !
73  DOUBLE PRECISION d_e(dimension_3, 0:dimension_m)
74 !
75 ! Average volume fraction at momentum cell centers
76  DOUBLE PRECISION EPGA !S. Dartevelle, LANL, Feb.2004
77 ! Usual Indices
78  INTEGER M, L, I, IJK, IJKE !S. Dartevelle, LANL, Feb.2004
79 ! Average solid volume fraction at momentum cell centers
80  DOUBLE PRECISION EPSA(dimension_m) !S. Dartevelle, LANL, Feb.2004
81 ! ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
82  DOUBLE PRECISION other_ratio_1, FOA1 !S. Dartevelle, LANL, Feb.2004
83 ! sum of All Solid-Gas VolxDrag
84  DOUBLE PRECISION SUM_VXF_GS !S. Dartevelle, LANL, Feb.2004
85 ! What phase momentum equation is activated?
86  LOGICAL Pass1, Pass2 !S. Dartevelle, LANL, Feb.2004
87 ! numerator needed for solid phase M time EPSA
88  DOUBLE PRECISION numeratorxEP(dimension_m) !S. Dartevelle, LANL, Feb.2004
89 ! denominator needed for solid phase M
90  DOUBLE PRECISION denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
91 ! denominator needed for solid phase M
92  DOUBLE PRECISION other_denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
93 ! total solids volume fraction
94  DOUBLE PRECISION EPStmp
95 !-----------------------------------------------
96 
97  pass1 = .false. !initialization
98  pass2 = .false.
99  m = mmax
100  if (momentum_x_eq(0) .AND. momentum_x_eq(m)) then
101  pass1 = .true. !we have at least one solid phase X-momentum equation
102  !with the gas phase X-momentum equation
103  elseif (momentum_x_eq(m)) then
104  pass2 = .true. !we have at least one solid phase X-momentum equation
105  !but the gas phase X-momentum is not solved
106  endif
107 
108 
109  IF (pass1) THEN !Gas and at least one solid phase X-momentum equation
110 !!!$omp parallel do private(I,IJK, IJKE, EPGA, EPSA,EPStmp, numeratorxEP, &
111 !!!$omp& M, L, Lp, LpL, LM, other_ratio_1, FOA1, &
112 !!!$omp& SUM_VXF_GS, other_denominator, denominator ),&
113 !!!$omp& schedule(static)
114  DO ijk = ijkstart3, ijkend3
115  IF (ip_at_e(ijk) .OR. mflow_at_e(ijk)) THEN !impermeable
116  DO m= 0, mmax
117  d_e(ijk,m) = zero
118  END DO
119  ELSE
120  i = i_of(ijk)
121  ijke = east_of(ijk)
122  epga = avg_x(ep_g(ijk),ep_g(ijke),i)
123 
124  sum_vxf_gs = zero
125  epsa(mmax) = zero
126  DO m= 1, smax
127  epsa(m) = avg_x(ep_s(ijk,m),ep_s(ijke,m),i)
128  epsa(mmax) = epsa(mmax) + epsa(m)
129  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m) !Gas - All Solids VolxDrag summation
130  END DO
131 
132  m = mmax
133  other_ratio_1 = ( vxf_gs(ijk,m)* (-a_m(ijk,0,m)) /&
134  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
135  )
136 
137  IF (model_b) THEN !Model B
138  !Linking velocity correction coefficient to pressure - GAS Phase
139  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
140  d_e(ijk,0) = p_scale*ayz(ijk)/( (-a_m(ijk,0,0))+other_ratio_1 )
141  else
142  d_e(ijk,0) = zero
143  endif
144  !Linking velocity correction coefficient to pressure - SOLID Phase
145  m = mmax
146  if ( momentum_x_eq(m) .AND. ( ((-a_m(ijk,0,m))>small_number) .OR. &
147  (vxf_gs(ijk,m)>small_number) ) ) then
148  d_e(ijk,m) = d_e(ijk,0)*(&
149  vxf_gs(ijk,m)/((-a_m(ijk,0,m))+vxf_gs(ijk,m))&
150  )
151  else
152  d_e(ijk,m) = zero
153  endif
154  ELSE !Model A
155  foa1 = zero
156  m = mmax
157  foa1 = foa1 + (epsa(m)*vxf_gs(ijk,m)/&
158  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
159  )
160  other_denominator(m) = vxf_gs(ijk,m)*( (-a_m(ijk,0,0))/&
161  ((-a_m(ijk,0,0))+sum_vxf_gs+small_number)&
162  )
163  numeratorxep(m) = zero
164  denominator(m) = zero
165  !Linking velocity correction coefficient to pressure - GAS Phase
166  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
167  d_e(ijk,0) = p_scale*ayz(ijk)*(epga+foa1)/( (-a_m(ijk,0,0))+other_ratio_1 )
168  else
169  d_e(ijk,0) = zero
170  endif
171  !Linking velocity correction coefficient to pressure - SOLID Phase
172  m = mmax
173  if ( momentum_x_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
174  (other_denominator(m)>small_number) ) ) then
175  d_e(ijk,m) = p_scale*ayz(ijk)*(&
176  ( epsa(m) + (vxf_gs(ijk,m)*epga/((-a_m(ijk,0,0))+sum_vxf_gs+small_number)) )/&
177  ( (-a_m(ijk,0,m))+other_denominator(m) )&
178  )
179  else
180  d_e(ijk,m) = zero
181  endif
182  ENDIF !end of Model_B/Model_A if then condition
183  ENDIF
184  ENDDO
185 
186  ELSE IF (momentum_x_eq(0)) THEN !the solid X-momentum equations are not solved
187  !only gas phase X-momentum
188 !!!$omp parallel do private(IJK, I, IJKE, EPGA, M, SUM_VXF_GS), &
189 !!!$omp& schedule(static)
190  DO ijk = ijkstart3, ijkend3
191  IF (ip_at_e(ijk) .OR. mflow_at_e(ijk)) THEN !impermeable
192  d_e(ijk,0) = zero
193  ELSE
194  i = i_of(ijk)
195  ijke = east_of(ijk)
196  epga = avg_x(ep_g(ijk),ep_g(ijke),i)
197 
198 
199  sum_vxf_gs = vxf_gs(ijk,mmax) !Gas - All Solids VolxDrag summation
200 
201  IF ( ((-a_m(ijk,0,0))>small_number) .OR. (sum_vxf_gs>small_number) ) THEN
202  IF (model_b) THEN
203  d_e(ijk,0) = p_scale*ayz(ijk)/((-a_m(ijk,0,0))+sum_vxf_gs)
204  ELSE
205  d_e(ijk,0) = p_scale*ayz(ijk)*epga/((-a_m(ijk,0,0))+sum_vxf_gs)
206  ENDIF
207  ELSE
208  d_e(ijk,0) = zero
209  ENDIF
210  ENDIF
211  END DO
212 
213  ELSE IF (pass2) THEN !at least one solid phase X-momentum equation is solved
214  !but the gas phase X-momentum is not solved
215 !!!$omp parallel do private(IJK, I, IJKE, EPSA,EPStmp, L, Lp, M, &
216 !!!$omp& numeratorxEP, denominator), &
217 !!!$omp& schedule(static)
218  DO ijk = ijkstart3, ijkend3
219  IF (ip_at_e(ijk) .OR. mflow_at_e(ijk) .OR. model_b) THEN
220  DO m= 1, mmax
221  d_e(ijk,m) = zero
222  END DO
223  ELSE
224  i = i_of(ijk)
225  ijke = east_of(ijk)
226 
227  m = mmax
228  epstmp = zero
229  DO l=1,smax
230  epstmp = epstmp+ avg_x(ep_s(ijk,l),ep_s(ijke,l),i)
231  ENDDO
232  epsa(m) = epstmp
233 
234  !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
235  m = mmax
236  if ( momentum_x_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
237  (vxf_gs(ijk,m)>small_number) ) ) then
238  d_e(ijk,m) = p_scale*ayz(ijk)*( epsa(m) )/&
239  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m) )
240  else
241  d_e(ijk,m) = zero
242  endif
243  ENDIF
244  END DO
245 
246  ENDIF
247 
248  RETURN
249  END SUBROUTINE calc_d_ghd_e
250 
251 
252 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
253 ! C
254 ! Module name: CALC_D_ghd_n(A_m, VxF_gs, d_n, IER) C
255 ! Purpose: calculte coefficients linking velocity correction to C
256 ! pressure correction -- North C
257 ! C
258 ! Author: M. Syamlal Date: 21-JUN-96 C
259 ! Reviewer: Date: C
260 ! C
261 ! Note: MFIX convention: center coeff is negative, hence C
262 ! (-A_M(IJK,0,M)) > or = 0 C
263 ! The NORTH face area is AXZ C
264 ! Modifed for GHD theory C
265 ! C
266 ! Revision Number: 2 C
267 ! Purpose: allow multiparticle in D_N calculation in accounting for C
268 ! the averaged Solid-Solid drag C
269 ! Author: S. Dartevelle, LANL Date: 28-FEB-04 C
270 ! Reviewer: Date: dd-mmm-yy C
271 ! C
272 ! Revision Number: 3 C
273 ! Purpose: more accurate formulas in the D_s(M) for Model_A C
274 ! C
275 ! Author: S. Dartevelle, LANL Date: 17-MAR-04 C
276 ! Reviewer: Date: dd-mmm-yy C
277 ! C
278 ! Literature/Document References: C
279 ! C
280 ! Variables referenced: C
281 ! Variables modified: C
282 ! C
283 ! Local variables: C
284 ! C
285 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
286 !
287  SUBROUTINE calc_d_ghd_n(A_M, VXF_GS, D_N)
288 !
289 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98
290 !...Switches: -xf
291 !
292 !-----------------------------------------------
293 ! M o d u l e s
294 !-----------------------------------------------
295  USE param
296  USE param1
297  USE parallel
298  USE fldvar
299  USE geometry
300  USE indices
301  USE physprop
302  USE run
303  USE scales
304  USE compar
305  USE sendrecv
306  USE fun_avg
307  USE functions
308  IMPLICIT NONE
309 !-----------------------------------------------
310 ! G l o b a l P a r a m e t e r s
311 !-----------------------------------------------
312 !-----------------------------------------------
313 ! D u m m y A r g u m e n t s
314 !-----------------------------------------------
315 !
316 ! Septadiagonal matrix A_m
317  DOUBLE PRECISION A_m(dimension_3, -3:3, 0:dimension_m)
318 !
319 ! Volume x average at momentum cell centers
320  DOUBLE PRECISION VxF_gs(dimension_3, dimension_m)
321 !
322  DOUBLE PRECISION d_n(dimension_3, 0:dimension_m)
323 !
324 ! Average volume fraction at momentum cell centers
325  DOUBLE PRECISION EPGA !S. Dartevelle, LANL, Feb.2004
326 ! Usual Indices
327  INTEGER M, L, J, IJK, IJKN !S. Dartevelle, LANL, Feb.2004
328 ! Average solid volume fraction at momentum cell centers
329  DOUBLE PRECISION EPSA(dimension_m) !S. Dartevelle, LANL, Feb.2004
330 ! ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
331  DOUBLE PRECISION other_ratio_1, FOA1 !S. Dartevelle, LANL, Feb.2004
332 ! sum of All Solid-Gas VolxDrag
333  DOUBLE PRECISION SUM_VXF_GS !S. Dartevelle, LANL, Feb.2004
334 ! What phase momentum equation is activated?
335  LOGICAL Pass1, Pass2 !S. Dartevelle, LANL, Feb.2004
336 ! numerator needed for solid phase M time EPSA
337  DOUBLE PRECISION numeratorxEP(dimension_m) !S. Dartevelle, LANL, Feb.2004
338 ! denominator needed for solid phase M
339  DOUBLE PRECISION denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
340 ! denominator needed for solid phase M
341  DOUBLE PRECISION other_denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
342 ! total solids volume fraction
343  DOUBLE PRECISION EPStmp
344 !-----------------------------------------------
345 
346  pass1 = .false. !initialization
347  pass2 = .false.
348  m = mmax
349  if (momentum_y_eq(0) .AND. momentum_y_eq(m)) then
350  pass1 = .true. !we have at least one solid phase Y-momentum equation
351  !with the gas phase Y-momentum equation
352  elseif (momentum_y_eq(m)) then
353  pass2 = .true. !we have at least one solid phase Y-momentum equation
354  !but the gas phase Y-momentum is not solved
355  endif
356 
357 
358  IF (pass1) THEN !Gas and at least one solid phases Y-momentum equation
359 !!!$omp parallel do private(J,IJK, IJKN, EPGA, EPSA, EPStmp, numeratorxEP, &
360 !!!$omp& M, L, Lp, LpL, LM, other_ratio_1, FOA1, &
361 !!!$omp& SUM_VXF_GS, other_denominator, denominator ),&
362 !!!$omp& schedule(static)
363  DO ijk = ijkstart3, ijkend3
364  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk)) THEN !impermeable
365  DO m= 0, mmax
366  d_n(ijk,m) = zero
367  END DO
368  ELSE
369  j = j_of(ijk)
370  ijkn = north_of(ijk)
371  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
372 
373  sum_vxf_gs = zero
374  epsa(mmax) = zero
375  DO m= 1, smax
376  epsa(m) = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
377  epsa(mmax) = epsa(mmax) + epsa(m)
378  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m) !Gas - All Solids VolxDrag summation
379  END DO
380 
381  m= mmax
382  other_ratio_1 = ( vxf_gs(ijk,m)* (-a_m(ijk,0,m)) /&
383  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
384  )
385 
386  IF (model_b) THEN !Model B
387  !Linking velocity correction coefficient to pressure - GAS Phase
388  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
389  d_n(ijk,0) = p_scale*axz(ijk)/( (-a_m(ijk,0,0))+other_ratio_1 )
390  else
391  d_n(ijk,0) = zero
392  endif
393  !Linking velocity correction coefficient to pressure - SOLID Phase
394  m = mmax
395  if ( momentum_y_eq(m) .AND. ( ((-a_m(ijk,0,m))>small_number) .OR. &
396  (vxf_gs(ijk,m)>small_number) ) ) then
397  d_n(ijk,m) = d_n(ijk,0)*(&
398  vxf_gs(ijk,m)/((-a_m(ijk,0,m))+vxf_gs(ijk,m))&
399  )
400  else
401  d_n(ijk,m) = zero
402  endif
403  ELSE !Model A
404  foa1 = zero
405  m = mmax
406  foa1 = foa1 + (epsa(m)*vxf_gs(ijk,m)/&
407  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
408  )
409  other_denominator(m) = vxf_gs(ijk,m)*( (-a_m(ijk,0,0))/&
410  ((-a_m(ijk,0,0))+sum_vxf_gs+small_number)&
411  )
412  numeratorxep(m) = zero
413  denominator(m) = zero
414  !Linking velocity correction coefficient to pressure - GAS Phase
415  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
416  d_n(ijk,0) = p_scale*axz(ijk)*(epga+foa1)/( (-a_m(ijk,0,0))+other_ratio_1 )
417  else
418  d_n(ijk,0) = zero
419  endif
420  !Linking velocity correction coefficient to pressure - SOLID Phase
421  m = mmax
422  if ( momentum_y_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
423  (other_denominator(m)>small_number) ) ) then
424  d_n(ijk,m) = p_scale*axz(ijk)*(&
425  ( epsa(m) + (vxf_gs(ijk,m)*epga/((-a_m(ijk,0,0))+sum_vxf_gs+small_number)) )/&
426  ( (-a_m(ijk,0,m))+other_denominator(m) )&
427  )
428  else
429  d_n(ijk,m) = zero
430  endif
431  ENDIF !end of Model_B/Model_A if then condition
432  ENDIF
433  END DO
434 
435  ELSE IF (momentum_y_eq(0)) THEN !the solid Y-momentum equations are not solved
436  !only the gas phase
437 !!!$omp parallel do private(IJK, J, IJKN, EPGA, M, SUM_VXF_GS), &
438 !!!$omp& schedule(static)
439  DO ijk = ijkstart3, ijkend3
440  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk)) THEN !impermeable
441  d_n(ijk,0) = zero
442  ELSE
443  j = j_of(ijk)
444  ijkn = north_of(ijk)
445  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
446 
447 
448  sum_vxf_gs = vxf_gs(ijk,mmax) !Gas - All Solids VolxDrag summation
449 
450  IF ( ((-a_m(ijk,0,0))>small_number) .OR. (sum_vxf_gs>small_number) ) THEN
451  IF (model_b) THEN
452  d_n(ijk,0) = p_scale*axz(ijk)/((-a_m(ijk,0,0))+sum_vxf_gs)
453  ELSE
454  d_n(ijk,0) = p_scale*axz(ijk)*epga/((-a_m(ijk,0,0))+sum_vxf_gs)
455  ENDIF
456  ELSE
457  d_n(ijk,0) = zero
458  ENDIF
459  ENDIF
460  END DO
461 
462  ELSE IF (pass2) THEN !at least one solid phase momentum Y-equation is solved
463  !but the gas phase Y-momentum is not solved
464 !!!$omp parallel do private(IJK, J, IJKN, EPSA, EPStmp, L, Lp, M, &
465 !!!$omp& numeratorxEP, denominator), &
466 !!!$omp& schedule(static)
467  DO ijk = ijkstart3, ijkend3
468  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk) .OR. model_b) THEN
469  DO m= 1, mmax
470  d_n(ijk,m) = zero
471  END DO
472  ELSE
473  j = j_of(ijk)
474  ijkn = north_of(ijk)
475 
476  m = mmax
477  epstmp = zero
478  DO l=1,smax
479  epstmp = epstmp+ avg_y(ep_s(ijk,l),ep_s(ijkn,l),j)
480  ENDDO
481  epsa(m) = epstmp
482 
483  !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
484  m = mmax
485  if ( momentum_y_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
486  (vxf_gs(ijk,m)>small_number) ) ) then
487  d_n(ijk,m) = p_scale*axz(ijk)*( epsa(m) )/&
488  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m) )
489  else
490  d_n(ijk,m) = zero
491  endif
492  ENDIF
493  END DO
494 
495  ENDIF
496 
497  RETURN
498  END SUBROUTINE calc_d_ghd_n
499 
500 
501 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
502 ! C
503 ! Module name: CALC_D_ghd_t(A_m, VxF_gs, d_t, IER) C
504 ! Purpose: calculte coefficients linking velocity correction to C
505 ! pressure correction -- Top C
506 ! C
507 ! Author: M. Syamlal Date: 21-JUN-96 C
508 ! Reviewer: Date: C
509 ! C
510 ! Note: MFIX convention: center coeff is negative, hence C
511 ! (-A_M(IJK,0,M)) > or = 0 C
512 ! The TOP face area is AXY C
513 ! Modifed for GHD theory C
514 ! C
515 ! Revision Number: 2 C
516 ! Purpose: allow multiparticle in D_T calculation in accounting for C
517 ! the averaged Solid-Solid drag C
518 ! Author: S. Dartevelle, LANL Date: 28-FEB-04 C
519 ! Reviewer: Date: dd-mmm-yy C
520 ! C
521 ! Revision Number: 3 C
522 ! Purpose: more accurate formulas in the D_s(M) for Model_A C
523 ! C
524 ! Author: S. Dartevelle, LANL Date: 17-MAR-04 C
525 ! Reviewer: Date: dd-mmm-yy C
526 ! C
527 ! Literature/Document References: C
528 ! C
529 ! Variables referenced: C
530 ! Variables modified: C
531 ! C
532 ! Local variables: C
533 ! C
534 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
535 !
536  SUBROUTINE calc_d_ghd_t(A_M, VXF_GS, D_T)
537 !
538 !...Translated by Pacific-Sierra Research VAST-90 2.06G5 12:17:31 12/09/98
539 !...Switches: -xf
540 !
541 !-----------------------------------------------
542 ! M o d u l e s
543 !-----------------------------------------------
544  USE param
545  USE param1
546  USE parallel
547  USE fldvar
548  USE geometry
549  USE indices
550  USE physprop
551  USE run
552  USE scales
553  USE compar
554  USE sendrecv
555  USE fun_avg
556  USE functions
557  IMPLICIT NONE
558 !-----------------------------------------------
559 ! G l o b a l P a r a m e t e r s
560 !-----------------------------------------------
561 !-----------------------------------------------
562 ! D u m m y A r g u m e n t s
563 !-----------------------------------------------
564 !
565 ! Septadiagonal matrix A_m
566  DOUBLE PRECISION A_m(dimension_3, -3:3, 0:dimension_m)
567 !
568 ! Volume x average at momentum cell centers
569  DOUBLE PRECISION VxF_gs(dimension_3, dimension_m)
570 !
571  DOUBLE PRECISION d_t(dimension_3, 0:dimension_m)
572 !
573 ! Average volume fraction at momentum cell centers
574  DOUBLE PRECISION EPGA !S. Dartevelle, LANL, Feb.2004
575 ! Usual Indices
576  INTEGER M, L, K, IJK, IJKT !S. Dartevelle, LANL, Feb.2004
577 ! Average solid volume fraction at momentum cell centers
578  DOUBLE PRECISION EPSA(dimension_m) !S. Dartevelle, LANL, Feb.2004
579 ! ratio of drag and A0 or A_solid and other sum of Solid-Solid drag
580  DOUBLE PRECISION other_ratio_1, FOA1 !S. Dartevelle, LANL, Feb.2004
581 ! sum of All Solid-Gas VolxDrag
582  DOUBLE PRECISION SUM_VXF_GS !S. Dartevelle, LANL, Feb.2004
583 ! What phase momentum equation is activated?
584  LOGICAL Pass1, Pass2 !S. Dartevelle, LANL, Feb.2004
585 ! numerator needed for solid phase M time EPSA
586  DOUBLE PRECISION numeratorxEP(dimension_m) !S. Dartevelle, LANL, Feb.2004
587 ! denominator needed for solid phase M
588  DOUBLE PRECISION denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
589 ! denominator needed for solid phase M
590  DOUBLE PRECISION other_denominator(dimension_m) !S. Dartevelle, LANL, Feb.2004
591 ! tmp variable for total solids volume fraction
592  DOUBLE PRECISION EPStmp
593 !-----------------------------------------------
594 
595  pass1 = .false. !initialization
596  pass2 = .false.
597  m = mmax
598  if (momentum_z_eq(0) .AND. momentum_z_eq(m)) then
599  pass1 = .true. !we have at least one solid phase Z-momentum equation
600  !with the gas phase Z-momentum equation
601  elseif (momentum_z_eq(m)) then
602  pass2 = .true. !we have at least one solid phase Z-momentum equation
603  !but the gas phase Z-momentum is not solved
604  endif
605 
606 
607  IF (pass1) THEN !Gas and at least one solid phases Z-momentum equation
608 !!!$omp parallel do private(K,IJK, IJKT, EPGA, EPSA, EPStmp, numeratorxEP, &
609 !!!$omp& M, L, Lp, LpL, LM, other_ratio_1, FOA1, &
610 !!!$omp& SUM_VXF_GS, other_denominator, denominator ),&
611 !!!$omp& schedule(static)
612  DO ijk = ijkstart3, ijkend3
613  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk)) THEN !impermeable
614  DO m= 0, mmax
615  d_t(ijk,m) = zero
616  END DO
617  ELSE
618  k = k_of(ijk)
619  ijkt = top_of(ijk)
620  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
621 
622  sum_vxf_gs = zero
623  epsa(mmax) = zero
624  DO m= 1, smax
625  epsa(m) = avg_z(ep_s(ijk,m),ep_s(ijkt,m),k)
626  epsa(mmax) = epsa(mmax) + epsa(m)
627  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m) !Gas - All Solids VolxDrag summation
628  END DO
629 
630  m = mmax
631  other_ratio_1 = ( vxf_gs(ijk,m)* (-a_m(ijk,0,m)) /&
632  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
633  )
634 
635  IF (model_b) THEN !Model B
636  !Linking velocity correction coefficient to pressure - GAS Phase
637  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
638  d_t(ijk,0) = p_scale*axy(ijk)/( (-a_m(ijk,0,0))+other_ratio_1 )
639  else
640  d_t(ijk,0) = zero
641  endif
642  !Linking velocity correction coefficient to pressure - SOLID Phase
643  m = mmax
644  if ( momentum_z_eq(m) .AND. ( ((-a_m(ijk,0,m))>small_number) .OR. &
645  (vxf_gs(ijk,m)>small_number) ) ) then
646  d_t(ijk,m) = d_t(ijk,0)*(&
647  vxf_gs(ijk,m)/((-a_m(ijk,0,m))+vxf_gs(ijk,m))&
648  )
649  else
650  d_t(ijk,m) = zero
651  endif
652  ELSE !Model A
653  foa1 = zero
654  m = mmax
655  foa1 = foa1 + (epsa(m)*vxf_gs(ijk,m)/&
656  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m)+small_number )&
657  )
658  other_denominator(m) = vxf_gs(ijk,m)*( (-a_m(ijk,0,0))/&
659  ((-a_m(ijk,0,0))+sum_vxf_gs+small_number)&
660  )
661  numeratorxep(m) = zero
662  denominator(m) = zero
663  !Linking velocity correction coefficient to pressure - GAS Phase
664  if ( ((-a_m(ijk,0,0))>small_number) .OR. (other_ratio_1>small_number) ) then
665  d_t(ijk,0) = p_scale*axy(ijk)*(epga+foa1)/( (-a_m(ijk,0,0))+other_ratio_1 )
666  else
667  d_t(ijk,0) = zero
668  endif
669  !Linking velocity correction coefficient to pressure - SOLID Phase
670  m = mmax
671  if ( momentum_z_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
672  (other_denominator(m)>small_number) ) ) then
673  d_t(ijk,m) = p_scale*axy(ijk)*(&
674  ( epsa(m) + (vxf_gs(ijk,m)*epga/((-a_m(ijk,0,0))+sum_vxf_gs+small_number)) )/&
675  ( (-a_m(ijk,0,m))+other_denominator(m) )&
676  )
677  else
678  d_t(ijk,m) = zero
679  endif
680  ENDIF !end of Model_B/Model_A if then condition
681  ENDIF
682  END DO
683 
684  ELSE IF (momentum_z_eq(0)) THEN !the solid Z-momentum equations are not solved
685  !only the gas phase Z-momentum is solved
686 !!!$omp parallel do private(IJK, K, IJKT, EPGA, M, SUM_VXF_GS), &
687 !!!$omp& schedule(static)
688  DO ijk = ijkstart3, ijkend3
689  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk)) THEN !impermeable
690  d_t(ijk,0) = zero
691  ELSE
692  k = k_of(ijk)
693  ijkt = top_of(ijk)
694  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
695 
696  sum_vxf_gs = vxf_gs(ijk,mmax) !Gas - All Solids VolxDrag summation
697 
698  IF ( ((-a_m(ijk,0,0))>small_number) .OR. (sum_vxf_gs>small_number) ) THEN
699  IF (model_b) THEN
700  d_t(ijk,0) = p_scale*axy(ijk)/((-a_m(ijk,0,0))+sum_vxf_gs)
701  ELSE
702  d_t(ijk,0) = p_scale*axy(ijk)*epga/((-a_m(ijk,0,0))+sum_vxf_gs)
703  ENDIF
704  ELSE
705  d_t(ijk,0) = zero
706  ENDIF
707  ENDIF
708  END DO
709 
710  ELSE IF (pass2) THEN !at least one solid phase momentum Z-equation is solved
711  !but the gas phase Z-momentum is not solved
712 !!!$omp parallel do private(IJK, K, IJKT, EPSA, EPStmp, L, Lp, M, &
713 !!!$omp& numeratorxEP, denominator), &
714 !!!$omp& schedule(static)
715  DO ijk = ijkstart3, ijkend3
716  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk) .OR. model_b) THEN
717  DO m= 1, mmax
718  d_t(ijk,m) = zero
719  END DO
720  ELSE
721  k = k_of(ijk)
722  ijkt = top_of(ijk)
723 
724  m = mmax
725  epstmp = zero
726  DO l=1,smax
727  epstmp = epstmp+ avg_z(ep_s(ijk,l),ep_s(ijkt,l),k)
728  ENDDO
729  epsa(m) = epstmp
730 
731  !Linking velocity correction coefficient to pressure - SOLID Phase (Model_A only)
732  m = mmax
733  if ( momentum_z_eq(m) .AND. ( (-a_m(ijk,0,m)>small_number) .OR. &
734  (vxf_gs(ijk,m)>small_number) ) ) then
735  d_t(ijk,m) = p_scale*axy(ijk)*( epsa(m) )/&
736  ( (-a_m(ijk,0,m))+vxf_gs(ijk,m) )
737  else
738  d_t(ijk,m) = zero
739  endif
740  ENDIF
741  END DO
742 
743  ENDIF
744 
745  RETURN
746  END SUBROUTINE calc_d_ghd_t
747 
748 !// Comments on the modifications for DMP version implementation
749 !// 001 Include header file and common declarations for parallelization
750 !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
run::momentum_y_eq
logical, dimension(0:dim_m) momentum_y_eq
Definition: run_mod.f:77
calc_d_ghd_e
subroutine calc_d_ghd_e(A_M, VXF_GS, D_E)
Definition: calc_d_ghd.f:39
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
run::momentum_x_eq
logical, dimension(0:dim_m) momentum_x_eq
Definition: run_mod.f:74
scales::p_scale
double precision p_scale
Definition: scales_mod.f:13
run::momentum_z_eq
logical, dimension(0:dim_m) momentum_z_eq
Definition: run_mod.f:80
geometry::ayz
double precision, dimension(:), allocatable ayz
Definition: geometry_mod.f:206
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
calc_d_ghd_n
subroutine calc_d_ghd_n(A_M, VXF_GS, D_N)
Definition: calc_d_ghd.f:288
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
run
Definition: run_mod.f:13
geometry::axz
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
param
Definition: param_mod.f:2
calc_d_ghd_t
subroutine calc_d_ghd_t(A_M, VXF_GS, D_T)
Definition: calc_d_ghd.f:537
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
fldvar::ep_s
double precision function ep_s(IJK, xxM)
Definition: fldvar_mod.f:178
run::model_b
logical model_b
Definition: run_mod.f:88
physprop::smax
integer smax
Definition: physprop_mod.f:22
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27