MFIX  2016-1
calc_d_n.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2 ! !
3 ! Subroutine: CALC_d_n !
4 ! Author: M. Syamlal Date: 21-JUN-96 !
5 ! !
6 ! Purpose: calculate coefficients linking velocity correction to !
7 ! pressure correction -- North !
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_n(A_M, VXF_GS, VXF_SS, D_N, IER)
16 
17 ! Global Variables:
18 !---------------------------------------------------------------------//
19 ! Number of solids phases.
20  use physprop, only: mmax
21 ! Flag: Solve the Y-Momentum Equations
22  use run, only: momentum_y_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(OUT) :: D_N(dimension_3, 0:dimension_m)
52 ! Integer error flag.
53  INTEGER, INTENT(OUT) :: 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_Y_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_y_momentum = any(momentum_y_eq(1:mmax))
82 
83  IF(momentum_y_eq(0))THEN
84  IF(any_solids_y_momentum) THEN
85  CALL calc_d_n_gas_and_solids(am0, vxf_gs, vxf_ss, d_n)
86  ELSE
87  CALL calc_d_n_gas_only(am0, vxf_gs, d_n)
88  ENDIF
89 
90  ELSEIF(any_solids_y_momentum) THEN
91  CALL calc_d_n_solids_only(am0, vxf_gs, vxf_ss, d_n)
92  ENDIF
93 
94  deallocate(am0)
95 
96  RETURN
97  END SUBROUTINE calc_d_n
98 
99 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
100 ! !
101 ! Subroutine: CALC_D_N_GAS_AND_SOLIDS !
102 ! Author: M. Syamlal Date: 21-JUN-96 !
103 ! !
104 ! Purpose: calculate coefficients linking velocity correction to !
105 ! pressure correction -- North !
106 ! !
107 ! CASE: The gas phase Y-momentum equation is solved and at least one !
108 ! solids phase (m>0) y-momentum equation is solved. !
109 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
110  SUBROUTINE calc_d_n_gas_and_solids(AM0, VXF_GS, VXF_SS, D_N)
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 Y-Momentum Equations
119  use run, only: momentum_y_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 XZ face.
127  use geometry, only: axz
128 ! Function to average across Y face.
129  use fun_avg, only: avg_y
130 ! Fluid grid loop bounds.
131  use compar, only: ijkstart3, ijkend3
132 ! Function to convert to phases to single array, IJK of cell to north
133  use functions, only: funlm, north_of
134 ! Flags: Impermeable surface and mass flow at north face
135  use functions, only: ip_at_n, mflow_at_n
136 ! Indices: J index of cell
137  use indices, only: j_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_N(dimension_3, 0:dimension_m)
161 
162 ! Local variables:
163 !---------------------------------------------------------------------//
164 ! Usual Indices
165  INTEGER :: LM, M, L, LPL, LP
166  INTEGER :: J, IJK, IJKN
167 ! Average volume fraction at momentum cell centers
168  DOUBLE PRECISION :: EPGA
169 ! Average solid volume fraction at momentum cell centers
170  DOUBLE PRECISION :: EPSA(dimension_m)
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 !$omp private(ijk,area_face,j,ijkn,epga,sum_vxf_gs,epsa,lm,sum_vxf_ss,den_mgas,num_mgas, &
205 !$omp tmpdp,num_msol_lgas,den_msol_lgas,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
206 !$omp shared(ijkstart3,ijkend3,d_n,mmax,cartesian_grid,axz,j_of,ep_g,vxf_gs,vxf_ss,momentum_y_eq,am0,model_b,p_scale)
207 !$omp do
208  DO ijk = ijkstart3, ijkend3
209  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk)) THEN
210  d_n(ijk,0:mmax) = zero
211  cycle
212  ENDIF
213 
214  area_face = merge(one, axz(ijk), cartesian_grid)
215  j = j_of(ijk)
216  ijkn = north_of(ijk)
217  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
218 
219  sum_vxf_gs = zero
220  DO m= 1, mmax
221  epsa(m) = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
222 ! Gas - All Solids VolxDrag summation
223  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m)
224  sum_vxf_ss(m) = zero
225  DO l = 1, mmax
226  IF (l .NE. m) THEN
227  lm = funlm(l,m)
228 ! Solids M - All other Solids VolxDrag summation
229  sum_vxf_ss(m) = sum_vxf_ss(m) + vxf_ss(ijk,lm)
230  ENDIF
231  ENDDO
232  ENDDO
233 
234 ! Calculating DEN_MGas and NUM_MGas
235  num_mgas = zero
236  den_mgas = zero
237  DO m= 1, mmax
238  IF(momentum_y_eq(m)) THEN
239  num_mgas = num_mgas + (epsa(m)*vxf_gs(ijk,m) / &
240  ((-am0(ijk,m)) + vxf_gs(ijk,m) + sum_vxf_ss(m) + &
241  small_number))
242  den_mgas = den_mgas + (vxf_gs(ijk,m)* &
243  ((-am0(ijk,m)) + sum_vxf_ss(m)) / &
244  ((-am0(ijk,m)) + vxf_gs(ijk,m) + sum_vxf_ss(m) + &
245  small_number))
246  ELSE
247  den_mgas = den_mgas + vxf_gs(ijk,m)
248  ENDIF
249  ENDDO
250 
251 ! Model B
252 ! -------------------------------------------------------------------
253  IF (model_b) THEN
254 ! Linking velocity correction coefficient to pressure - GAS Phase
255  tmpdp = -am0(ijk,0) + den_mgas
256  IF(abs(tmpdp) > small_number ) THEN
257  d_n(ijk,0) = p_scale*area_face/tmpdp
258  ELSE
259  d_n(ijk,0) = zero
260  ENDIF
261 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
262  DO m = 1, mmax
263  IF(momentum_y_eq(m)) THEN
264  tmpdp = -am0(ijk,m) + vxf_gs(ijk,m)
265  IF(abs(tmpdp) > small_number) THEN
266  d_n(ijk,m) = d_n(ijk,0)*vxf_gs(ijk,m) / tmpdp
267  ELSE
268  d_n(ijk,m) = zero
269  ENDIF
270  ELSE
271  d_n(ijk,m) = zero
272  ENDIF
273  ENDDO
274 
275 ! Model A
276 ! -------------------------------------------------------------------
277  ELSE
278 
279 ! Linking velocity correction coefficient to pressure - GAS Phase
280  tmpdp = -am0(ijk,0) + den_mgas
281  IF(abs(tmpdp) > small_number) THEN
282  d_n(ijk,0) = p_scale*area_face*(epga+num_mgas) / tmpdp
283  ELSE
284  d_n(ijk,0) = zero
285  ENDIF
286 
287  DO m= 1, mmax
288 ! calculating NUM_MSol_LGas and DEN_MSol_LGas
289  num_msol_lgas = vxf_gs(ijk,m)*epga/ &
290  ((-am0(ijk,0))+sum_vxf_gs+small_number)
291  den_msol_lgas = vxf_gs(ijk,m)*( &
292  ((-am0(ijk,0))+(sum_vxf_gs - vxf_gs(ijk,m))) / &
293  ((-am0(ijk,0))+sum_vxf_gs + small_number) )
294 
295 ! calculating NUM_MSol_LSol and DEN_MSol_LSol
296  num_msol_lsol(m) = zero
297  den_msol_lsol(m) = zero
298  DO l = 1, mmax
299  IF (l .NE. m) THEN
300  lm = funlm(l,m)
301  sum_vxf_ss_wt_m = zero
302  DO lp = 1, mmax
303  IF((lp /= l) .AND. (lp /= m)) THEN
304  lpl = funlm(lp,l)
305 ! Solids L - All other Solids VolxDrag but M summation
306  sum_vxf_ss_wt_m = sum_vxf_ss_wt_m + &
307  vxf_ss(ijk,lpl)
308  ENDIF
309  ENDDO
310 
311  IF(momentum_y_eq(l)) THEN
312  num_msol_lsol(m) = num_msol_lsol(m) + &
313  (vxf_ss(ijk,lm)*epsa(l) / &
314  ((-am0(ijk,l))+vxf_gs(ijk,l)+sum_vxf_ss(l)+ &
315  small_number ))
316 
317  den_msol_lsol(m) = den_msol_lsol(m) + &
318  vxf_ss(ijk,lm)*(((-am0(ijk,l)) + &
319  vxf_gs(ijk,l) + sum_vxf_ss_wt_m ) / &
320  ((-am0(ijk,l)) + vxf_gs(ijk,l)+ &
321  sum_vxf_ss(l) + small_number) )
322  ELSE
323  den_msol_lsol(m) = den_msol_lsol(m) + &
324  vxf_ss(ijk,lm)
325  ENDIF
326  ENDIF ! end if (l.ne.m)
327  ENDDO ! end do (l=1,mmax)
328 
329 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
330  IF(momentum_y_eq(m)) THEN
331  tmpdp = -am0(ijk,m) + den_msol_lgas + den_msol_lsol(m)
332 
333  IF(abs(tmpdp) > small_number) THEN
334  d_n(ijk,m) = p_scale*area_face*(epsa(m) + &
335  num_msol_lsol(m) + num_msol_lgas) / tmpdp
336  ELSE
337  d_n(ijk,m) = zero
338  ENDIF
339  ELSE
340  d_n(ijk,m) = zero
341  ENDIF
342  ENDDO ! end do (m=1,mmax)
343 
344  ENDIF !end if/else branch Model_B/Model_A
345  ENDDO ! end do loop (ijk=ijkstart3,ijkend3)
346 !$omp end parallel
347 
348  RETURN
349  END SUBROUTINE calc_d_n_gas_and_solids
350 
351 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
352 ! !
353 ! Subroutine: CALC_D_N_GAS_ONLY !
354 ! Author: M. Syamlal Date: 21-JUN-96 !
355 ! !
356 ! Purpose: Calculate coefficients linking velocity correction to !
357 ! pressure correction -- North !
358 ! !
359 ! CASE: Only the gas phase (M=0) Y-momentnum equation is solved; the !
360 ! solids Y-momentum equations are NOT solved. This routine is where a !
361 ! coupled DEM simulation should be directed for proper evaluation of !
362 ! pressure correction terms. !
363 ! !
364 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
365  SUBROUTINE calc_d_n_gas_only(AM0, VXF_GS, D_N)
366 
367 ! Global Variables:
368 !---------------------------------------------------------------------//
369 ! Volume fractions of gas and solids phases.
370  use fldvar, only: ep_g
371 ! Number of solids phases.
372  use physprop, only: mmax
373 ! Flag: Use MODEL_B
374  use run, only: model_b
375 ! Pressure scale factor
376  use scales, only: p_scale
377 ! Flag: Cartesian grid simulation
378  use cutcell, only: cartesian_grid
379 ! Area of cell's XZ face.
380  use geometry, only: axz
381 ! Volume of V-momentum cell.
382  use geometry, only: vol_v
383 ! Function to average across Y face.
384  use fun_avg, only: avg_y
385 ! Fluid grid loop bounds.
386  use compar, only: ijkstart3, ijkend3
387 ! Flags: Impermeable surface and mass flow at north face, IJK of cell to north
388  use functions, only: ip_at_n, mflow_at_n, north_of
389 ! Indices: J index of cell
390  use indices, only: j_of
391 ! Flag and variables for QMOM implementation.
392  use qmom_kinetic_equation, only: qmomk, qmomk_nn
393  use qmom_kinetic_equation, only: qmomk_f_gs, qmomk_f_gs
394 
395 ! Global Parameters:
396 !---------------------------------------------------------------------//
397 ! Size of IJK arrays and size of solids phase arrays.
398  use param, only: dimension_3, dimension_m
399 ! Double precision parameters.
400  use param1, only: zero, small_number, one
401 
402  IMPLICIT NONE
403 
404 ! Dummy arguments
405 !---------------------------------------------------------------------//
406 ! Temporary variable to store A matrix for local manipulation
407  DOUBLE PRECISION, INTENT(IN) :: AM0(dimension_3, 0:dimension_m)
408 ! Volume x average at momentum cell centers
409  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
410 ! Coefficients for pressure correction equation. These coefficients are
411 ! initialized to zero in the subroutine time_march before the time loop
412  DOUBLE PRECISION, INTENT(INOUT) :: d_n(dimension_3, 0:dimension_m)
413 
414 ! Local variables:
415 !---------------------------------------------------------------------//
416 ! Usual Indices
417  INTEGER :: M
418  INTEGER :: INN
419  INTEGER :: J, IJK, IJKN
420 ! Average volume fraction at momentum cell centers
421  DOUBLE PRECISION :: EPGA
422 ! local alias for face area
423  DOUBLE PRECISION :: AREA_FACE
424 ! sum of All Solid-Gas VolxDrag
425  DOUBLE PRECISION :: SUM_VXF_GS
426 ! Temp variable for double precision values.
427  DOUBLE PRECISION :: TMPdp
428 
429 !......................................................................!
430 
431 !$omp parallel default(none) &
432 !$omp private(ijk,area_face,j,ijkn,epga,sum_vxf_gs,tmpdp) &
433 !$omp shared(ijkstart3,ijkend3,d_n,mmax,cartesian_grid,axz,j_of,ep_g,vxf_gs,am0,model_b,p_scale,qmomk,vol_v,qmomk_f_gs)
434 !$omp do
435  DO ijk = ijkstart3, ijkend3
436 
437  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk)) THEN
438  d_n(ijk,0) = zero
439  cycle
440  ENDIF
441 
442  area_face = merge(one, axz(ijk), cartesian_grid)
443  j = j_of(ijk)
444  ijkn = north_of(ijk)
445  epga = avg_y(ep_g(ijk),ep_g(ijkn),j)
446 
447  sum_vxf_gs = zero
448  IF (.NOT. qmomk) THEN
449  DO m= 1, mmax
450 ! Gas - All Solids VolxDrag summation
451  sum_vxf_gs = sum_vxf_gs + vxf_gs(ijk,m)
452  ENDDO
453  ELSE
454  DO inn = 1, qmomk_nn
455  DO m = 1, mmax
456  sum_vxf_gs = sum_vxf_gs + vol_v(ijk)*avg_y( &
457  qmomk_f_gs(inn,ijk,m),qmomk_f_gs(inn,ijkn,m),j)
458  ENDDO
459  ENDDO
460  ENDIF
461 
462  tmpdp = -am0(ijk,0) + sum_vxf_gs
463  IF(abs(tmpdp) > small_number) THEN
464  IF (model_b) THEN
465  d_n(ijk,0) = p_scale*area_face/tmpdp
466  ELSE
467  d_n(ijk,0) = p_scale*area_face*epga/tmpdp
468  ENDIF !end if/else branch Model_B/Model_A
469  ELSE
470  d_n(ijk,0) = zero
471  ENDIF
472 
473  ENDDO ! end do (ijk=ijkstart3,ijkend3)
474 !$omp end parallel
475 
476  RETURN
477  END SUBROUTINE calc_d_n_gas_only
478 
479 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
480 ! !
481 ! Subroutine: CALC_D_N_SOLIDS_ONLY !
482 ! Author: M. Syamlal Date: 21-JUN-96 !
483 ! !
484 ! Purpose: calculate coefficients linking velocity correction to !
485 ! pressure correction -- North !
486 ! !
487 ! CASE: At least one solids phase Y-momentum equation is solved but !
488 ! the gas phase (m=0) Y-momentum equation is NOT solved. !
489 ! !
490 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
491  SUBROUTINE calc_d_n_solids_only(AM0, VXF_GS, VXF_SS, D_N)
492 
493 ! Global Variables:
494 !---------------------------------------------------------------------//
495 ! Volume fractions of gas and solids phases.
496  use fldvar, only: ep_s
497 ! Number of solids phases.
498  use physprop, only: mmax
499 ! Flag: Solve the Y-Momentum Equations
500  use run, only: momentum_y_eq
501 ! Flag: Use MODEL_B
502  use run, only: model_b
503 ! Pressure scale factor
504  use scales, only: p_scale
505 ! Flag: Cartesian grid simulation
506  use cutcell, only: cartesian_grid
507 ! Area of cell's XZ face.
508  use geometry, only: axz
509 ! Function to average across Y face.
510  use fun_avg, only: avg_y
511 ! Fluid grid loop bounds.
512  use compar, only: ijkstart3, ijkend3
513 ! Function to convert to phases to single array, IJK of cell to north
514  use functions, only: funlm, north_of
515 ! Flags: Impermeable surface and mass flow at north face
516  use functions, only: ip_at_n, mflow_at_n
517 ! Indices: J index of cell
518  use indices, only: j_of
519 
520 ! Global Parameters:
521 !---------------------------------------------------------------------//
522 ! Size of IJK arrays and size of solids phase arrays.
523  use param, only: dimension_3, dimension_m
524 ! Size of MxM upper triangular matrix.
525  use param1, only: dimension_lm
526 ! Double precision parameters.
527  use param1, only: zero, small_number, one
528 
529  IMPLICIT NONE
530 
531 ! Dummy arguments
532 !---------------------------------------------------------------------//
533 ! Temporary variable to store A matrix for local manipulation
534  DOUBLE PRECISION, INTENT(IN) :: AM0(dimension_3, 0:dimension_m)
535 ! Volume x average at momentum cell centers
536  DOUBLE PRECISION, INTENT(IN) :: VxF_gs(dimension_3, dimension_m)
537 ! Volume x average at momentum cell centers Solid-Solid Drag
538  DOUBLE PRECISION, INTENT(IN) :: VxF_ss(dimension_3, dimension_lm)
539 ! Coefficients for pressure correction equation. These coefficients are
540 ! initialized to zero in the subroutine time_march before the time loop
541  DOUBLE PRECISION, INTENT(INOUT) :: D_N(dimension_3, 0:dimension_m)
542 
543 ! Local variables:
544 !---------------------------------------------------------------------//
545 ! Usual Indices
546  INTEGER :: LM, M, L, LPL, LP
547  INTEGER :: J, IJK, IJKN
548 ! Average solid volume fraction at momentum cell centers
549  DOUBLE PRECISION :: EPSA(dimension_m)
550 ! local alias for face area
551  DOUBLE PRECISION :: AREA_FACE
552 ! sum of Solid M - All other Solid drag
553  DOUBLE PRECISION :: SUM_VXF_SS(dimension_m)
554 ! sum of Solid L - All other Solid VolxDrag but M
555  DOUBLE PRECISION :: SUM_VXF_SS_wt_M
556 ! component of the total numerator needed for evaluation of pressure
557 ! correction coefficient for SOLIDS phase M. specifically, for when the
558 ! sum L is over solids phase only
559  DOUBLE PRECISION :: NUM_MSol_LSol(dimension_m)
560 ! component of the total denominator needed for evaluation of pressure
561 ! correction coefficient for SOLIDS phase M. specifically, for when the
562 ! sum L is over solids phases only
563  DOUBLE PRECISION :: DEN_MSol_LSol(dimension_m)
564 ! Temp variable for double precision values.
565  DOUBLE PRECISION :: TMPdp
566 
567 !......................................................................!
568 
569 !$omp parallel default(none) &
570 !$omp private(ijk,area_face,j,ijkn,epsa,lm,sum_vxf_ss,tmpdp,den_msol_lsol,num_msol_lsol,sum_vxf_ss_wt_m,lpl) &
571 !$omp shared(ijkstart3,ijkend3,d_n,mmax,cartesian_grid,axz,j_of,vxf_gs,vxf_ss,momentum_y_eq,am0,model_b,p_scale)
572 !$omp do
573  DO ijk = ijkstart3, ijkend3
574  IF (ip_at_n(ijk) .OR. mflow_at_n(ijk) .OR. model_b) THEN
575  d_n(ijk,1:mmax) = zero
576  cycle
577  ENDIF
578 
579  area_face = merge(one, axz(ijk), cartesian_grid)
580  j = j_of(ijk)
581  ijkn = north_of(ijk)
582 
583  DO m= 1, mmax
584  epsa(m) = avg_y(ep_s(ijk,m),ep_s(ijkn,m),j)
585  sum_vxf_ss(m) = zero
586  DO l = 1, mmax
587  IF (l .NE. m) THEN
588  lm = funlm(l,m)
589 ! Solids M - All other Solids VolxDrag summation
590  sum_vxf_ss(m) = sum_vxf_ss(m) + vxf_ss(ijk,lm)
591  ENDIF
592  ENDDO
593  ENDDO
594 
595  DO m= 1, mmax
596 ! calculating NUM_MSol_LSol and DEN_MSol_LSol
597  num_msol_lsol(m) = zero
598  den_msol_lsol(m) = zero
599  DO l = 1, mmax
600  IF (l .NE. m) THEN
601  lm = funlm(l,m)
602  sum_vxf_ss_wt_m = zero
603  DO lp = 1, mmax
604  IF((lp /= l) .AND. (lp /= m)) THEN
605  lpl = funlm(lp,l)
606 ! Solids L - All other Solids VolxDrag but M summation
607  sum_vxf_ss_wt_m=sum_vxf_ss_wt_m+vxf_ss(ijk,lpl)
608  ENDIF
609  ENDDO
610  IF(momentum_y_eq(l)) THEN
611  num_msol_lsol(m) = num_msol_lsol(m) + &
612  (vxf_ss(ijk,lm)*epsa(l)/((-am0(ijk,l)) + &
613  vxf_gs(ijk,l)+sum_vxf_ss(l)+small_number ))
614  den_msol_lsol(m) = den_msol_lsol(m) + &
615  vxf_ss(ijk,lm)*(((-am0(ijk,l)) + vxf_gs(ijk,l)+&
616  sum_vxf_ss_wt_m )/((-am0(ijk,l))+vxf_gs(ijk,l)+&
617  sum_vxf_ss(l) + small_number))
618  ELSE
619  den_msol_lsol(m) = den_msol_lsol(m)+vxf_ss(ijk,lm)
620  ENDIF
621  ENDIF ! end if (l.ne.m)
622  ENDDO ! end do (l=1,mmax)
623 
624 ! Linking velocity correction coefficient to pressure - SOLIDs Phase
625  IF (momentum_y_eq(m)) THEN
626  tmpdp = -am0(ijk,m) + vxf_gs(ijk,m) + den_msol_lsol(m)
627  IF(abs(tmpdp) > small_number) THEN
628  d_n(ijk,m) = p_scale*area_face* (epsa(m) + &
629  num_msol_lsol(m))/tmpdp
630  ELSE
631  d_n(ijk,m) = zero
632  ENDIF
633  ELSE
634  d_n(ijk,m) = zero
635  ENDIF
636  ENDDO ! end do (m=1,mmax)
637  ENDDO ! end do (ijk=ijkstart3,ijkend3)
638 !$omp end parallel
639 
640  RETURN
641  END SUBROUTINE calc_d_n_solids_only
run::momentum_y_eq
logical, dimension(0:dim_m) momentum_y_eq
Definition: run_mod.f:77
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
scales::p_scale
double precision p_scale
Definition: scales_mod.f:13
calc_d_n
subroutine calc_d_n(A_M, VXF_GS, VXF_SS, D_N, IER)
Definition: calc_d_n.f:16
calc_d_n_gas_and_solids
subroutine calc_d_n_gas_and_solids(AM0, VXF_GS, VXF_SS, D_N)
Definition: calc_d_n.f:111
fun_avg
Definition: fun_avg_mod.f:1
calc_d_n_gas_only
subroutine calc_d_n_gas_only(AM0, VXF_GS, D_N)
Definition: calc_d_n.f:366
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
calc_d_n_solids_only
subroutine calc_d_n_solids_only(AM0, VXF_GS, VXF_SS, D_N)
Definition: calc_d_n.f:492
run
Definition: run_mod.f:13
geometry::axz
double precision, dimension(:), allocatable axz
Definition: geometry_mod.f:208
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
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
geometry::vol_v
double precision, dimension(:), allocatable vol_v
Definition: geometry_mod.f:233