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

1     ! TO DO:
2     ! 1. check the formulation based on MCp.
3     
4     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
5     !                                                                      C
6     !  Subroutine: CORRECT_1                                               C
7     !  Purpose: Correct the solids volume fraction                         C
8     !                                                                      C
9     !  Notes: MCP must be defined to call this routine.                    C
10     !         The solids volume fraction correction routines are           C
11     !         explicitly setup for a single solids phases.  Therefore,     C
12     !         even if multiple solids phases are present, only one of      C
13     !         the solids phases (M=MCP) can employ the correction          C
14     !         routines                                                     C
15     !                                                                      C
16     !  Author: M. Syamlal                                 Date: 25-SEP-96  C
17     !  Reviewer:                                          Date:            C
18     !                                                                      C
19     !  Revision Number: 1                                                  C
20     !  Purpose: To incorporate Cartesian grid modifications                C
21     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
22     !                                                                      C
23     !  Literature/Document References:                                     C
24     !  Variables referenced:                                               C
25     !  Variables modified:                                                 C
26     !  Local variables:                                                    C
27     !                                                                      C
28     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
29     
30           SUBROUTINE CORRECT_1()
31     
32     !-----------------------------------------------
33     ! Modules
34     !-----------------------------------------------
35           USE param1, only: one, zero, undefined_i
36           USE fldvar, only: u_s, v_s, w_s, rop_s, ro_s, ep_s
37           USE physprop, only: close_packed
38           USE pscor, only: mcp, epp, k_cp, e_e, e_n, e_t
39           USE ur_facs, only: ur_fac
40           USE indices
41           USE geometry
42           USE compar
43           USE sendrecv
44           USE cutcell, only: cartesian_grid, cut_cell_at, cg_ur_fac
45           USE visc_s, only: ep_star_array
46           USE functions
47           IMPLICIT NONE
48     !-----------------------------------------------
49     ! Local variables
50     !-----------------------------------------------
51     ! corrected solids volume fraction
52           DOUBLE PRECISION :: EPcor
53     ! dPodEP_s(EP_s(IJK, M)) * EPp(IJK)
54           DOUBLE PRECISION :: Pp_P
55     ! indices
56           INTEGER :: IJK, IJKE, IJKN, IJKT
57     ! solids index
58           INTEGER :: M
59     ! solids volume fraction at maximum packing
60           DOUBLE PRECISION :: EP_S_CP
61     !-----------------------------------------------
62     
63           IF (MCP == UNDEFINED_I) THEN
64     ! this error should be caught earlier in the routines so that this
65     ! branch should never be entered
66              RETURN
67           ELSE
68     ! the lowest solids phase index of those solids phases that can close
69     ! pack (i.e. close_packed=T) and the index of the solids phase that is
70     ! used to form the solids correction equation.
71              M = MCP
72           ENDIF
73     
74     
75     ! by definition M must be close_packed (this is a redundant check)
76           IF (CLOSE_PACKED(M)) THEN
77     
78     ! Correct solids volume fraction
79     ! ---------------------------------------------------------------->>>
80     !!$omp    parallel do &
81     !!$omp&   private( IJK, EPCOR )
82              DO IJK = ijkstart3, ijkend3
83                 IF (FLUID_AT(IJK)) THEN
84     
85                    EPCOR = EP_S(IJK,M) + EPP(IJK)
86                    EP_S_CP = ONE - EP_STAR_ARRAY(IJK)
87                    IF (EPCOR>EP_S_CP .AND. EPP(IJK)>ZERO) THEN
88                       EPP(IJK) = UR_FAC(2)*EPP(IJK)
89     
90                       IF(CARTESIAN_GRID) THEN
91     ! JFD: Using a low value of CG_UR_FAC(2) for the cut cell tends to
92     ! stabilize code (limits occurrence of negative void fraction)
93                          IF(CUT_CELL_AT(IJK)) EPP(IJK) = CG_UR_FAC(2)*EPP(IJK)
94                       ENDIF
95     
96                       EPCOR = EP_S(IJK,M) + EPP(IJK)
97                    ENDIF
98                    ROP_S(IJK,M) = MAX(ZERO,RO_S(IJK,M)*EPCOR)
99                 ENDIF
100              ENDDO
101     ! end correct solids volume fraction
102     ! ----------------------------------------------------------------<<<
103     
104     ! Correct solids velocities
105     ! ---------------------------------------------------------------->>>
106     !!$omp    parallel do &
107     !!$omp&   private( IJK, PP_P, IJKE, IJKN, IJKT )
108              DO IJK = ijkstart3, ijkend3
109                 IF (FLUID_AT(IJK)) THEN
110     
111                    PP_P = K_CP(IJK)*EPP(IJK)
112                    IF (FLOW_AT_E(IJK)) THEN
113                       IJKE = EAST_OF(IJK)
114                       U_S(IJK,M)=U_S(IJK,M)-E_E(IJK)*&
115                          (K_CP(IJKE)*EPP(IJKE)-PP_P)
116                    ENDIF
117                    IF (FLOW_AT_N(IJK)) THEN
118                       IJKN = NORTH_OF(IJK)
119                       V_S(IJK,M)=V_S(IJK,M)-E_N(IJK)*&
120                          (K_CP(IJKN)*EPP(IJKN)-PP_P)
121                    ENDIF
122                    IF (DO_K) THEN
123                       IF (FLOW_AT_T(IJK)) THEN
124                          IJKT = TOP_OF(IJK)
125                          W_S(IJK,M) = W_S(IJK,M) - E_T(IJK)*&
126                             (K_CP(IJKT)*EPP(IJKT)-PP_P)
127                       ENDIF
128                    ENDIF
129                 ENDIF
130              ENDDO
131     ! end correct solids velocities
132     ! ----------------------------------------------------------------<<<
133     
134           ENDIF   ! end if (close_packed)
135     
136           RETURN
137           END SUBROUTINE CORRECT_1
138     
139     
140