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