File: /nfs/home/0/users/jenkins/mfix.git/model/bc_theta.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 SUBROUTINE BC_THETA(M, A_m, B_m)
31
32
33
34
35 USE param
36 USE param1
37 USE parallel
38 USE matrix
39 USE scales
40 USE constant
41 USE toleranc
42 USE run
43 USE physprop
44 USE fldvar
45 USE visc_s
46 USE geometry
47 USE output
48 USE indices
49 USE bc
50 USE compar
51 USE mpi_utility
52 USE fun_avg
53 USE functions
54 IMPLICIT NONE
55
56
57
58
59 INTEGER :: M
60
61 DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
62
63 DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
64
65
66
67
68 INTEGER :: L
69
70 INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK,&
71 IM, JM, KM
72
73 DOUBLE PRECISION :: Gw, Hw, Cw
74
75
76
77 DO L = 1, DIMENSION_BC
78 IF( BC_DEFINED(L) ) THEN
79 IF(BC_TYPE(L) .EQ. 'NO_SLIP_WALL' .OR.&
80 BC_TYPE(L) .EQ. 'FREE_SLIP_WALL' .OR.&
81 BC_TYPE(L) .EQ. 'PAR_SLIP_WALL' ) THEN
82 I1 = BC_I_w(L)
83 I2 = BC_I_e(L)
84 J1 = BC_J_s(L)
85 J2 = BC_J_n(L)
86 K1 = BC_K_b(L)
87 K2 = BC_K_t(L)
88
89 IF(BC_JJ_PS(L).GT.0) THEN
90 DO K = K1, K2
91 DO J = J1, J2
92 DO I = I1, I2
93
94 IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
95 IF (DEAD_CELL_AT(I,J,K)) CYCLE
96 = FUNIJK(I, J, K)
97 IF (FLOW_AT(IJK)) CYCLE
98 = Im1(I)
99 JM = Jm1(J)
100 KM = Km1(K)
101
102
103
104 (IJK, e, M) = ZERO
105 A_m(IJK, w, M) = ZERO
106 A_m(IJK, n, M) = ZERO
107 A_m(IJK, s, M) = ZERO
108 A_m(IJK, t, M) = ZERO
109 A_m(IJK, b, M) = ZERO
110 A_m(IJK, 0, M) = -ONE
111 b_m(IJK, M) = ZERO
112
113
114 IF(FLUID_AT(EAST_OF(IJK)))THEN
115 IF(EP_s(EAST_OF(IJK),M).LE.DIL_EP_s)THEN
116 A_m(IJK, e, M) = ONE
117 ELSE
118 IF (BC_JJ_PS(L).EQ.3) THEN
119
120
121 = 1d0
122 Hw = 0d0
123 Cw = 0d0
124 ELSE
125
126
127 CALL CALC_THETA_BC(IJK,EAST_OF(IJK),'E',&
128 M,L,gw,hw,cw)
129
130 IF (BC_JJ_PS(L).EQ.2) cw=0d0
131 ENDIF
132
133 A_m(IJK, e, M) = -(HALF*hw - oDX_E(I)*gw)
134 A_m(IJK, 0, M) = -(HALF*hw + oDX_E(I)*gw)
135
136 IF (BC_JJ_PS(L) .EQ. 1) THEN
137 b_m(IJK, M) = -cw
138 ELSE
139 b_m(IJK,M) = ZERO
140 ENDIF
141 ENDIF
142
143 ELSEIF(FLUID_AT(WEST_OF(IJK)))THEN
144 IF(EP_s(WEST_OF(IJK),M).LE.DIL_EP_s)THEN
145 A_m(IJK, w, M) = ONE
146 ELSE
147 IF (BC_JJ_PS(L).EQ.3) THEN
148 Gw = 1d0
149 Hw = 0d0
150 Cw = 0d0
151 ELSE
152 CALL CALC_THETA_BC(IJK,WEST_OF(IJK),'W',&
153 M,L,gw,hw,cw)
154 IF (BC_JJ_PS(L).EQ.2) cw=0d0
155 ENDIF
156
157 A_m(IJK, w, M) = -(HALF*hw - oDX_E(IM)*gw)
158 A_m(IJK, 0, M) = -(HALF*hw + oDX_E(IM)*gw)
159
160 IF (BC_JJ_PS(L) .EQ. 1) THEN
161 b_m(IJK, M) = -cw
162 ELSE
163 b_m(IJK,M) = ZERO
164 ENDIF
165 ENDIF
166
167 ELSEIF(FLUID_AT(NORTH_OF(IJK)))THEN
168 IF(EP_s(NORTH_OF(IJK),M).LE.DIL_EP_s)THEN
169 A_m(IJK, n, M) = ONE
170 ELSE
171 IF (BC_JJ_PS(L).EQ.3) THEN
172 Gw = 1d0
173 Hw = 0d0
174 Cw = 0d0
175 ELSE
176 CALL CALC_THETA_BC(IJK,NORTH_OF(IJK),'N',&
177 M,L,gw,hw,cw)
178 IF (BC_JJ_PS(L).EQ.2) cw=0d0
179 ENDIF
180
181 A_m(IJK, n, M) = -(HALF*hw - oDY_N(J)*gw)
182 A_m(IJK, 0, M) = -(HALF*hw + oDY_N(J)*gw)
183
184 IF (BC_JJ_PS(L) .EQ. 1) THEN
185 b_m(IJK, M) = -cw
186 ELSE
187 b_m(IJK,M) = ZERO
188 ENDIF
189 ENDIF
190
191 ELSEIF(FLUID_AT(SOUTH_OF(IJK)))THEN
192 IF(EP_s(SOUTH_OF(IJK),M).LE.DIL_EP_s)THEN
193 A_m(IJK, s, M) = ONE
194 ELSE
195 IF (BC_JJ_PS(L).EQ.3) THEN
196 Gw = 1d0
197 Hw = 0d0
198 Cw = 0d0
199 ELSE
200 CALL CALC_THETA_BC(IJK,SOUTH_OF(IJK),'S',&
201 M,L,gw,hw,cw)
202 IF (BC_JJ_PS(L).EQ.2) cw=0d0
203 ENDIF
204
205 A_m(IJK, s, M) = -(HALF*hw - oDY_N(JM)*gw)
206 A_m(IJK, 0, M) = -(HALF*hw + oDY_N(JM)*gw)
207
208 IF (BC_JJ_PS(L) .EQ. 1) THEN
209 b_m(IJK, M) = -cw
210 ELSE
211 b_m(IJK,M) = ZERO
212 ENDIF
213 ENDIF
214
215 ELSEIF(FLUID_AT(TOP_OF(IJK)))THEN
216 IF(EP_s(TOP_OF(IJK),M).LE.DIL_EP_s)THEN
217 A_m(IJK, t, M) = ONE
218 ELSE
219 IF (BC_JJ_PS(L).EQ.3) THEN
220 Gw = 1d0
221 Hw = 0d0
222 Cw = 0d0
223 ELSE
224 CALL CALC_THETA_BC(IJK,TOP_OF(IJK),'T',&
225 M,L,gw,hw,cw)
226 IF (BC_JJ_PS(L).EQ.2) cw=0d0
227 ENDIF
228
229 A_m(IJK, t, M) = -(HALF*hw - oX(I)*oDZ_T(K)*gw)
230 A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(K)*gw)
231
232 IF (BC_JJ_PS(L) .EQ. 1) THEN
233 b_m(IJK, M) = -cw
234 ELSE
235 b_m(IJK,M) = ZERO
236 ENDIF
237 ENDIF
238
239 ELSEIF(FLUID_AT(BOTTOM_OF(IJK)))THEN
240 IF(EP_s(BOTTOM_OF(IJK),M).LE.DIL_EP_s)THEN
241 A_m(IJK, b , M) = ONE
242 ELSE
243 IF (BC_JJ_PS(L).EQ.3) THEN
244 Gw = 1d0
245 Hw = 0d0
246 Cw = 0d0
247 ELSE
248 CALL CALC_THETA_BC(IJK,BOTTOM_OF(IJK),'B',&
249 M,L,gw,hw,cw)
250 IF (BC_JJ_PS(L).EQ.2) cw=0d0
251 ENDIF
252
253 A_m(IJK, b, M) = -(HALF*hw - oX(I)*oDZ_T(KM)*gw)
254 A_m(IJK, 0, M) = -(HALF*hw + oX(I)*oDZ_T(KM)*gw)
255
256 IF (BC_JJ_PS(L) .EQ. 1) THEN
257 b_m(IJK, M) = -cw
258 ELSE
259 b_m(IJK,M) = ZERO
260 ENDIF
261 ENDIF
262 ENDIF
263
264 ENDDO
265 ENDDO
266 ENDDO
267
268 ENDIF
269 ENDIF
270 ENDIF
271 ENDDO
272
273 RETURN
274 END SUBROUTINE BC_THETA
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292 SUBROUTINE CALC_THETA_BC(IJK1,IJK2,FCELL,M,L,Gw,Hw,Cw)
293
294
295
296
297 USE bc
298 USE compar
299 USE constant
300 USE fldvar
301 USE fun_avg
302 USE functions
303 USE geometry
304 USE indices
305 USE mpi_utility
306 USE param
307 USE param1
308 USE physprop
309 USE rdf
310 USE run
311 USE rxns
312 USE toleranc
313 USE turb
314 USE visc_s
315 IMPLICIT NONE
316
317
318
319
320 INTEGER, INTENT(IN) :: IJK1, IJK2
321
322 CHARACTER, INTENT(IN) :: FCELL
323
324 INTEGER, INTENT(IN) :: M
325
326 INTEGER, INTENT(IN) :: L
327
328 DOUBLE PRECISION, INTENT(INOUT) :: Gw, Hw, Cw
329
330
331
332
333 INTEGER :: IJK3
334
335 INTEGER :: MM
336
337 DOUBLE PRECISION :: EPg_avg, Mu_g_avg, RO_g_avg, smallTheta
338
339 DOUBLE PRECISION :: ep_star_avg
340
341 DOUBLE PRECISION :: EPs_avg(DIMENSION_M), DP_avg(DIMENSION_M), &
342 TH_avg(DIMENSION_M)
343
344 DOUBLE PRECISION ROs_avg(DIMENSION_M)
345
346 DOUBLE PRECISION :: K_12_avg, Tau_12_avg
347
348 DOUBLE PRECISION :: USCM, VSCM, WSCM
349
350 DOUBLE PRECISION :: UGC, VGC, WGC
351
352 DOUBLE PRECISION :: VREL
353
354 DOUBLE PRECISION :: VSLIPSQ
355
356 DOUBLE PRECISION :: VWDOTN (DIMENSION_M)
357
358
359 DOUBLE PRECISION :: GNUWDOTN (DIMENSION_M)
360
361
362 DOUBLE PRECISION :: GTWDOTN (DIMENSION_M)
363
364 CHARACTER(LEN=80) LINE(1)
365
366 DOUBLE PRECISION :: g0(DIMENSION_M)
367
368 DOUBLE PRECISION :: g0EPs_avg
369
370
371
372
373 DOUBLE PRECISION :: PHIP_JJ
374
375
376
377
378
379 = (to_SI)**4 * ZERO_EP_S
380
381
382 = EP_g(IJK2)
383 ep_star_avg = EP_star_array(IJK2)
384 Mu_g_avg = Mu_g(IJK2)
385 RO_g_avg = RO_g(IJK2)
386 g0EPs_avg = ZERO
387
388
389 IF(KT_TYPE_ENUM == SIMONIN_1996 .OR. &
390 KT_TYPE_ENUM == AHMADI_1995) THEN
391 K_12_avg = K_12(IJK2)
392 Tau_12_avg = Tau_12(IJK2)
393 ELSE
394 K_12_avg = ZERO
395 Tau_12_avg = ZERO
396 ENDIF
397
398 DO MM = 1, SMAX
399 g0(MM) = G_0(IJK2, M, MM)
400 EPs_avg(MM) = EP_s(IJK2,MM)
401 DP_avg(MM) = D_P(IJK2,MM)
402 g0EPs_avg = g0EPs_avg + G_0(IJK2, M, MM)*EP_s(IJK2,MM)
403 ROs_avg(MM) = RO_S(IJK2,MM)
404 ENDDO
405
406
407 IF(FCELL .EQ. 'N')THEN
408 DO MM = 1, SMAX
409 TH_avg(MM) = AVG_Y(Theta_m(IJK1,MM),Theta_m(IJK2,MM),J_OF(IJK1))
410 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
411
412
413
414
415
416 (MM) = -1.d0*V_S(IJK1,MM)
417
418
419
420
421
422 = NORTH_OF(IJK2)
423 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
424 (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDY_N(J_OF(IJK2))
425
426
427
428 (MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
429 oDY_N(J_OF(IJK2))
430 ENDDO
431
432
433 = AVG_Y(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
434 AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
435 J_OF(IJK1))
436 VGC = V_g(IJK1)
437 WGC = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
438 AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
439 J_OF(IJK1))
440 USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
441 AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
442 J_OF(IJK1))
443 VSCM = V_s(IJK1,M)
444 WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
445 AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
446 J_OF(IJK1))
447
448
449 =(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
450
451
452 ELSEIF(FCELL .EQ. 'S')THEN
453 DO MM = 1, SMAX
454 TH_avg(MM) = AVG_Y(Theta_m(IJK2, MM),Theta_m(IJK1, MM),J_OF(IJK2))
455 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
456
457
458
459
460
461 (MM) = 1.d0*V_S(IJK2,MM)
462
463
464
465
466
467 = SOUTH_OF(IJK2)
468 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
469 (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDY_N(J_OF(IJK3))
470
471
472
473 (MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
474 oDY_N(J_OF(IJK3))
475 ENDDO
476
477
478 = AVG_Y(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
479 AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
480 J_OF(IJK2))
481 VGC = V_g(IJK2)
482 WGC = AVG_Y(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
483 AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
484 J_OF(IJK2))
485 USCM = AVG_Y(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
486 AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
487 J_OF(IJK2))
488 VSCM = V_s(IJK2,M)
489 WSCM = AVG_Y(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
490 AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
491 J_OF(IJK2))
492
493
494 =(WSCM-BC_Ww_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
495
496
497 ELSEIF(FCELL== 'E')THEN
498 DO MM = 1, SMAX
499 TH_avg(MM) = AVG_X(Theta_m(IJK1,MM),Theta_m(IJK2,MM),I_OF(IJK1))
500 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
501
502
503
504
505
506 (MM) = -1.d0*U_S(IJK1,MM)
507
508
509
510
511
512 = EAST_OF(IJK2)
513 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
514 (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oDX_E(I_OF(IJK2))
515
516
517
518 (MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
519 oDX_E(I_OF(IJK2))
520 ENDDO
521
522
523 = U_g(IJK1)
524 VGC = AVG_X(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
525 AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
526 I_OF(IJK1))
527 WGC = AVG_X(AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
528 AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
529 I_OF(IJK1))
530 USCM = U_s(IJK1,M)
531 VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
532 AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
533 I_OF(IJK1))
534 WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
535 AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
536 I_OF(IJK1))
537
538
539 =(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
540
541
542 ELSEIF(FCELL== 'W')THEN
543 DO MM = 1, SMAX
544 TH_avg(MM) = AVG_X(Theta_m(IJK2,MM),Theta_m(IJK1,MM),I_OF(IJK2))
545 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
546
547
548
549
550
551 (MM) = 1.d0*U_S(IJK2,MM)
552
553
554
555
556
557 = WEST_OF(IJK2)
558 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
559 (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oDX_E(I_OF(IJK3))
560
561
562
563 (MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
564 oDX_E(I_OF(IJK3))
565 ENDDO
566
567
568 = U_g(IJK2)
569 VGC = AVG_X(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
570 AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
571 I_OF(IJK2))
572 WGC = AVG_X(AVG_Z_T(W_g(KM_OF(IJK2)), W_g(IJK2)),&
573 AVG_Z_T(W_g(KM_OF(IJK1)), W_g(IJK1)),&
574 I_OF(IJK2))
575 USCM = U_s(IJK2,M)
576 VSCM = AVG_X(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
577 AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
578 I_OF(IJK2))
579 WSCM = AVG_X(AVG_Z_T(W_s(KM_OF(IJK2),M), W_s(IJK2,M)),&
580 AVG_Z_T(W_s(KM_OF(IJK1),M), W_s(IJK1,M)),&
581 I_OF(IJK2))
582
583
584 =(WSCM-BC_Ww_s(L, M))**2 + (VSCM-BC_Vw_s(L, M))**2
585
586
587 ELSEIF(FCELL== 'T')THEN
588 DO MM = 1, SMAX
589 TH_avg(MM) = AVG_Z(Theta_m(IJK1,MM),Theta_m(IJK2,MM),K_OF(IJK1))
590 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
591
592
593
594
595
596 (MM) = -1.d0*W_S(IJK1,MM)
597
598
599
600
601
602 = TOP_OF(IJK2)
603 GNUWDOTN(MM) = -1.d0*(6.d0/(PI*DP_avg(MM)))*&
604 (EP_s(IJK3,MM)-EP_s(IJK2,MM))*oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
605
606
607
608 (MM) = -1.d0*(Theta_m(IJK3,MM)-Theta_m(IJK2,MM))*&
609 oX(I_of(IJK2))*oDZ_T(K_OF(IJK2))
610 ENDDO
611
612
613 = AVG_Z(AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
614 AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
615 K_OF(IJK1))
616 VGC = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
617 AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
618 K_OF(IJK1))
619 WGC = W_g(IJK1)
620 USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
621 AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
622 K_OF(IJK1))
623 VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
624 AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
625 K_OF(IJK1))
626 WSCM = W_s(IJK1,M)
627
628
629 =(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
630
631
632 ELSEIF(FCELL== 'B')THEN
633 DO MM = 1, SMAX
634 TH_avg(MM) = AVG_Z(Theta_m(IJK2,MM),Theta_m(IJK1,MM),K_OF(IJK2))
635 IF(TH_avg(MM) < ZERO) TH_avg(MM) = smallTheta
636
637
638
639
640
641 (MM) = 1.d0*W_S(IJK2,MM)
642
643
644
645
646
647 = BOTTOM_OF(IJK2)
648 GNUWDOTN(MM) = 1.d0*(6.d0/(PI*DP_avg(MM)))*&
649 (EP_s(IJK2,MM)-EP_s(IJK3,MM))*oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
650
651
652
653 (MM) = 1.d0*(Theta_m(IJK2,MM)-Theta_m(IJK3,MM))*&
654 oX(I_of(IJK3))*oDZ_T(K_OF(IJK3))
655 ENDDO
656
657
658 = AVG_Z(AVG_X_E(U_g(IM_OF(IJK2)),U_g(IJK2),I_OF(IJK2)),&
659 AVG_X_E(U_g(IM_OF(IJK1)),U_g(IJK1),I_OF(IJK1)),&
660 K_OF(IJK2))
661 VGC = AVG_Z(AVG_Y_N(V_g(JM_OF(IJK2)), V_g(IJK2)),&
662 AVG_Y_N(V_g(JM_OF(IJK1)), V_g(IJK1)),&
663 K_OF(IJK2))
664 WGC = W_g(IJK2)
665 USCM = AVG_Z(AVG_X_E(U_s(IM_OF(IJK2),M),U_s(IJK2,M),I_OF(IJK2)),&
666 AVG_X_E(U_s(IM_OF(IJK1),M),U_s(IJK1,M),I_OF(IJK1)),&
667 K_OF(IJK2))
668 VSCM = AVG_Z(AVG_Y_N(V_s(JM_OF(IJK2),M), V_s(IJK2,M)),&
669 AVG_Y_N(V_s(JM_OF(IJK1),M), V_s(IJK1,M)),&
670 K_OF(IJK2))
671 WSCM = W_s(IJK2,M)
672
673
674 =(VSCM-BC_Vw_s(L, M))**2 + (USCM-BC_Uw_s(L, M))**2
675
676 ELSE
677 WRITE(LINE,'(A, A)') 'Error: Unknown FCELL'
678 CALL WRITE_ERROR('CALC_THETA_BC', LINE, 1)
679 call exitMPI(myPE)
680 ENDIF
681
682
683 = DSQRT( (UGC-USCM)**2 + (VGC-VSCM)**2 + (WGC-WSCM)**2 )
684
685 CALL THETA_Hw_Cw(g0, EPs_avg, EPg_avg, ep_star_avg, &
686 g0EPs_avg, TH_avg,Mu_g_avg,RO_g_avg, ROs_avg, &
687 DP_avg, K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
688 GNUWDOTN,GTWDOTN,M,Gw,Hw,Cw,L)
689
690
691
692 IF(BC_JJ_M .and. PHIP_OUT_JJ) THEN
693 IF(PHIP_OUT_ITER.eq.1) THEN
694 ReactionRates(IJK1,1)= PHIP_JJ(dsqrt(VSLIPSQ),TH_avg(1))
695 ENDIF
696 ENDIF
697
698 RETURN
699 END SUBROUTINE CALC_THETA_BC
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743 SUBROUTINE THETA_HW_CW(g0,EPS,EPG, ep_star_avg, &
744 g0EPs_avg,TH,Mu_g_avg,RO_g_avg, ROs_avg, DP_avg, &
745 K_12_avg,Tau_12_avg,VREL,VSLIPSQ,VWDOTN,&
746 GNUWDOTN,GTWDOTN,M,GW,HW,CW,L)
747
748
749
750
751 USE param
752 USE param1
753 USE physprop
754 USE run
755 USE constant
756 USE fldvar
757 USE toleranc
758 USE bc
759 USE compar
760 IMPLICIT NONE
761
762
763
764
765
766 DOUBLE PRECISION, INTENT(IN) :: g0(DIMENSION_M)
767
768 DOUBLE PRECISION, INTENT(IN) :: EPS(DIMENSION_M)
769
770 DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
771
772 DOUBLE PRECISION, INTENT(INOUT) :: g0EPs_avg
773
774 DOUBLE PRECISION, INTENT(INOUT) :: TH (DIMENSION_M)
775
776 DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
777
778 DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
779
780 DOUBLE PRECISION, INTENT(IN) :: ROS_avg(DIMENSION_M)
781
782 DOUBLE PRECISION, INTENT(IN) :: DP_avg(DIMENSION_M)
783
784 DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg
785
786 DOUBLE PRECISION, INTENT(IN) :: VREL
787
788 DOUBLE PRECISION, INTENT(IN) :: VSLIPSQ
789
790 DOUBLE PRECISION, INTENT(IN) :: VWDOTN (DIMENSION_M)
791
792
793 DOUBLE PRECISION, INTENT(IN) :: GNUWDOTN (DIMENSION_M)
794
795
796 DOUBLE PRECISION, INTENT(IN) :: GTWDOTN (DIMENSION_M)
797
798 INTEGER, INTENT(IN) :: M
799
800
801
802
803 DOUBLE PRECISION, INTENT(INOUT) :: GW, HW, CW
804
805 INTEGER, INTENT(IN) :: L
806
807
808
809
810 INTEGER :: LL
811
812 DOUBLE PRECISION :: K_1
813
814 DOUBLE PRECISION :: Kgran
815
816 DOUBLE PRECISION :: Kgran_star
817
818 DOUBLE PRECISION :: Beta, DgA
819
820 DOUBLE PRECISION :: Re_g
821
822 DOUBLE PRECISION :: C_d
823
824 DOUBLE PRECISION :: D_PM, M_PM, NU_PM
825
826 DOUBLE PRECISION :: D_PL, M_PL, NU_PL, MPSUM, DPSUMo2
827 DOUBLE PRECISION :: Ap_lm, Dp_lm, R0p_lm, R1p_lm, R8p_lm, R9p_lm, &
828 Bp_lm, R5p_lm, R6p_lm, R7p_lm
829 DOUBLE PRECISION :: K_s_sum, K_s_MM, K_sM_LM
830 DOUBLE PRECISION :: Kvel_s_sum, Kth_s_sum
831 DOUBLE PRECISION :: Knu_s_sum, K_common_term
832
833 DOUBLE PRECISION :: c_star, zeta0_star, press_star, eta0, &
834 kappa0, nu_kappa_star, kappa_k_star, &
835 kappa_star
836
837 DOUBLE PRECISION :: Re_T, Chi, vfrac, mu2_0, mu4_0, mu4_1, &
838 A2_gtshW, zeta_star, nu0, NuK, Kth0, KthK, EDT_s
839
840 DOUBLE PRECISION :: Zeta_c, Omega_c, Tau_2_c, Kappa_kin, &
841 Kappa_Col, Tau_12_st
842
843 DOUBLE PRECISION :: Ps, PsoTheta
844
845
846
847
848 DOUBLE PRECISION :: PHIP_JJ
849 DOUBLE PRECISION :: S_star
850
851
852 IF(TH(M) .LE. ZERO)THEN
853 TH(M) = 1D-8
854 IF (myPE.eq.PE_IO) THEN
855 WRITE(*,*) &
856 'Warning: Negative granular temp at wall set to 1e-8'
857
858 ENDIF
859 ENDIF
860
861
862 = DP_avg(M)
863 M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
864 NU_PM = (EPS(M)*ROS_avg(M))/M_PM
865
866
867 = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
868 IF (Re_g .lt. 1000.d0) THEN
869 C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
870 ELSE
871 C_d = 0.44d0
872 ENDIF
873 DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(DP_avg(M)*EPG**(2.65d0))
874 IF(VREL == ZERO) DgA = LARGE_NUMBER
875 Beta = EPS(M)*DgA
876
877
878
879
880
881 SELECT CASE (KT_TYPE_ENUM)
882
883 CASE (LUN_1984)
884 Kgran = 75d0*ROs_avg(M)*DP_avg(M)*DSQRT(Pi*TH(M))/&
885 (48*Eta*(41d0-33d0*Eta))
886
887 IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
888 Kgran_star = Kgran
889 ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
890 Kgran_star = ZERO
891 ELSE
892 Kgran_star = ROs_avg(M)*EPS(M)* g0(M)*TH(M)* Kgran/ &
893 (ROs_avg(M)*g0EPs_avg*TH(M) + &
894 1.2d0*DgA/ROs_avg(M)* Kgran)
895 ENDIF
896
897
898 = Kgran_star/g0(M)*( &
899 ( ONE + (12d0/5.d0)*Eta*g0EPs_avg ) * &
900 ( ONE + (12d0/5.d0)*Eta*Eta*(4d0*Eta-3d0) *g0EPs_avg ) + &
901 (64d0/(25d0*Pi)) * (41d0-33d0*Eta) *(Eta*g0EPs_avg)**2 )
902
903
904 = ROs_avg(M)*EPS(M)*(ONE+4.D0*Eta*g0EPs_avg)
905 Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
906
907
908 CASE (SIMONIN_1996)
909
910 = ROs_avg(M)/(DgA+small_number)
911 Zeta_c = (ONE+ C_e)*(49.d0-33.d0*C_e)/100.d0
912 Omega_c = 3.d0*(ONE+ C_e)**2 *(2.d0*C_e-ONE)/5.d0
913 Tau_2_c = DP_avg(M)/(6.d0*EPS(M)*g0(M) &
914 *DSQRT(16.d0*(TH(M)+Small_number)/PI))
915
916
917 = (9.d0/10.d0*K_12_avg*(Tau_12_avg/Tau_12_st) &
918 + 3.0D0/2.0D0 * TH(M)*(ONE+ Omega_c*EPS(M)*g0(M)))/ &
919 (9.d0/(5.d0*Tau_12_st) + zeta_c/Tau_2_c)
920
921 Kappa_Col = 18.d0/5.d0*EPS(M)*g0(M)*Eta* (Kappa_kin+ &
922 5.d0/9.d0*DP_avg(M)*DSQRT(TH(M)/PI))
923
924
925 = EPS(M)*ROs_avg(M)*(Kappa_kin + Kappa_Col)
926
927
928 = ROs_avg(M)*EPS(M)*(ONE+4.D0*Eta*g0EPs_avg)
929 Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+4.D0*Eta*g0EPs_avg)
930
931
932 CASE (AHMADI_1995)
933
934 = 0.1306D0*ROs_avg(M)*DP_avg(M)*(ONE+C_e**2)* ( &
935 ONE/g0(M)+4.8D0*EPS(M)+12.1184D0 *EPS(M)*EPS(M)*g0(M) )*&
936 DSQRT(TH(M))
937
938
939 = ROs_avg(M)*EPS(M)* &
940 ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
941 Ps = ROs_avg(M)*EPS(M)*TH(M)*&
942 ((ONE + 4.0D0*g0EPs_avg)+HALF*(ONE -C_e*C_e))
943
944
945 CASE (IA_2005)
946
947 IF(.NOT. SWITCH_IA) g0EPs_avg = EPS(M)*ROS_avg(M)
948
949 K_s_sum = ZERO
950 Kvel_s_sum = ZERO
951 Knu_s_sum = ZERO
952 Kth_s_sum = ZERO
953
954 Kgran = (75.d0/384.d0)*ROs_avg(M)*D_PM*DSQRT(Pi*TH(M)/M_PM)
955
956 IF(SWITCH == ZERO .OR. RO_g_avg == ZERO)THEN
957 Kgran_star = Kgran
958 ELSEIF(TH(M)/M_PM .LT. SMALL_NUMBER)THEN
959 Kgran_star = ZERO
960 ELSE
961 Kgran_star = Kgran*g0(M)*EPS(M)/ &
962 (g0EPs_avg+ 1.2d0*DgA*Kgran / (ROS_avg(M)**2 *(TH(M)/M_PM)))
963 ENDIF
964
965 K_s_MM = (Kgran_star/(M_PM*g0(M)))*&
966 (1.d0+(3.d0/5.d0)*(1.d0+C_E)*(1.d0+C_E)*g0EPs_avg)**2
967
968 DO LL = 1, SMAX
969 D_PL = DP_avg(LL)
970 M_PL = (PI/6.d0)*(D_PL**3.)*ROS_avg(LL)
971 MPSUM = M_PM + M_PL
972 DPSUMo2 = (D_PM+D_PL)/2.d0
973 NU_PL = (EPS(LL)*ROS_avg(LL))/M_PL
974
975 IF ( LL .eq. M) THEN
976 K_s_sum = K_s_sum + K_s_MM
977
978
979
980 ELSE
981 Ap_lm = (M_PM*TH(LL)+M_PL*TH(M))/(2.d0)
982 Bp_lm = (M_PM*M_PL*(TH(LL)-TH(M) ))/(2.d0*MPSUM)
983 Dp_lm = (M_PL*M_PM*(M_PM*TH(M)+M_PL*TH(LL) ))/&
984 (2.d0*MPSUM*MPSUM)
985 R0p_lm = ( 1.d0/( Ap_lm**1.5 * Dp_lm**2.5 ) )+ &
986 ( (15.d0*Bp_lm*Bp_lm)/( 2.d0* Ap_lm**2.5 * Dp_lm**3.5 ) )+&
987 ( (175.d0*(Bp_lm**4))/( 8.d0*Ap_lm**3.5 * Dp_lm**4.5 ) )
988 R1p_lm = ( 1.d0/( (Ap_lm**1.5)*(Dp_lm**3) ) )+ &
989 ( (9.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**4 ) )+&
990 ( (30.d0*Bp_lm**4) /( 2.d0*Ap_lm**3.5 * Dp_lm**5 ) )
991 R5p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**3 ) )+ &
992 ( (5.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**4 ) )+&
993 ( (14.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5 ) )
994 R6p_lm = ( 1.d0/( Ap_lm**3.5 * Dp_lm**3 ) )+ &
995 ( (7.d0*Bp_lm*Bp_lm)/( Ap_lm**4.5 * Dp_lm**4 ) )+&
996 ( (126.d0*Bp_lm**4)/( 5.d0*Ap_lm**5.5 * Dp_lm**5 ) )
997 R7p_lm = ( 3.d0/( 2.d0*Ap_lm**2.5 * Dp_lm**4 ) )+ &
998 ( (10.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**5 ) )+&
999 ( (35.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**6 ) )
1000 R8p_lm = ( 1.d0/( 2.d0*Ap_lm**1.5 * Dp_lm**4 ) )+ &
1001 ( (6.d0*Bp_lm*Bp_lm)/( Ap_lm**2.5 * Dp_lm**5 ) )+&
1002 ( (25.d0*Bp_lm**4)/( Ap_lm**3.5 * Dp_lm**6 ) )
1003 R9p_lm = ( 1.d0/( Ap_lm**2.5 * Dp_lm**3 ) )+ &
1004 ( (15.d0*Bp_lm*Bp_lm)/( Ap_lm**3.5 * Dp_lm**4 ) )+&
1005 ( (70.d0*Bp_lm**4)/( Ap_lm**4.5 * Dp_lm**5 ) )
1006 K_common_term = DPSUMo2**3 * M_PL*M_PM/(2.d0*MPSUM)*&
1007 (1.d0+C_E)*g0(LL) * (M_PM*M_PL)**1.5
1008
1009
1010
1011 = - K_common_term*NU_PM*NU_PL*(&
1012 ((DPSUMo2*DSQRT(PI)/16.d0)*(3.d0/2.d0)*Bp_lm*R5p_lm)+&
1013 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*PI/6.d0)*&
1014 (3.d0/2.d0)*R1p_lm)-(&
1015 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM/8.d0)*Bp_lm*R6p_lm)+&
1016 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1017 8.d0)*M_PM*R9p_lm)+&
1018 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PL*M_PM/(MPSUM*MPSUM))*&
1019 M_PL*Bp_lm*R7p_lm)+&
1020 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*(DPSUMo2*DSQRT(PI)/&
1021 2.d0)*(M_PM/MPSUM)**2 * M_PL*R8p_lm)+&
1022 ((DPSUMo2*DSQRT(PI)/16.d0)*(M_PM*M_PL/(2.d0*MPSUM))*&
1023 R9p_lm)-&
1024 ((M_PL/(8.d0*MPSUM))*(1.d0-C_E)*DPSUMo2*DSQRT(PI)*&
1025 (M_PM*M_PL/MPSUM)*Bp_lm*R7p_lm) )*TH(LL) )*&
1026 (TH(M)**2 * TH(LL)**3)
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062 = K_s_sum + K_sM_LM
1063
1064
1065
1066
1067 ENDIF
1068 ENDDO
1069
1070
1071 = K_s_sum
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089 CASE (GD_1999)
1090
1091 = DP_avg(M)
1092 M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
1093 NU_PM = (EPS(M)*ROS_avg(M))/M_PM
1094
1095
1096 = ONE + 2.d0*(1.d0+C_E)*EPS(M)*G0(M)
1097
1098
1099 = 5.0d0*M_PM*DSQRT(TH(M)/PI) / (16.d0*D_PM*D_PM)
1100
1101 c_star = 32.0d0*(1.0d0 - C_E)*(1.d0 - 2.0d0*C_E*C_E) &
1102 / (81.d0 - 17.d0*C_E + 30.d0*C_E*C_E*(1.0d0-C_E))
1103
1104 zeta0_star = (5.d0/12.d0)*G0(M)*(1.d0 - C_E*C_E) &
1105 * (1.d0 + (3.d0/32.d0)*c_star)
1106
1107 kappa0 = (15.d0/4.d0)*eta0
1108
1109 nu_kappa_star = (G0(M)/3.d0)*(1.d0+C_E) * ( 1.d0 + &
1110 (33.d0/16.d0)*(1.d0-C_E) + ((19.d0-3.d0*C_E)/1024.d0)*c_star)
1111
1112
1113 = (2.d0/3.d0)*(1.d0 + 0.5d0*(1.d0+press_star)*c_star + &
1114 (3.d0/5.d0)*EPS(M)*G0(M)*(1.d0+C_E)*(1.d0+C_E) * &
1115 (2.d0*C_E - 1.d0 + ( 0.5d0*(1.d0+C_E) - 5.d0/(3*(1.d0+C_E))) * &
1116 c_star ) ) / (nu_kappa_star - 2.d0*zeta0_star)
1117
1118 kappa_star = kappa_k_star * (1.d0 + (6.d0/5.d0)*EPS(M)* &
1119 G0(M)*(1.d0+C_E) ) + (256.d0/25.d0)*(EPS(M)* &
1120 EPS(M)/PI)*G0(M)*(1.d0+C_E)*(1.d0+(7.d0/32.d0)* &
1121 c_star)
1122
1123 IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
1124 Kgran_star = kappa0
1125 ELSEIF(TH(M) .LT. SMALL_NUMBER)THEN
1126 Kgran_star = ZERO
1127 ELSE
1128 Kgran_star = ROS_avg(M)*EPS(M)* G0(M)*TH(M)* kappa0/ &
1129 (ROS_avg(M)*G0(M)*EPS(M)*TH(M) + &
1130 1.2d0*DgA*kappa0/ROs_avg(M))
1131 ENDIF
1132
1133
1134
1135 = Kgran_star * kappa_star
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151 = ROs_avg(M)*EPS(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1152
1153 = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1154
1155
1156 CASE (GTSH_2012)
1157 D_PM = DP_avg(M)
1158 M_PM = (PI/6.d0)*(D_PM**3.)*ROS_avg(M)
1159 NU_PM = (EPS(M)*ROS_avg(M))/M_PM
1160 Chi = g0(M)
1161 vfrac = EPS(M)
1162
1163 Re_g = EPG*RO_g_avg*D_PM*VREL/Mu_g_avg
1164 Re_T = RO_g_avg*D_PM*dsqrt(TH(M)) / Mu_g_avg
1165 mu2_0 = dsqrt(2d0*pi) * Chi * (one-C_E**2)
1166 = (4.5d0+C_E**2) * mu2_0
1167 = (6.46875d0+0.9375d0*C_E**2)*mu2_0 + 2d0*dsqrt(2d0*pi)* &
1168 Chi*(one+C_E)
1169 = zero
1170 if((vfrac> small_number) .and. (TH(M) > small_number)) then
1171 zeta_star = 4.5d0*dsqrt(2d0*Pi)*(RO_g_avg/ROs_avg(M))**2*Re_g**2 * &
1172 S_star(vfrac,Chi) / (vfrac*(one-vfrac)**2 * Re_T**4)
1173 A2_gtshW = (5d0*mu2_0 - mu4_0) / (mu4_1 - 5d0* &
1174 (19d0/16d0*mu2_0 - 1.5d0*zeta_star))
1175 endif
1176 eta0 = 0.3125d0/(dsqrt(pi)*D_PM**2)*M_pm*dsqrt(TH(M))
1177 nu0 = (96.d0/5.d0)*(vfrac/D_PM)*DSQRT(TH(M)/PI)
1178 NuK = nu0*(one+C_E)/3d0*Chi*( one+2.0625d0*(one-C_E)+ &
1179 ((947d0-579*C_E)/256d0*A2_gtshW) )
1180 Kth0 = 3.75d0*eta0/M_pm
1181 EDT_s = 4d0/3d0*dsqrt(pi)*(one-C_E**2)*Chi* &
1182 (one+0.1875d0*A2_gtshW)*NU_PM*D_PM**2*dsqrt(TH(M))
1183 KthK = zero
1184 if(vfrac > small_number) KthK = 2d0/3d0*Kth0*nu0/(NuK - 2d0*EDT_s) * &
1185 (one+2d0*A2_gtshW+0.6d0*vfrac*Chi* &
1186 (one+C_E)**2*(2*C_E-one+A2_gtshW*(one+C_E)))
1187
1188
1189
1190 = KthK*(one+1.2d0*vfrac*Chi*(one+C_E)) + (10.24d0/pi* &
1191 vfrac**2*Chi*(one+C_E)*(one+0.4375d0*A2_gtshW)*Kth0)
1192
1193
1194 = ROs_avg(M)*EPS(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1195 Ps = ROs_avg(M)*EPS(M)*TH(M)*(ONE+2.d0*(ONE+C_E)*g0EPs_avg)
1196
1197
1198 CASE DEFAULT
1199
1200 WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1201 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1202 call mfix_exit(myPE)
1203 END SELECT
1204
1205
1206
1207
1208
1209
1210
1211 SELECT CASE (KT_TYPE_ENUM)
1212 CASE (LUN_1984, SIMONIN_1996, AHMADI_1995, GD_1999)
1213
1214
1215
1216 = K_1
1217
1218
1219 IF(JENKINS) THEN
1220 HW = 3.D0/8.D0*DSQRT(3.D0*TH(M))*((1d0-e_w))*&
1221 PsoTheta
1222
1223
1224
1225
1226 = tan_Phi_w*tan_Phi_w*(ONE+e_w)*21.d0/16.d0*&
1227 DSQRT(3.D0*TH(M)) * Ps
1228
1229 ELSE
1230 HW = (Pi*DSQRT(3d0)/(4.D0*(ONE-ep_star_avg)))*&
1231 (1d0-e_w*e_w)*&
1232 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))
1233
1234 IF(.NOT. BC_JJ_M) THEN
1235 CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*PHIP*&
1236 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1237 ELSE
1238 CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*&
1239 PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1240 EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1241 ENDIF
1242 ENDIF
1243
1244
1245 CASE (GTSH_2012)
1246
1247
1248
1249 = M_PM * K_1
1250
1251 HW = (Pi*DSQRT(3d0)/(4.D0*(ONE-ep_star_avg)))*&
1252 (1d0-e_w*e_w)*&
1253 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))
1254
1255 IF(.NOT. BC_JJ_M) THEN
1256 CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*PHIP*&
1257 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1258 ELSE
1259 CW = (Pi*DSQRT(3d0)/(6.D0*(ONE-ep_star_avg)))*&
1260 PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1261 EPS(M)*g0(M)*DSQRT(TH(M))*VSLIPSQ
1262 ENDIF
1263
1264
1265 CASE (IA_2005)
1266
1267
1268
1269 = M_PM * K_1
1270
1271 HW = (PI*DSQRT(3.d0)/(4.d0*(ONE-ep_star_avg)))*&
1272 (1.d0-e_w*e_w)*&
1273 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M)/M_PM)
1274
1275 IF(.NOT. BC_JJ_M) THEN
1276 CW = (PI*DSQRT(3.d0)/(6.d0*(ONE-ep_star_avg)))*PHIP*&
1277 ROs_avg(M)*EPS(M)*g0(M)*DSQRT(TH(M)*M_PM)*VSLIPSQ
1278 ELSE
1279 CW = (PI*DSQRT(3.d0)/(6.d0*(ONE-ep_star_avg)))*&
1280 PHIP_JJ(DSQRT(vslipsq),TH(M))*ROs_avg(M)*&
1281 EPS(M)*g0(M)*DSQRT(TH(M)*M_PM)*VSLIPSQ
1282 ENDIF
1283
1284
1285
1286 CASE DEFAULT
1287
1288 WRITE (*, '(A)') 'BC_THETA => THETA_HW_CW'
1289 WRITE (*, '(A,A)') 'Unknown KT_TYPE: ', KT_TYPE
1290 call mfix_exit(myPE)
1291 END SELECT
1292
1293
1294 RETURN
1295 END SUBROUTINE THETA_HW_CW
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315 DOUBLE PRECISION FUNCTION phip_jj(uslip,g_theta)
1316
1317
1318
1319 use constant, only : e_w, PI, k4phi, phip0
1320 implicit none
1321
1322
1323
1324 double precision, INTENT(IN) :: uslip
1325 double precision, INTENT(IN) :: g_theta
1326
1327
1328
1329 double precision :: r4phi
1330
1331
1332
1333 =uslip/dsqrt(3.d0*g_theta)
1334
1335 IF(r4phi .le. 4.d0*k4phi/7.d0/dsqrt(6.d0*PI)/phip0) THEN
1336 phip_jj=-7.d0*dsqrt(6.d0*PI)*phip0**2/8.d0/k4phi*r4phi+phip0
1337 ELSE
1338 phip_jj=2.d0/7.d0*k4phi/r4phi/dsqrt(6.d0*PI)
1339 ENDIF
1340
1341 RETURN
1342
1343 END FUNCTION PHIP_JJ
1344