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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SOLVE_K_Epsilon_EQ(IER)                                C
4     !  Purpose: Solve K & Epsilon equations for a turbulent flow           C
5     !                                                                      C
6     !                                                                      C
7     !  Author: S. Benyahia                                Date: MAY-13-04  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !                                                                      C
11     !  Literature/Document References:                                     C
12     !                                                                      C
13     !  Variables referenced:                                               C
14     !  Variables modified:                                                 C
15     !                                                                      C
16     !  Local variables:                                                    C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     !
20           SUBROUTINE SOLVE_K_Epsilon_EQ(IER)
21     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
22     !...Switches: -xf
23     !
24     !  Include param.inc file to specify parameter values
25     !
26     !-----------------------------------------------
27     !   M o d u l e s
28     !-----------------------------------------------
29           !USE matrix
30           USE bc
31           USE compar
32           USE constant
33           USE cutcell
34           USE drag
35           USE energy
36           USE fldvar
37           USE fun_avg
38           USE functions
39           USE geometry
40           USE indices
41           USE leqsol
42           USE mflux
43           USE output
44           USE param
45           USE param1
46           USE pgcor
47           USE physprop
48           USE pscor
49           USE residual
50           USE run
51           USE rxns
52           USE toleranc
53           USE turb
54           USE ur_facs
55           USE usr
56           Use ambm
57           Use tmp_array, S_p => Array1, S_c => Array2, EPs => Array3, VxGama => Array4
58     
59           IMPLICIT NONE
60     !-----------------------------------------------
61     !   G l o b a l   P a r a m e t e r s
62     !-----------------------------------------------
63     !-----------------------------------------------
64     !   D u m m y   A r g u m e n t s
65     !-----------------------------------------------
66     !
67     !                      Error index
68           INTEGER          IER
69     !-----------------------------------------------
70     !   L o c a l   P a r a m e t e r s
71     !-----------------------------------------------
72     !-----------------------------------------------
73     !   L o c a l   V a r i a b l e s
74     !-----------------------------------------------
75     !
76     !
77     !                      phase index
78           INTEGER          m, I, J, K, LC
79     !
80           DOUBLE PRECISION apo
81     
82     !                      temporary variables in residual computation
83           DOUBLE PRECISION res1, mres1, num_res, den_res
84           INTEGER          ires1
85     !
86     !                      Indices
87           INTEGER          IJK
88     !
89     !                      linear equation solver method and iterations
90           INTEGER          LEQM, LEQI
91     
92     !                      A default zero flux will be defined for both K & Epsilon at walls
93           DOUBLE PRECISION BC_hw_K_Turb_G (DIMENSION_BC),  BC_hw_E_Turb_G (DIMENSION_BC)
94           DOUBLE PRECISION BC_K_Turb_GW (DIMENSION_BC),   BC_E_Turb_GW (DIMENSION_BC)
95           DOUBLE PRECISION BC_C_K_Turb_G (DIMENSION_BC),  BC_C_E_Turb_G (DIMENSION_BC)
96     !
97     !                      small value of K or E, 1 cm2/s2 = 1e-4 m2/s2 = 1e-4 m2/s3
98           DOUBLE PRECISION smallTheta
99     !
100           character(LEN=8) :: Vname
101     !-----------------------------------------------
102     
103           IF( .NOT. K_Epsilon) RETURN
104     
105           call lock_ambm
106           call lock_tmp_array
107     
108     !
109           smallTheta = (to_SI)**4 * ZERO_EP_S
110     
111           RESID(RESID_ke,0) = ZERO
112           NUM_RESID(RESID_ke,0) = ZERO
113           DEN_RESID(RESID_ke,0) = ZERO
114           MAX_RESID(RESID_ke,0) = ZERO
115           IJK_RESID(RESID_ke,0) = 0
116     
117     ! Setting default zero flux for K & Epsilon since we use wall functions.
118     ! If an expert user want to use Low Re K-Epilon model and needs to set
119     ! the turbulence quatities to zero at walls, then set the hw's to UNDEFINE will
120     ! do it. All the variables below can be changed in the same way as in the
121     ! MFIX data file in the boundary conditions section.
122     !
123           DO LC = 1, DIMENSION_BC
124             BC_hw_K_Turb_G (LC) = ZERO
125             BC_hw_E_Turb_G (LC) = ZERO
126             BC_K_Turb_GW (LC) = ZERO
127             BC_E_Turb_GW (LC) = ZERO
128             BC_C_K_Turb_G (LC) = ZERO
129             BC_C_E_Turb_G (LC) = ZERO
130           ENDDO
131     ! End of setting default zero flux for K & Epsilon wall boundary conditions
132     !
133     ! Equations solved for gas phase, thus M = 0
134           M = 0
135               CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
136     
137     ! Solve first fot the K_Turb_G Equation
138     
139               DO IJK = IJKSTART3, IJKEND3
140     !
141                  I = I_OF(IJK)
142                  J = J_OF(IJK)
143                  K = K_OF(IJK)
144                    IF (FLUID_AT(IJK)) THEN
145     
146                       APO = ROP_G(IJK)*VOL(IJK)*ODT
147                       S_P(IJK) = APO + (  ZMAX(SUM_R_G(IJK)) &
148                                         + K_Turb_G_p(IJK) )*VOL(IJK)
149                       S_C(IJK) =   APO*K_Turb_GO(IJK) &
150                                 + K_Turb_G(IJK)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) &
151                                 + K_Turb_G_c(IJK) *VOL(IJK)
152                    ELSE
153     !
154                       S_P(IJK) = ZERO
155                       S_C(IJK) = ZERO
156     !
157                    ENDIF
158     
159                 END DO
160     
161                 IF(.NOT.ADDED_MASS) THEN
162                    CALL CONV_DIF_PHI (K_Turb_G, DIF_K_Turb_G, DISCRETIZE(9), &
163                                    U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, M, A_M, B_M, IER)
164                 ELSE
165                    CALL CONV_DIF_PHI (K_Turb_G, DIF_K_Turb_G, DISCRETIZE(9), &
166                                    U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, M, A_M, B_M, IER)
167                 ENDIF
168     !
169     !
170                 CALL BC_PHI (K_Turb_G, BC_K_Turb_G, BC_K_Turb_GW, BC_HW_K_Turb_G, &
171                              BC_C_K_Turb_G, M, A_M, B_M, IER)
172     !
173     !
174                 CALL SOURCE_PHI (S_P, S_C, EP_G, K_Turb_G, M, A_M, B_M)
175     !
176                 CALL CALC_RESID_S (K_Turb_G, A_M, B_M, M, num_res, den_res, res1, &
177                                    mres1, ires1, ZERO, IER)
178     
179                 RESID(RESID_ke,0) = RESID(RESID_ke,0)+res1
180                 NUM_RESID(RESID_ke,0) = NUM_RESID(RESID_ke,0)+num_res
181                 DEN_RESID(RESID_ke,0) = DEN_RESID(RESID_ke,0)+den_res
182                 if(mres1 .gt. MAX_RESID(RESID_ke,0))then
183                   MAX_RESID(RESID_ke,0) = mres1
184                   IJK_RESID(RESID_ke,0) = ires1
185                 endif
186     !
187                 CALL UNDER_RELAX_S (K_Turb_G, A_M, B_M, M, UR_FAC(9), IER)
188     !
189     !          call check_ab_m(a_m, b_m, m, .false., ier)
190     !          call write_ab_m(a_m, b_m, ijkmax2, m, ier)
191     !          write(*,*) res1, mres1, &
192     !           ires1
193     !
194     !          call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, DO_K, &
195     !          ier)
196     !
197                 CALL ADJUST_LEQ (RESID(RESID_ke,0), LEQ_IT(9), LEQ_METHOD(9), &
198                    LEQI, LEQM, IER)
199     !
200                 write(Vname, '(A,I2)')'K_Turb_G'
201                 CALL SOLVE_LIN_EQ (Vname, 9, K_Turb_G, A_M, B_M, M, LEQI, LEQM, &
202                                  LEQ_SWEEP(9), LEQ_TOL(9),  LEQ_PC(9), IER)
203     !          call out_array(K_Turb_G, Vname)
204     !
205     ! remove small negative K values generated by linear solver
206     ! same as adjust_theta.f
207     !
208                DO IJK = IJKSTART3, IJKEND3
209                 IF (FLUID_AT(IJK)) THEN
210                  IF(K_Turb_G(IJK) < smallTheta) K_Turb_G(IJK) = smallTheta
211                 ENDIF
212                END DO
213     
214     !!!!!!!!!!!!!
215     ! Now, we'll solve for the E_Turb_G (dissipation) Equation.
216     !!!!!!!!!!!!!
217     !
218     ! Initiate (again) the Am Bm matrix. This has to be done for every scalar equation.
219               CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
220     
221               DO IJK = IJKSTART3, IJKEND3
222     !
223                  I = I_OF(IJK)
224                  J = J_OF(IJK)
225                  K = K_OF(IJK)
226                    IF (FLUID_AT(IJK)) THEN
227     
228                       APO = ROP_G(IJK)*VOL(IJK)*ODT
229                       S_P(IJK) = APO + (  ZMAX(SUM_R_G(IJK)) &
230                                         + E_Turb_G_p(IJK)  )*VOL(IJK)
231                       S_C(IJK) =   APO*E_Turb_GO(IJK) &
232                                 + E_Turb_G(IJK)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) &
233                                 + E_Turb_G_c(IJK) *VOL(IJK)
234                    ELSE
235     !
236                       S_P(IJK) = ZERO
237                       S_C(IJK) = ZERO
238     !
239                    ENDIF
240     
241                 END DO
242     
243                 IF(.NOT.ADDED_MASS) THEN
244                    CALL CONV_DIF_PHI (E_Turb_G, DIF_E_Turb_G, DISCRETIZE(9), &
245                                    U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, M, A_M, B_M, IER)
246                 ELSE
247                    CALL CONV_DIF_PHI (E_Turb_G, DIF_E_Turb_G, DISCRETIZE(9), &
248                                    U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, M, A_M, B_M, IER)
249                 ENDIF
250     !
251     !
252                 CALL BC_PHI (E_Turb_G, BC_E_Turb_G, BC_E_Turb_GW, BC_HW_E_Turb_G, &
253                              BC_C_E_Turb_G, M, A_M, B_M, IER)
254     !
255     !
256                 CALL SOURCE_PHI (S_P, S_C, EP_G, E_Turb_G, M, A_M, B_M)
257     !
258     ! When implementing the wall functions, The Epsilon (dissipation) value at the fluid cell
259     ! near the walls needs to be set.
260     
261                 DO IJK = IJKSTART3, IJKEND3
262                  I = I_OF(IJK)
263                  J = J_OF(IJK)
264                  K = K_OF(IJK)
265     !
266                    IF (FLUID_AT(IJK)) THEN
267     !
268                       IF(WALL_AT(JP_OF(IJK)).OR.WALL_AT(JM_OF(IJK))) THEN
269                          A_M(IJK,1,M) = ZERO
270                          A_M(IJK,-1,M) = ZERO
271                          A_M(IJK,2,M) = ZERO
272                          A_M(IJK,-2,M) = ZERO
273                          A_M(IJK,3,M) = ZERO
274                          A_M(IJK,-3,M) = ZERO
275                          A_M(IJK,0,M) = -ONE
276                          B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/DY(J) &
277                                      *2.0D+0/0.42D+0
278     
279                       ELSE IF(WALL_AT(KP_OF(IJK)).OR.WALL_AT(KM_OF(IJK))) THEN
280                          A_M(IJK,1,M) = ZERO
281                          A_M(IJK,-1,M) = ZERO
282                          A_M(IJK,2,M) = ZERO
283                          A_M(IJK,-2,M) = ZERO
284                          A_M(IJK,3,M) = ZERO
285                          A_M(IJK,-3,M) = ZERO
286                          A_M(IJK,0,M) = -ONE
287                          B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)* &
288                                       (ODZ(K)*OX(I)*2.0D+0)/0.42D+0
289                       ENDIF  !for identifying wall cells in J or K direction
290     
291                       IF(CYLINDRICAL) THEN
292                          IF (WALL_AT(IP_OF(IJK)))  THEN
293                             A_M(IJK,1,M) = ZERO
294                             A_M(IJK,-1,M) = ZERO
295                             A_M(IJK,2,M) = ZERO
296                             A_M(IJK,-2,M) = ZERO
297                             A_M(IJK,3,M) = ZERO
298                             A_M(IJK,-3,M) = ZERO
299                             A_M(IJK,0,M) = -ONE
300                             B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/DX(I) &
301                                         *2.0D+0/0.42D+0
302     
303                          ENDIF! for wall cells in I direction
304     
305                       ELSE IF (WALL_AT(IP_OF(IJK)).OR.WALL_AT(IM_OF(IJK))) THEN
306                          A_M(IJK,1,M) = ZERO
307                          A_M(IJK,-1,M) = ZERO
308                          A_M(IJK,2,M) = ZERO
309                          A_M(IJK,-2,M) = ZERO
310                          A_M(IJK,3,M) = ZERO
311                          A_M(IJK,-3,M) = ZERO
312                          A_M(IJK,0,M) = -ONE
313                          B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/DX(I) &
314                                        *2.0D+0/0.42D+0
315     
316                       ENDIF ! for cylindrical
317     
318                       IF(CUT_CELL_AT(IJK)) THEN
319                          A_M(IJK,1,M) = ZERO
320                          A_M(IJK,-1,M) = ZERO
321                          A_M(IJK,2,M) = ZERO
322                          A_M(IJK,-2,M) = ZERO
323                          A_M(IJK,3,M) = ZERO
324                          A_M(IJK,-3,M) = ZERO
325                          A_M(IJK,0,M) = -ONE
326     
327                          B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/(0.42D+0*DELH_Scalar(IJK))
328                       ENDIF
329     
330     !
331                    ENDIF  !for fluid at ijk
332                 ENDDO
333     
334                 CALL CALC_RESID_S (E_Turb_G, A_M, B_M, M, num_res, den_res, res1, &
335                                    mres1, ires1, ZERO, IER)
336     
337                 RESID(RESID_ke,0) = RESID(RESID_ke,0)+res1
338                 NUM_RESID(RESID_ke,0) = NUM_RESID(RESID_ke,0)+num_res
339                 DEN_RESID(RESID_ke,0) = DEN_RESID(RESID_ke,0)+den_res
340                 if(mres1 .gt. MAX_RESID(RESID_ke,0))then
341                   MAX_RESID(RESID_ke,0) = mres1
342                   IJK_RESID(RESID_ke,0) = ires1
343                 endif
344     !
345                 CALL UNDER_RELAX_S (E_Turb_G, A_M, B_M, M, UR_FAC(9), IER)
346     !
347     !          call check_ab_m(a_m, b_m, m, .false., ier)
348     !          call write_ab_m(a_m, b_m, ijkmax2, m, ier)
349     !          write(*,*) res1, mres1, &
350     !           ires1
351     !
352     !          call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, DO_K, &
353     !          ier)
354     !
355                 CALL ADJUST_LEQ (RESID(RESID_ke,0), LEQ_IT(9), LEQ_METHOD(9), &
356                    LEQI, LEQM, IER)
357     !
358                 write(Vname, '(A,I2)')'E_Turb_G'
359                 CALL SOLVE_LIN_EQ (Vname, 9, E_Turb_G, A_M, B_M, M, LEQI, LEQM, &
360                                  LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
361     !          call out_array(E_Turb_G, Vname)
362     !
363     ! remove small negative Epsilon values generated by linear solver
364     ! same as adjust_theta.f
365     !
366                DO IJK = IJKSTART3, IJKEND3
367                 IF (FLUID_AT(IJK)) THEN
368                  IF(E_Turb_G(IJK) < smallTheta) E_Turb_G(IJK) = smallTheta
369                 ENDIF
370                END DO
371     
372           call unlock_ambm
373           call unlock_tmp_array
374     
375           RETURN
376           END SUBROUTINE SOLVE_K_Epsilon_EQ
377     
378     
379     !// Comments on the modifications for DMP version implementation
380     !// 001 Include header file and common declarations for parallelization
381     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
382