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