File: /nfs/home/0/users/jenkins/mfix.git/model/calc_p_star.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 SUBROUTINE CALC_P_STAR(EP_G, P_STAR)
27
28
29
30
31 USE param
32 USE param1
33 USE parallel
34 USE geometry
35 USE indices
36 USE physprop
37 USE constant
38 USE pgcor
39 USE pscor
40 USE ur_facs
41 USE residual
42 USE compar
43 USE run
44 USE visc_s
45 USE solids_pressure
46 USE functions
47 IMPLICIT NONE
48
49
50
51
52 DOUBLE PRECISION, INTENT(IN) :: EP_g(DIMENSION_3)
53
54 DOUBLE PRECISION, INTENT(INOUT) :: P_star(DIMENSION_3)
55
56
57
58
59
60
61
62 INTEGER :: IJK
63
64 DOUBLE PRECISION :: blend
65
66
67
68 DOUBLE PRECISION :: CALC_EP_STAR
69 DOUBLE PRECISION, EXTERNAL :: BLEND_FUNCTION
70
71
72
73
74
75 DO IJK = ijkstart3, ijkend3
76 IF (FLUID_AT(IJK)) THEN
77
78
79 IF (YU_STANDISH .OR. FEDORS_LANDEL) THEN
80
81
82
83 (ijk) = calc_ep_star(ijk)
84
85
86
87 IF(BLENDING_STRESS.AND.TANH_BLEND) THEN
88 ep_g_blend_start(ijk) = ep_star_array(ijk) * 0.99d0
89 ep_g_blend_end(ijk) = ep_star_array(ijk) * 1.01d0
90 ELSEIF(BLENDING_STRESS.AND.SIGM_BLEND) THEN
91 ep_g_blend_start(ijk) = EP_star_array(ijk) * 0.97d0
92 ep_g_blend_end(ijk) = EP_star_array(ijk) * 1.01d0
93 ELSE
94 ep_g_blend_start(ijk) = ep_star_array(ijk)
95 ep_g_blend_end(ijk) = ep_star_array(ijk)
96 ENDIF
97 ENDIF
98
99 IF (EP_G(IJK) < EP_g_blend_end(ijk)) THEN
100 P_STAR(IJK) = NEG_H(EP_G(IJK),EP_g_blend_end(ijk))
101 IF(BLENDING_STRESS) THEN
102 blend = blend_function(IJK)
103 P_STAR (IJK) = (1.0d0-blend) * P_STAR (IJK)
104 ENDIF
105 ELSE
106 P_STAR(IJK) = ZERO
107 ENDIF
108 ENDIF
109 ENDDO
110
111 RETURN
112 END SUBROUTINE CALC_P_STAR
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134 DOUBLE PRECISION FUNCTION CALC_ep_star(IJK)
135
136
137
138
139 USE param
140 USE param1
141 USE fldvar
142 USE geometry
143 USE indices
144 USE physprop
145 USE constant
146 USE toleranc
147 USE compar
148 USE run
149 USE fun_avg
150 USE functions
151 IMPLICIT NONE
152
153
154
155
156 INTEGER, INTENT(IN) :: IJK
157
158
159
160
161 INTEGER :: I, J
162
163
164
165
166
167 DOUBLE PRECISION :: P_IT(MMAX)
168
169 DOUBLE PRECISION :: EPs_max_local
170
171 DOUBLE PRECISION :: P_IJ(MMAX, MMAX)
172
173 DOUBLE PRECISION :: R_IJ(MMAX, MMAX)
174
175 DOUBLE PRECISION :: X_IJ(MMAX, MMAX)
176
177
178 DOUBLE PRECISION :: COMP_X_I(MMAX), SUM_LOCAL
179
180
181
182 DOUBLE PRECISION :: DP_TMP(MMAX), EPs_TMP(MMAX), &
183 EPs_max_TMP(MMAX), old_value
184
185
186
187 IF (CALL_DQMOM) THEN
188
189
190 DO I = 1, SMAX
191 DP_TMP(I) = D_P(IJK,I)
192 EPs_TMP(I) = EP_s(IJK,I)
193 EPs_max_TMP(I) = ep_s_max(I)
194 ENDDO
195
196 DO I = 1, SMAX
197 DO J = I , SMAX
198
199 IF(DP_TMP(I) < DP_TMP(J)) THEN
200
201 = DP_TMP(I)
202
203 (I) = DP_TMP(J)
204
205 (J) = old_value
206
207 old_value = EPs_TMP(I)
208 EPs_TMP(I) = EPs_TMP(J)
209 EPs_TMP(J) = old_value
210
211 old_value = EPs_max_TMP(I)
212 EPs_max_TMP(I) = EPs_max_TMP(J)
213 EPs_max_TMP(J) = old_value
214 ENDIF
215 ENDDO
216 ENDDO
217
218 ELSE
219
220
221 DO I = 1, SMAX
222 DP_TMP(I) = D_P(IJK,M_MAX(I))
223 EPs_TMP(I) = EP_s(IJK,M_MAX(I))
224 EPs_max_TMP(I) = ep_s_max(M_MAX(I))
225 ENDDO
226 ENDIF
227
228
229
230
231
232 DO I = 1, SMAX
233 SUM_LOCAL = ZERO
234 DO J = 1, SMAX
235 IF(I .GE. J) THEN
236 R_IJ(I,J) = DP_TMP(I)/DP_TMP(J)
237 ELSE
238 R_IJ(I,J) = DP_TMP(J)/DP_TMP(I)
239 ENDIF
240 SUM_LOCAL = SUM_LOCAL + EPs_TMP(J)
241 ENDDO
242
243 IF(SUM_LOCAL > DIL_EP_s) THEN
244
245 (I) = EPs_TMP(I)/SUM_LOCAL
246 ELSE
247
248 = ONE - EPs_max_TMP(1)
249 RETURN
250 ENDIF
251 ENDDO
252
253
254
255 IF(YU_STANDISH) THEN
256
257 DO I = 1, SMAX
258 DO J = 1, SMAX
259 IF(R_IJ(I,J) .LE. 0.741d0) THEN
260 IF(J .LT. I) THEN
261 X_IJ(I,J) = (ONE - R_IJ(I,J)*R_IJ(I,J))/&
262 (2.0d0 - EPs_max_TMP(I))
263 ELSE
264 X_IJ(I,J) = ONE - (ONE - R_IJ(I,J)*R_IJ(I,J))/&
265 (2.0d0 - EPs_max_TMP(I))
266 ENDIF
267 P_IJ(I, J) = EPs_max_TMP(I) + EPs_max_TMP(I)*&
268 (ONE-EPs_max_TMP(I)) * (ONE - 2.35d0*R_IJ(I,J) + &
269 1.35d0*R_IJ(I,J)*R_IJ(I,J) )
270 ELSE
271 P_IJ(I, J) = EPs_max_TMP(I)
272 ENDIF
273 ENDDO
274 ENDDO
275
276
277 = ONE
278 DO I = 1, SMAX
279 SUM_LOCAL = ZERO
280
281 IF(I .GE. 2) THEN
282 DO J = 1, (I-1)
283 IF(P_IJ(I,J) == EPs_max_TMP(I)) THEN
284 SUM_LOCAL = SUM_LOCAL
285 ELSE
286 SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
287 P_IJ(I,J))*COMP_X_I(J)/X_IJ(I,J)
288 ENDIF
289 ENDDO
290 ENDIF
291
292 IF((I+1) .LE. SMAX) THEN
293 DO J = (I+1), SMAX
294 IF( P_IJ(I, J) == EPs_max_TMP(I) ) THEN
295 SUM_LOCAL = SUM_LOCAL
296 ELSE
297 SUM_LOCAL = SUM_LOCAL + (ONE - EPs_max_TMP(I)/&
298 P_IJ(I, J))*COMP_X_I(J)/X_IJ(I, J)
299 ENDIF
300 ENDDO
301 ENDIF
302
303 IF (SUM_LOCAL .NE. ZERO) THEN
304 P_IT(I) = EPs_max_TMP(I)/(ONE - SUM_LOCAL)
305 ELSE
306
307 (I) = ONE
308 ENDIF
309
310 EPs_max_local = MIN(P_IT(I), EPs_max_local)
311 ENDDO
312
313
314
315 IF (EPs_max_local == ONE) EPs_max_local = EPs_max_TMP(1)
316 CALC_EP_star = ONE - EPs_max_local
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343 ELSEIF(FEDORS_LANDEL) THEN
344
345 IF(COMP_X_I(1) .LE. (EPs_max_TMP(1)/(EPs_max_TMP(1)+ &
346 (ONE - EPs_max_TMP(1))*EPs_max_TMP(2))) ) THEN
347
348 CALC_EP_star = (EPs_max_TMP(1) - EPs_max_TMP(2) + &
349 (1 - sqrt(R_IJ(2,1))) * (ONE - EPs_max_TMP(1)) * &
350 EPs_max_TMP(2) )*&
351 (EPs_max_TMP(1) + (ONE - EPs_max_TMP(1)) * &
352 EPs_max_TMP(2)) * COMP_X_I(1)/EPs_max_TMP(1) + &
353 EPs_max_TMP(2)
354 ELSE
355 CALC_EP_star = (ONE-sqrt(R_IJ(2,1))) * (EPs_max_TMP(1)+&
356 (ONE-EPs_max_TMP(1)) * EPs_max_TMP(2)) * &
357 (ONE - COMP_X_I(1)) + EPs_max_TMP(1)
358 ENDIF
359
360 = ONE - CALC_EP_star
361 ENDIF
362
363
364
365 RETURN
366 END FUNCTION CALC_ep_star
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387 DOUBLE PRECISION FUNCTION blend_function(IJK)
388
389
390
391
392 USE param
393 USE param1
394 USE fldvar
395 USE geometry
396 USE indices
397 USE physprop
398 USE constant
399 USE toleranc
400 USE compar
401 USE run
402 USE visc_s
403 USE fun_avg
404 USE functions
405 IMPLICIT NONE
406
407
408
409
410 INTEGER, INTENT(IN) :: IJK
411
412
413
414
415 LOGICAL,SAVE:: FIRST_PASS = .TRUE.
416
417 Double Precision:: blend, blend_right
418
419 Double Precision, Save:: scale
420
421 Double Precision, Save:: ep_mid_point
422
423
424
425 IF(TANH_BLEND) THEN
426 IF(EP_g(IJK) .LT. ep_g_blend_end(ijk).AND. EP_g(IJK) .GT. ep_g_blend_start(ijk)) THEN
427 ep_mid_point = (ep_g_blend_end(IJK)+ep_g_blend_start(IJK))/2.0d0
428 blend = tanh(2.0d0*pi*(ep_g(IJK)-ep_mid_point)/ &
429 (ep_g_blend_end(IJK)-ep_g_blend_start(IJK)))
430 blend = (blend+1.0d0)/2.0d0
431 ELSEIF(EP_g(IJK) .GE. ep_g_blend_end(ijk)) THEN
432 blend = 1.0d0
433 ELSEIF(EP_g(IJK) .LE. ep_g_blend_start(ijk)) THEN
434 blend = 0.0d0
435 ENDIF
436
437
438 ELSEIF(SIGM_BLEND) THEN
439 IF(FIRST_PASS) THEN
440 blend_right = 1.0d0/(1+0.01d0**((ep_g_blend_end(IJK)-ep_star_array(IJK))&
441 /(ep_g_blend_end(IJK)-ep_g_blend_start(IJK))))
442 blend_right = (blend_right+1.0d0)/2.0d0
443 scale = 1.0d0/blend_right
444 write(*,*) 'Blending value at end and scaling factor', blend_right, scale
445 FIRST_PASS = .FALSE.
446 ENDIF
447 IF(EP_g(IJK) .LT. ep_g_blend_end(ijk)) THEN
448 blend = scale/(1+0.01d0**((ep_g(IJK)-ep_star_array(IJK))&
449 /(ep_g_blend_end(IJK)-ep_g_blend_start(IJK))))
450 ELSEIF(EP_g(IJK) .GE. ep_g_blend_end(ijk)) THEN
451 blend = 1.0d0
452 ENDIF
453
454 ENDIF
455
456 blend_function = blend
457
458 RETURN
459 END FUNCTION blend_function
460
461