File: /nfs/home/0/users/jenkins/mfix.git/model/source_granular_energy.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42 SUBROUTINE SOURCE_GRANULAR_ENERGY(SOURCELHS, &
43 SOURCERHS, IJK, M)
44
45
46
47
48 USE compar
49 USE constant
50 USE drag
51 USE fldvar
52 USE fun_avg
53 USE functions
54 USE geometry
55 USE indices
56 USE kintheory
57 USE mms
58 USE parallel
59 USE param
60 USE param1
61 USE physprop
62 USE rdf
63 USE run
64 USE solids_pressure
65 USE toleranc
66 USE trace
67 USE turb
68 USE visc_g
69 USE visc_s
70
71 IMPLICIT NONE
72
73
74
75
76 DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
77
78 DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
79
80 INTEGER, INTENT(IN) :: M
81
82 INTEGER, INTENT(IN) :: IJK
83
84
85
86
87 INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM, &
88 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
89
90 INTEGER :: MM
91
92 DOUBLE PRECISION :: Tau_12_st
93
94 DOUBLE PRECISION :: SUM_EpsGo
95
96 DOUBLE PRECISION :: VSLIP
97
98 DOUBLE PRECISION :: Knu_e, Knu_w, Knu_n, Knu_s, &
99 Knu_t, Knu_b
100
101 DOUBLE PRECISION :: S1_rhs, S2_rhs, S3_rhs, &
102 S5a_rhs, S5b_rhs, S5c_rhs, S5_rhs, &
103 S7_rhs
104
105 DOUBLE PRECISION :: S1_lhs, S3_lhs, S4_lhs, S6_lhs
106
107
108 = I_OF(IJK)
109 J = J_OF(IJK)
110 K = K_OF(IJK)
111 IM = Im1(I)
112 JM = Jm1(J)
113 KM = Km1(K)
114 IMJK = IM_OF(IJK)
115 IJMK = JM_OF(IJK)
116 IJKM = KM_OF(IJK)
117 IJKE = EAST_OF(IJK)
118 IJKW = WEST_OF(IJK)
119 IJKN = NORTH_OF(IJK)
120 IJKS = SOUTH_OF(IJK)
121 IJKT = TOP_OF(IJK)
122 IJKB = BOTTOM_OF(IJK)
123
124
125 = ZERO
126 S1_lhs = ZERO
127 S2_rhs = ZERO
128 S3_rhs = ZERO
129 S3_lhs = ZERO
130 S4_lhs = ZERO
131 S6_lhs = ZERO
132 S7_rhs = ZERO
133
134
135
136 = ZERO
137 DO MM = 1, MMAX
138 SUM_EpsGo = SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,MM,MM)
139 ENDDO
140
141
142
143 = P_S_C(IJK,M)*ZMAX(( -TRD_S_C(IJK,M) ))
144 S1_lhs = P_S_C(IJK,M)*ZMAX(( TRD_S_C(IJK,M) ))
145
146
147
148 = 2.d0*MU_S_C(IJK,M)*TRD_S2(IJK,M)
149
150
151
152 = (TRD_S_C(IJK,M)**2)*ZMAX( LAMBDA_S_C(IJK,M) )
153 S3_lhs = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_S_C(IJK,M) )
154
155
156 = (48.d0/DSQRT(PI))*ETA*(ONE-ETA)*ROP_S(IJK,M)*&
157 SUM_EpsGo*DSQRT(THETA_M(IJK,M))/D_P(IJK,M)
158
159
160 IF (USE_MMS) S4_lhs = ZERO
161
162
163
164 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
165 KT_TYPE_ENUM == AHMADI_1995) THEN
166 S6_lhs = 3.d0 * F_GS(IJK,M)
167 ELSE IF(SWITCH > ZERO .AND. RO_g0 /= ZERO) THEN
168 S6_lhs = SWITCH * 3.d0 * F_GS(IJK,M)
169 ENDIF
170
171
172
173 IF(KT_TYPE_ENUM == SIMONIN_1996) THEN
174 S7_rhs = F_GS(IJK,M)*K_12(IJK)
175 ELSEIF(KT_TYPE_ENUM == AHMADI_1995) THEN
176
177 IF(Ep_s(IJK,M) > DIL_EP_S .AND. F_GS(IJK,1) > small_number) THEN
178 Tau_12_st = Ep_s(IJK,M)*RO_S(IJK,M)/F_GS(IJK,1)
179 S7_rhs = 2.d0*F_GS(IJK,M) * (ONE/(ONE+Tau_12_st/ &
180 (Tau_1(ijk)+small_number)))*K_Turb_G(IJK)
181 ELSE
182 S7_rhs = ZERO
183 ENDIF
184 ELSEIF(SWITCH > ZERO .AND. RO_g0 /= ZERO) THEN
185 VSLIP = DSQRT( (U_S(IJK,M)-U_G(IJK))**2 + &
186 (V_S(IJK,M)-V_G(IJK))**2 + (W_S(IJK,M)-W_G(IJK))**2 )
187 S7_rhs = SWITCH*81.d0*EP_S(IJK,M)*(MU_G(IJK)*VSLIP)**2 /&
188 (G_0(IJK,M,M)*D_P(IJK,M)**3 * RO_S(IJK,M)*&
189 DSQRT(PI)*DSQRT( THETA_M(IJK,M)+SMALL_NUMBER ) )
190
191 IF (USE_MMS) S7_rhs = ZERO
192
193 ENDIF
194
195
196
197
198
199
200
201 IF (.FALSE.) THEN
202 Knu_e = AVG_X_H(Kphi_s(IJK,M), Kphi_s(IJKE,M),I)
203 Knu_w = AVG_X_H(Kphi_s(IJKW,M),Kphi_s(IJK,M),IM)
204 Knu_n = AVG_Y_H(Kphi_s(IJK,M), Kphi_s(IJKN,M),J)
205 Knu_s = AVG_Y_H(Kphi_s(IJKS,M),Kphi_s(IJK,M),JM)
206 Knu_t = AVG_Z_H(Kphi_s(IJK,M), Kphi_s(IJKT,M),K)
207 Knu_b = AVG_Z_H(Kphi_s(IJKB,M),Kphi_s(IJK,M),KM)
208
209 S5a_rhs = Knu_e*(EP_s(IJKE,M)-EP_s(IJK,M))*&
210 oDX_E(I)*AYZ(IJK) - &
211 Knu_w*(EP_s(IJK,M)-EP_s(IJKW,M))*&
212 oDX_E(IM)*AYZ(IMJK)
213 S5b_rhs = Knu_n*(EP_s(IJKN,M)-EP_s(IJK,M))*&
214 oDY_N(J)*AXZ(IJK) - &
215 Knu_s*(EP_s(IJK,M)-EP_s(IJKS,M))*&
216 oDY_N(JM)*AXZ(IJMK)
217 S5c_rhs = Knu_t*(EP_s(IJKT,M)-EP_s(IJK,M))*&
218 oX(I)*oDZ_T(K)*AXY(IJK) - &
219 Knu_b*(EP_s(IJK,M)-EP_s(IJKB,M))*&
220 oX(I)*oDZ_T(KM)*AXY(IJKM)
221 S5_rhs = S5a_rhs + S5b_rhs + S5c_rhs
222 ENDIF
223 S5_rhs = ZERO
224
225
226 SOURCERHS = (S1_rhs + S2_rhs + S3_rhs + S7_rhs)*VOL(IJK) + &
227 S5_rhs
228
229 SOURCELHS = ( ((S1_lhs + S3_lhs)/(THETA_M(IJK,M)+SMALL_NUMBER)) + &
230 S4_lhs + S6_lhs) * VOL(IJK)
231
232 RETURN
233 END SUBROUTINE SOURCE_GRANULAR_ENERGY
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256 SUBROUTINE SOURCE_GRANULAR_ENERGY_IA(SOURCELHS, &
257 SOURCERHS, IJK, M)
258
259
260
261
262 USE param
263 USE param1
264 USE parallel
265 USE physprop
266 USE run
267 USE drag
268 USE geometry
269 USE fldvar
270 USE visc_g
271 USE visc_s
272 USE trace
273 USE turb
274 USE indices
275 USE constant
276 USE toleranc
277 USE residual
278 use kintheory
279 USE compar
280 USE fun_avg
281 USE functions
282 IMPLICIT NONE
283
284
285
286
287 DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
288
289 DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
290
291 INTEGER, INTENT(IN) :: M
292
293 INTEGER, INTENT(IN) :: IJK
294
295
296
297
298 INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
299 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
300
301 INTEGER :: L, LM
302
303 DOUBLE PRECISION :: UsM_e, UsM_w, VsM_n, VsM_s, WsM_t, WsM_b,&
304 UsL_e, UsL_w, VsL_n, VsL_s, WsL_t, WsL_b,&
305 UsM_p, VsM_p, WsM_P, UsL_p, VsL_p, WsL_p
306
307 DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
308 NU_PM_B, NU_PM_p,&
309 NU_PL_E, NU_PL_W, NU_PL_N, NU_PL_S, NU_PL_T,&
310 NU_PL_B, NU_PL_p
311
312 DOUBLE PRECISION :: T_PL_E, T_PL_W, T_PL_N, T_PL_S, T_PL_T,&
313 T_PL_B, T_PL_p
314
315 DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL
316
317 DOUBLE PRECISION :: Knu_sL_e, Knu_sL_w, Knu_sL_n, Knu_sL_s, Knu_sL_t,&
318 Knu_sL_b, Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
319 Knu_sM_t, Knu_sM_b,&
320 Kvel_s_e, Kvel_s_w, Kvel_s_n, Kvel_s_s, Kvel_s_t,&
321 Kvel_s_b,&
322 Kth_sL_e, Kth_sL_w, Kth_sL_n, Kth_sL_s, Kth_sL_t,&
323 Kth_sL_b
324
325 DOUBLE PRECISION :: S10_rhs, S15_rhs, S16_rhs,&
326 S11_sum_rhs, S12_sum_rhs,&
327 S13_sum_rhs, S17_sum_rhs, S18_sum_rhs
328
329 DOUBLE PRECISION :: S14a_sum, S14b_sum, S14c_sum, &
330 S9_sum, s21a_sum, s21b_sum, s21c_sum
331
332 DOUBLE PRECISION :: S10_lhs, S16_lhs,&
333 S11_sum_lhs, S12_sum_lhs, S13_sum_lhs,&
334 S17_sum_lhs, S18_sum_lhs, S20_sum_lhs
335
336
337 = I_OF(IJK)
338 J = J_OF(IJK)
339 K = K_OF(IJK)
340 IM = Im1(I)
341 JM = Jm1(J)
342 KM = Km1(K)
343 IMJK = IM_OF(IJK)
344 IJMK = JM_OF(IJK)
345 IJKM = KM_OF(IJK)
346 IJKE = EAST_OF(IJK)
347 IJKW = WEST_OF(IJK)
348 IJKN = NORTH_OF(IJK)
349 IJKS = SOUTH_OF(IJK)
350 IJKT = TOP_OF(IJK)
351 IJKB = BOTTOM_OF(IJK)
352
353
354 = ZERO
355 S14a_sum = ZERO
356 S14b_sum = ZERO
357 S14c_sum = ZERO
358 S21a_sum = ZERO
359 S21b_sum = ZERO
360 S21c_sum = ZERO
361
362 S11_sum_rhs = ZERO
363 S12_sum_rhs = ZERO
364 S13_sum_rhs = ZERO
365 S17_sum_rhs = ZERO
366 S18_sum_rhs = ZERO
367
368 S11_sum_lhs = ZERO
369 S12_sum_lhs = ZERO
370 S13_sum_lhs = ZERO
371 S17_sum_lhs = ZERO
372 S18_sum_lhs = ZERO
373 S20_sum_lhs = ZERO
374
375 UsM_e = U_S(IJK,M)
376 UsM_w = U_S(IMJK,M)
377 VsM_n = V_S(IJK,M)
378 VsM_s = V_S(IJMK,M)
379 WsM_t = W_S(IJK,M)
380 WsM_b = W_S(IJKM,M)
381 UsM_p = AVG_X_E(U_S(IMJK,M),U_S(IJK,M),I)
382 VsM_p = AVG_Y_N(V_S(IJMK,M),V_S(IJK,M) )
383 WsM_p = AVG_Z_T(W_S(IJKM,M),W_S(IJK,M) )
384
385 D_PM = D_P(IJK,M)
386 M_PM = (Pi/6.d0)*D_PM**3 * RO_S(IJK,M)
387 NU_PM_p = ROP_S(IJK,M)/M_PM
388 NU_PM_E = ROP_S(IJKE,M)/M_PM
389 NU_PM_W = ROP_S(IJKW,M)/M_PM
390 NU_PM_N = ROP_S(IJKN,M)/M_PM
391 NU_PM_S = ROP_S(IJKS,M)/M_PM
392 NU_PM_T = ROP_S(IJKT,M)/M_PM
393 NU_PM_B = ROP_S(IJKB,M)/M_PM
394
395
396
397 = P_S_C(IJK,M) * ZMAX(TRD_S_C(IJK,M))
398 S10_rhs = P_S_C(IJK,M) * ZMAX(-TRD_S_C(IJK,M))
399
400
401
402
403 = 2.d0*Mu_s_c(IJK,M)*TRD_S2(IJK,M)
404
405
406
407 = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_s_c(IJK,M) )
408 S16_rhs = (TRD_S_C(IJK,M)**2)*ZMAX( LAMBDA_s_C(IJK,M) )
409
410 DO L = 1, MMAX
411 D_PL = D_P(IJK,L)
412 M_PL = (Pi/6.d0)*D_PL**3 * RO_S(IJK,L)
413 NU_PL_p = ROP_S(IJK,L)/M_PL
414 NU_PL_E = ROP_S(IJKE,L)/M_PL
415 NU_PL_W = ROP_S(IJKW,L)/M_PL
416 NU_PL_N = ROP_S(IJKN,L)/M_PL
417 NU_PL_S = ROP_S(IJKS,L)/M_PL
418 NU_PL_T = ROP_S(IJKT,L)/M_PL
419 NU_PL_B = ROP_S(IJKB,L)/M_PL
420
421 T_PL_p = Theta_m(IJK,L)
422 T_PL_E = Theta_m(IJKE,L)
423 T_PL_W = Theta_m(IJKW,L)
424 T_PL_N = Theta_m(IJKN,L)
425 T_PL_S = Theta_m(IJKS,L)
426 T_PL_T = Theta_m(IJKT,L)
427 T_PL_B = Theta_m(IJKB,L)
428
429 UsL_e = U_S(IJK,L)
430 UsL_w = U_S(IMJK,L)
431 VsL_n = V_S(IJK,L)
432 VsL_s = V_S(IJMK,L)
433 WsL_t = W_S(IJK,L)
434 WsL_b = W_S(IJKM,L)
435 UsL_p = AVG_X_E(U_S(IMJK,L),U_S(IJK,L),I)
436 VsL_p = AVG_Y_N(V_S(IJMK,L),V_S(IJK,L) )
437 WsL_p = AVG_Z_T(W_S(IJKM,L),W_S(IJK,L) )
438
439
440
441 = S20_sum_lhs + EDT_s_ip(IJK,M,L)
442
443
444
445 = S11_sum_lhs + ZMAX(-EDvel_sL_ip(IJK,M,L)* &
446 TRD_S_C(IJK,L) ) * VOL(IJK)
447 S11_sum_rhs = S11_sum_rhs + ZMAX( EDvel_sL_ip(IJK,M,L)* &
448 TRD_S_C(IJK,L) ) * VOL(IJK)
449
450
451
452 = S12_sum_lhs + ZMAX(-EDvel_sM_ip(IJK,M,L)*&
453 TRD_S_C(IJK,M) ) * VOL(IJK)
454 S12_sum_rhs = S12_sum_rhs + ZMAX( EDvel_sM_ip(IJK,M,L)*&
455 TRD_S_C(IJK,M) ) * VOL(IJK)
456
457 IF (M .NE. L) THEN
458 LM = FUNLM(L,M)
459
460
461
462 = S17_sum_lhs + 2.d0*MU_sL_ip(IJK,M,L)*&
463 ZMAX( - TRD_s2_ip(IJK,M,L) )
464 S17_sum_rhs = S17_sum_rhs + 2.d0*MU_sL_ip(IJK,M,L)*&
465 ZMAX( TRD_s2_ip(IJK,M,L) )
466
467
468
469
470
471
472
473
474
475
476 = S18_sum_lhs + ZMAX(-1.d0* &
477 (Xi_sL_ip(IJK,M,L)-(2.d0/3.d0)*Mu_sL_ip(IJK,M,L))*&
478 TRD_S_C(IJK,M)*TRD_S_C(IJK,L) )
479 S18_sum_rhs = S18_sum_rhs + ZMAX(&
480 (Xi_sL_ip(IJK,M,L)-(2.d0/3.d0)*Mu_sL_ip(IJK,M,L))*&
481 TRD_S_C(IJK,M)*TRD_S_C(IJK,L) )
482
483
484
485
486 = AVG_X_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKE,M,L),I)
487 Kth_sL_w = AVG_X_S(Kth_sL_ip(IJKW,M,L),Kth_sL_ip(IJK,M,L), IM)
488 Kth_sL_n = AVG_Y_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKN,M,L),J)
489 Kth_sL_s = AVG_Y_S(Kth_sL_ip(IJKS,M,L),Kth_sL_ip(IJK,M,L), JM)
490 Kth_sL_t = AVG_Z_S(Kth_sL_ip(IJK,M,L), Kth_sL_ip(IJKT,M,L),K)
491 Kth_sL_b = AVG_Z_S(Kth_sL_ip(IJKB,M,L),Kth_sL_ip(IJK,M,L), KM)
492
493 S21a_sum = S21a_sum + ( (Kth_sL_e*(T_PL_E-T_PL_p) )*&
494 ODX_E(I)*AYZ(IJK) - (Kth_sL_w*(T_PL_p-T_PL_W) )*ODX_E(IM)*&
495 AYZ(IMJK) )
496 S21b_sum = S21b_sum + ( (Kth_sL_n*(T_PL_N-T_PL_p) )*&
497 ODY_N(J)*AXZ(IJK) - (Kth_sL_s*(T_PL_p-T_PL_S) )*ODY_N(JM)*&
498 AXZ(IJMK) )
499 S21c_sum = S21c_sum + ( (Kth_sL_t*(T_PL_T-T_PL_p) )*&
500 ODZ_T(K)*OX(I)*AXY(IJK) - (Kth_sL_b*(T_PL_p-T_PL_B) )*&
501 ODZ_T(KM)*OX(I)*AXY(IJKM) )
502
503
504
505
506 = AVG_X_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKE,M,L),I)
507 Knu_sM_e = AVG_X_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKE,M,L),I)
508 Knu_sL_w = AVG_X_S(Knu_sL_ip(IJKW,M,L),Knu_sL_ip(IJK,M,L), IM)
509 Knu_sM_w = AVG_X_S(Knu_sM_ip(IJKW,M,L),Knu_sM_ip(IJK,M,L), IM)
510 Knu_sL_n = AVG_Y_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKN,M,L),J)
511 Knu_sM_n = AVG_Y_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKN,M,L),J)
512 Knu_sL_s = AVG_Y_S(Knu_sL_ip(IJKS,M,L),Knu_sL_ip(IJK,M,L), JM)
513 Knu_sM_s = AVG_Y_S(Knu_sM_ip(IJKS,M,L),Knu_sM_ip(IJK,M,L), JM)
514 Knu_sL_t = AVG_Z_S(Knu_sL_ip(IJK,M,L), Knu_sL_ip(IJKT,M,L),K)
515 Knu_sM_t = AVG_Z_S(Knu_sM_ip(IJK,M,L), Knu_sM_ip(IJKT,M,L),K)
516 Knu_sL_b = AVG_Z_S(Knu_sL_ip(IJKB,M,L),Knu_sL_ip(IJK,M,L), KM)
517 Knu_sM_b = AVG_Z_S(Knu_sM_ip(IJKB,M,L),Knu_sM_ip(IJK,M,L), KM)
518
519 S14a_sum = S14a_sum + ( (Knu_sL_e*(NU_PL_E-NU_PL_p) - &
520 Knu_sM_e*(NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK) - (Knu_sL_w*&
521 (NU_PL_p-NU_PL_W) - Knu_sM_w*(NU_PM_p-NU_PM_W) )*ODX_E(IM)*&
522 AYZ(IMJK) )
523 S14b_sum = S14b_sum + ( (Knu_sL_n*(NU_PL_N-NU_PL_p) - &
524 Knu_sM_n*(NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK) - (Knu_sL_s*&
525 (NU_PL_p-NU_PL_S) - Knu_sM_s*(NU_PM_p-NU_PM_S) )*ODY_N(JM)*&
526 AXZ(IJMK) )
527 S14c_sum = S14c_sum + ( (Knu_sL_t*(NU_PL_T-NU_PL_p) - &
528 Knu_sM_t*(NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK) - &
529 (Knu_sL_b*(NU_PL_p-NU_PL_B) - Knu_sM_b*(NU_PM_p-NU_PM_B) )*&
530 ODZ_T(KM)*OX(I)*AXY(IJKM) )
531
532
533
534
535 = AVG_X_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKE,M,L),I)
536 Kvel_s_w = AVG_X_H(Kvel_s_ip(IJKW,M,L),Kvel_s_ip(IJK,M,L), IM)
537 Kvel_s_n = AVG_Y_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKN,M,L),J)
538 Kvel_s_s = AVG_Y_H(Kvel_s_ip(IJKS,M,L),Kvel_s_ip(IJK,M,L), JM)
539 Kvel_s_t = AVG_Z_H(Kvel_s_ip(IJK,M,L), Kvel_s_ip(IJKT,M,L),K)
540 Kvel_s_b = AVG_Z_H(Kvel_s_ip(IJKB,M,L),Kvel_s_ip(IJK,M,L), KM)
541
542 S9_sum = S9_sum + ( Kvel_s_e*(UsM_e-UsL_e)*AYZ(IJK) - &
543 Kvel_s_w*(UsM_w-UsL_w)*AYZ(IMJK) + Kvel_s_n*(VsM_n-VsL_n)*AXZ(IJK)-&
544 Kvel_s_s*(VsM_s-VsL_s)*AXZ(IJMK) + Kvel_s_t*(WsM_t-WsL_t)*AXY(IJK)-&
545 Kvel_s_b*(WsM_b-WsL_b)*AXY(IJKM) )
546
547 ENDIF
548
549 ENDDO
550
551
552
553
554 = ( (S11_sum_lhs+S12_sum_lhs)+&
555 (S10_lhs+S16_lhs+S17_sum_lhs+&
556 S18_sum_lhs-S20_sum_lhs+S13_sum_lhs)*VOL(IJK) + &
557 ZMAX(S21a_sum+S21b_sum+S21c_sum)+ &
558 ZMAX(S14a_sum+S14b_sum+S14c_sum)+ ZMAX(S9_sum) ) / &
559 Theta_m(IJK,M)
560
561 SOURCERHS = ( S10_rhs+S15_rhs+S16_rhs+S17_sum_rhs+S18_sum_rhs+&
562 S13_sum_rhs) * VOL(IJK) + S11_sum_rhs+S12_sum_rhs+ &
563 ZMAX(- (S14a_sum+S14b_sum+S14c_sum) ) + ZMAX(-S9_sum) + &
564 ZMAX(- (S21a_sum+S21b_sum+S21c_sum) )
565
566
567 RETURN
568 END SUBROUTINE SOURCE_GRANULAR_ENERGY_IA
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588 SUBROUTINE SOURCE_GRANULAR_ENERGY_GD(SOURCELHS, &
589 SOURCERHS, IJK, M)
590
591
592
593
594 USE compar
595 USE constant
596 USE drag
597 USE fldvar
598 USE fun_avg
599 USE functions
600 USE geometry
601 USE indices
602 USE parallel
603 USE param
604 USE param1
605 USE physprop
606 USE rdf
607 USE residual
608 USE run
609 USE toleranc
610 USE trace
611 USE turb
612 USE visc_g
613 USE visc_s
614 use kintheory
615 IMPLICIT NONE
616
617
618
619
620 DOUBLE PRECISION, INTENT(INOUT) :: sourcerhs
621
622 DOUBLE PRECISION, INTENT(INOUT) :: sourcelhs
623
624 INTEGER, INTENT(IN) :: M
625
626 INTEGER, INTENT(IN) :: IJK
627
628
629
630
631 INTEGER :: I, J, K, IM, JM, KM, IMJK, IJMK, IJKM,&
632 IJKE, IJKW, IJKN, IJKS, IJKT, IJKB
633
634 DOUBLE PRECISION :: NU_PM_E, NU_PM_W, NU_PM_N, NU_PM_S, NU_PM_T,&
635 NU_PM_B, NU_PM_p
636
637 DOUBLE PRECISION :: M_PM, D_PM
638
639 DOUBLE PRECISION :: Knu_sM_e, Knu_sM_w, Knu_sM_n, Knu_sM_s,&
640 Knu_sM_t, Knu_sM_b
641
642 DOUBLE PRECISION :: S8_rhs, S9_rhs, S10_rhs, S11_rhs, S12_rhs,&
643 S13a_rhs, S13b_rhs, S14a_rhs, S14b_rhs, S15a_rhs,&
644 S15b_rhs
645
646 DOUBLE PRECISION :: S8_lhs, S10_lhs, S11_lhs, S12_lhs, S13a_lhs,&
647 S13b_lhs, S14a_lhs, S14b_lhs, S15a_lhs, S15b_lhs
648
649 DOUBLE PRECISION :: VSLIP, Kslip, Tslip_rhs, Tslip_lhs, Tvis_lhs
650 DOUBLE PRECISION :: ep_sM
651
652 DOUBLE PRECISION :: chi, Sgama_lhs, Spsi_rhs
653
654
655
656 DOUBLE PRECISION, EXTERNAL :: G_gtsh
657
658
659 = I_OF(IJK)
660 J = J_OF(IJK)
661 K = K_OF(IJK)
662 IM = Im1(I)
663 JM = Jm1(J)
664 KM = Km1(K)
665 IMJK = IM_OF(IJK)
666 IJMK = JM_OF(IJK)
667 IJKM = KM_OF(IJK)
668 IJKE = EAST_OF(IJK)
669 IJKW = WEST_OF(IJK)
670 IJKN = NORTH_OF(IJK)
671 IJKS = SOUTH_OF(IJK)
672 IJKT = TOP_OF(IJK)
673 IJKB = BOTTOM_OF(IJK)
674
675 S8_lhs = ZERO
676 S10_lhs = ZERO
677 S11_lhs = ZERO
678 S12_lhs = ZERO
679 S13a_lhs = ZERO
680 S13b_lhs = ZERO
681 S14a_lhs = ZERO
682 S14b_lhs = ZERO
683 S15a_lhs = ZERO
684 S15b_lhs = ZERO
685 S8_rhs = ZERO
686 S9_rhs = ZERO
687 S10_rhs = ZERO
688 S11_rhs = ZERO
689 S12_rhs = ZERO
690 S13a_rhs = ZERO
691 S13b_rhs = ZERO
692 S14a_rhs = ZERO
693 S14b_rhs = ZERO
694 S15a_rhs = ZERO
695 S15b_rhs = ZERO
696 Sgama_lhs = ZERO
697 = ZERO
698
699 = D_P(IJK,M)
700 M_PM = (Pi/6.d0)*D_PM**3 * RO_S(IJK,M)
701 EP_SM = EP_s(IJK,M)
702 CHI = G_0(IJK,M,M)
703
704 NU_PM_p = ROP_S(IJK,M)/M_PM
705 NU_PM_E = ROP_S(IJKE,M)/M_PM
706 NU_PM_W = ROP_S(IJKW,M)/M_PM
707 NU_PM_N = ROP_S(IJKN,M)/M_PM
708 NU_PM_S = ROP_S(IJKS,M)/M_PM
709 NU_PM_T = ROP_S(IJKT,M)/M_PM
710 NU_PM_B = ROP_S(IJKB,M)/M_PM
711
712
713
714 = P_S_C(IJK,M) * ZMAX(TRD_S_C(IJK,M))
715 S8_rhs = P_S_C(IJK,M) * ZMAX(-TRD_S_C(IJK,M))
716
717
718
719 = 2.d0*Mu_s_c(IJK,M)*TRD_S2(IJK,M)
720
721
722
723 = (TRD_S_C(IJK,M)**2)*ZMAX( -LAMBDA_s_C(IJK,M) )
724 S10_rhs = (TRD_S_C(IJK,M)**2)*ZMAX( LAMBDA_s_C(IJK,M) )
725
726
727
728 = (3.d0/2.d0)*EDT_s_ip(IJK,M,M)
729 S11_rhs = (1.d0/2.d0)*EDT_s_ip(IJK,M,M)*Theta_m(IJK,M)
730
731
732
733 = ZMAX( EDvel_sM_ip(IJK,M,M) * TRD_S_C(IJK,M) )
734 S12_rhs = ZMAX( -EDvel_sM_ip(IJK,M,M) * TRD_S_C(IJK,M) )*Theta_m(IJK,M)
735
736
737
738
739
740 IF(KT_TYPE_ENUM == GTSH_2012) THEN
741 S11_lhs = 1.5d0*ROP_s(IJK,M) * S11_lhs
742 S11_rhs = 1.5d0*ROP_s(IJK,M) * S11_rhs
743 S12_lhs = 1.5d0*ROP_s(IJK,M) * S12_lhs
744 S12_rhs = 1.5d0*ROP_s(IJK,M) * S12_rhs
745 Sgama_lhs = 3d0*NU_PM_p * G_gtsh(EP_SM, chi, IJK, M)
746 Spsi_rhs = 1.5d0*ROP_s(IJK,M) * xsi_gtsh(ijk)
747 ENDIF
748
749
750
751 = AVG_X_S(Kphi_s(IJK,M), Kphi_s(IJKE,M),I)
752 Knu_sM_w = AVG_X_S(Kphi_s(IJKW,M),Kphi_s(IJK,M), IM)
753 Knu_sM_n = AVG_Y_S(Kphi_s(IJK,M), Kphi_s(IJKN,M),J)
754 Knu_sM_s = AVG_Y_S(Kphi_s(IJKS,M),Kphi_s(IJK,M), JM)
755 Knu_sM_t = AVG_Z_S(Kphi_s(IJK,M), Kphi_s(IJKT,M),K)
756 Knu_sM_b = AVG_Z_S(Kphi_s(IJKB,M),Kphi_s(IJK,M), KM)
757
758 S13a_rhs = Knu_sM_e*ZMAX( (NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK)
759 S13a_lhs = Knu_sM_e*ZMAX( -(NU_PM_E-NU_PM_p) )*ODX_E(I)*AYZ(IJK)
760
761 S13b_rhs = Knu_sM_w*ZMAX( -(NU_PM_p-NU_PM_W) )*ODX_E(IM)*AYZ(IMJK)
762 S13b_lhs = Knu_sM_w*ZMAX( (NU_PM_p-NU_PM_W) )*ODX_E(IM)*AYZ(IMJK)
763
764 S14a_rhs = Knu_sM_n*ZMAX( (NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK)
765 S14a_lhs = Knu_sM_n*ZMAX( -(NU_PM_N-NU_PM_p) )*ODY_N(J)*AXZ(IJK)
766
767 S14b_rhs = Knu_sM_s*ZMAX( -(NU_PM_p-NU_PM_S) )*ODY_N(JM)*AXZ(IJMK)
768 S14b_lhs = Knu_sM_s*ZMAX( (NU_PM_p-NU_PM_S) )*ODY_N(JM)*AXZ(IJMK)
769
770 S15a_rhs = Knu_sM_t*ZMAX( (NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK)
771 S15a_lhs = Knu_sM_t*ZMAX( -(NU_PM_T-NU_PM_p) )*ODZ_T(K)*OX(I)*AXY(IJK)
772
773 S15b_rhs = Knu_sM_b*ZMAX( -(NU_PM_p-NU_PM_B) )*ODZ_T(KM)*OX(I)*AXY(IJKM)
774 S15b_lhs = Knu_sM_b*ZMAX( (NU_PM_p-NU_PM_B) )*ODZ_T(KM)*OX(I)*AXY(IJKM)
775
776 SOURCELHS = ( (S8_lhs+S10_lhs)*VOL(IJK) + &
777 S13a_lhs+S13b_lhs+S14a_lhs+S14b_lhs+S15a_lhs+S15b_lhs)/Theta_m(IJK,M)&
778 + (S11_lhs + S12_lhs + Sgama_lhs)*VOL(IJK)
779
780
781 SOURCERHS = ( S8_rhs+S9_rhs+S10_rhs+S11_rhs+S12_rhs+Spsi_rhs) * VOL(IJK) + &
782 S13a_rhs+S13b_rhs+S14a_rhs+S14b_rhs+S15a_rhs+S15b_rhs
783
784
785
786 IF(SWITCH > ZERO .AND. RO_g0 /= ZERO .AND. &
787 (KT_TYPE_ENUM == GD_1999)) THEN
788
789 VSLIP = (U_S(IJK,M)-U_G(IJK))**2 + (V_S(IJK,M)-V_G(IJK))**2 +&
790 (W_S(IJK,M)-W_G(IJK))**2
791 VSLIP = DSQRT(VSLIP)
792
793
794 = SWITCH*81.d0*EP_SM*(MU_G(IJK)*VSLIP)**2.d0 / &
795 (chi*D_P(IJK,M)**3.D0*RO_S(IJK,M)*DSQRT(PI))
796
797 Tslip_rhs = 1.5d0*Kslip/( (THETA_M(IJK,M)+SMALL_NUMBER)**0.50)*VOL(IJK)
798 Tslip_lhs = 0.5d0*Kslip/( (THETA_M(IJK,M)+SMALL_NUMBER)**1.50)*VOL(IJK)
799
800
801 = SWITCH*3d0*F_GS(IJK,M)*VOL(IJK)
802
803 SOURCELHS = SOURCELHS + Tslip_lhs + Tvis_lhs
804 SOURCERHS = SOURCERHS + Tslip_rhs
805 ENDIF
806
807 RETURN
808 END SUBROUTINE SOURCE_GRANULAR_ENERGY_GD
809
810
811