File: N:\mfix\model\k_epsilon_prop.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: K_Epsilon_PROP                                          C
4     !  Purpose: Calculate diffusion coefficeint and sources for K &        C
5     !           Epsilon equations                                          C
6     !                                                                      C
7     !  Author:                                       Date:                 C
8     !  Modified: S. Benyahia                         Date:May-13-04        C
9     !                                                                      C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !  Wilcox, D.C., Turbulence Modeling for CFD. DCW Industries, Inc.     C
13     !     La Canada, Ca. 1994.                                             C
14     !                                                                      C
15     !                                                                      C
16     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
17           SUBROUTINE K_Epsilon_PROP()
18     
19     ! Modules
20     !---------------------------------------------------------------------//
21           USE param
22           USE param1
23           USE parallel
24           USE physprop
25           USE drag
26           USE run
27           USE output
28           USE geometry
29           USE fldvar
30           USE visc_g
31           USE visc_s
32           USE trace
33           USE indices
34           USE constant
35           Use vshear
36           USE turb
37           USE toleranc
38           USE compar
39           USE TAU_G
40           USE sendrecv
41     
42           USE cutcell
43           USE fun_avg
44           USE functions
45     
46           IMPLICIT NONE
47     
48     ! Dummy arguments
49     !---------------------------------------------------------------------//
50     
51     ! Local parameters
52     !---------------------------------------------------------------------//
53           DOUBLE PRECISION, PARAMETER :: F2O3 = 2.D0/3.D0
54     
55     ! Local variables
56     !---------------------------------------------------------------------//
57     ! Indices
58           INTEGER :: I, J, K, IJK
59           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
60     ! Loop indices
61           INTEGER :: I1, J1, K1
62     ! Solids phase index
63           INTEGER :: M
64     ! U_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
65           DOUBLE PRECISION :: UgE, UgW
66     ! U_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
67           DOUBLE PRECISION :: UgN, UgS
68     ! U_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
69           DOUBLE PRECISION :: UgT, UgB
70     ! U_g at the center of a scalar cell (i, j, k)
71     ! Calculated for Cylindrical coordinates only.
72           DOUBLE PRECISION :: UgcC
73     ! Cell center value of U_g
74           DOUBLE PRECISION UGC
75     ! V_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
76           DOUBLE PRECISION :: VgE, VgW
77     ! V_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
78           DOUBLE PRECISION :: VgN, VgS
79     ! V_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
80           DOUBLE PRECISION :: VgT, VgB
81     ! Cell center value of V_g
82           DOUBLE PRECISION VGC
83     ! W_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
84           DOUBLE PRECISION :: WgE, WgW
85     ! W_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
86           DOUBLE PRECISION :: WgN, WgS
87     ! W_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
88           DOUBLE PRECISION :: WgT, WgB
89     ! W_g at the center of a scalar cell (i, j, k).
90     ! Calculated for Cylindrical coordinates only.
91           DOUBLE PRECISION :: WgcC
92     ! Cell center value of W_g
93           DOUBLE PRECISION WGC
94     ! trace_g and eddy viscosity
95           DOUBLE PRECISION :: Trace_G, Mu_gas_t
96           DOUBLE PRECISION :: C_Eps_Pope, Xsi_Pope
97     ! gradient of velocity
98           DOUBLE PRECISION :: DelV_G(3,3)
99     ! rate of strain tensor
100           DOUBLE PRECISION :: D_g(3,3)
101     
102     !  Production of Turb. Due to shear, Turb Visc, and Constants.
103     !  See Wilcox PP. 89
104           DOUBLE PRECISION Tauij_gDUi_gODxj, C_MU, Sigma_k, Sigma_E, Kappa
105           DOUBLE PRECISION Pos_Tauij_gDUi_gODxj, Neg_Tauij_gDUi_gODxj
106           DOUBLE PRECISION Ceps_1, Ceps_2, C_Eps_3, Check_Log
107           DOUBLE PRECISION Pos_PI_kq_2, Neg_PI_kq_2
108     ! Modif. for Sof Local Var.
109           INTEGER :: P,Q
110     !---------------------------------------------------------------------//
111     
112           IF( .NOT. K_Epsilon) RETURN
113     
114     ! M should be forced to be equal to one to get some info.
115     ! from solids phase.
116           M = 1
117     
118     ! Add constants. Most of these constants have the same names and values
119     ! as the ones defined in Wilcox book (turbulence modeling for CFD).
120     ! Some are necessary only for Simonin turbulence model
121     
122           C_MU = 9.0D-02
123           Kappa = 0.42D+0
124           Sigma_k = 1.0D0
125           Sigma_E = 1.3D0
126           Ceps_1 = 1.44D0 !should be 1.6 for axisymmetric cases with no Pope Correction.
127           Ceps_2 = 1.92D0
128           C_Eps_3 = 1.2D0 ! for Simonin model
129           C_Eps_Pope = 0.79D0 ! for Pope Correction.
130     
131     !!!$omp  parallel do private(ijk, L)
132           DO IJK = IJKSTART3, IJKEND3
133              IF (FLUID_AT(IJK)) THEN
134                 I = I_OF(IJK)
135                 J = J_OF(IJK)
136                 K = K_OF(IJK)
137     
138     ! Velocity components at the scalar faces
139                 CALL GET_FACE_VEL_GAS(IJK, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
140                                            vge, vgw, vgn, vgs, vgt, vgb, &
141                                            wge, wgw, wgn, wgs, wgt, wgb, wgcc)
142     
143     ! Velocity derivatives (gradient and rate of strain tensor)
144                 CALL CALC_DERIV_VEL_GAS(IJK, DelV_g, D_g)
145     
146     ! Value of trace is also available via trd_g but the latter is
147     ! calculated only once per time step
148                 Trace_g = D_G(1,1) + D_G(2,2) + D_G(3,3)
149     
150     ! Pope's correction in 2-D to the Epsilon Equation for a round-Jet
151     ! from Wilcox book, Page 103.
152                 Xsi_Pope = ZERO
153                 DO I1 = 1,3
154                    DO J1 = 1,3
155                       DO K1 = 1,3
156                          Xsi_Pope = Xsi_Pope + (DelV_g(I1,J1) - DelV_g(J1,I1))* &
157                                     (DelV_g(J1,K1) - DelV_g(K1,J1))*&
158                                     (DelV_g(K1,I1) + DelV_g(I1,K1))
159                       ENDDO
160                    ENDDO
161                 ENDDO
162                 Xsi_Pope = Xsi_Pope/6. * K_Turb_G(IJK)**2/(E_Turb_G(IJK)+Small_number)
163     
164     ! This IF statment is to ensure that we are using the turbulent viscosity
165     ! and NOT the effective viscosity.
166                 IF (MU_GT(IJK) .GE. MU_g(IJK)) THEN
167                    Mu_gas_t = MU_GT(IJK) - MU_g(IJK)
168                 ELSE
169                    Mu_gas_t = ZERO
170                 ENDIF
171     
172     ! Calculate Tau(i,j)*dUi/dXj (production term in the K Equation
173                 IF(.NOT.CUT_CELL_AT(IJK)) THEN
174                    Tauij_gDUi_gODxj = 2D0*Mu_gas_t*(             &
175                       D_G(1,1) * D_G(1,1) + &
176                       D_G(1,2) * (UgN-UgS)*ODY(J) + &
177                       D_G(1,3) * ((UgT-UgB)*(OX(I)*ODZ(K))-WgcC*OX(I)) + &
178                       D_G(2,1) * (VgE-VgW)*ODX(I) + &
179                       D_G(2,2) * D_G(2,2) + &
180                       D_G(2,3) * (VgT-VgB)*(OX(I)*ODZ(K)) + &
181                       D_G(3,1) * (WgE-WgW)*ODX(I) + &
182                       D_G(3,2) * (WgN-WgS)*ODY(J) + &
183                       D_G(3,3) * D_G(3,3)) - &
184                       F2O3 * RO_G(IJK) * K_Turb_G(IJK)*Trace_g - &
185                       F2O3 * Mu_gas_t * Trace_g**2
186     
187                 ELSE  ! CUT_CELL
188     ! This is actually not used because of wall functions in cut cells
189     !               Tauij_gDUi_gODxj = 2D0*Mu_gas_t*(                             &
190     !                  D_G(1,1) * UG(1,1)  +                      &
191     !                  D_G(1,2) * UG(1,2)  +                      &
192     !                  D_G(1,3) * UG(1,3)  +                      &
193     !                  D_G(2,1) * UG(2,1)  +                      &
194     !                  D_G(2,2) * UG(2,2)  +                      &
195     !                  D_G(2,3) * UG(2,3)  +                      &
196     !                  D_G(3,1) * UG(3,1)  +                      &
197     !                  D_G(3,2) * UG(3,2)  +                      &
198     !                  D_G(3,3) * UG(3,3)) -                      &
199     !                  F2O3 * RO_G(IJK) * K_Turb_G(IJK)*Trace_g-  &
200     !                  F2O3 * Mu_gas_t * Trace_g**2
201                 ENDIF
202     
203     
204     ! To avoid very small negative numbers
205                 IF (Tauij_gDUi_gODxj .GE. ZERO) THEN
206                    Pos_Tauij_gDUi_gODxj = Tauij_gDUi_gODxj
207                    Neg_Tauij_gDUi_gODxj = ZERO
208                 ELSE
209                    Pos_Tauij_gDUi_gODxj = ZERO
210                    Neg_Tauij_gDUi_gODxj = Tauij_gDUi_gODxj
211                 ENDIF
212     
213     
214     ! Interaction terms in the K-Epsilon equations FOR USE Simonin and Ahmadi models
215                 IF(SIMONIN) THEN
216                    Pos_PI_kq_2 = F_GS(IJK,1)*K_12(IJK)
217                    Neg_PI_kq_2 = F_GS(IJK,1)*2.0D0* K_Turb_G(IJK)
218                 ELSE IF(AHMADI) THEN
219                    Pos_PI_kq_2 = F_GS(IJK,1)*3.0D0*Theta_m(IJK,M)
220                    Neg_PI_kq_2 = F_GS(IJK,1)*2.0D0* K_Turb_G(IJK)
221                    C_Eps_3 = zero ! no extra terms in epsilon equation for Ahmadi model !
222                 ELSE
223                    Pos_PI_kq_2 = ZERO
224                    Neg_PI_kq_2 = ZERO
225                 ENDIF
226     
227     
228                 IF(K_Turb_G(IJK) > Small_number .AND. &
229                    E_Turb_G(IJK) > Small_number) THEN
230     
231     ! Start Adding source terms to the K equation
232                    K_Turb_G_c (IJK) = (EP_g(IJK)*Pos_Tauij_gDUi_gODxj +  &
233                                        Pos_PI_kq_2 )
234                    K_Turb_G_p (IJK) =(-EP_g(IJK)*Neg_Tauij_gDUi_gODxj +  &
235                       EP_g(IJK)*RO_G(IJK)*E_Turb_G(IJK)+Neg_PI_kq_2)/ &
236                       K_Turb_G(IJK)
237     
238     
239     ! Calculate velocity components at i, j, k to be used in wall functions
240                 IMJK = IM_OF(IJK)
241                 IPJK = IP_OF(IJK)
242                 IJMK = JM_OF(IJK)
243                 IJPK = JP_OF(IJK)
244                 IJKM = KM_OF(IJK)
245                 IJKP = KP_OF(IJK)
246                 UGC = AVG_X_E(U_G(IMJK),U_G(IJK),I)
247                 VGC = AVG_Y_N(V_G(IJMK),V_G(IJK))
248                 WGC = AVG_Z_T(W_G(IJKM),W_G(IJK))
249     
250     ! Implementing wall functions targeted to fluid cells next to walls...
251     ! Setting Source and sink terms in the K equation since the production
252     ! in the K eq. due to shear should include the LOG law of the wall.
253     ! North/South wall
254                    IF(WALL_AT(IJPK).OR.WALL_AT(IJMK)) THEN
255                       Check_Log = LOG(9.81*C_mu**0.25*                     &
256                          RO_G(IJK)*SQRT(K_Turb_G(IJK))*DY(J)/2.0D+0/       &
257                          Mu_g(IJK))
258                       IF(Check_Log .LE. ZERO) THEN
259                          K_Turb_G_c (IJK) = zero
260                          K_Turb_G_p (IJK) = zero
261                       ELSE
262                          K_Turb_G_c (IJK) = SQRT(C_mu) * 2.D+0/DY(J)*      &
263                             MAX(ABS(UGC),ABS(WGC))*EP_g(IJK)*RO_G(IJK)*    &
264                             K_Turb_G(IJK)/Check_Log
265                          K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)*         &
266                             C_mu**0.75*K_Turb_G(IJK)**1.5/DY(J)*           &
267                             2.0D+0/Kappa))/K_Turb_G(IJK)
268                       ENDIF !for check_log less than zero
269     ! Top/Bottom wall
270                    ELSEIF(WALL_AT(IJKP).OR.WALL_AT(IJKM)) THEN
271                       Check_Log = LOG(9.81*C_mu**0.25*                     &
272                          RO_G(IJK)*SQRT(K_Turb_G(IJK))/                    &
273                          (ODZ(K)*OX(I)*2.0D+0)/Mu_g(IJK))
274                       IF(Check_Log .LE. ZERO) THEN
275                          K_Turb_G_c (IJK) = zero
276                          K_Turb_G_p (IJK) = zero
277                       ELSE
278                          K_Turb_G_c (IJK) = SQRT(C_mu)*(ODZ(K)*OX(I)*      &
279                             2.0D+0)* MAX(ABS(UGC),ABS(VGC))*EP_g(IJK)*     &
280                             RO_G(IJK)*K_Turb_G(IJK)/Check_Log
281                          K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)*         &
282                             C_mu**0.75*K_Turb_G(IJK)**1.5/Kappa*           &
283                             (ODZ(K)*OX(I)*2.0D+0)))/K_Turb_G(IJK)
284                       ENDIF !for check_log less than zero
285                    ENDIF !For walls
286     
287     ! For Cylindrical cases, wall_at (IP) is a wall cell, but wall_at (IM) is
288     ! the axis of symmetry where wall functions obviously don't apply.
289                    IF(CYLINDRICAL) THEN
290                       IF (WALL_AT(IPJK))  THEN
291                          Check_Log = LOG(9.81*C_mu**0.25* RO_G(IJK)*       &
292                          SQRT(K_Turb_G(IJK))*DX(I)/2.0D+0/Mu_g(IJK))
293                          IF(Check_Log .LE. ZERO) THEN
294                             K_Turb_G_c (IJK) = zero
295                             K_Turb_G_p (IJK) = zero
296                          ELSE
297                             K_Turb_G_c (IJK) = SQRT(C_mu)*2.D+0/DX(I)*     &
298                                MAX(ABS(VGC),ABS(WGC)) *EP_g(IJK)*RO_G(IJK)*&
299                                K_Turb_G(IJK)/Check_Log
300                             K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)*      &
301                                C_mu**0.75*K_Turb_G(IJK)**1.5/DX(I)*        &
302                                2.0D+0/Kappa))/ K_Turb_G(IJK)
303                          ENDIF !for check_log less than zero
304                       ENDIF  ! for wall cells in I direction
305     ! East/West wall
306                    ELSEIF (WALL_AT(IPJK).OR.WALL_AT(IMJK)) THEN
307                       Check_Log = LOG(9.81*C_mu**0.25*RO_G(IJK)*           &
308                          SQRT(K_Turb_G(IJK))*DX(I)/2.0D+0/Mu_g(IJK))
309                       IF(Check_Log .LE. ZERO) THEN
310                          K_Turb_G_c (IJK) = zero
311                          K_Turb_G_p (IJK) = zero
312                       ELSE
313                          K_Turb_G_c (IJK) = SQRT(C_mu)*2.D+0/DX(I)*        &
314                             MAX(ABS(VGC),ABS(WGC))*EP_g(IJK)*RO_G(IJK)*    &
315                             K_Turb_G(IJK)/Check_Log
316                          K_Turb_G_p (IJK) = ((EP_g(IJK)*RO_G(IJK)*         &
317                             C_mu**0.75*K_Turb_G(IJK)**1.5/DX(I)*           &
318                             2.0D+0/Kappa))/K_Turb_G(IJK)
319                       ENDIF !for check_log less than zero
320                    ENDIF ! for cylindrical
321     
322                    IF(CUT_CELL_AT(IJK)) THEN
323                       Check_Log = LOG(9.81*C_mu**0.25* RO_G(IJK)*          &
324                       SQRT(K_Turb_G(IJK))*DELH_Scalar(IJK)/Mu_g(IJK))
325                       IF(Check_Log .LE. ZERO) THEN
326                          K_Turb_G_c (IJK) = zero
327                          K_Turb_G_p (IJK) = zero
328                       ELSE
329     !                     K_Turb_G_c (IJK) = SQRT(C_mu)/DELH_Scalar(IJK)*   &
330     !                        MAX(ABS(UGC),ABS(VGC),ABS(WGC))*               &
331     !                        EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK) /Check_Log
332                          K_Turb_G_c (IJK) = SQRT(C_mu)/DELH_Scalar(IJK)*   &
333                             DSQRT(UGC**2 + VGC**2 + WGC**2) *              &
334                             EP_g(IJK)*RO_G(IJK)*K_Turb_G(IJK)/Check_Log
335     
336                          K_Turb_G_p (IJK) = (EP_g(IJK)*RO_G(IJK)*          &
337                             C_mu**0.75*K_Turb_G(IJK)**1.5)/                &
338                            (DELH_Scalar(IJK)*Kappa)/K_Turb_G(IJK)
339                       ENDIF !for check_log less than zero
340                    ENDIF
341     
342     
343     ! Diffusion coefficient for turbulent kinetic energy (K)
344                    Dif_K_Turb_G(IJK) = EP_g(IJK)* &
345                       (MU_G(IJK) + Mu_gas_t /Sigma_k)
346     
347     ! Add here Dissipation of Turbulence source terms
348                    E_Turb_G_c (IJK) = (Ceps_1 *&
349                       EP_g(IJK)*Pos_Tauij_gDUi_gODxj*                      &
350                       E_Turb_G(IJK)/K_Turb_G(IJK) +                        &
351                       C_Eps_3*Pos_PI_kq_2*E_Turb_G(IJK)/K_Turb_G(IJK))
352     
353     ! Pope Correction in Xsi_Pope, Add it to E_Turb_G_c to use this option.
354                       ! + C_Eps_Pope*RO_G(IJK)*EP_g(IJK)*&
355                        !  ZMAX(Xsi_Pope)
356     
357                    E_Turb_G_p (IJK) = -Ceps_1 * &
358                       EP_g(IJK)*Neg_Tauij_gDUi_gODxj /K_Turb_G(IJK) +      &
359                       Ceps_2 * EP_g(IJK) *RO_G(IJK)*                       &
360                       E_Turb_G(IJK)/K_Turb_G(IJK) +                        &
361                       C_Eps_3*(Neg_PI_kq_2)/K_Turb_G(IJK)
362     
363     ! Diffusion coefficient for Dissipation of turbulent energy (Epsilon)
364                    Dif_E_Turb_G(IJK) =EP_g(IJK)* &
365                       (MU_G(IJK) + Mu_gas_t /Sigma_E)
366     
367                 ELSE
368                    K_Turb_G_c (IJK) = zero
369                    K_Turb_G_p (IJK) = zero
370                    E_Turb_G_c (IJK) = zero
371                    E_Turb_G_p (IJK) = zero
372                    Dif_K_Turb_G(IJK) = zero
373                    Dif_E_Turb_G(IJK) = zero
374                 ENDIF !for K_turb_G and E_Turb_G having very small numbers
375              ENDIF
376           ENDDO
377     
378           RETURN
379           END SUBROUTINE K_Epsilon_PROP
380