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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SOLVE_K_Epsilon_EQ                                      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     !                                                                      C
9     !                                                                      C
10     !  Literature/Document References:                                     C
11     !                                                                      C
12     !                                                                      C
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14           SUBROUTINE SOLVE_K_Epsilon_EQ(IER)
15     
16     ! Modules
17     !---------------------------------------------------------------------//
18           use ambm, only: a_m, b_m, lock_ambm, unlock_ambm
19           use bc, only: bc_k_turb_g, bc_e_turb_g
20           use compar, only: ijkstart3, ijkend3
21           use constant, only: to_si
22           use fldvar, only: rop_g, ep_g, u_g, v_g, w_g
23           use fldvar, only: k_turb_g, k_turb_go
24           use fldvar, only: e_turb_g, e_turb_go
25           use functions, only: fluid_at, zmax
26           use indices, only: i_of, j_of, k_of
27           use geometry, only: vol, ijkmax2
28           use leqsol, only: leq_it, leq_method, leq_sweep, leq_pc, leq_tol
29           use mflux, only: flux_ge, flux_gn, flux_gt
30           use mflux, only: flux_gse, flux_gsn, flux_gst
31           use param, only: dimension_bc, dimension_3
32           use param1, only: zero
33           use residual, only: resid, max_resid, ijk_resid
34           use residual, only: num_resid, den_resid
35           use residual, only: resid_ke
36           use run, only: discretize, k_epsilon, added_mass, odt
37           use rxns, only: sum_r_g
38           use toleranc, only: zero_ep_s
39           use turb, only: dif_k_turb_g, k_turb_g_c, k_turb_g_p
40           use turb, only: dif_e_turb_g, e_turb_g_c, e_turb_g_p
41           use ur_facs, only: ur_fac
42           use usr_src, only: call_usr_source, calc_usr_source
43           use usr_src, only: k_epsilon_e, k_epsilon_k
44           IMPLICIT NONE
45     
46     ! Dummy arguments
47     !---------------------------------------------------------------------//
48     ! Error index
49           INTEGER, INTENT(INOUT) :: IER
50     
51     ! Local variables
52     !---------------------------------------------------------------------//
53     ! phase index
54           INTEGER :: M
55     ! indices
56           INTEGER :: IJK, I, J, K
57     ! loop counter
58           INTEGER :: LC
59     !
60           DOUBLE PRECISION :: apo
61     ! linear equation solver method and iterations
62           INTEGER :: LEQM, LEQI
63     ! temporary variables in residual computation
64           DOUBLE PRECISION :: res1, mres1, num_res, den_res
65           INTEGER :: ires1
66     ! A default zero flux will be defined for both K & Epsilon at walls
67           DOUBLE PRECISION :: BC_hw_K_Turb_G (DIMENSION_BC)
68           DOUBLE PRECISION :: BC_hw_E_Turb_G (DIMENSION_BC)
69           DOUBLE PRECISION :: BC_K_Turb_GW (DIMENSION_BC)
70           DOUBLE PRECISION :: BC_E_Turb_GW (DIMENSION_BC)
71           DOUBLE PRECISION :: BC_C_K_Turb_G (DIMENSION_BC)
72           DOUBLE PRECISION :: BC_C_E_Turb_G (DIMENSION_BC)
73     
74     ! small value of K or E, 1 cm2/s2 = 1e-4 m2/s2 = 1e-4 m2/s3
75           DOUBLE PRECISION smallTheta
76     
77     ! source vector: coefficient of dependent variable
78     ! becomes part of a_m matrix; must be positive
79           DOUBLE PRECISION :: S_P(DIMENSION_3)
80     ! source vector: constant part becomes part of b_m vector
81           DOUBLE PRECISION :: S_C(DIMENSION_3)
82     !
83           character(LEN=8) :: Vname
84     !---------------------------------------------------------------------//
85     
86           IF( .NOT. K_Epsilon) RETURN
87     
88           call lock_ambm
89     
90           smallTheta = (to_SI)**4 * ZERO_EP_S
91           RESID(RESID_ke,0) = ZERO
92           NUM_RESID(RESID_ke,0) = ZERO
93           DEN_RESID(RESID_ke,0) = ZERO
94           MAX_RESID(RESID_ke,0) = ZERO
95           IJK_RESID(RESID_ke,0) = 0
96     
97     ! Setting default zero flux for K & Epsilon since we use wall functions.
98     ! If an expert user want to use Low Re K-Epilon model and needs to set
99     ! the turbulence quantities to zero at walls, then set the hw's to UNDEFINE will
100     ! do it. All the variables below can be changed in the same way as in the
101     ! MFIX data file in the boundary conditions section.
102           DO LC = 1, DIMENSION_BC
103             BC_hw_K_Turb_G (LC) = ZERO
104             BC_K_Turb_GW (LC) = ZERO
105             BC_C_K_Turb_G (LC) = ZERO
106             BC_hw_E_Turb_G (LC) = ZERO
107             BC_E_Turb_GW (LC) = ZERO
108             BC_C_E_Turb_G (LC) = ZERO
109           ENDDO
110     ! End of setting default zero flux for K & Epsilon wall boundary conditions
111     
112     ! Equations solved for gas phase, thus M = 0
113           M = 0
114           CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
115     
116     ! Solve the K_Turb_G equation first
117     ! ---------------------------------------------------------------->>>
118           DO IJK = IJKSTART3, IJKEND3
119              IF (FLUID_AT(IJK)) THEN
120                 APO = ROP_G(IJK)*VOL(IJK)*ODT
121                 S_P(IJK) = APO + (ZMAX(SUM_R_G(IJK)) + K_Turb_G_p(IJK))*&
122                            VOL(IJK)
123                 S_C(IJK) = APO*K_Turb_GO(IJK) + &
124                            K_Turb_G(IJK)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) + &
125                            K_Turb_G_c(IJK)*VOL(IJK)
126              ELSE
127                 S_P(IJK) = ZERO
128                 S_C(IJK) = ZERO
129              ENDIF
130           ENDDO
131     
132           IF(.NOT.ADDED_MASS) THEN
133              CALL CONV_DIF_PHI (K_Turb_G, DIF_K_Turb_G, DISCRETIZE(9), &
134                                 U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, &
135                                 M, A_M, B_M)
136           ELSE
137              CALL CONV_DIF_PHI (K_Turb_G, DIF_K_Turb_G, DISCRETIZE(9), &
138                                 U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, &
139                                 M, A_M, B_M)
140           ENDIF
141     
142           CALL BC_PHI (K_Turb_G, BC_K_Turb_G, BC_K_Turb_GW, BC_HW_K_Turb_G, &
143                        BC_C_K_Turb_G, M, A_M, B_M)
144           CALL SOURCE_PHI (S_P, S_C, EP_G, K_Turb_G, M, A_M, B_M)
145     
146     ! calculate any usr source terms
147           IF (CALL_USR_SOURCE(9)) CALL CALC_USR_SOURCE(K_EPSILON_K,&
148                                 A_M, B_M, lM=0)
149     
150           CALL CALC_RESID_S (K_Turb_G, A_M, B_M, M, num_res, den_res, res1, &
151                              mres1, ires1, ZERO)
152     
153           RESID(RESID_ke,0) = RESID(RESID_ke,0)+res1
154           NUM_RESID(RESID_ke,0) = NUM_RESID(RESID_ke,0)+num_res
155           DEN_RESID(RESID_ke,0) = DEN_RESID(RESID_ke,0)+den_res
156           if(mres1 .gt. MAX_RESID(RESID_ke,0)) then
157              MAX_RESID(RESID_ke,0) = mres1
158              IJK_RESID(RESID_ke,0) = ires1
159           endif
160     
161           CALL UNDER_RELAX_S (K_Turb_G, A_M, B_M, M, UR_FAC(9))
162     
163     !      call check_ab_m(a_m, b_m, m, .false., ier)
164     !      call write_ab_m(a_m, b_m, ijkmax2, m, ier)
165     !      write(*,*) res1, mres1, ires1
166     !      call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
167     !                       DO_K, ier)
168     
169           CALL ADJUST_LEQ (RESID(RESID_ke,0), LEQ_IT(9), LEQ_METHOD(9), &
170                    LEQI, LEQM)
171     
172           write(Vname, '(A,I2)')'K_Turb_G'
173           CALL SOLVE_LIN_EQ (Vname, 9, K_Turb_G, A_M, B_M, M, LEQI, LEQM, &
174                              LEQ_SWEEP(9), LEQ_TOL(9),  LEQ_PC(9), IER)
175     
176     !      call out_array(K_Turb_G, Vname)
177     
178     ! remove small negative K values generated by linear solver
179     ! same as adjust_theta.f
180           DO IJK = IJKSTART3, IJKEND3
181              IF (FLUID_AT(IJK)) THEN
182                 IF(K_Turb_G(IJK) < smallTheta) K_Turb_G(IJK) = smallTheta
183              ENDIF
184           ENDDO
185     
186     
187     ! Now solve the E_Turb_G (dissipation) equation.
188     ! ---------------------------------------------------------------->>>
189     ! Initiate (again) the Am Bm matrix.
190     ! This has to be done for every scalar equation.
191           CALL INIT_AB_M (A_M, B_M, IJKMAX2, M)
192     
193           DO IJK = IJKSTART3, IJKEND3
194              I = I_OF(IJK)
195              J = J_OF(IJK)
196              K = K_OF(IJK)
197              IF (FLUID_AT(IJK)) THEN
198                 APO = ROP_G(IJK)*VOL(IJK)*ODT
199                 S_P(IJK) = APO + (ZMAX(SUM_R_G(IJK)) + E_Turb_G_p(IJK))*&
200                            VOL(IJK)
201                 S_C(IJK) = APO*E_Turb_GO(IJK) + &
202                            E_Turb_G(IJK)*ZMAX((-SUM_R_G(IJK)))*VOL(IJK) + &
203                            E_Turb_G_c(IJK)*VOL(IJK)
204              ELSE
205                 S_P(IJK) = ZERO
206                 S_C(IJK) = ZERO
207              ENDIF
208           ENDDO
209     
210           IF(.NOT.ADDED_MASS) THEN
211              CALL CONV_DIF_PHI (E_Turb_G, DIF_E_Turb_G, DISCRETIZE(9), &
212                                 U_G, V_G, W_G, Flux_gE, Flux_gN, Flux_gT, &
213                                 M, A_M, B_M)
214           ELSE
215              CALL CONV_DIF_PHI (E_Turb_G, DIF_E_Turb_G, DISCRETIZE(9), &
216                                 U_G, V_G, W_G, Flux_gSE, Flux_gSN, Flux_gST, &
217                                 M, A_M, B_M)
218           ENDIF
219     
220           CALL BC_PHI (E_Turb_G, BC_E_Turb_G, BC_E_Turb_GW, BC_HW_E_Turb_G, &
221                        BC_C_E_Turb_G, M, A_M, B_M)
222     
223           CALL SOURCE_PHI (S_P, S_C, EP_G, E_Turb_G, M, A_M, B_M)
224     
225     ! calculate any usr source terms
226           IF (CALL_USR_SOURCE(9)) CALL CALC_USR_SOURCE(K_EPSILON_E,&
227                                 A_M, B_M, lM=0)
228     
229     ! set epsilon in fluid cells adjacent to wall cells
230           CALL SOURCE_K_EPSILON_BC(A_M, B_M, M)
231     
232           CALL CALC_RESID_S (E_Turb_G, A_M, B_M, M, num_res, den_res, res1, &
233                              mres1, ires1, ZERO)
234     
235           RESID(RESID_ke,0) = RESID(RESID_ke,0)+res1
236           NUM_RESID(RESID_ke,0) = NUM_RESID(RESID_ke,0)+num_res
237           DEN_RESID(RESID_ke,0) = DEN_RESID(RESID_ke,0)+den_res
238           if(mres1 .gt. MAX_RESID(RESID_ke,0))then
239             MAX_RESID(RESID_ke,0) = mres1
240             IJK_RESID(RESID_ke,0) = ires1
241           endif
242     
243           CALL UNDER_RELAX_S (E_Turb_G, A_M, B_M, M, UR_FAC(9))
244     
245     !      call check_ab_m(a_m, b_m, m, .false., ier)
246     !      call write_ab_m(a_m, b_m, ijkmax2, m, ier)
247     !      write(*,*) res1, mres1, ires1
248     !      call test_lin_eq(ijkmax2, ijmax2, imax2, a_m(1, -3, 0), 1, &
249     !                       DO_K, ier)
250     
251            CALL ADJUST_LEQ (RESID(RESID_ke,0), LEQ_IT(9), LEQ_METHOD(9), &
252                             LEQI, LEQM)
253     
254           write(Vname, '(A,I2)')'E_Turb_G'
255           CALL SOLVE_LIN_EQ (Vname, 9, E_Turb_G, A_M, B_M, M, LEQI, LEQM, &
256                              LEQ_SWEEP(9), LEQ_TOL(9), LEQ_PC(9), IER)
257     
258     !      call out_array(E_Turb_G, Vname)
259     
260     ! remove small negative Epsilon values generated by linear solver
261     ! same as adjust_theta.f
262           DO IJK = IJKSTART3, IJKEND3
263              IF (FLUID_AT(IJK)) THEN
264                 IF(E_Turb_G(IJK) < smallTheta) E_Turb_G(IJK) = smallTheta
265              ENDIF
266           ENDDO
267     
268           call unlock_ambm
269     
270           RETURN
271           END SUBROUTINE SOLVE_K_Epsilon_EQ
272     
273     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
274     !                                                                      C
275     !  Subroutine: SOURCE_K_EPSILON_BC                                     C
276     !  Purpose: When implementing the wall functions, the epsilon          C
277     !  (dissipation) value at the fluid cell near the walls needs to be    C
278     !  set.                                                                C
279     !                                                                      C
280     !                                                                      C
281     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
282           SUBROUTINE SOURCE_K_EPSILON_BC(A_M, B_M, M)
283     
284     ! Modules
285     !---------------------------------------------------------------------//
286           use compar, only: ijkstart3, ijkend3
287           use cutcell, only: cut_cell_at, delh_scalar
288           use fldvar, only: k_turb_G, e_turb_g
289           use functions, only: fluid_at, wall_at
290           use functions, only: im_of, jm_of, km_of
291           use functions, only: ip_of, jp_of, kp_of
292           use geometry, only: dy, odz, ox, dx, cylindrical
293           use indices, only: i_of, j_of, k_of
294           use param, only: dimension_3, dimension_m
295           use param1, only: zero, one
296           IMPLICIT NONE
297     
298     ! Dummy Arguments
299     !---------------------------------------------------------------------//
300     ! Septadiagonal matrix A_m
301           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
302     ! Vector b_m
303           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
304     ! phase index
305           INTEGER, INTENT(IN) :: M
306     ! Local variables
307     !---------------------------------------------------------------------//
308           INTEGER :: IJK, I, J, K
309     !---------------------------------------------------------------------//
310     
311           DO IJK = IJKSTART3, IJKEND3
312              I = I_OF(IJK)
313              J = J_OF(IJK)
314              K = K_OF(IJK)
315              IF (FLUID_AT(IJK)) THEN
316                 IF(WALL_AT(JP_OF(IJK)).OR.WALL_AT(JM_OF(IJK))) THEN
317                    A_M(IJK,1,M) = ZERO
318                    A_M(IJK,-1,M) = ZERO
319                    A_M(IJK,2,M) = ZERO
320                    A_M(IJK,-2,M) = ZERO
321                    A_M(IJK,3,M) = ZERO
322                    A_M(IJK,-3,M) = ZERO
323                    A_M(IJK,0,M) = -ONE
324                    B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/DY(J) &
325                                *2.0D+0/0.42D+0
326                 ELSEIF(WALL_AT(KP_OF(IJK)).OR.WALL_AT(KM_OF(IJK))) THEN
327                    A_M(IJK,1,M) = ZERO
328                    A_M(IJK,-1,M) = ZERO
329                    A_M(IJK,2,M) = ZERO
330                    A_M(IJK,-2,M) = ZERO
331                    A_M(IJK,3,M) = ZERO
332                    A_M(IJK,-3,M) = ZERO
333                    A_M(IJK,0,M) = -ONE
334                    B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)* &
335                                 (ODZ(K)*OX(I)*2.0D+0)/0.42D+0
336                 ENDIF  !for identifying wall cells in J or K direction
337     
338                 IF(CYLINDRICAL) THEN
339                    IF (WALL_AT(IP_OF(IJK)))  THEN
340                       A_M(IJK,1,M) = ZERO
341                       A_M(IJK,-1,M) = ZERO
342                       A_M(IJK,2,M) = ZERO
343                       A_M(IJK,-2,M) = ZERO
344                       A_M(IJK,3,M) = ZERO
345                       A_M(IJK,-3,M) = ZERO
346                       A_M(IJK,0,M) = -ONE
347                       B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/&
348                                    DX(I)*2.0D+0/0.42D+0
349     
350                    ENDIF! for wall cells in I direction
351                 ELSEIF (WALL_AT(IP_OF(IJK)).OR.WALL_AT(IM_OF(IJK))) THEN
352                    A_M(IJK,1,M) = ZERO
353                    A_M(IJK,-1,M) = ZERO
354                    A_M(IJK,2,M) = ZERO
355                    A_M(IJK,-2,M) = ZERO
356                    A_M(IJK,3,M) = ZERO
357                    A_M(IJK,-3,M) = ZERO
358                    A_M(IJK,0,M) = -ONE
359                    B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/DX(I) &
360                                  *2.0D+0/0.42D+0
361                 ENDIF ! for cylindrical
362     
363                 IF(CUT_CELL_AT(IJK)) THEN
364                    A_M(IJK,1,M) = ZERO
365                    A_M(IJK,-1,M) = ZERO
366                    A_M(IJK,2,M) = ZERO
367                    A_M(IJK,-2,M) = ZERO
368                    A_M(IJK,3,M) = ZERO
369                    A_M(IJK,-3,M) = ZERO
370                    A_M(IJK,0,M) = -ONE
371                    B_M(IJK,M) =-((0.09D+0)**0.75*K_Turb_G(IJK)**1.5)/&
372                                 (0.42D+0*DELH_Scalar(IJK))
373                 ENDIF
374              ENDIF  !for fluid at ijk
375           ENDDO   ! enddo ijk
376           RETURN
377           END SUBROUTINE SOURCE_K_EPSILON_BC
378