File: /nfs/home/0/users/jenkins/mfix.git/model/calc_k_cp.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 SUBROUTINE CALC_K_cp(Kcp)
26
27
28
29
30 USE compar
31 USE constant
32 USE fldvar
33 USE functions
34 USE geometry
35 USE indices
36 USE param
37 USE param1
38 USE physprop
39 USE pscor
40 USE rdf
41 USE run
42 USE sendrecv
43 USE solids_pressure
44 USE trace
45 USE visc_s
46 IMPLICIT NONE
47
48
49
50
51 DOUBLE PRECISION :: Kcp(DIMENSION_3)
52
53
54
55
56 INTEGER :: IJK, M
57
58 DOUBLE PRECISION :: Pc, DPcoDEPs, Mu, Mu_b, Mu_zeta, ZETA
59 DOUBLE PRECISION :: F2, DF2oDEPs, Pf, Pfmax, N_Pff
60
61 Double Precision :: blend
62
63
64
65 DOUBLE PRECISION :: DZETAoDEPs
66 DOUBLE PRECISION, EXTERNAL :: BLEND_FUNCTION
67
68
69
70 (:) = ZERO
71
72 IF (MCP == UNDEFINED_I) THEN
73
74
75 RETURN
76 ELSE
77
78
79
80 = MCP
81 ENDIF
82
83
84
85 IF(CLOSE_PACKED(M)) THEN
86
87
88
89
90
91
92 DO IJK = ijkstart3, ijkend3
93 IF(.NOT.WALL_AT(IJK))THEN
94
95 IF (FRICTION) THEN
96
97 IF ((ONE-EP_G(IJK)).GT.EPS_f_min) THEN
98
99
100
101 IF ((ONE-EP_G(IJK)).GT.(ONE-ep_star_array(ijk))) THEN
102
103
104
105 = (to_SI*Fr)*((delta**5)*&
106 (2d0*(ONE-ep_star_array(IJK)-delta) - &
107 2d0*eps_f_min)+&
108 ((ONE-ep_star_array(ijk)-delta)-eps_f_min)*&
109 (5*delta**4))/(delta**10)
110
111 Pc = (to_SI*Fr)*&
112 ( ((ONE-ep_star_array(IJK)-delta)-&
113 EPS_f_min)**N_Pc )/(delta**D_Pc)
114 Pc = Pc + DPcoDEPS*( (ONE-EP_G(IJK))+delta-&
115 (ONE-ep_star_array(IJK)))
116
117
118
119
120
121 ELSE
122 Pc = Fr*(((ONE-EP_G(IJK)) - EPS_f_min)**N_Pc)/&
123 (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK))&
124 + SMALL_NUMBER)**D_Pc)
125 DPcoDEPs =&
126 Fr*(((ONE-EP_G(IJK)) - EPS_f_min)**(N_Pc - ONE))&
127 *(N_Pc*((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)))&
128 +D_Pc*((ONE-EP_G(IJK)) - EPS_f_min))&
129 / (((ONE-ep_star_array(ijk)) - (ONE-EP_G(IJK)) + &
130 SMALL_NUMBER)**(D_Pc + ONE))
131 ENDIF
132
133 Mu = (5d0*DSQRT(Pi*Theta_m(IJK,M))&
134
135 *D_p(IJK,M)*RO_S(IJK,M))/96d0
136
137 Mu_b = (256d0*Mu*EP_s(IJK,M)*EP_s(IJK,M)&
138 *G_0(IJK,M,M))/(5d0*Pi)
139
140 IF (SAVAGE.EQ.1) THEN
141 Mu_zeta =&
142 ((2d0+ALPHA)/3d0)*((Mu/(Eta*(2d0-Eta)*&
143 G_0(IJK,M,M)))*(1d0+1.6d0*Eta*EP_s(IJK,M)*&
144 G_0(IJK,M,M))*(1d0+1.6d0*Eta*(3d0*Eta-2d0)*&
145 EP_s(IJK,M)*G_0(IJK,M,M))+(0.6d0*Mu_b*Eta))
146
147 = ((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*EP_s(IJK,M)*&
148 EP_s(IJK,M)*G_0(IJK,M,M)*&
149 (Theta_m(IJK,M)**1.5d0))/&
150 (SQRT_Pi*D_p(IJK,M)*2d0*Mu_zeta))**0.5d0
151
152 ELSEIF (SAVAGE.EQ.0) THEN
153 ZETA = (SMALL_NUMBER +&
154 trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
155 trD_s_C(IJK,M))/3.d0))**0.5d0
156
157 ELSE
158 ZETA = ((Theta_m(IJK,M)/(D_p(IJK,M)*D_p(IJK,M))) +&
159 (trD_s2(IJK,M) - ((trD_s_C(IJK,M)*&
160 trD_s_C(IJK,M))/3.d0)))**0.5d0
161
162 ENDIF
163
164 IF (trD_s_C(IJK,M) .GE. ZERO) THEN
165 N_Pff = DSQRT(3d0)/(2d0*Sin_Phi)
166 ELSE
167 N_Pff = N_Pf
168 ENDIF
169
170 IF ((trD_s_C(IJK,M)/(ZETA*N_Pff*DSQRT(2d0)*&
171 Sin_Phi)) .GT. 1d0) THEN
172 F2 = 0d0
173 DF2oDEPs = ZERO
174 ELSEIF(trD_s_C(IJK,M) == ZERO) THEN
175 F2 = ONE
176 DF2oDEPs = ZERO
177 ELSE
178 F2 = (1d0 - (trD_s_C(IJK,M)/(ZETA*N_Pff*&
179 DSQRT(2d0)*Sin_Phi)))**(N_Pff-1d0)
180 IF (SAVAGE.EQ.1) THEN
181 DF2oDEPs = (N_Pff-1d0)*(F2**(N_Pff-2d0))*&
182 trD_s_C(IJK,M)&
183 *DZETAoDEPs(EP_s(IJK,M), IJK, M)&
184 / (ZETA*ZETA*&
185 N_Pff*DSQRT(2d0)*Sin_Phi)
186 ELSE
187 DF2oDEPs=ZERO
188 ENDIF
189
190 Pf = Pc*F2
191 Pfmax = Pc*((N_Pf/(N_Pf-1d0))**(N_Pf-1d0))
192
193 IF (Pf> Pfmax) THEN
194 F2 = (N_Pf/(N_Pf-1d0))**(N_Pf-1d0)
195 DF2oDEPS = ZERO
196 ENDIF
197 ENDIF
198
199
200
201
202 (IJK) = F2*DPcoDEPS + Pc*DF2oDEPS
203
204 ELSE
205
206 (IJK) = ZERO
207
208 ENDIF
209
210
211
212 ELSE
213
214
215 IF(EP_g(IJK) .LT. ep_g_blend_end(ijk)) THEN
216
217
218
219
220 (IJK) = dPodEP_s(EP_s(IJK, M),ep_g_blend_end(ijk))
221
222 IF(BLENDING_STRESS) THEN
223 blend = blend_function(IJK)
224 Kcp(IJK) = (1.0d0-blend) * Kcp(IJK)
225 ENDIF
226 ELSE
227
228
229 (IJK) = ZERO
230 ENDIF
231
232
233
234 ENDIF
235
236 ELSE
237 (IJK) = ZERO
238 ENDIF
239
240 ENDDO
241 ENDIF
242
243 CALL send_recv(Kcp, 2)
244
245 RETURN
246 END
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267 DOUBLE PRECISION FUNCTION DZETAoDEPs(EPs, IJK, M)
268
269
270
271
272 USE compar
273 USE constant
274 USE fldvar
275 USE functions
276 USE geometry
277 USE indices
278 USE param
279 USE param1
280 USE physprop
281 USE pscor
282 USE rdf
283 USE run
284 USE trace
285 USE visc_s
286 IMPLICIT NONE
287
288
289
290
291 DOUBLE PRECISION, INTENT(IN) :: EPs
292
293 INTEGER, INTENT(IN) :: IJK,M
294
295
296
297
298 DOUBLE PRECISION :: g0
299
300 DOUBLE PRECISION :: Mu, Mu_b, DEPs2G_0oDEPs, F1, DF1oDEPs
301
302
303 = G_0(IJK, M, M)
304
305 = (5d0*DSQRT(Pi*Theta_m(IJK,M))*D_p(IJK,M)*RO_S(IJK,M))/96d0
306
307 Mu_b = (256d0*Mu*EPs*EPs*g0/(5d0*Pi))
308
309 DEPs2G_0oDEPs = EPs*EPs*DG_0DNU(EPs) + 2d0*EPs*g0
310
311 F1 = ((2d0+ALPHA)/3d0)*((2*Mu/(Eta*(2d0-Eta)*&
312 g0))*(1d0+1.6d0*Eta*EPs*g0)*(1d0+1.6d0*Eta*(3d0*Eta-2d0)&
313 *EPs*g0)+(1.2d0*Mu_b*Eta))
314
315 DF1oDEPs = ((2d0+ALPHA)/3d0)*((2*Mu/(Eta*(2d0-Eta))*&
316 ((-DG_0DNU(EPs)/(g0*g0)) + (1.6d0*Eta*(3d0*Eta-1d0))&
317 + (64d0*Eta*Eta*(3d0*Eta-2d0)*DEPs2G_0oDEPs/25d0))) +&
318 3.2d0*Eta*RO_S(IJK,M)*D_p(IJK,M)*((Theta_m(IJK,M)/Pi)**0.5d0)&
319 *DEPs2G_0oDEPs)
320
321 DZETAoDEPs = 0.5d0*((48d0*Eta*(1d0-Eta)*RO_S(IJK,M)*F1*&
322 (Theta_m(IJK,M)**1.5d0)/&
323 (SQRT_Pi*D_p(IJK,M)*EPs*EPs*g0))**0.5d0)*&
324 (F1*DEPs2G_0oDEPs - EPs*Eps*g0*DF1oDEPs)/(F1*F1)
325
326
327 RETURN
328 END
329
330