File: N:\mfix\model\source_pp_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOURCE_Pp_g                                             C
4     !  Purpose: Determine source terms for pressure correction equation.   C
5     !                                                                      C
6     !  Notes: The off-diagonal coefficients are positive. The center       C
7     !         coefficient and the source vector are negative. See          C
8     !         conv_Pp_g                                                    C
9     !                                                                      C
10     !  Author: M. Syamlal                                 Date: 21-JUN-96  C
11     !                                                                      C
12     !  Revision Number: 1                                                  C
13     !  Purpose: To incorporate Cartesian grid modifications                C
14     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
15     !                                                                      C
16     !                                                                      C
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
18           SUBROUTINE SOURCE_PP_G(A_M, B_M, B_MMAX)
19     
20     ! Modules
21     !---------------------------------------------------------------------//
22           USE bc, ONLY: IJK_P_G
23           USE compar, ONLY: IJKSTART3, IJKEND3
24           USE cutcell, ONLY: CARTESIAN_GRID, A_UPG_E, A_VPG_N, A_WPG_T
25           USE fldvar, ONLY: U_G, V_G, W_G, U_S, V_S, W_S, ROP_G, ROP_S
26           USE fldvar, only: ROP_GO, ROP_SO, RO_G
27           USE geometry, ONLY: VOL
28           USE param, only: dimension_3, dimension_m
29           USE param, only: east, west, south, north, top, bottom
30           USE param1, ONLY: SMALL_NUMBER, ONE, ZERO, UNDEFINED
31           USE pgcor, ONLY: D_E, D_N, D_T
32           USE physprop, ONLY: CLOSE_PACKED, MMAX, RO_G0
33           USE run, ONLY: SHEAR, ODT, UNDEFINED_I
34           USE rxns, ONLY: SUM_R_G, SUM_R_S
35           USE vshear, ONLY: VSH
36           IMPLICIT NONE
37     
38     ! Dummy arguments
39     !---------------------------------------------------------------------//
40     ! Septadiagonal matrix A_m
41           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
42     ! Vector b_m
43           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
44     ! maximum term in b_m expression
45           DOUBLE PRECISION, INTENT(INOUT) :: B_mmax(DIMENSION_3, 0:DIMENSION_M)
46     
47     ! Local Variables
48     !---------------------------------------------------------------------//
49     ! solids phase index
50           INTEGER :: M
51     ! Indices
52           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
53     ! terms of bm expression
54           DOUBLE PRECISION bma, bme, bmw, bmn, bms, bmt, bmb, bmr
55     ! error message
56           CHARACTER(LEN=80) :: LINE(1)
57     !---------------------------------------------------------------------//
58     ! loezos
59     ! update to true velocity
60           IF (SHEAR) THEN
61     !!$omp parallel do private(IJK)
62              DO IJK = IJKSTART3, IJKEND3
63                 IF (FLUID_AT(IJK)) THEN
64                    V_G(IJK)=V_G(IJK)+VSH(IJK)
65                 ENDIF
66              ENDDO
67           ENDIF
68     
69     ! Calculate convection-diffusion fluxes through each of the faces
70     
71     !$omp parallel default(none) &
72     !$omp          private(IJK, IMJK, IJMK, IJKM, M, bma,bme,bmw,bmn,&
73     !$omp                  bms,bmt,bmb,bmr,line)  &
74     !$omp          shared(ijkstart3,ijkend3,cartesian_grid,rop_g,rop_go,&
75     !$omp                 rop_s,rop_so,vol,odt,u_g,v_g,w_g,u_s,v_s,w_s,b_m, &
76     !$omp                 b_mmax,d_e,d_n,d_t,a_m,a_upg_e,a_vpg_n,a_wpg_t,&
77     !$omp                 mmax,close_packed,sum_r_s,sum_r_g,ro_g0)
78     !$omp do
79           DO IJK = ijkstart3, ijkend3
80              IF (FLUID_AT(IJK)) THEN
81                 IMJK = IM_OF(IJK)
82                 IJMK = JM_OF(IJK)
83                 IJKM = KM_OF(IJK)
84     
85                 bma = (ROP_G(IJK)-ROP_GO(IJK))*VOL(IJK)*ODT
86                 bme = A_M(IJK,east,0)*U_G(IJK)
87                 bmw = A_M(IJK,west,0)*U_G(IMJK)
88                 bmn = A_M(IJK,north,0)*V_G(IJK)
89                 bms = A_M(IJK,south,0)*V_G(IJMK)
90                 bmt = A_M(IJK,top,0)*W_G(IJK)
91                 bmb = A_M(IJK,bottom,0)*W_G(IJKM)
92                 bmr = SUM_R_G(IJK)*VOL(IJK)
93                 B_M(IJK,0) = -((-(bma + bme - bmw + &
94                                   bmn - bms + bmt - bmb ))+ bmr )
95                 B_MMAX(IJK,0) = max(abs(bma),abs(bme),abs(bmw),&
96                    abs(bmn),abs(bms),abs(bmt),abs(bmb),abs(bmr) )
97     
98                 A_M(IJK,east,0) = A_M(IJK,east,0)*D_E(IJK,0)
99                 A_M(IJK,west,0) = A_M(IJK,west,0)*D_E(IMJK,0)
100                 A_M(IJK,north,0) = A_M(IJK,north,0)*D_N(IJK,0)
101                 A_M(IJK,south,0) = A_M(IJK,south,0)*D_N(IJMK,0)
102                 A_M(IJK,top,0) = A_M(IJK,top,0)*D_T(IJK,0)
103                 A_M(IJK,bottom,0) = A_M(IJK,bottom,0)*D_T(IJKM,0)
104     
105                 IF(CARTESIAN_GRID) THEN
106                    A_M(IJK,east,0) = A_M(IJK,east,0) * A_UPG_E(IJK)
107                    A_M(IJK,west,0) = A_M(IJK,west,0) * A_UPG_E(IMJK)
108                    A_M(IJK,north,0) = A_M(IJK,north,0) * A_VPG_N(IJK)
109                    A_M(IJK,south,0) = A_M(IJK,south,0) * A_VPG_N(IJMK)
110                    A_M(IJK,top,0) = A_M(IJK,top,0) * A_WPG_T(IJK)
111                    A_M(IJK,bottom,0) = A_M(IJK,bottom,0) * A_WPG_T(IJKM)
112                 ENDIF
113     
114     ! solids phase terms are needed for those solids phases that do not
115     ! become close packed. these terms are used to essentially make the
116     ! matrix equation for gas pressure correction a mixture pressure
117     ! correction equation by adding solids phases that do not become
118     ! close packed
119                 DO M = 1, MMAX
120                    IF (.NOT.CLOSE_PACKED(M)) THEN
121                       B_M(IJK,0) = B_M(IJK,0) - &
122                          ((-((ROP_S(IJK,M)-ROP_SO(IJK,M))*VOL(IJK)*ODT+&
123                          A_M(IJK,east,M)*U_S(IJK,M)-A_M(IJK,west,M)*U_S(IMJK,M)+&
124                          A_M(IJK,north,M)*V_S(IJK,M)-A_M(IJK,south,M)*V_S(IJMK,M)+&
125                          A_M(IJK,top,M)*W_S(IJK,M)-A_M(IJK,bottom,M)*W_S(IJKM,M)))+&
126                          SUM_R_S(IJK,M)*VOL(IJK))
127     
128                       IF(.NOT.CARTESIAN_GRID) THEN
129                          A_M(IJK,east,0) = A_M(IJK,east,0) + A_M(IJK,east,M)*D_E(IJK,M)
130                          A_M(IJK,west,0) = A_M(IJK,west,0) + A_M(IJK,west,M)*D_E(IMJK,M)
131                          A_M(IJK,north,0) = A_M(IJK,north,0) + A_M(IJK,north,M)*D_N(IJK,M)
132                          A_M(IJK,south,0) = A_M(IJK,south,0) + A_M(IJK,south,M)*D_N(IJMK,M)
133                          A_M(IJK,top,0) = A_M(IJK,top,0) + A_M(IJK,top,M)*D_T(IJK,M)
134                          A_M(IJK,bottom,0) = A_M(IJK,bottom,0) + A_M(IJK,bottom,M)*D_T(IJKM,M)
135                       ELSE
136                          A_M(IJK,east,0) = A_M(IJK,east,0) + A_M(IJK,east,M)*D_E(IJK,M)  * A_UPG_E(IJK)
137                          A_M(IJK,west,0) = A_M(IJK,west,0) + A_M(IJK,west,M)*D_E(IMJK,M) * A_UPG_E(IMJK)
138                          A_M(IJK,north,0) = A_M(IJK,north,0) + A_M(IJK,north,M)*D_N(IJK,M)  * A_VPG_N(IJK)
139                          A_M(IJK,south,0) = A_M(IJK,south,0) + A_M(IJK,south,M)*D_N(IJMK,M) * A_VPG_N(IJMK)
140                          A_M(IJK,top,0) = A_M(IJK,top,0) + A_M(IJK,top,M)*D_T(IJK,M)  * A_WPG_T(IJK)
141                          A_M(IJK,bottom,0) = A_M(IJK,bottom,0) + A_M(IJK,bottom,M)*D_T(IJKM,M) * A_WPG_T(IJKM)
142                       ENDIF
143                    ENDIF
144                 ENDDO
145     
146                 A_M(IJK,0,0) = -(A_M(IJK,east,0)+A_M(IJK,west,0)+&
147                                  A_M(IJK,north,0)+A_M(IJK,south,0)+&
148                                  A_M(IJK,top,0)+A_M(IJK,bottom,0))
149     
150                 IF (ABS(A_M(IJK,0,0)) < SMALL_NUMBER) THEN
151                    IF (ABS(B_M(IJK,0)) < SMALL_NUMBER) THEN
152                       A_M(IJK,0,0) = -ONE
153                       B_M(IJK,0) = ZERO
154                    ELSEIF (RO_G0 .NE. UNDEFINED) THEN !This is an error only in incompressible flow
155     !!$omp             critical
156                       WRITE (LINE, '(A,I6,A,I1,A,G12.5)') 'Error: At IJK = ', IJK, &
157                          ' M = ', 0, ' A = 0 and b = ', B_M(IJK,0)
158                       CALL WRITE_ERROR ('SOURCE_Pp_g', LINE, 1)
159     !!$omp             end critical
160                    ENDIF
161                 ENDIF
162     
163     
164              ELSE   ! if/else branch .not.fluid_at(ijk)
165     ! set the value (correction) in all wall and flow boundary cells to zero
166     ! note the matrix coefficients and source vector should already be zero
167     ! from the initialization of A and B but the following ensures the zero
168     ! value.
169                 A_M(IJK,east,0) = ZERO
170                 A_M(IJK,west,0) = ZERO
171                 A_M(IJK,north,0) = ZERO
172                 A_M(IJK,south,0) = ZERO
173                 A_M(IJK,top,0) = ZERO
174                 A_M(IJK,bottom,0) = ZERO
175                 A_M(IJK,0,0) = -ONE
176                 B_M(IJK,0) = ZERO
177              ENDIF   ! end if/else branch fluid_at(ijk)
178           ENDDO    ! end do loop (ijk=ijkstart3,ijkend3)
179     !$omp end parallel
180     
181     
182     ! loezos
183           IF (SHEAR) THEN
184     !!$omp parallel do private(IJK)
185              DO IJK = IJKSTART3, IJKEND3
186                 IF (FLUID_AT(IJK)) THEN
187                    V_G(IJK)=V_G(IJK)-VSH(IJK)
188                 ENDIF
189              ENDDO
190           ENDIF
191     
192     ! Modification for compressibility 
193           IF (RO_G0 == UNDEFINED) CALL COMPRESSIBLE_PP_G(A_M)
194     
195     ! Remove the asymmetry in matrix caused by the pressure outlet or inlet
196     ! boundaries.  Because the P' at such boundaries is zero we may set the
197     ! coefficient in the neighboring fluid cell to zero without affecting
198     ! the linear equation set.
199     !!$omp    parallel do                                                     &
200     !!$omp&   private(IJK,IMJK, IPJK, IJMK, IJPK, IJKM, IJKP)
201           DO IJK = ijkstart3, ijkend3
202              IF (FLUID_AT(IJK)) THEN
203                 IMJK = IM_OF(IJK)
204                 IPJK = IP_OF(IJK)
205                 IJMK = JM_OF(IJK)
206                 IJPK = JP_OF(IJK)
207                 IJKM = KM_OF(IJK)
208                 IJKP = KP_OF(IJK)
209     ! Cutting the neighbor link between fluid cell and adjacent p_flow_at cell
210                 if(p_flow_at(imjk)) A_m(IJK, west, 0) = ZERO
211                 if(p_flow_at(ipjk)) A_m(IJK, east, 0) = ZERO
212                 if(p_flow_at(ijmk)) A_m(IJK, south, 0) = ZERO
213                 if(p_flow_at(ijpk)) A_m(IJK, north, 0) = ZERO
214                 if(p_flow_at(ijkm)) A_m(IJK, bottom, 0) = ZERO
215                 if(p_flow_at(ijkp)) A_m(IJK, top, 0) = ZERO
216              ENDIF
217           ENDDO
218     
219     ! Specify P' to zero for incompressible flows. Check set_bc0
220     ! for details on selection of IJK_P_g.
221           IF (IJK_P_G /= UNDEFINED_I) THEN
222              B_M(IJK_P_G,0) = ZERO
223              A_M(IJK_P_G,:,0) = ZERO
224              A_M(IJK_P_G,0,0) = -ONE
225           ENDIF
226     
227           RETURN
228     
229     CONTAINS
230     
231           INCLUDE 'functions.inc'
232     
233     END SUBROUTINE SOURCE_PP_G
234     
235     
236     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
237     !                                                                      C
238     !  Subroutine: COMPRESSIBLE_PP_G                                       C
239     !  Purpose: Correction for incompressible flows.                       C
240     !                                                                      C
241     !  Notes: A HS_CORRECT option is available as a better approximation   C
242     !  for high speed flows because it considers density changes in the    C
243     !  neighboring C cells. However, the code runs faster without it for   C
244     !  low speed flows. The gas phase mass balance cannot be maintained    C
245     !  to machine precision with this approximation.                       C
246     !                                                                      C
247     !                                                                      C
248     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
249           SUBROUTINE COMPRESSIBLE_PP_G(A_M)
250     
251     ! Modules
252     !---------------------------------------------------------------------//
253           use compar, only: ijkstart3, ijkend3
254           use eos, only: droodp_g
255           use fldvar, only: ep_g, u_g, v_g, w_g, rop_g, ro_g, p_g
256           use functions, only: fluid_at, im_of, jm_of, km_of
257           use functions, only: east_of, west_of, north_of, south_of
258           use functions, only: top_of, bottom_of
259           use geometry, only: vol, ayz, axz, axy, do_k
260           use param, only: dimension_3, dimension_m
261           use param, only: east, west, south, north, top, bottom
262           use param1, only: one
263           use run, only: discretize, odt
264           use ur_facs, only: ur_fac
265           use xsi
266           IMPLICIT NONE
267     
268     ! Parameters
269     !---------------------------------------------------------------------//
270           LOGICAL, PARAMETER :: HS_CORRECT = .FALSE.
271     
272     ! Dummy arguments
273     !---------------------------------------------------------------------//
274     ! Septadiagonal matrix A_m
275           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
276     
277     ! Local variables
278     !---------------------------------------------------------------------//
279     ! loezos: used for including shearing
280           integer :: incr
281     ! temporary use of global arrays:
282     ! xsi_array: convection weighting factors
283           DOUBLE PRECISION :: XSI_e(DIMENSION_3)
284           DOUBLE PRECISION :: XSI_n(DIMENSION_3)
285           DOUBLE PRECISION :: XSI_t(DIMENSION_3)
286     ! under relaxation factor for pressure
287           double precision :: fac
288     ! Indices
289           INTEGER :: IJK, IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
290           INTEGER :: IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
291     !---------------------------------------------------------------------//
292     
293     ! since p_g = p_g* + ur_fac * pp_g
294           fac = UR_FAC(1)  
295     
296          IF (.NOT.HS_CORRECT) THEN
297     !!$omp    parallel do  &
298     !!$omp&   private(IJK)
299           DO IJK = ijkstart3, ijkend3
300              IF (FLUID_AT(IJK)) THEN
301     ! account for change in density w.r.t. pressure
302                 A_M(IJK,0,0) = A_M(IJK,0,0) - ur_fac(1)*&
303                    DROODP_G(RO_G(IJK),P_G(IJK))*EP_G(IJK)*VOL(IJK)*ODT
304              ENDIF
305           ENDDO
306           ENDIF
307     
308           IF (HS_CORRECT) THEN
309     ! since p_g = p_g* + ur_fac * pp_g
310     ! loezos
311           incr=0
312           CALL CALC_XSI(DISCRETIZE(1),ROP_G,U_G,V_G,W_G,XSI_E,XSI_N,XSI_T,incr)
313     
314     !!$omp    parallel do  &
315     !!$omp&   private(IJK,I,J,K, &
316     !!$omp&            IMJK,IJMK,IJKM,IJKE,IJKW,IJKN,IJKS,IJKT,IJKB)
317           DO IJK = ijkstart3, ijkend3
318              IF (FLUID_AT(IJK)) THEN
319                 IMJK = IM_OF(IJK)
320                 IJMK = JM_OF(IJK)
321                 IJKM = KM_OF(IJK)
322                 IJKE = EAST_OF(IJK)
323                 IJKW = WEST_OF(IJK)
324                 IJKN = NORTH_OF(IJK)
325                 IJKS = SOUTH_OF(IJK)
326                 IJKT = TOP_OF(IJK)
327                 IJKB = BOTTOM_OF(IJK)
328                 A_M(IJK,0,0) = A_M(IJK,0,0) - fac*DROODP_G(RO_G(IJK),P_G(IJK))*&
329                    EP_G(IJK)*((ONE - XSI_E(IJK))*U_G(IJK)*AYZ(IJK)-&
330                                      XSI_E(IMJK)*U_G(IMJK)*AYZ(IMJK)+&
331                                  (ONE-XSI_N(IJK))*V_G(IJK)*AXZ(IJK)-&
332                                       XSI_N(IJMK)*V_G(IJMK)*AXZ(IJMK))
333     
334                 A_M(IJK,east,0) = A_M(IJK,east,0) - EP_G(IJKE)*fac*&
335                    DROODP_G(RO_G(IJKE),P_G(IJKE))*XSI_E(IJK)*U_G(IJK)*AYZ(IJK)
336                 A_M(IJK,west,0) = A_M(IJK,west,0) + EP_G(IJKW)*fac*&
337                    DROODP_G(RO_G(IJKW),P_G(IJKW))*(ONE - XSI_E(IMJK))*U_G(IMJK)*AYZ(IMJK)
338                 A_M(IJK,north,0) = A_M(IJK,north,0) - EP_G(IJKN)*fac*&
339                    DROODP_G(RO_G(IJKN),P_G(IJKN))*XSI_N(IJK)*V_G(IJK)*AXZ(IJK)
340                 A_M(IJK,south,0) = A_M(IJK,south,0) + EP_G(IJKS)*fac*&
341                    DROODP_G(RO_G(IJKS),P_G(IJKS))*(ONE - XSI_N(IJMK))*V_G(IJMK)*AXZ(IJMK)
342                 IF (DO_K) THEN
343                    A_M(IJK,0,0) = A_M(IJK,0,0) - fac*DROODP_G(RO_G(IJK),P_G(IJK))*&
344                       EP_G(IJK)*((ONE - XSI_T(IJK))*W_G(IJK)*AXY(IJK)-&
345                                         XSI_T(IJKM)*W_G(IJKM)*AXY(IJKM))
346                    A_M(IJK,top,0) = A_M(IJK,top,0) - EP_G(IJKT)*fac*&
347                       DROODP_G(RO_G(IJKT),P_G(IJKT))*XSI_T(IJK)*W_G(IJK)*AXY(IJK)
348                    A_M(IJK,bottom,0) = A_M(IJK,bottom,0) + EP_G(IJKB)*fac*&
349                       DROODP_G(RO_G(IJKB),P_G(IJKB))*(ONE - XSI_T(IJKM))*W_G(IJKM)*AXY(IJKM)
350                 ENDIF
351     
352              ENDIF   !end if (fluid_at(ijk))
353           ENDDO    ! end do (ijk=ijkstart3,ijkend3)
354           ENDIF   ! end if (hs_correct)
355     
356           RETURN
357           END SUBROUTINE COMPRESSIBLE_PP_G
358     
359     
360