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