File: /nfs/home/0/users/jenkins/mfix.git/model/calc_u_friction.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 SUBROUTINE CALC_Gw_Hw_Cw(g0, EPS, EPG, ep_star_avg, &
29 g0EP_avg, TH, Mu_g_avg, RO_g_avg, ROs_avg, &
30 DP_avg, K_12_avg, Tau_12_avg, Tau_1_avg, VREL, VSLIP,&
31 DEL_U, S_S, S_dd, VEL, W_VEL, M, gw, hw, cw)
32
33
34
35
36 USE param
37 USE param1
38 USE constant
39 USE physprop
40 USE fldvar
41 USE bc
42 USE run
43 USE mpi_utility
44 IMPLICIT NONE
45
46
47
48
49
50 DOUBLE PRECISION, INTENT(IN) :: g0
51
52 DOUBLE PRECISION, INTENT(IN) :: EPS
53
54 DOUBLE PRECISION, INTENT(IN) :: EPG, ep_star_avg
55
56 DOUBLE PRECISION, INTENT(IN) :: g0EP_avg
57
58 DOUBLE PRECISION, INTENT(INOUT) :: TH
59
60 DOUBLE PRECISION, INTENT(IN) :: Mu_g_avg
61
62 DOUBLE PRECISION, INTENT(IN) :: RO_g_avg
63
64 DOUBLE PRECISION, INTENT(IN) :: ROS_avg
65
66 DOUBLE PRECISION, INTENT(IN) :: DP_avg
67
68 DOUBLE PRECISION, INTENT(IN) :: K_12_avg, Tau_12_avg, Tau_1_avg
69
70 DOUBLE PRECISION, INTENT(IN) :: VREL
71
72 DOUBLE PRECISION, INTENT(IN) :: VSLIP
73
74 DOUBLE PRECISION, INTENT(IN) :: VEL
75
76 DOUBLE PRECISION, INTENT(IN) :: W_VEL
77
78 DOUBLE PRECISION, INTENT(IN) :: DEL_U
79
80 DOUBLE PRECISION, INTENT(IN) :: S_S
81
82 DOUBLE PRECISION, INTENT(IN) :: S_dd
83
84 INTEGER, INTENT(IN) :: M
85
86
87 DOUBLE PRECISION, INTENT(INOUT) :: Gw
88
89 DOUBLE PRECISION, INTENT(INOUT) :: Hw
90
91 DOUBLE PRECISION, INTENT(INOUT) :: Cw
92
93
94
95
96 DOUBLE PRECISION :: Mu_s
97
98 DOUBLE PRECISION :: F_2
99
100 DOUBLE PRECISION :: Pf
101
102 DOUBLE PRECISION :: Pc
103
104 DOUBLE PRECISION :: Chi, N_Pff
105
106 DOUBLE PRECISION :: Mu
107
108 DOUBLE PRECISION :: Mu_b
109
110 DOUBLE PRECISION :: Mu_star
111
112 DOUBLE PRECISION :: Re_g
113
114 DOUBLE PRECISION :: C_d
115
116 DOUBLE PRECISION :: Beta, DgA
117
118 DOUBLE PRECISION :: ZETA
119
120 DOUBLE PRECISION :: Sigma_c, Tau_2_c, Tau_12_st, Nu_t
121 DOUBLE PRECISION :: Tau_2, zeta_c_2, MU_2_T_Kin, Mu_2_Col
122 DOUBLE PRECISION :: Tmp_Ahmadi_Const
123
124 DOUBLE PRECISION :: dpc_dphi
125
126 DOUBLE PRECISION :: g0_loc
127
128
129
130
131 IF(TH .LE. ZERO)THEN
132 TH = 1D-8
133 IF (myPE.eq.PE_IO) THEN
134 WRITE(*,*) &
135 'Warning: Negative granular temp at wall set to 1e-8'
136
137 ENDIF
138 ENDIF
139
140 g0_loc = g0
141
142
143 IF(JENKINS) THEN
144
145 IF (VSLIP == ZERO) THEN
146
147 = zero
148
149 ELSE
150 IF(AHMADI) THEN
151
152
153
154
155
156
157 = tan_Phi_w*ROS_avg*EPS* &
158 ((ONE + 4.0D0*g0EP_avg) + HALF*(ONE -C_e*C_e))*TH/VSLIP
159
160 ELSE
161
162 = tan_Phi_w*ROS_avg*EPS*(1d0+ 4.D0 * Eta *g0EP_avg)*TH/VSLIP
163 ENDIF
164
165 ENDIF
166
167 ELSE
168
169 = (PHIP*DSQRT(3d0*TH)*Pi*ROS_avg*EPS*g0_loc)&
170 /(6d0*(ONE-ep_star_avg))
171
172 ENDIF
173
174
175 = (5d0*DSQRT(Pi*TH)*Dp_avg*ROS_avg)/96d0
176 Mu_b = (256d0*Mu*EPS*g0EP_avg)/(5d0*Pi)
177
178
179 = EPG*RO_g_avg*Dp_avg*VREL/Mu_g_avg
180 IF (Re_g.lt.1000d0) THEN
181 C_d = (24.d0/(Re_g+SMALL_NUMBER))*(ONE + 0.15d0 * Re_g**0.687d0)
182 ELSE
183 C_d = 0.44d0
184 ENDIF
185 DgA = 0.75d0*C_d*Ro_g_avg*EPG*VREL/(Dp_avg*EPG**(2.65d0))
186 IF(VREL == ZERO) DgA = LARGE_NUMBER
187 Beta = SWITCH*EPS*DgA
188
189
190
191
192
193
194 IF(SWITCH == ZERO .OR. Ro_g_avg == ZERO)THEN
195 Mu_star = Mu
196 ELSEIF(TH .LT. SMALL_NUMBER)THEN
197 MU_star = ZERO
198 ELSE
199 Mu_star = ROS_avg*EPS* g0_loc*TH* Mu/ &
200 (ROS_avg*g0EP_avg*TH + 2.0d0*SWITCH*DgA/ROS_avg* Mu)
201 ENDIF
202
203 Mu_s = ((2d0+ALPHA)/3d0)*((Mu_star/(Eta*(2d0-Eta)*&
204 g0_loc))*(ONE+1.6d0*Eta*g0EP_avg&
205 )*(ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
206 g0EP_avg)+(0.6d0*Mu_b*Eta))
207
208
209 = ROS_avg/(DgA+small_number)
210
211 IF(SIMONIN) THEN
212
213 = (ONE+ C_e)*(3.d0-C_e)/5.d0
214 Tau_2_c = DP_avg/(6.d0*EPS*g0_loc*DSQRT(16.d0*(TH+Small_number)/PI))
215 zeta_c_2= 2.D0/5.D0*(ONE+ C_e)*(3.d0*C_e-ONE)
216 Nu_t = Tau_12_avg/Tau_12_st
217 Tau_2 = ONE/(2.D0/Tau_12_st+Sigma_c/Tau_2_c)
218 MU_2_T_Kin = (2.0D0/3.0D0*K_12_avg*Nu_t + TH * &
219 (ONE+ zeta_c_2*EPS*g0_loc))*Tau_2
220 Mu_2_Col = 8.D0/5.D0*EPS*g0_loc*Eta* (MU_2_T_Kin+ &
221 Dp_avg*DSQRT(TH/PI))
222 Mu_s = EPS*ROS_avg*(MU_2_T_Kin + Mu_2_Col)
223 ELSEIF(AHMADI) THEN
224 IF(EPS < (ONE-ep_star_avg)) THEN
225 Tmp_Ahmadi_Const = &
226 ONE/(ONE+ Tau_1_avg/Tau_12_st * (ONE-EPS/(ONE-ep_star_avg))**3)
227 ELSE
228 Tmp_Ahmadi_Const = ONE
229 ENDIF
230 Mu_s = Tmp_Ahmadi_Const &
231 *0.1045D0*(ONE/g0_loc+3.2D0*EPS+12.1824D0*g0_loc*EPS*EPS) &
232 *Dp_avg*ROS_avg* DSQRT(TH)
233 ENDIF
234
235
236 IF ((ONE-EPG)<= EPS_f_min) THEN
237 Pf = ZERO
238 Chi = ZERO
239 ZETA = 1d0
240 ELSE
241 IF (SAVAGE.EQ.1) THEN
242 = ((48d0*Eta*(1d0-Eta)*ROS_avg*EPS*EPS*g0_loc*&
243 (TH**1.5d0))/&
244 (SQRT_Pi*Dp_avg*2d0*Mu_s))**0.5d0
245 ELSEIF (SAVAGE.EQ.0) THEN
246 = DSQRT(S_S)
247 ELSE
248 ZETA = DSQRT(S_S + (TH/(Dp_avg*Dp_avg)))
249 ENDIF
250
251 IF (EPG < ep_star_avg) THEN
252 dpc_dphi = (to_SI*Fr)*((delta**5)*(2d0*(ONE-ep_star_avg-delta) - &
253 2d0*eps_f_min)+((ONE-ep_star_avg-delta)-eps_f_min)&
254 *(5*delta**4))/(delta**10)
255
256 Pc = (to_SI*Fr)*(((ONE-ep_star_avg-delta) - EPS_f_min)**N_Pc)/(delta**D_Pc)
257
258 ELSE
259 Pc = Fr*(((ONE-EPG) - EPS_f_min)**N_Pc)/ &
260 (((ONE-ep_star_avg) - (ONE-EPG) + SMALL_NUMBER)**D_Pc)
261 ENDIF
262
263 IF (DEL_U .GE. ZERO) THEN
264 N_Pff = DSQRT(3d0)/(2d0*Sin_Phi)
265 ELSE
266 N_Pff = N_Pf
267 ENDIF
268
269 IF ((DEL_U/(ZETA*N_Pff*DSQRT(2d0)*Sin_Phi)) .GT. 1d0) THEN
270 Pf = ZERO
271 ELSEIF( DEL_U == ZERO ) THEN
272 Pf = Pc
273 ELSE
274 Pf = Pc*(1d0 - (DEL_U/(ZETA*N_Pff*DSQRT(2d0)*Sin_Phi)))**&
275 (N_Pff-1d0)
276 ENDIF
277
278 Chi = DSQRT(2d0)*Pf*Sin_Phi*(N_Pff - (N_Pff-1d0)*&
279 ((Pf/Pc)**(1d0/(N_Pff-1d0))))
280
281 IF (Chi< ZERO) THEN
282 Pf = Pc*((N_Pff/(N_Pff-1d0))**(N_Pff-1d0))
283 Chi = ZERO
284 ENDIF
285
286
287
288
289 = Pf * EPS/(ONE-EPG)
290 Chi = Chi * EPS/(ONE-EPG)
291
292 ENDIF
293
294
295
296 = ONE
297
298 IF(VSLIP == ZERO .OR. ZETA == ZERO) THEN
299 Hw = ZERO
300 ELSE
301 Hw = F_2 + Pf*tan_Phi_w - Chi*S_dd*tan_Phi_w/ZETA
302 Hw = Hw / (MU_s + Chi/(2d0*ZETA))
303 ENDIF
304 Cw = hw * W_VEL
305
306 RETURN
307 END SUBROUTINE CALC_Gw_Hw_Cw
308