File: /nfs/home/0/users/jenkins/mfix.git/model/solve_k_epsilon_eq.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE SOLVE_K_Epsilon_EQ(IER)
21
22
23
24
25
26
27
28
29
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
62
63
64
65
66
67
68 INTEGER IER
69
70
71
72
73
74
75
76
77
78 INTEGER m, I, J, K, LC
79
80 DOUBLE PRECISION apo
81
82
83 DOUBLE PRECISION res1, mres1, num_res, den_res
84 INTEGER ires1
85
86
87 INTEGER IJK
88
89
90 INTEGER LEQM, LEQI
91
92
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
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 = (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
118
119
120
121
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
132
133
134 = 0
135 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
136
137
138
139 DO IJK = IJKSTART3, IJKEND3
140
141 = 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 (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
190
191
192
193
194
195
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
204
205
206
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
216
217
218
219 CALL INIT_AB_M (A_M, B_M, IJKMAX2, M, IER)
220
221 DO IJK = IJKSTART3, IJKEND3
222
223 = 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 (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
259
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
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
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
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
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
348
349
350
351
352
353
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
362
363
364
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
380
381
382