File: /nfs/home/0/users/jenkins/mfix.git/model/set_fluidbed_p.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SET_FLUIDBED_P                                          C
4     !  Purpose: Set the pressure field inside the bed assuming a fluidized C
5     !           bed with gravity acting the -ve y-direction                C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 21-JAN-92  C
8     !  Reviewer:M. Syamlal, S. Venkatesan, P. Nicoletti,  Date: 29-JAN-92  C
9     !           W. Rogers                                                  C
10     !                                                                      C
11     !  Revision Number: 1                                                  C
12     !  Purpose: Modifications for including cylindrical geometry           C
13     !  Author: M. Syamlal                                 Date:  6-MAR-92  C
14     !  Reviewer: M. Syamlal                               Date: 11-DEC-92  C
15     !  Revision Number: 2                                                  C
16     !  Purpose: Set pressure drop for cyclic boundary condition w/         C
17     !           pressure drop                                              C
18     !  Author: M. Syamlal                                 Date: 29-APR-94  C
19     !                                                                      C
20     !  Literature/Document References:                                     C
21     !                                                                      C
22     !  Variables referenced: BC_DEFINED, BC_TYPE, IC_P_g, BC_P_g           C
23     !                        EP_g, MW_MIX_G, RO_g0, T_g,                   C
24     !                        SMAX, ROP_s,                                  C
25     !                        DX, DY, DZ, BFY_G, DELP_X, DELP_Y, DELP_Z,    C
26     !                        DO_I, DO_J, DO_K, IMIN1, KMIN1, JMIN1, IMAX1, C
27     !                        IMAX2, JMAX1, JMAX2, KMAX1, KMAX2             C
28     !  Variables modified: P_g                                             C
29     !  Local variables:                                                    C
30     !                                                                      C
31     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
32     
33           SUBROUTINE SET_FLUIDBED_P
34     
35     !-----------------------------------------------
36     ! Modules
37     !-----------------------------------------------
38           USE bc
39           USE bodyforce
40           USE compar
41           USE constant
42           USE discretelement
43           USE eos, ONLY: EOSG
44           USE fldvar
45           USE functions
46           USE funits
47           USE geometry
48           USE ic
49           USE indices
50           USE mpi_utility
51           USE param
52           USE param1
53           USE physprop
54           USE scales
55           USE sendrecv
56           IMPLICIT NONE
57     !-----------------------------------------------
58     ! Local variables
59     !-----------------------------------------------
60     ! indices
61           INTEGER :: I, J, K, IJK, M
62     ! Local loop counter
63           INTEGER :: L
64     ! Gas pressure at the axial location j
65           DOUBLE PRECISION :: PJ
66     ! Bed weight per unit area
67           DOUBLE PRECISION :: BED_WEIGHT
68     ! Total area of a x-z plane
69           DOUBLE PRECISION :: AREA
70     ! x-z plane area of one cell
71           DOUBLE PRECISION :: dAREA
72     ! Average pressure drop per unit length
73           DOUBLE PRECISION :: DPoDX, DPoDY, DPoDZ
74     !-----------------------------------------------
75     
76     ! If any initial pressures are unspecified skip next section
77     ! calculations.
78           DO L = 1, DIMENSION_IC
79              IF (IC_DEFINED(L)) THEN
80                 IF (IC_P_G(L) == UNDEFINED) GOTO 60
81                 PJ = IC_P_G(L)
82              ENDIF
83           ENDDO
84     
85     ! Here the pressure in each cell is determined from a specified pressure
86     ! drop across the domain length. This section requires that the pressure
87     ! is already defined in all initial condition regions (otherwise this
88     ! section would be skipped)
89     ! ---------------------------------------------------------------->>>
90           IF (DO_I .AND. DELP_X/=UNDEFINED) THEN
91              DPODX = DELP_X/XLENGTH
92              PJ = PJ - DPODX*HALF*(DX(IMAX1)+DX(IMAX2))
93              DO I = IMAX1, IMIN1, -1
94                 PJ = PJ + DPODX*HALF*(DX(I)+DX(I+1))
95                 DO K = KMIN1, KMAX1
96                    DO J = JMIN1, JMAX1
97     ! Bound Checking
98                       IF(.NOT.IS_ON_MYPE_OWNS(I,J,K)) CYCLE
99                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
100                       IJK = FUNIJK(I,J,K)
101                       IF (FLUID_AT(IJK)) P_G(IJK) = SCALE(PJ)
102                    ENDDO
103                 ENDDO
104              ENDDO
105           ENDIF
106     
107           IF (DO_J .AND. DELP_Y/=UNDEFINED) THEN
108              DPODY = DELP_Y/YLENGTH
109              PJ = PJ - DPODY*HALF*(DY(JMAX1)+DY(JMAX2))
110              DO J = JMAX1, JMIN1, -1
111                 PJ = PJ + DPODY*HALF*(DY(J)+DY(J+1))
112                 DO K = KMIN1, KMAX1
113                    DO I = IMIN1, IMAX1
114     ! Bound Checking
115                       IF(.NOT.IS_ON_MYPE_OWNS(I,J,K)) CYCLE
116                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
117                       IJK = FUNIJK(I,J,K)
118                       IF (FLUID_AT(IJK)) P_G(IJK) = SCALE(PJ)
119                    ENDDO
120                 ENDDO
121              ENDDO
122           ENDIF
123     
124           IF (DO_K .AND. DELP_Z/=UNDEFINED) THEN
125              DPODZ = DELP_Z/ZLENGTH
126              PJ = PJ - DPODZ*HALF*(DZ(KMAX1)+DZ(KMAX2))
127              DO K = KMAX1, KMIN1, -1
128                 PJ = PJ + DPODZ*HALF*(DZ(K)+DZ(K+1))
129                 DO J = JMIN1, JMAX1
130                    DO I = IMIN1, IMAX1
131     ! Bound Checking
132                       IF(.NOT.IS_ON_MYPE_OWNS(I,J,K)) CYCLE
133                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
134                       IJK = FUNIJK(I,J,K)
135                       IF (FLUID_AT(IJK)) P_G(IJK) = SCALE(PJ)
136                    ENDDO
137                 ENDDO
138              ENDDO
139           ENDIF
140     ! ----------------------------------------------------------------<<<
141           GOTO 100   ! pressure in all intial condition region cells was defined
142     
143        60 CONTINUE   ! pressure in an initial condition region cell was undefined
144     
145     
146     ! ---------------------------------------------------------------->>>
147     ! Search for an outflow boundary condition where pressure is specified
148           PJ = UNDEFINED
149           DO L = 1, DIMENSION_BC
150              IF (BC_DEFINED(L) .AND. BC_TYPE(L)=='P_OUTFLOW') PJ = BC_P_G(L)
151           ENDDO
152     
153           IF (PJ == UNDEFINED) THEN
154     ! either a PO was not specified and/or a PO was specified but not the
155     ! pressure at the outlet
156              IF (RO_G0 /= UNDEFINED) THEN
157     ! If incompressible flow set P_g to zero
158                 DO IJK = IJKSTART3, IJKEND3
159                    IF (FLUID_AT(IJK)) P_G(IJK) = ZERO
160                 ENDDO
161                 GOTO 100
162     
163              ELSE   ! compressible case
164     
165     ! Error condition -- no pressure outflow boundary condition is specified
166     ! if a case is compressible and pressure in any of the initial
167     ! conditions regions is unspecified, then a PO is effectively required
168     ! (i.e., is specifies a bc_p_g).
169                 CALL START_LOG
170                 IF(DMP_LOG)WRITE (UNIT_LOG, 1000)
171                 CALL MFIX_EXIT(myPE)
172              ENDIF
173           ENDIF
174     
175     
176     ! Set an approximate pressure field assuming that the pressure drop
177     ! balances the weight of the bed, if the initial pressure-field is not
178     ! specified
179           DO J = JMAX2, JMIN1, -1
180     
181     ! Find the average weight per unit area over an x-z slice
182              BED_WEIGHT = 0.0
183              AREA = 0.0
184              DO K = KMIN1, KMAX1
185                 DO I = IMIN1, IMAX1
186                    IF(.NOT.IS_ON_MYPE_OWNS(I,J,K)) CYCLE
187                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
188                    IJK = FUNIJK(I,J,K)
189                    IF (FLUID_AT(IJK)) THEN
190                       IF (COORDINATES == 'CARTESIAN') THEN
191                          DAREA = DX(I)*DZ(K)
192                       ELSE IF (CYLINDRICAL) THEN
193                          DAREA = DX(I)*X(I)*DZ(K)
194                       ENDIF
195                       AREA = AREA + DAREA
196                       IF (RO_G0 == UNDEFINED) THEN
197                          BED_WEIGHT = BED_WEIGHT - DY(J)*BFY_G(IJK)*EP_G(IJK)*EOSG(&
198                             MW_MIX_G(IJK),PJ,T_G(IJK))*DAREA
199                       ELSE
200                          BED_WEIGHT = BED_WEIGHT - DY(J)*BFY_G(IJK)*EP_G(IJK)*RO_G0&
201                             *DAREA
202                       ENDIF
203     ! This code is turned off for DEM runs until the value of rop_s can be
204     ! ensured valid values for a DEM run at this point in the code.
205                       IF (.NOT.DISCRETE_ELEMENT) THEN
206                          DO M = 1, SMAX
207                             BED_WEIGHT = BED_WEIGHT - DY(J)*BFY_S(IJK,M)*ROP_S(IJK,M)*&
208                                DAREA
209                          ENDDO
210                       ENDIF  ! end if (.not.discrete_element)
211                    ENDIF  ! end if (fluid_at(ijk))
212                 ENDDO    ! end do loop (i=imin1,imax1)
213              ENDDO    ! end do loop (k=kmin1,kmax1)
214     
215     ! Global Sum
216              call global_all_sum(bed_weight)
217              call global_all_sum(area)
218              IF (AREA /= 0.0) BED_WEIGHT = BED_WEIGHT/AREA
219     
220              PJ = PJ + BED_WEIGHT
221              DO K = KMIN1, KMAX1
222                 DO I = IMIN1, IMAX1
223                    IF(.NOT.IS_ON_MYPE_OWNS(I,J,K)) CYCLE
224                    IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
225                    IJK = FUNIJK(I,J,K)
226                    IF(FLUID_AT(IJK).AND.P_G(IJK)==UNDEFINED)P_G(IJK)=SCALE(PJ)
227                 ENDDO    ! end do (i=imin1,imax1)
228              ENDDO   ! end do (k = kmin1,kmax1)
229           ENDDO   ! end do (j=jmax2,jimn1, -1)
230     ! end setting an undefined pressure in an initial condition region
231     ! ----------------------------------------------------------------<<<
232     
233       100 CONTINUE
234     
235           call send_recv(P_G,2)
236     
237           RETURN
238     
239      1000 FORMAT(/1X,70('*')//' From: SET_FLUIDBED_P'/' Message: Outflow ',&
240              'pressure boundary condition (P_OUTFLOW) not found.',/&
241              'All the initial pressures (IC_P_g) or at least one P_OUTFLOW',/&
242              'condition need to be specified',/1X,70('*')/)
243     
244           END SUBROUTINE SET_FLUIDBED_P
245     
246     
247