MFIX  2016-1
calc_d_t.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2 ! !
3 ! Subroutine: CALC_D_T !
4 ! Author: M. Syamlal Date: 21-JUN-96 !
5 ! !
6 ! Purpose: Calculate coefficients linking velocity correction to !
7 ! pressure correction -- Top !
8 ! !
9 ! Note that the D_T coefficients for phases M>0 are generally not !
10 ! used unless the solids phase has close_packed=F, in which case the !
11 ! D_E coefficients for that phase are employed in a mixture pressure !
12 ! correction equation and for correcting velocities. !
13 ! !
14 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
15  SUBROUTINE calc_d_t(A_M, VXF_GS, VXF_SS, D_T, IER)
16 
17 ! Global Variables:
18 !---------------------------------------------------------------------//
19 ! Number of solids phases.
20  use physprop, only: mmax
21 ! Flag: Solve the Z-Momentum Equations
22  use run, only: momentum_z_eq
23 ! Flag: Coupled DEM, PIC, or Hybrid simulation
24  use discretelement, only: des_continuum_coupled
25 ! Flag: TMF-DEM solids hybrid model
26  use discretelement, only: des_continuum_hybrid
27 ! Volume x average at momentum cell center drag for DEM/PIC
28  use discretelement, only: vxf_gds, vxf_sds
29 ! Number of solids phases.
30  use physprop, only: mmax
31 
32 ! Global Parameters:
33 !---------------------------------------------------------------------//
34 ! Size of IJK arrays and size of solids phase arrays.
35  use param, only: dimension_3, dimension_m
36 ! Size of MxM upper triangular matrix.
37  use param1, only: dimension_lm
38 
39  IMPLICIT NONE
40 
41 ! Dummy arguments
42 !---------------------------------------------------------------------//
43 ! Septadiagonal matrix A_m
44  DOUBLE PRECISION, INTENT(IN):: A_m(dimension_3,-3:3,0:dimension_m)
45 ! Volume x average at momentum cell centers
46  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
47 ! Volume x average at momentum cell centers Solid-Solid Drag
48  DOUBLE PRECISION, INTENT(IN) :: VxF_ss(dimension_3, dimension_lm)
49 ! Coefficients for pressure correction equation. These coefficients are
50 ! initialized to zero in the subroutine time_march before the time loop.
51  DOUBLE PRECISION, INTENT(INOUT) :: d_t(dimension_3, 0:dimension_m)
52 ! Error index
53  INTEGER, INTENT(INOUT) :: IER
54 
55 ! Local variables:
56 !---------------------------------------------------------------------//
57 ! Numerator and demoninator needed for evaluation of pressure correction
58 ! coefficient for GAS phase. A temporary variable is used for local
59 ! manipulations.
60  DOUBLE PRECISION, dimension(:,:), allocatable :: AM0
61 ! Flag: One or more solids momentum equations are solved.
62  LOGICAL :: ANY_SOLIDS_Z_MOMENTUM
63 
64 !......................................................................!
65 
66  allocate(am0(dimension_3, 0:dimension_m))
67 
68 ! Initialize the error flag.
69  ier = 0
70 
71 ! Copy the gas phase A_M coefficients into temp array.
72  am0(:,:) = a_m(:,0,:)
73 
74 ! Add DEM temp A_M so they are included in the pressure correciton eq.
75  IF(des_continuum_coupled) THEN
76  am0(:,0) = am0(:,0) - vxf_gds(:)
77  IF (des_continuum_hybrid) &
78  am0(:,1:mmax) = am0(:,1:mmax) - vxf_sds(:,1:mmax)
79  ENDIF
80 
81  any_solids_z_momentum = any(momentum_z_eq(1:mmax))
82 
83 ! Determine which calculations are needed
84  IF (momentum_z_eq(0)) THEN
85  IF(any_solids_z_momentum)THEN
86  CALL calc_d_t_gas_and_solids(am0, vxf_gs, vxf_ss, d_t)
87  ELSE
88  CALL calc_d_t_gas_only(am0, vxf_gs, d_t)
89  ENDIF
90  ELSEIF (any_solids_z_momentum) THEN
91  CALL calc_d_t_solids_only(am0, vxf_gs, vxf_ss, d_t)
92  ENDIF
93 
94  deallocate(am0)
95 
96  RETURN
97  END SUBROUTINE calc_d_t
98 
99 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
100 ! !
101 ! Subroutine: CALC_D_T_GAS_AND_SOLIDS !
102 ! Author: M. Syamlal Date: 21-JUN-96 !
103 ! !
104 ! Purpose: Calculate coefficients linking velocity correction to !
105 ! pressure correction -- Top !
106 ! !
107 ! CASE: The gas phase Z-momentum equation is solved and at least one !
108 ! solids phase (m>0) Z-momentum equation is solved. !
109 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
110  SUBROUTINE calc_d_t_gas_and_solids(AM0, VXF_GS, VXF_SS, D_T)
111 
112 ! Global Variables:
113 !---------------------------------------------------------------------//
114 ! Volume fractions of gas and solids phases.
115  use fldvar, only: ep_g, ep_s
116 ! Number of solids phases.
117  use physprop, only: mmax
118 ! Flag: Solve the Z-Momentum Equations
119  use run, only: momentum_z_eq
120 ! Flag: Use MODEL_B
121  use run, only: model_b
122 ! Pressure scale factor
123  use scales, only: p_scale
124 ! Flag: Cartesian grid simulation
125  use cutcell, only: cartesian_grid
126 ! Area of cell's XY face.
127  use geometry, only: axy
128 ! Function to average across Z face.
129  use fun_avg, only: avg_z
130 ! Fluid grid loop bounds.
131  use compar, only: ijkstart3, ijkend3
132 ! Function to convert to phases to single array, IJK of cell to top
133  use functions, only: funlm, top_of
134 ! Flags: Impermeable surface and mass flow at top face
135  use functions, only: ip_at_t, mflow_at_t
136 ! Indices: K index of cell
137  use indices, only: k_of
138 
139 ! Global Parameters:
140 !---------------------------------------------------------------------//
141 ! Size of IJK arrays and size of solids phase arrays.
142  use param, only: dimension_3, dimension_m
143 ! Size of MxM upper triangular matrix.
144  use param1, only: dimension_lm
145 ! Double precision parameters.
146  use param1, only: zero, small_number, one
147 
148  IMPLICIT NONE
149 
150 ! Dummy arguments
151 !---------------------------------------------------------------------//
152 ! Temporary variable to store A matrix for local manipulation
153  DOUBLE PRECISION, INTENT(IN) :: AM0(dimension_3, 0:dimension_m)
154 ! Volume x average at momentum cell centers
155  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
156 ! Volume x average at momentum cell centers Solid-Solid Drag
157  DOUBLE PRECISION, INTENT(IN) :: VxF_ss(dimension_3, dimension_lm)
158 ! Coefficients for pressure correction equation. These coefficients are
159 ! initialized to zero in the subroutine time_march before the time loop.
160  DOUBLE PRECISION, INTENT(INOUT) :: D_T(dimension_3, 0:dimension_m)
161 
162 ! Local variables:
163 !---------------------------------------------------------------------//
164 ! Usual Indices
165  INTEGER :: LM, M, L, LPL, LP
166  INTEGER :: IJK, IJKT, K
167 ! Average solid volume fraction at momentum cell centers
168  DOUBLE PRECISION :: EPSA(dimension_m)
169 ! Average volume fraction at momentum cell centers
170  DOUBLE PRECISION :: EPGA
171 ! local alias for face area
172  DOUBLE PRECISION :: AREA_FACE
173 ! sum of All Solid-Gas VolxDrag
174  DOUBLE PRECISION :: SUM_VXF_GS
175 ! sum of Solid M - All other Solid drag
176  DOUBLE PRECISION :: SUM_VXF_SS(dimension_m)
177 ! sum of Solid L - All other Solid VolxDrag but M
178  DOUBLE PRECISION :: SUM_VXF_SS_wt_M
179 ! numerator and demoninator needed for evaluation of pressure correction
180 ! coefficient for GAS phase
181  DOUBLE PRECISION :: DEN_MGas, NUM_MGas
182 ! component of the total numerator needed for evaluation of pressure
183 ! correction coefficient for SOLIDS phase M. specifically, for when the
184 ! sum L is over solids phase only
185  DOUBLE PRECISION :: NUM_MSol_LSol(dimension_m)
186 ! component of the total denominator needed for evaluation of pressure
187 ! correction coefficient for SOLIDS phase M. specifically, for when the
188 ! sum L is over solids phases only
189  DOUBLE PRECISION :: DEN_MSol_LSol(dimension_m)
190 ! component of the total numerator needed for evaluation of pressure
191 ! correction coefficient for SOLIDS phase M. specifically, for when the
192 ! sum L is over gas phase only
193  DOUBLE PRECISION :: NUM_MSol_LGas
194 ! component of the total denominator needed for evaluation of pressure
195 ! correction coefficient for SOLIDS phase M. specifically, for when the
196 ! sum L is over gas phase only
197  DOUBLE PRECISION :: DEN_MSol_LGas
198 ! Temp variable for double precision values.
199  DOUBLE PRECISION :: TMPdp
200 
201 !......................................................................!
202 
203 !$omp parallel default(none) &
204 
205 !$omp private(ijk,area_face,k,ijkt,epga,sum_vxf_gs,epsa,lm,sum_vxf_ss,den_mgas,num_mgas, &
206 !$omp tmpdp,num_msol_lgas,den_msol_lgas,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
207 !$omp shared(ijkstart3,ijkend3,d_t,mmax,cartesian_grid,axy,k_of,ep_g,vxf_gs,vxf_ss,momentum_z_eq,am0,model_b,p_scale)
208 !$omp do
209  DO ijk = ijkstart3, ijkend3
210 
211  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk)) THEN !impermeable
212  DO m= 0, mmax
213  d_t(ijk,m) = zero
214  ENDDO
215  cycle
216  ENDIF
217 
218  area_face = merge(one, axy(ijk), cartesian_grid)
219  k = k_of(ijk)
220  ijkt = top_of(ijk)
221  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
222 
223  sum_vxf_gs = zero
224  DO m= 1, mmax
225  epsa(m) = avg_z(ep_s(ijk,m),ep_s(ijkt,m),k)
226 ! Gas - All Solids VolxDrag summation
227  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m)
228  sum_vxf_ss(m) = zero
229  DO l = 1, mmax
230  IF (l .NE. m) THEN
231  lm = funlm(l,m)
232 ! Solid M - All other Solids VolxDrag summation
233  sum_vxf_ss(m) = sum_vxf_ss(m) + vxf_ss(ijk,lm)
234  ENDIF
235  ENDDO
236  ENDDO
237 
238 ! calculating DEN_MGas and NUM_MGas
239  den_mgas = zero
240  num_mgas = zero
241  DO m= 1, mmax
242  IF(momentum_z_eq(m)) THEN
243  num_mgas = num_mgas + (epsa(m)*vxf_gs(ijk,m)/ &
244  ((-am0(ijk,m))+vxf_gs(ijk,m)+sum_vxf_ss(m) + &
245  small_number))
246  den_mgas = den_mgas + (vxf_gs(ijk,m)* &
247  ((-am0(ijk,m))+sum_vxf_ss(m))/ &
248  ((-am0(ijk,m))+vxf_gs(ijk,m)+sum_vxf_ss(m)+ &
249  small_number))
250  ELSE
251  den_mgas = den_mgas + vxf_gs(ijk,m)
252  ENDIF
253  ENDDO
254 
255 ! Model B
256 ! -------------------------------------------------------------------
257  IF (model_b) THEN
258 ! Linking velocity correction coefficient to pressure - GAS Phase
259  tmpdp = -am0(ijk,0) + den_mgas
260  IF(abs(tmpdp) > small_number) THEN
261  d_t(ijk,0) = p_scale*area_face/tmpdp
262  ELSE
263  d_t(ijk,0) = zero
264  ENDIF
265 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
266  DO m = 1, mmax
267  IF(momentum_z_eq(m)) THEN
268  tmpdp = -am0(ijk,m) + vxf_gs(ijk,m)
269  IF(abs(tmpdp) > small_number) THEN
270  d_t(ijk,m) = d_t(ijk,0)*vxf_gs(ijk,m)/tmpdp
271  ELSE
272  d_t(ijk,m) = zero
273  ENDIF
274  ELSE
275  d_t(ijk,m) = zero
276  ENDIF
277  ENDDO
278 
279 ! Model A
280 ! -------------------------------------------------------------------
281  ELSE
282 
283 ! Linking velocity correction coefficient to pressure - GAS Phase
284  tmpdp = -am0(ijk,0)+den_mgas
285  IF(abs(tmpdp) > small_number) THEN
286  d_t(ijk,0) = p_scale*area_face*(epga+num_mgas)/tmpdp
287  ELSE
288  d_t(ijk,0) = zero
289  ENDIF
290 
291  DO m= 1, mmax
292 ! calculating NUM_MSol_LGas and DEN_MSol_LGas
293  num_msol_lgas = vxf_gs(ijk,m)*epga/ &
294  ((-am0(ijk,0))+sum_vxf_gs+small_number)
295  den_msol_lgas = vxf_gs(ijk,m)*( &
296  ((-am0(ijk,0)) + (sum_vxf_gs - vxf_gs(ijk,m)))/ &
297  ((-am0(ijk,0)) + sum_vxf_gs+small_number))
298 ! calculating NUM_MSol_LSol and DEN_MSol_LSol
299  num_msol_lsol(m) = zero
300  den_msol_lsol(m) = zero
301  DO l = 1, mmax
302  IF (l .NE. m) THEN
303  lm = funlm(l,m)
304  sum_vxf_ss_wt_m = zero
305  DO lp = 1, mmax
306  IF( (lp /= l) .AND. (lp /= m) ) THEN
307  lpl = funlm(lp,l)
308 ! Solids L - All other Solids VolxDrag but M summation
309  sum_vxf_ss_wt_m = &
310  sum_vxf_ss_wt_m + vxf_ss(ijk,lpl)
311  ENDIF
312  ENDDO
313 
314  IF(momentum_z_eq(l)) THEN
315  num_msol_lsol(m) = num_msol_lsol(m) + &
316  (vxf_ss(ijk,lm)*epsa(l)/((-am0(ijk,l)) + &
317  vxf_gs(ijk,l)+sum_vxf_ss(l)+ small_number))
318 
319  den_msol_lsol(m) = den_msol_lsol(m) + &
320  vxf_ss(ijk,lm)*(((-am0(ijk,l)) + &
321  vxf_gs(ijk,l)+sum_vxf_ss_wt_m )/ &
322  ((-am0(ijk,l))+vxf_gs(ijk,l)+sum_vxf_ss(l)+ &
323  small_number))
324  ELSE
325  den_msol_lsol(m) = &
326  den_msol_lsol(m) + vxf_ss(ijk,lm)
327  ENDIF
328  ENDIF ! end if (l.ne.m)
329  ENDDO ! end do (l=1,mmax)
330 
331 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
332  IF(momentum_z_eq(m)) THEN
333  tmpdp = -am0(ijk,m) + den_msol_lgas + den_msol_lsol(m)
334  IF(abs(tmpdp) > small_number) THEN
335  d_t(ijk,m) = p_scale*area_face*(epsa(m) + &
336  num_msol_lsol(m) + num_msol_lgas)/tmpdp
337  ELSE
338  d_t(ijk,m) = zero
339  ENDIF
340  ELSE
341  d_t(ijk,m) = zero
342  ENDIF
343  ENDDO ! end do (m=1,mmax)
344 
345  ENDIF !end if/else branch Model_B/Model_A
346  ENDDO ! end do loop (ijk=ijkstart3,ijkend3)
347 !$omp end parallel
348 
349  RETURN
350  END SUBROUTINE calc_d_t_gas_and_solids
351 
352 
353 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
354 ! !
355 ! Subroutine: CALC_D_T_GAS_ONLY !
356 ! Author: M. Syamlal Date: 21-JUN-96 !
357 ! !
358 ! Purpose: Calculate coefficients linking velocity correction to !
359 ! pressure correction -- Top !
360 ! !
361 ! CASE: Only the gas phase (M=0) Z-momentnum equation is solved; the !
362 ! solids Z-momentum equations are NOT solved. This routine is where a !
363 ! coupled DEM simulation should be directed for proper evaluation of !
364 ! pressure correction terms. !
365 ! !
366 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
367  SUBROUTINE calc_d_t_gas_only(AM0, VXF_GS, D_T)
368 
369 ! Global Variables:
370 !---------------------------------------------------------------------//
371 ! Volume fractions of gas and solids phases.
372  use fldvar, only: ep_g
373 ! Number of solids phases.
374  use physprop, only: mmax
375 ! Flag: Use MODEL_B
376  use run, only: model_b
377 ! Pressure scale factor
378  use scales, only: p_scale
379 ! Flag: Cartesian grid simulation
380  use cutcell, only: cartesian_grid
381 ! Area of cell's XY face.
382  use geometry, only: axy
383 ! Volume of W-momentum cell.
384  use geometry, only: vol_w
385 ! Function to average across Z face.
386  use fun_avg, only: avg_z
387 ! Fluid grid loop bounds.
388  use compar, only: ijkstart3, ijkend3
389 ! Flags: Impermeable surface and mass flow at top face, IJK of cell to top
390  use functions, only: ip_at_t, mflow_at_t, top_of
391 ! Indices: K index of cell
392  use indices, only: k_of
393 ! Flag and variables for QMOM implementation.
394  use qmom_kinetic_equation, only: qmomk, qmomk_nn
395  use qmom_kinetic_equation, only: qmomk_f_gs, qmomk_f_gs
396 
397 ! Global Parameters:
398 !---------------------------------------------------------------------//
399 ! Size of IJK arrays and size of solids phase arrays.
400  use param, only: dimension_3, dimension_m
401 ! Double precision parameters.
402  use param1, only: zero, small_number, one
403 
404  IMPLICIT NONE
405 
406 ! Dummy arguments
407 !---------------------------------------------------------------------//
408 ! Temporary variable to store A matrix for local manipulation
409  DOUBLE PRECISION, INTENT(IN) :: AM0(dimension_3, 0:dimension_m)
410 ! Volume x average at momentum cell centers
411  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
412 ! Coefficients for pressure correction equation. These coefficients are
413 ! initialized to zero in the subroutine time_march before the time loop.
414  DOUBLE PRECISION, INTENT(INOUT) :: d_t(dimension_3, 0:dimension_m)
415 
416 
417 ! Local variables:
418 !---------------------------------------------------------------------//
419 ! Usual Indices
420  INTEGER :: M
421  INTEGER :: INN
422  INTEGER :: IJK, IJKT, K
423 ! Average volume fraction at momentum cell centers
424  DOUBLE PRECISION :: EPGA
425 ! local alias for face area
426  DOUBLE PRECISION :: AREA_FACE
427 ! sum of All Solid-Gas VolxDrag
428  DOUBLE PRECISION :: SUM_VXF_GS
429 ! Temp variable for double precision values.
430  DOUBLE PRECISION :: TMPdp
431 
432 !......................................................................!
433 
434 !$omp parallel default(none) &
435 !$omp private(ijk,area_face,k,ijkt,epga,sum_vxf_gs,tmpdp) &
436 !$omp shared(ijkstart3,ijkend3,d_t,mmax,cartesian_grid,axy,k_of,ep_g,vxf_gs,am0,model_b,p_scale,qmomk,vol_w,qmomk_f_gs)
437 !$omp do
438  DO ijk = ijkstart3, ijkend3
439  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk)) THEN
440  d_t(ijk,0) = zero
441  cycle
442  ENDIF
443 
444  area_face = merge(one, axy(ijk), cartesian_grid)
445  k = k_of(ijk)
446  ijkt = top_of(ijk)
447  epga = avg_z(ep_g(ijk),ep_g(ijkt),k)
448 
449  sum_vxf_gs = zero
450  IF (.NOT. qmomk) THEN
451  DO m= 1, mmax
452  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m)
453  ENDDO
454  ELSE
455  DO inn = 1, qmomk_nn
456  DO m = 1, mmax
457  sum_vxf_gs = sum_vxf_gs + vol_w(ijk)*avg_z( &
458  qmomk_f_gs(inn,ijk,m), qmomk_f_gs(inn,ijkt,m),k)
459  ENDDO
460  ENDDO
461  ENDIF
462 
463  tmpdp = -am0(ijk,0)+sum_vxf_gs
464  IF(abs(tmpdp) > small_number) THEN
465  IF (model_b) THEN
466  d_t(ijk,0) = p_scale*area_face/tmpdp
467  ELSE
468  d_t(ijk,0) = p_scale*area_face*epga/tmpdp
469  ENDIF !end if/else branch Model_B/Model_A
470  ELSE
471  d_t(ijk,0) = zero
472  ENDIF
473  ENDDO ! end do (ijk=ijkstart3,ijkend3)
474 !$omp end parallel
475 
476  RETURN
477  END SUBROUTINE calc_d_t_gas_only
478 
479 
480 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
481 ! !
482 ! Subroutine: CALC_D_T_SOLIDS_ONLY !
483 ! Author: M. Syamlal Date: 21-JUN-96 !
484 ! !
485 ! Purpose: Calculate coefficients linking velocity correction to !
486 ! pressure correction -- Top !
487 ! !
488 ! CASE: At least one solids phase Z-momentum equation is solved but !
489 ! the gas phase (m=0) Z-momentum equation is NOT solved. !
490 ! !
491 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
492  SUBROUTINE calc_d_t_solids_only(AM0, VXF_GS, VXF_SS, D_T)
493 
494 ! Global Variables:
495 !---------------------------------------------------------------------//
496 ! Volume fractions of gas and solids phases.
497  use fldvar, only: ep_s
498 ! Number of solids phases.
499  use physprop, only: mmax
500 ! Flag: Solve the Z-Momentum Equations
501  use run, only: momentum_z_eq
502 ! Flag: Use MODEL_B
503  use run, only: model_b
504 ! Pressure scale factor
505  use scales, only: p_scale
506 ! Flag: Cartesian grid simulation
507  use cutcell, only: cartesian_grid
508 ! Area of cell's XY face.
509  use geometry, only: axy
510 ! Function to average across Z face.
511  use fun_avg, only: avg_z
512 ! Fluid grid loop bounds.
513  use compar, only: ijkstart3, ijkend3
514 ! Function to convert to phases to single array, IJK of cell to top
515  use functions, only: funlm, top_of
516 ! Flags: Impermeable surface and mass flow at top face
517  use functions, only: ip_at_t, mflow_at_t
518 ! Indices: K index of cell
519  use indices, only: k_of
520 
521 ! Global Parameters:
522 !---------------------------------------------------------------------//
523 ! Size of IJK arrays and size of solids phase arrays.
524  use param, only: dimension_3, dimension_m
525 ! Size of MxM upper triangular matrix.
526  use param1, only: dimension_lm
527 ! Double precision parameters.
528  use param1, only: zero, small_number, one
529 
530  IMPLICIT NONE
531 
532 ! Dummy arguments
533 !---------------------------------------------------------------------//
534 ! Temporary variable to store A matrix for local manipulation
535  DOUBLE PRECISION, INTENT(IN) :: AM0(dimension_3, 0:dimension_m)
536 ! Volume x average at momentum cell centers
537  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
538 ! Volume x average at momentum cell centers Solid-Solid Drag
539  DOUBLE PRECISION, INTENT(IN) :: VxF_ss(dimension_3, dimension_lm)
540 ! Coefficients for pressure correction equation. These coefficients are
541 ! initialized to zero in the subroutine time_march before the time loop.
542  DOUBLE PRECISION, INTENT(INOUT) :: d_t(dimension_3, 0:dimension_m)
543 
544 ! Local variables:
545 !---------------------------------------------------------------------//
546 ! Usual Indices
547  INTEGER :: LM, M, L, LPL, LP
548  INTEGER :: IJK, IJKT, K
549 ! Average solid volume fraction at momentum cell centers
550  DOUBLE PRECISION :: EPSA(dimension_m)
551 ! local alias for face area
552  DOUBLE PRECISION :: AREA_FACE
553 ! sum of Solid M - All other Solid drag
554  DOUBLE PRECISION :: SUM_VXF_SS(dimension_m)
555 ! sum of Solid L - All other Solid VolxDrag but M
556  DOUBLE PRECISION :: SUM_VXF_SS_wt_M
557 ! component of the total numerator needed for evaluation of pressure
558 ! correction coefficient for SOLIDS phase M. specifically, for when the
559 ! sum L is over solids phase only
560  DOUBLE PRECISION :: NUM_MSol_LSol(dimension_m)
561 ! component of the total denominator needed for evaluation of pressure
562 ! correction coefficient for SOLIDS phase M. specifically, for when the
563 ! sum L is over solids phases only
564  DOUBLE PRECISION :: DEN_MSol_LSol(dimension_m)
565 ! Temp variable for double precision values.
566  DOUBLE PRECISION :: TMPdp
567 
568 !......................................................................!
569 
570 !$omp parallel default(none) &
571 !$omp private(ijk,area_face,k,ijkt,epsa,lm,sum_vxf_ss,tmpdp,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
572 !$omp shared(ijkstart3,ijkend3,d_t,mmax,cartesian_grid,axy,k_of,vxf_gs,vxf_ss,momentum_z_eq,am0,model_b,p_scale)
573 !$omp do
574  DO ijk = ijkstart3, ijkend3
575  IF (ip_at_t(ijk) .OR. mflow_at_t(ijk) .OR. model_b) THEN
576  DO m= 1, mmax
577  d_t(ijk,m) = zero
578  ENDDO
579  cycle
580  ENDIF
581 
582  area_face = merge(one, axy(ijk), cartesian_grid)
583  k = k_of(ijk)
584  ijkt = top_of(ijk)
585 
586  DO m= 1, mmax
587  epsa(m) = avg_z(ep_s(ijk,m),ep_s(ijkt,m),k)
588  sum_vxf_ss(m) = zero
589  DO l = 1, mmax
590  IF (l .NE. m) THEN
591  lm = funlm(l,m)
592 ! Solids M - All other Solids VolxDrag summation
593  sum_vxf_ss(m) = sum_vxf_ss(m) + vxf_ss(ijk,lm)
594  ENDIF
595  ENDDO
596  ENDDO
597 
598  DO m= 1, mmax
599 ! calculating NUM_MSol_LSol and DEN_MSol_LSol
600  num_msol_lsol(m) = zero
601  den_msol_lsol(m) = zero
602  DO l = 1, mmax
603  IF (l .NE. m) THEN
604  lm = funlm(l,m)
605  sum_vxf_ss_wt_m = zero
606  DO lp = 1, mmax
607  IF ( (lp .NE. l) .AND. (lp .NE. m) ) THEN
608  lpl = funlm(lp,l)
609 ! Solids L - All other Solids VolxDrag but M summation
610  sum_vxf_ss_wt_m = sum_vxf_ss_wt_m + vxf_ss(ijk,lpl)
611  ENDIF
612  ENDDO
613  IF (momentum_z_eq(l)) THEN
614  num_msol_lsol(m) = num_msol_lsol(m) + &
615  (vxf_ss(ijk,lm)*epsa(l)/&
616  ((-am0(ijk,l))+vxf_gs(ijk,l)+sum_vxf_ss(l)+&
617  small_number))
618  den_msol_lsol(m) = den_msol_lsol(m) +vxf_ss(ijk,lm)*(&
619  ((-am0(ijk,l))+vxf_gs(ijk,l)+sum_vxf_ss_wt_m )/&
620  ((-am0(ijk,l))+vxf_gs(ijk,l)+sum_vxf_ss(l)+&
621  small_number ))
622  ELSE
623  den_msol_lsol(m) = den_msol_lsol(m) + vxf_ss(ijk,lm)
624  ENDIF
625  ENDIF ! end if (l.ne.m)
626  ENDDO ! end do (l=1,mmax)
627 
628 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
629  IF(momentum_z_eq(m)) THEN
630  tmpdp = (-am0(ijk,m))+vxf_gs(ijk,m)+den_msol_lsol(m)
631 
632  IF(abs(tmpdp) > small_number) THEN
633  d_t(ijk,m) = p_scale*area_face*&
634  (epsa(m) + num_msol_lsol(m))/tmpdp
635  ELSE
636  d_t(ijk,m) = zero
637  ENDIF
638  ELSE
639  d_t(ijk,m) = zero
640  ENDIF
641  ENDDO ! end do (m=1,mmax)
642  ENDDO ! end do (ijk=ijkstart3,ijkend3)
643 !$omp end parallel
644 
645  RETURN
646  END SUBROUTINE calc_d_t_solids_only
geometry::vol_w
double precision, dimension(:), allocatable vol_w
Definition: geometry_mod.f:242
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
fldvar::ep_g
double precision, dimension(:), allocatable ep_g
Definition: fldvar_mod.f:15
param1::dimension_lm
integer dimension_lm
Definition: param1_mod.f:13
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
geometry::axy
double precision, dimension(:), allocatable axy
Definition: geometry_mod.f:210
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
calc_d_t_solids_only
subroutine calc_d_t_solids_only(AM0, VXF_GS, VXF_SS, D_T)
Definition: calc_d_t.f:493
calc_d_t
subroutine calc_d_t(A_M, VXF_GS, VXF_SS, D_T, IER)
Definition: calc_d_t.f:16
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
fun_avg
Definition: fun_avg_mod.f:1
calc_d_t_gas_and_solids
subroutine calc_d_t_gas_and_solids(AM0, VXF_GS, VXF_SS, D_T)
Definition: calc_d_t.f:111
physprop::mmax
integer mmax
Definition: physprop_mod.f:19
param1::small_number
double precision, parameter small_number
Definition: param1_mod.f:24
run
Definition: run_mod.f:13
param
Definition: param_mod.f:2
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
qmom_kinetic_equation
Definition: qmom_kinetic_equation_mod.f:12
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
calc_d_t_gas_only
subroutine calc_d_t_gas_only(AM0, VXF_GS, D_T)
Definition: calc_d_t.f:368
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
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
qmom_kinetic_equation::qmomk_f_gs
double precision, dimension(:,:,:), allocatable qmomk_f_gs
Definition: qmom_kinetic_equation_mod.f:41
qmom_kinetic_equation::qmomk
logical qmomk
Definition: qmom_kinetic_equation_mod.f:48