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