File: RELATIVE:/../../../mfix.git/model/set_bc1.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SET_BC1                                                 C
4     !  Purpose: Set transient flow boundary conditions                     C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 29-JAN-92  C
7     !                                                                      C
8     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
9           SUBROUTINE SET_BC1
10     
11     ! Modules
12     !---------------------------------------------------------------------//
13           use bc, only: bc_defined, bc_type
14           USE param, only: dimension_bc
15           IMPLICIT NONE
16     
17     ! Local variables
18     !---------------------------------------------------------------------//
19     ! index for boundary condition
20           INTEGER :: L
21     
22     
23     ! Set the boundary conditions
24           DO L = 1, DIMENSION_BC
25              IF (BC_DEFINED(L)) THEN
26     
27                 SELECT CASE(TRIM(BC_TYPE(L)))
28                 CASE ('P_OUTFLOW')
29                    CALL SET_OUTFLOW(L)
30                    CALL SET_BC1_REPORT_OUTFLOW(L)
31                 CASE ('MASS_OUTFLOW')
32                    CALL SET_OUTFLOW(L)
33                    CALL SET_BC1_ADJUST_OUTFLOW(L)
34                 CASE ('MASS_INFLOW')
35                    CALL SET_BC1_JET(L)
36                 CASE ('P_INFLOW')
37                    CALL SET_OUTFLOW(L)
38                 CASE ('OUTFLOW')
39                    CALL SET_OUTFLOW(L)
40                    CALL SET_BC1_REPORT_OUTFLOW(L)
41                 END SELECT
42              ENDIF   ! end if (bc_defined(l))
43           ENDDO    ! end do loop (l=1,dimension_bc)
44     
45           RETURN
46           END SUBROUTINE SET_BC1
47     
48     
49     
50     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
51     !                                                                      C
52     !  Subroutine: set_bc1_jet                                             C
53     !  Purpose: update transient jet conditions                            C
54     !                                                                      C
55     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
56           SUBROUTINE SET_BC1_JET(BCV)
57     
58     ! Modules
59     !---------------------------------------------------------------------//
60           use bc, only: bc_plane
61           use bc, only: bc_k_b, bc_k_t
62           use bc, only: bc_j_s, bc_j_n
63           use bc, only: bc_i_w, bc_i_e
64           use bc, only: bc_jet_g
65           use bc, only: bc_time
66           use bc, only: bc_jet_gh, bc_dt_h
67           use bc, only: bc_jet_gl, bc_dt_l
68     
69           use fldvar, only: u_g, v_g, w_g
70     
71           use run, only: time, dt
72     
73           use param1, only: undefined
74     
75           use functions, only: im_of, jm_of, km_of
76           use functions, only: is_on_mype_plus2layers
77           use functions, only: funijk
78           use compar, only: dead_cell_at
79     
80           IMPLICIT NONE
81     ! Dummy arguments
82     !---------------------------------------------------------------------//
83     ! index for boundary condition
84           INTEGER, INTENT(IN) :: BCV
85     
86     ! Local variables
87     !---------------------------------------------------------------------//
88     ! indices
89           INTEGER :: I, J, K, IJK
90     ! IJK index for setting velocity bc
91           INTEGER :: IJK2
92     !---------------------------------------------------------------------//
93     
94           IF (TIME + 0.1d0*DT>=BC_TIME(BCV) .AND. &
95               BC_JET_G(BCV)/=UNDEFINED) THEN
96     
97              IF (BC_JET_G(BCV) == BC_JET_GH(BCV)) THEN
98                 BC_JET_G(BCV) = BC_JET_GL(BCV)
99                 BC_TIME(BCV) = TIME + BC_DT_L(BCV)
100              ELSEIF (BC_JET_G(BCV) == BC_JET_GL(BCV)) THEN
101                 BC_JET_G(BCV) = BC_JET_GH(BCV)
102                 BC_TIME(BCV) = TIME + BC_DT_H(BCV)
103              ELSE
104                 BC_JET_G(BCV) = BC_JET_GH(BCV)
105                 BC_TIME(BCV) = TIME + BC_DT_H(BCV)
106              ENDIF
107     
108              DO K = BC_K_B(BCV), BC_K_T(BCV)
109              DO J = BC_J_S(BCV), BC_J_N(BCV)
110              DO I = BC_I_W(BCV), BC_I_E(BCV)
111                 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
112                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
113                 IJK = FUNIJK(I,J,K)
114     ! Why is the velocity of the boundary cell not always set (in case of
115     ! w, s, or b plane then velocity of the adjacent fluid cell is set)?
116     ! It should not really matter for MI...
117                 SELECT CASE (TRIM(BC_PLANE(BCV)))
118                 CASE ('W')
119                    IJK2 = IM_OF(IJK)
120                    U_G(IJK2) = BC_JET_G(BCV)
121                 CASE ('E')
122                    U_G(IJK) = BC_JET_G(BCV)
123                 CASE ('S')
124                    IJK2 = JM_OF(IJK)
125                    V_G(IJK2) = BC_JET_G(BCV)
126                 CASE ('N')
127                    V_G(IJK) = BC_JET_G(BCV)
128                 CASE ('B')
129                    IJK2 = KM_OF(IJK)
130                    W_G(IJK2) = BC_JET_G(BCV)
131                 CASE ('T')
132                    W_G(IJK) = BC_JET_G(BCV)
133                 END SELECT
134              ENDDO
135              ENDDO
136              ENDDO
137           ENDIF
138           RETURN
139           END SUBROUTINE SET_BC1_JET
140     
141     
142     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
143     !                                                                      C
144     !  Subroutine: set_bc1_report_outflow                                  C
145     !  Purpose: print out outflow conditions                               C
146     !                                                                      C
147     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
148           SUBROUTINE SET_BC1_REPORT_OUTFLOW(BCV)
149     
150     ! Modules
151     !---------------------------------------------------------------------//
152           use bc, only: bc_dt_0, bc_time
153           use bc, only: bc_out_n
154           use bc, only: bc_mout_g, bc_mout_s
155           use bc, only: bc_vout_g, bc_vout_s
156           use funits, only: dmp_log, unit_log
157           use param1, only: undefined, zero
158           use physprop, only: smax
159           use machine, only: start_log, end_log
160           use run, only: time, dt, tstop
161           IMPLICIT NONE
162     
163     ! Dummy arguments
164     !---------------------------------------------------------------------//
165     ! index for boundary condition
166           INTEGER, INTENT(IN) :: BCV
167     
168     ! Local variables
169     !---------------------------------------------------------------------//
170     ! Solids phase index
171           INTEGER :: M
172     
173           IF (BC_DT_0(BCV) == UNDEFINED) RETURN
174     
175           CALL CALC_OUTFLOW(BCV)
176     
177     ! Calculate and accumulate the actual mass and volume outflow
178           IF (TIME + 0.1d0*DT>=BC_TIME(BCV) .OR. &
179               TIME+0.1d0*DT>=TSTOP) THEN
180              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
181     
182     ! Average and print out the flow rates
183              BC_MOUT_G(BCV) = ABS(BC_MOUT_G(BCV))/BC_OUT_N(BCV)
184              BC_VOUT_G(BCV) = ABS(BC_VOUT_G(BCV))/BC_OUT_N(BCV)
185              CALL START_LOG
186              IF(DMP_LOG)WRITE (UNIT_LOG, 1000) BCV, TIME
187              IF(DMP_LOG)WRITE (UNIT_LOG, 1100) BC_MOUT_G(BCV), &
188                 BC_VOUT_G(BCV)
189              BC_MOUT_G(BCV) = ZERO
190              BC_VOUT_G(BCV) = ZERO
191              DO M = 1, SMAX
192                 BC_MOUT_S(BCV,M) = ABS(BC_MOUT_S(BCV,M))/BC_OUT_N(BCV)
193                 BC_VOUT_S(BCV,M) = ABS(BC_VOUT_S(BCV,M))/BC_OUT_N(BCV)
194                 IF(DMP_LOG)WRITE (UNIT_LOG, 1200) M, &
195                    BC_MOUT_S(BCV,M), BC_VOUT_S(BCV,M)
196                 BC_MOUT_S(BCV,M) = ZERO
197                 BC_VOUT_S(BCV,M) = ZERO
198              ENDDO
199              CALL END_LOG
200              BC_OUT_N(BCV) = 0
201           ENDIF
202     
203      1000 FORMAT(/,1X,'Average outflow rates at BC No. ',I2,'  At Time = ',G12.5)
204      1100 FORMAT(3X,'Gas : Mass flow = ',G12.5,'     Volumetric flow = ',G12.5)
205      1200 FORMAT(3X,'Solids-',I1,' : Mass flow = ',G12.5,&
206              '     Volumetric flow = ',G12.5)
207     
208           RETURN
209           END SUBROUTINE SET_BC1_REPORT_OUTFLOW
210     
211     
212     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
213     !                                                                      C
214     !  Subroutine: set_bc1_adjust_outflow                                  C
215     !  Purpose: Adjust velocities to get specified mass or volumetric      C
216     !  flow rate based on average outflow rate                             C
217     !                                                                      C
218     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
219           SUBROUTINE SET_BC1_ADJUST_OUTFLOW(BCV)
220     
221     ! Modules
222     !---------------------------------------------------------------------//
223           use bc, only: bc_dt_0, bc_time
224           use bc, only: bc_i_w, bc_i_e
225           use bc, only: bc_j_s, bc_j_n
226           use bc, only: bc_k_b, bc_k_t
227           use bc, only: bc_massflow_g, bc_massflow_s
228           use bc, only: bc_mout_g, bc_mout_s
229           use bc, only: bc_out_n
230           use bc, only: bc_plane
231           use bc, only: bc_u_g, bc_v_g, bc_w_g
232           use bc, only: bc_u_s, bc_v_s, bc_w_s
233           use bc, only: bc_volflow_g, bc_volflow_s
234           use bc, only: bc_vout_g, bc_vout_s
235           use compar, only: dead_cell_at
236           use fldvar, only: u_g, v_g, w_g
237           use fldvar, only: u_s, v_s, w_s, rop_s
238           use functions, only: funijk
239           use functions, only: im_of, jm_of, km_of
240           use functions, only: is_on_mype_plus2layers
241           use funits, only: dmp_log, unit_log
242           use machine, only: start_log, end_log
243           use param1, only: undefined, zero, small_number
244           use physprop, only: smax, mmax
245           use run, only: kt_type_enum, ghd_2007
246           use run, only: time, dt, tstop
247     
248           IMPLICIT NONE
249     ! Dummy arguments
250     !---------------------------------------------------------------------//
251     ! index for boundary condition
252           INTEGER, INTENT(IN) :: BCV
253     
254     ! Local variables
255     !---------------------------------------------------------------------//
256     ! indices
257           INTEGER :: I, J, K, IJK
258     ! IJK index for setting velocity bc
259           INTEGER :: IJK2
260     ! Solids phase index
261           INTEGER :: M
262     ! local solids velocity for mixture (for ghd)
263           DOUBLE PRECISION :: lvel_s
264     !---------------------------------------------------------------------//
265     
266           CALL CALC_OUTFLOW(BCV)
267     
268     ! Calculate and accumulate the actual mass and volume outflow
269           IF (TIME + 0.1d0*DT>=BC_TIME(BCV) .OR. &
270               TIME+0.1d0*DT>=TSTOP) THEN
271              BC_TIME(BCV) = TIME + BC_DT_0(BCV)
272     
273     ! Average and print out the flow rates
274              BC_MOUT_G(BCV) = ABS(BC_MOUT_G(BCV))/BC_OUT_N(BCV)
275              BC_VOUT_G(BCV) = ABS(BC_VOUT_G(BCV))/BC_OUT_N(BCV)
276              CALL START_LOG
277              IF(DMP_LOG)WRITE (UNIT_LOG, 1000) BCV, TIME
278              IF(DMP_LOG)WRITE (UNIT_LOG, 1100) BC_MOUT_G(BCV), &
279                 BC_VOUT_G(BCV)
280              DO M = 1, SMAX
281                 BC_MOUT_S(BCV,M) = ABS(BC_MOUT_S(BCV,M))/BC_OUT_N(BCV)
282                 BC_VOUT_S(BCV,M) = ABS(BC_VOUT_S(BCV,M))/BC_OUT_N(BCV)
283                 IF(DMP_LOG)WRITE (UNIT_LOG, 1200) M, &
284                    BC_MOUT_S(BCV,M), BC_VOUT_S(BCV,M)
285              ENDDO
286              CALL END_LOG
287              BC_OUT_N(BCV) = 0
288     
289     ! Now that we know the mass and volume outflow update the bc velocities
290     ! (gas phase)
291              IF (BC_MASSFLOW_G(BCV) /= UNDEFINED) THEN
292                 IF (BC_MOUT_G(BCV) > SMALL_NUMBER) THEN
293                    SELECT CASE (TRIM(BC_PLANE(BCV)))
294                    CASE ('W', 'E')
295                       BC_U_G(BCV) = BC_U_G(BCV)*BC_MASSFLOW_G(BCV)/&
296                          BC_MOUT_G(BCV)
297                    CASE ('S', 'N')
298                       BC_V_G(BCV) = BC_V_G(BCV)*BC_MASSFLOW_G(BCV)/&
299                          BC_MOUT_G(BCV)
300                    CASE ('B', 'T')
301                       BC_W_G(BCV) = BC_W_G(BCV)*BC_MASSFLOW_G(BCV)/&
302                          BC_MOUT_G(BCV)
303                    END SELECT
304                 ENDIF
305              ELSEIF (BC_VOLFLOW_G(BCV) /= UNDEFINED) THEN
306                 IF (BC_VOUT_G(BCV) > SMALL_NUMBER) THEN
307                    SELECT CASE (TRIM(BC_PLANE(BCV)))
308                    CASE ('W', 'E')
309                       BC_U_G(BCV) = BC_U_G(BCV)*BC_VOLFLOW_G(BCV)/&
310                          BC_VOUT_G(BCV)
311                    CASE ('S', 'N')
312                       BC_V_G(BCV) = BC_V_G(BCV)*BC_VOLFLOW_G(BCV)/&
313                          BC_VOUT_G(BCV)
314                    CASE ('B', 'T')
315                       BC_W_G(BCV) = BC_W_G(BCV)*BC_VOLFLOW_G(BCV)/&
316                          BC_VOUT_G(BCV)
317                    END SELECT
318                 ENDIF
319              ENDIF
320     ! zero out counter for new cycle
321              BC_MOUT_G(BCV) = zero
322              BC_VOUT_G(BCV) = zero
323     
324     ! Now that we know the mass and volume outflow update the bc velocities
325     ! (solids phase)
326              DO M = 1, SMAX
327                 IF (BC_MASSFLOW_S(BCV,M) /= UNDEFINED) THEN
328                    IF (BC_MOUT_S(BCV,M) > SMALL_NUMBER) THEN
329                       SELECT CASE (TRIM(BC_PLANE(BCV)))
330                       CASE ('W', 'E')
331                          BC_U_S(BCV,M) = BC_U_S(BCV,M)*&
332                             BC_MASSFLOW_S(BCV,M)/BC_MOUT_S(BCV,M)
333                       CASE ('S', 'N')
334                          BC_V_S(BCV,M) = BC_V_S(BCV,M)*&
335                             BC_MASSFLOW_S(BCV,M)/BC_MOUT_S(BCV,M)
336                       CASE ('B', 'T')
337                          BC_W_S(BCV,M) = BC_W_S(BCV,M)*&
338                             BC_MASSFLOW_S(BCV,M)/BC_MOUT_S(BCV,M)
339                       END SELECT
340                    ENDIF
341                 ELSEIF (BC_VOLFLOW_S(BCV,M) /= UNDEFINED) THEN
342                    IF (BC_VOUT_S(BCV,M) > SMALL_NUMBER) THEN
343                       SELECT CASE (TRIM(BC_PLANE(BCV)))
344                       CASE ('W', 'E')
345                          BC_U_S(BCV,M) = BC_U_S(BCV,M)*&
346                             BC_VOLFLOW_S(BCV,M)/BC_VOUT_S(BCV,M)
347                       CASE ('S', 'N')
348                          BC_V_S(BCV,M) = BC_V_S(BCV,M)*&
349                             BC_VOLFLOW_S(BCV,M)/BC_VOUT_S(BCV,M)
350                       CASE ('B', 'T')
351                          BC_W_S(BCV,M) = BC_W_S(BCV,M)*&
352                             BC_VOLFLOW_S(BCV,M)/BC_VOUT_S(BCV,M)
353                      END SELECT
354                    ENDIF
355                 ENDIF
356     ! zero out counter for new cycle
357                 BC_MOUT_S(BCV,M) = zero
358                 BC_VOUT_S(BCV,M) = zero
359              ENDDO
360     
361     ! Apply updated boundary velocities - Define the field variables at the
362     ! boundaries according to user specifications with modifications from
363     ! the above calculations.
364     ! If the boundary plane is W, S, or B (i.e., the fluid cell is on the
365     ! west, south or bottom of the boundary cell) then define the velocity
366     ! of the adjacent fluid cell according to the boundary velocity rather
367     ! than the velocity of the boundary cell.
368     ! Why not set the velocity in the boundary cell itself?  Based on the
369     ! momentum bc routine it should not really matter for MO.
370              DO K = BC_K_B(BCV), BC_K_T(BCV)
371              DO J = BC_J_S(BCV), BC_J_N(BCV)
372              DO I = BC_I_W(BCV), BC_I_E(BCV)
373                 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
374                 IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
375                 IJK = FUNIJK(I,J,K)
376                 SELECT CASE (TRIM(BC_PLANE(BCV)))
377                 CASE ('W')
378                    IJK2 = IM_OF(IJK)
379                    U_G(IJK2) = BC_U_G(BCV)
380                 CASE ('E')
381                    U_G(IJK) = BC_U_G(BCV)
382                 CASE ('S')
383                    IJK2 = JM_OF(IJK)
384                    V_G(IJK2) = BC_V_G(BCV)
385                 CASE ('N')
386                    V_G(IJK) = BC_V_G(BCV)
387                 CASE ('B')
388                    IJK2 = KM_OF(IJK)
389                    W_G(IJK2) = BC_W_G(BCV)
390                 CASE ('T')
391                    W_G(IJK) = BC_W_G(BCV)
392                 END SELECT
393                 DO M = 1, SMAX
394                    SELECT CASE (TRIM(BC_PLANE(BCV)))
395                    CASE ('W')
396                       IJK2 = IM_OF(IJK)
397                       U_S(IJK2,M) = BC_U_S(BCV,M)
398                    CASE ('E')
399                       U_S(IJK,M) = BC_U_S(BCV,M)
400                    CASE ('S')
401                       IJK2 = JM_OF(IJK)
402                       V_S(IJK2,M) = BC_V_S(BCV,M)
403                    CASE ('N')
404                       V_S(IJK,M) = BC_V_S(BCV,M)
405                    CASE ('B')
406                       IJK2 = KM_OF(IJK)
407                       W_S(IJK2,M) = BC_W_S(BCV,M)
408                    CASE ('T')
409                       W_S(IJK,M) = BC_W_S(BCV,M)
410                    END SELECT
411                 ENDDO
412     
413     ! compute mixutre velocity BC for GHD theory
414                 IF(KT_TYPE_ENUM == GHD_2007) THEN
415                    lvel_s = zero
416     
417     ! bulk density is set by set_outflow and is set in bc according to
418     ! neighboring fluid appropriately. so we don't need to reference
419     ! bulk density with ijk vs ijk2 index (ijk value will = ijk2 value).
420                    DO M = 1, SMAX
421                       SELECT CASE (TRIM(BC_PLANE(BCV)))
422                       CASE ('W', 'E')
423                          lvel_s = lvel_s + &
424                             BC_U_S(BCV,M)*ROP_S(IJK,M)
425                       CASE ('S', 'N')
426                          lvel_s = lvel_s + &
427                             BC_V_S(BCV,M)*ROP_S(IJK,M)
428                       CASE ('B', 'T')
429                          lvel_s = lvel_s + &
430                             BC_W_S(BCV,M)*ROP_S(IJK,M)
431                       END SELECT
432                    ENDDO
433     
434                    IF (ROP_S(IJK,MMAX) > 0) THEN
435                       lvel_s = lvel_s /ROP_S(IJK,MMAX)
436                    ELSE
437                       lvel_s = ZERO
438                    ENDIF
439     
440                    SELECT CASE (TRIM(BC_PLANE(BCV)))
441                    CASE ('W')
442                       IJK2 = IM_OF(IJK)
443                       U_S(IJK2,MMAX) =  lvel_s
444                    CASE ('E')
445                       U_S(IJK,MMAX) = lvel_s
446                    CASE ('S')
447                       IJK2 = JM_OF(IJK)
448                       V_S(IJK2,MMAX) = lvel_s
449                    CASE ('N')
450                       V_S(IJK,MMAX) = lvel_s
451                    CASE ('B')
452                       IJK2 = KM_OF(IJK)
453                       W_S(IJK2,MMAX) = lvel_s
454                    CASE ('T')
455                       W_S(IJK,MMAX) = lvel_s
456                    END SELECT
457     
458                 ENDIF   ! end if (kt_type_enum==ghd_2007)
459              ENDDO
460              ENDDO
461              ENDDO
462           ENDIF   ! if time to update outflow condition
463     
464      1000 FORMAT(/,1X,'Average outflow rates at BC No. ',I2,'  At Time = ',G12.5)
465      1100 FORMAT(3X,'Gas : Mass flow = ',G12.5,'     Volumetric flow = ',G12.5)
466      1200 FORMAT(3X,'Solids-',I1,' : Mass flow = ',G12.5,&
467              '     Volumetric flow = ',G12.5)
468     
469           RETURN
470           END SUBROUTINE SET_BC1_ADJUST_OUTFLOW
471