File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/ghdmassflux.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20 SUBROUTINE GHDMASSFLUX ()
21
22
23
24
25 USE param
26 USE param1
27 USE geometry
28 USE compar
29 USE fldvar
30 USE indices
31 USE visc_s
32 USE ghdtheory
33 USE physprop
34 USE run
35 USE constant
36 USE toleranc
37 USE drag
38 USE fun_avg
39 USE functions
40 IMPLICIT NONE
41
42
43
44
45 INTEGER :: IJK, I, J, K
46
47 INTEGER :: IJKE, IJKN, IJKT
48
49 INTEGER :: M, L
50
51 DOUBLE PRECISION NjC, NjE, NjN, NjT
52
53
54 DOUBLE PRECISION :: Mi, Mj
55
56 DOUBLE PRECISION :: Ni
57
58 DOUBLE PRECISION :: ropsE, ropsN, ropsT, ThetaE, ThetaN, ThetaT
59 DOUBLE PRECISION :: EPSA
60
61 DOUBLE PRECISION :: DiTE, DiTN, DiTT
62 DOUBLE PRECISION :: DijE, DijN, DijT
63 DOUBLE PRECISION :: DijFE, DijFN, DijFT
64
65 DOUBLE PRECISION :: ordinDiffTermX, ordinDiffTermY, ordinDiffTermZ
66 DOUBLE PRECISION :: massMobilityTermX, massMobilityTermY, &
67 massMobilityTermZ
68 DOUBLE PRECISION :: massMobilityTermXvelupdate, &
69 massMobilityTermYvelupdate, &
70 massMobilityTermZvelupdate
71 DOUBLE PRECISION :: thermalDiffTermX, thermalDiffTermY, &
72 thermalDiffTermZ
73 DOUBLE PRECISION :: ropsme, ropsmn, ropsmt
74
75 DOUBLE PRECISION :: addtermx, &
76 addtermy, addtermz
77 DOUBLE PRECISION :: massMobilityTermNoDragX, &
78 massMobilityTermNoDragY, &
79 massMobilityTermNoDragZ
80 DOUBLE PRECISION :: DiTE_H, DiTE_A, DiTN_H, DiTN_A, DiTT_H, DiTT_A
81 DOUBLE PRECISION :: DijE_H, DijE_A, DijN_H, DijN_A, DijT_H, DijT_A
82 DOUBLE PRECISION :: DijFE_H, DijFE_A, DijFN_H, DijFN_A, DijFT_H, DijFT_A
83
84
85
86
87 DO M = 1, SMAX
88 DO 200 IJK = ijkstart3, ijkend3
89 I = I_OF(IJK)
90 J = J_OF(IJK)
91 K = K_OF(IJK)
92
93 IF ( FLUID_AT(IJK) ) THEN
94 Mi = (PI/6.d0)*D_P(IJK,M)**3 * RO_S(IJK,M)
95 Ni = ROP_s(IJK,M) / Mi
96
97 IJKE = EAST_OF(IJK)
98 IJKN = NORTH_OF(IJK)
99 IJKT = TOP_OF(IJK)
100
101
102 = AVG_X(ROP_S(IJK,MMAX),ROP_S(IJKE,MMAX),I)
103 ropsN = AVG_Y(ROP_S(IJK,MMAX),ROP_S(IJKN,MMAX),J)
104 ropsT = AVG_Z(ROP_S(IJK,MMAX),ROP_S(IJKT,MMAX),K)
105
106 ThetaE = AVG_X(THETA_M(IJK,MMAX),THETA_M(IJKE,MMAX),I)
107 ThetaN = AVG_Y(THETA_M(IJK,MMAX),THETA_M(IJKN,MMAX),J)
108 ThetaT = AVG_Z(THETA_M(IJK,MMAX),THETA_M(IJKT,MMAX),K)
109
110
111
112 = AVG_X_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
113 DiT(IJKE,M)*ROP_S(IJKE,MMAX)/Theta_m(IJKE,MMAX),I)
114 DiTE_A = AVG_X(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
115 DiT(IJKE,M)*ROP_S(IJKE,MMAX)/Theta_m(IJKE,MMAX),I)
116
117 DiTN_H = AVG_Y_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
118 DiT(IJKN,M)*ROP_S(IJKN,MMAX)/Theta_m(IJKN,MMAX),J)
119 DiTN_A = AVG_Y(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
120 DiT(IJKN,M)*ROP_S(IJKN,MMAX)/Theta_m(IJKN,MMAX),J)
121
122 DiTT_H = AVG_Z_S(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
123 DiT(IJKT,M)*ROP_S(IJKT,MMAX)/Theta_m(IJKT,MMAX),K)
124 DiTT_A = AVG_Z(DiT(IJK,M)*ROP_S(IJK,MMAX)/Theta_m(IJK,MMAX),&
125 DiT(IJKT,M)*ROP_S(IJKT,MMAX)/Theta_m(IJKT,MMAX),K)
126
127 IF(M .eq. 1)THEN
128 IF(MIN(ABS(DiTE_H),ABS(DiTE_A)) .eq. ABS(DiTE_H))THEN
129 DiTE = DiTE_H
130 DiT_HarmE(IJK) = .TRUE.
131 ELSE
132 DiTE = DiTE_A
133 DiT_HarmE(IJK) = .FALSE.
134 ENDIF
135
136 IF(MIN(ABS(DiTN_H),ABS(DiTN_A)) .eq. ABS(DiTN_H))THEN
137 DiTN = DiTN_H
138 DiT_HarmN(IJK) = .TRUE.
139 ELSE
140 DiTN = DiTN_A
141 DiT_HarmN(IJK) = .FALSE.
142 ENDIF
143
144 IF(MIN(ABS(DiTT_H),ABS(DiTT_A)) .eq. ABS(DiTT_H))THEN
145 DiTT = DiTT_H
146 DiT_HarmT(IJK) = .TRUE.
147 ELSE
148 DiTT = DiTT_A
149 DiT_HarmT(IJK) = .FALSE.
150 ENDIF
151 ELSE
152 IF(DiT_HarmE(IJK))THEN
153 DiTE = DiTE_H
154 ELSE
155 DiTE = DiTE_A
156 ENDIF
157
158 IF(DiT_HarmN(IJK))THEN
159 DiTN = DiTN_H
160 ELSE
161 DiTN = DiTN_A
162 ENDIF
163
164 IF(DiT_HarmT(IJK))THEN
165 DiTT = DiTT_H
166 ELSE
167 DiTT = DiTT_A
168 ENDIF
169 ENDIF
170
171
172 = ZERO
173 ordinDiffTermY = ZERO
174 ordinDiffTermZ = ZERO
175 massMobilityTermX = ZERO
176 massMobilityTermY = ZERO
177 massMobilityTermZ = ZERO
178 massMobilityTermXvelUpdate = ZERO
179 massMobilityTermYvelUpdate = ZERO
180 massMobilityTermZvelUpdate = ZERO
181 addtermx = ZERO
182 addtermy = ZERO
183 addtermz = ZERO
184 massMobilityTermNoDragX = ZERO
185 massMobilityTermNoDragY = ZERO
186 massMobilityTermNoDragZ = ZERO
187
188 DO L = 1, SMAX
189 Mj = (PI/6.d0)*D_P(IJK,L)**3 * RO_S(IJK,L)
190
191 NjC = ROP_s(IJK,L) / Mj
192 NjE = ROP_S(IJKE,L) / Mj
193 NjN = ROP_S(IJKN,L) / Mj
194 NjT = ROP_S(IJKT,L) / Mj
195
196 IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
197 (ROP_S(IJKE,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
198 DijE_H = AVG_X_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
199 Dij(IJKE,M,L)*Mi*Mj/ROP_S(IJKE,MMAX),I)
200 DijE_A = AVG_X(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
201 Dij(IJKE,M,L)*Mi*Mj/ROP_S(IJKE,MMAX),I)
202
203 IF(M .eq. 1)THEN
204 IF(MIN(ABS(DijE_H),ABS(DijE_A)) .eq. ABS(DijE_H))THEN
205 DijE = DijE_H
206 Dij_HarmE(IJK,L) = .TRUE.
207 ELSE
208 DijE = DijE_A
209 Dij_HarmE(IJK,L) = .FALSE.
210 ENDIF
211 ELSE
212 IF(Dij_HarmE(IJK,L))THEN
213 DijE = DijE_H
214 ELSE
215 DijE = DijE_A
216 ENDIF
217 ENDIF
218 ELSE
219 DijE = ZERO
220 ENDIF
221
222 IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
223 (ROP_S(IJKN,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
224 DijN_H = AVG_Y_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
225 Dij(IJKN,M,L)*Mi*Mj/ROP_S(IJKN,MMAX),J)
226 DijN_A = AVG_Y(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
227 Dij(IJKN,M,L)*Mi*Mj/ROP_S(IJKN,MMAX),J)
228 IF(M .eq. 1)THEN
229 IF(MIN(ABS(DijN_H),ABS(DijN_A)) .eq. ABS(DijN_H))THEN
230 DijN = DijN_H
231 Dij_HarmN(IJK,L) = .TRUE.
232 ELSE
233 DijN = DijN_A
234 Dij_HarmN(IJK,L) = .FALSE.
235 ENDIF
236 ELSE
237 IF(Dij_HarmN(IJK,L))THEN
238 DijN = DijN_H
239 ELSE
240 DijN = DijN_A
241 ENDIF
242 ENDIF
243 ELSE
244 DijN = ZERO
245 ENDIF
246
247 IF((ROP_S(IJK,MMAX)/RO_S(IJK,M) > DIL_EP_S) .and. &
248 (ROP_S(IJKT,MMAX)/RO_S(IJK,M) > DIL_EP_S))THEN
249 DijT_H = AVG_Z_S(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
250 Dij(IJKT,M,L)*Mi*Mj/ROP_S(IJKT,MMAX),K)
251 DijT_A = AVG_Z(Dij(IJK,M,L)*Mi*Mj/ROP_S(IJK,MMAX),&
252 Dij(IJKT,M,L)*Mi*Mj/ROP_S(IJKT,MMAX),K)
253 IF(M .eq. 1)THEN
254 IF(MIN(ABS(DijT_H),ABS(DijT_A)) .eq. ABS(DijT_H))THEN
255 DijT = DijT_H
256 Dij_HarmT(IJK,L) = .TRUE.
257 ELSE
258 DijT = DijT_A
259 Dij_HarmT(IJK,L) = .FALSE.
260 ENDIF
261 ELSE
262 IF(Dij_HarmT(IJK,L))THEN
263 DijT = DijT_H
264 ELSE
265 DijT = DijT_A
266 ENDIF
267 ENDIF
268 ELSE
269 DijT = ZERO
270 ENDIF
271
272
273 DijFE_H = AVG_X_S(DijF(IJK,M,L),DijF(IJKE,M,L),I)
274 DijFE_A = AVG_X(DijF(IJK,M,L),DijF(IJKE,M,L),I)
275 DijFN_H = AVG_Y_S(DijF(IJK,M,L),DijF(IJKN,M,L),J)
276 DijFN_A = AVG_Y(DijF(IJK,M,L),DijF(IJKN,M,L),J)
277 DijFT_H = AVG_Z_S(DijF(IJK,M,L),DijF(IJKT,M,L),K)
278 DijFT_A = AVG_Z(DijF(IJK,M,L),DijF(IJKT,M,L),K)
279
280 IF(M .eq. 1)THEN
281 IF(MIN(ABS(DijFE_H),ABS(DijFE_A)) .eq. ABS(DijFE_H))THEN
282 DijFE = DijFE_H
283 DijF_HarmE(IJK,L) = .TRUE.
284 ELSE
285 DijFE = DijFE_A
286 DijF_HarmE(IJK,L) = .FALSE.
287 ENDIF
288
289 IF(MIN(ABS(DijFN_H),ABS(DijFN_A)) .eq. ABS(DijFN_H))THEN
290 DijFN = DijFN_H
291 DijF_HarmN(IJK,L) = .TRUE.
292 ELSE
293 DijFN = DijFN_A
294 DijF_HarmN(IJK,L) = .FALSE.
295 ENDIF
296
297 IF(MIN(ABS(DijFT_H),ABS(DijFT_A)) .eq. ABS(DijFT_H))THEN
298 DijFT = DijFT_H
299 DijF_HarmT(IJK,L) = .TRUE.
300 ELSE
301 DijFT = DijFT_A
302 DijF_HarmT(IJK,L) = .FALSE.
303 ENDIF
304 ELSE
305 IF(DijF_HarmE(IJK,L))THEN
306 DijFE = DijFE_H
307 ELSE
308 DijFE = DijFE_A
309 ENDIF
310
311 IF(DijF_HarmN(IJK,L))THEN
312 DijFN = DijFN_H
313 ELSE
314 DijFN = DijFN_A
315 ENDIF
316
317 IF(DijF_HarmT(IJK,L))THEN
318 DijFT = DijFT_H
319 ELSE
320 DijFT = DijFT_A
321 ENDIF
322 ENDIF
323
324 = ordinDiffTermX + DijE * (NjE - NjC) * oDX_E(I)
325 ordinDiffTermY = ordinDiffTermY + DijN * (NjN - NjC) * oDY_N(J)
326 ordinDiffTermZ = ordinDiffTermZ + DijT * (NjT - NjC) * (oX_E(I)*oDZ_T(K))
327
328 massMobilityTermX = massMobilityTermX + DijFE * FiX(IJK,L)
329 massMobilityTermY = massMobilityTermY + DijFN * FiY(IJK,L)
330 massMobilityTermZ = massMobilityTermZ + DijFT * FiZ(IJK,L)
331
332 massMobilityTermXvelUpdate = massMobilityTermXvelUpdate + DijFE * FiXvel(IJK,L)
333 massMobilityTermYvelUpdate = massMobilityTermYvelUpdate + DijFN * FiYvel(IJK,L)
334 massMobilityTermZvelUpdate = massMobilityTermZvelUpdate + DijFT * FiZvel(IJK,L)
335
336 massMobilityTermNoDragX = massMobilityTermNoDragX + DijFE * FiMinusDragX(IJK,L)
337 massMobilityTermNoDragY = massMobilityTermNoDragY + DijFN * FiMinusDragY(IJK,L)
338 massMobilityTermNoDragZ = massMobilityTermNoDragZ + DijFT * FiMinusDragZ(IJK,L)
339
340 ENDDO
341
342
343 = DiTE * ( THETA_M(IJKE,MMAX) - THETA_M(IJK,MMAX) ) * oDX_E(I)
344 thermalDiffTermY = DiTN * ( THETA_M(IJKN,MMAX) - THETA_M(IJK,MMAX) ) * oDY_N(J)
345 thermalDiffTermZ = DiTT * ( THETA_M(IJKT,MMAX) - THETA_M(IJK,MMAX) ) * (oX_E(I)*oDZ_T(K))
346
347 JoiX(IJK,M) = -(ordinDiffTermX + thermalDiffTermX + massMobilityTermX)
348 JoiY(IJK,M) = -(ordinDiffTermY + thermalDiffTermY + massMobilityTermY)
349 JoiZ(IJK,M) = -(ordinDiffTermZ + thermalDiffTermZ + massMobilityTermZ)
350
351 JoiMinusDragX(IJK,M) = (ordinDiffTermX + thermalDiffTermX + massMobilityTermNoDragX)
352 JoiMinusDragY(IJK,M) = (ordinDiffTermY + thermalDiffTermY + massMobilityTermNoDragY)
353 JoiMinusDragZ(IJK,M) = (ordinDiffTermZ + thermalDiffTermZ + massMobilityTermNoDragZ)
354
355 ropsME=AVG_X(ROP_S(IJK,M),ROP_S(IJKE,M),I)
356 ropsMN=AVG_Y(ROP_S(IJK,M),ROP_S(IJKN,M),J)
357 ropsMT=AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K)
358
359 DELTAU(IJK,M) = -(ordinDiffTermX+thermalDiffTermX+massMobilityTermXvelupdate)
360 DELTAV(IJK,M) = -(ordinDiffTermY+thermalDiffTermY+massMobilityTermYvelupdate)
361 DELTAW(IJK,M) = -(ordinDiffTermZ+thermalDiffTermz+massMobilityTermZvelUpdate)
362
363
364 = AVG_X(ROP_S(IJK,M),ROP_S(IJKE,M),I) / RO_S(IJK,M)
365 IF(EPSA <= ZERO_EP_S) JoiX(IJK,M) = ZERO
366
367 EPSA = AVG_Y(ROP_S(IJK,M),ROP_S(IJKN,M),J) / RO_S(IJK,M)
368 IF(EPSA <= ZERO_EP_S) JoiY(IJK,M) = ZERO
369
370 EPSA = AVG_Z(ROP_S(IJK,M),ROP_S(IJKT,M),K) / RO_S(IJK,M)
371 IF(EPSA <= ZERO_EP_S) JoiZ(IJK,M) = ZERO
372
373
374 IF (IP_AT_E(IJK)) JoiX(IJK,M) = ZERO
375 IF (IP_AT_N(IJK)) JoiY(IJK,M) = ZERO
376 IF (IP_AT_T(IJK)) JoiZ(IJK,M) = ZERO
377
378 ELSE
379 JoiX(IJK,M) = ZERO
380 JoiY(IJK,M) = ZERO
381 JoiZ(IJK,M) = ZERO
382 ENDIF
383
384 CONTINUE
385 ENDDO
386
387 RETURN
388 END SUBROUTINE GHDMASSFLUX
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409 SUBROUTINE UpdateSpeciesVelocities ()
410
411
412
413
414 USE param
415 USE param1
416 USE geometry
417 USE compar
418 USE fldvar
419 USE indices
420 USE is
421 USE drag
422 USE visc_s
423 USE ghdtheory
424 USE physprop
425 USE run
426 USE constant
427 USE toleranc
428 USE fun_avg
429 USE functions
430 IMPLICIT NONE
431
432
433
434
435 INTEGER :: IJK, I, J, K
436
437 INTEGER :: IJKE, IJKN, IJKT, IJKW, IJKS, IJKB, IMJK, IJMK, &
438 IJKM
439
440 INTEGER :: S
441
442 integer :: kk, maxFluxS
443 double precision :: epgN, rogN, mugN, Vg
444 double precision :: Ur(smax), vrelSq(smax), vel, velup(smax)
445 double precision :: maxFlux
446 double precision :: rosN(smax), dp(smax)
447 double precision :: DijN(smax,smax), JoiM(smax), &
448 DijN_H(smax,smax), DijN_A(smax,smax)
449 double precision :: beta_cell(smax), beta_ij_cell(smax,smax)
450
451 integer :: ntrial
452 double precision tolx, tolf
453
454
455
456
457
458 = 100
459 tolx = 1d-14
460 tolf = 1d-14
461
462 DO 200 IJK = ijkstart3, ijkend3
463 I = I_OF(IJK)
464 J = J_OF(IJK)
465 K = K_OF(IJK)
466
467 IF ( FLUID_AT(IJK) ) THEN
468
469 IJKE = EAST_OF(IJK)
470 IJKW = WEST_OF(IJK)
471 IJKN = NORTH_OF(IJK)
472 IJKS = SOUTH_OF(IJK)
473 IJKT = TOP_OF(IJK)
474 IJKB = BOTTOM_OF(IJK)
475 IMJK = IM_OF(IJK)
476 IJMK = JM_OF(IJK)
477 IJKM = KM_OF(IJK)
478
479
480
481 IF(.NOT.IP_AT_E(IJK) .OR. .NOT.SIP_AT_E(IJK)) THEN
482 IF(RO_g0 == ZERO) THEN
483 do s = 1, smax
484 rosN(s) = AVG_X(ROP_S(IJK,s),ROP_S(IJKE,s),I)/ RO_S(IJK,s)
485 if(rosN(s) > zero_ep_s) then
486 u_s(ijk,s) = u_s(ijk,mmax) + JoiX(IJK,s)/(rosN(s)*ro_s(IJK,s))
487 else
488 u_s(ijk,s) = u_s(ijk,mmax)
489 endif
490 enddo
491 ELSE
492 Vg = U_g(ijk) - u_s(ijk,mmax)
493 epgN = AVG_X(EP_g(IJK),EP_g(IJKE),I)
494 rogN = AVG_X(ROP_g(IJK),ROP_g(IJKE),I)
495 mugN = AVG_X(MU_g(IJK),MU_g(IJKE),I)
496
497 do s = 1, smax
498 Ur(s) = u_g(ijk)-u_s(ijk,s)
499
500 (s) = (v_g(ijk)-v_s(ijk,s))**2 + (w_g(ijk)-w_s(ijk,s))**2
501 rosN(s) = AVG_X(ROP_S(IJK,s),ROP_S(IJKE,s),I)
502 velup(s) = 0.d0
503 beta_cell(s) = beta_cell_X(IJK,s)
504 dp(s) = D_P(IJK,s)
505
506 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
507 JoiM(s) = DELTAU(IJK,s)
508 ELSE
509 JoiM(s) = JoiMinusDragX(ijk,s)
510 ENDIF
511
512 do kk = 1, smax
513 DijN_H(s,kk) = AVG_X_S(DijF(IJK,s,kk),DijF(IJKE,s,kk),I)
514 DijN_A(s,kk) = AVG_X(DijF(IJK,s,kk),DijF(IJKE,s,kk),I)
515 if(DijF_HarmE(IJK,kk))THEN
516 DijN(s,kk) = DijN_H(s,kk)
517 ELSE
518 DijN(s,kk) = DijN_A(s,kk)
519 ENDIF
520 if(s .eq. kk)then
521 beta_ij_cell(s,kk)=0.d0
522 else
523 beta_ij_cell(s,kk)=beta_ij_cell_X(IJK,s,kk)
524 endif
525 enddo
526 enddo
527
528 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
529 vel=U_S(IJK,MMAX)
530 CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, &
531 JoiM, beta_cell, beta_ij_cell,vel)
532 ELSE
533 CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
534 epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
535 ENDIF
536
537
538 do s = 1, smax
539 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
540 U_S(IJK,s)=velup(s)
541 JoiX(IJK,s) = rosN(s) * (u_s(ijk,s)-u_s(ijk,mmax))
542 ELSE
543 u_s(ijk,s) = u_g(ijk) - Ur(s)
544 JoiX(IJK,s) = rosN(s) * (u_s(ijk,s)-u_s(ijk,mmax))
545 ENDIF
546 enddo
547 ENDIF
548 ENDIF
549
550 if(smax==2) then
551 (IJK,2)=-JoiX(IJK,1)
552 elseif(smax > 2) then
553 maxFlux = JoiX(IJK,1)
554 maxFluxS = 1
555 do s = 2, smax
556 if( abs(JoiX(IJK,s)) > abs(maxFlux) ) then
557 maxFlux = JoiX(IJK,s)
558 maxFluxS = s
559 endif
560 enddo
561 JoiX(IJK,maxFluxS) = 0d0
562 do s = 1, smax
563 if(s /= maxFluxS) JoiX(IJK,maxFluxS) = JoiX(IJK,maxFluxS) - JoiX(IJK,s)
564 enddo
565 endif
566
567
568
569
570
571
572
573 IF (.NOT.IP_AT_N(IJK) .OR. .NOT.SIP_AT_N(IJK)) THEN
574 IF(RO_g0 == ZERO) THEN
575 do s = 1, smax
576 rosN(s) = AVG_Y(ROP_S(IJK,s),ROP_S(IJKN,s),J)/ RO_S(IJK,s)
577 if(rosN(s) > zero_ep_s) then
578 v_s(ijk,s) = v_s(ijk,mmax) + JoiY(IJK,s)/(rosN(s)*ro_s(IJK,s))
579 else
580 v_s(ijk,s) = v_s(ijk,mmax)
581 endif
582 enddo
583 ELSE
584 Vg = V_g(ijk) - v_s(ijk,mmax)
585 epgN = AVG_Y(EP_g(IJK),EP_g(IJKN),J)
586 rogN = AVG_Y(ROP_g(IJK),ROP_g(IJKN),J)
587 mugN = AVG_Y(MU_g(IJK),MU_g(IJKN),J)
588
589 do s = 1, smax
590 Ur(s) = v_g(ijk)-v_s(ijk,s)
591
592 (s) = (u_g(ijk)-u_s(ijk,s))**2 + (w_g(ijk)-w_s(ijk,s))**2
593 rosN(s) = AVG_Y(ROP_S(IJK,s),ROP_S(IJKN,s),J)
594 velup(s) = 0.d0
595 beta_cell(s) = beta_cell_Y(IJK,s)
596 dp(s) = D_P(IJK,s)
597
598 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
599 JoiM(s) = DELTAV(IJK,s)
600 ELSE
601 JoiM(s) = JoiMinusDragY(ijk,s)
602 ENDIF
603
604 do kk = 1, smax
605
606 DijN_H(s,kk) = AVG_Y_S(DijF(IJK,s,kk),DijF(IJKN,s,kk),J)
607 DijN_A(s,kk) = AVG_Y(DijF(IJK,s,kk),DijF(IJKN,s,kk),J)
608
609 if(DijF_HarmN(IJK,kk))THEN
610 DijN(s,kk) = DijN_H(s,kk)
611 ELSE
612 DijN(s,kk) = DijN_A(s,kk)
613 ENDIF
614
615 if(s .eq. kk)then
616 beta_ij_cell(s,kk)=0.d0
617 else
618 beta_ij_cell(s,kk)=beta_ij_cell_Y(IJK,s,kk)
619 endif
620 enddo
621 enddo
622
623 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
624 vel=V_S(IJK,MMAX)
625 CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, JoiM, &
626 beta_cell, beta_ij_cell,vel)
627 ELSE
628 CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
629 epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
630 ENDIF
631
632
633 do s = 1, smax
634 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
635 V_S(IJK,s)=velup(s)
636 JoiY(IJK,s) = rosN(s) * (v_s(ijk,s)-v_s(ijk,mmax))
637 ELSE
638 v_s(ijk,s) = v_g(ijk) - Ur(s)
639 JoiY(IJK,s) = rosN(s) * (v_s(ijk,s)-v_s(ijk,mmax))
640 ENDIF
641 enddo
642 ENDIF
643 ENDIF
644
645 if(smax==2) then
646 (IJK,2)=-JoiY(IJK,1)
647 elseif(smax > 2) then
648 maxFlux = JoiY(IJK,1)
649 maxFluxS = 1
650 do s = 2, smax
651 if( abs(JoiY(IJK,s)) > abs(maxFlux) ) then
652 maxFlux = JoiY(IJK,s)
653 maxFluxS = s
654 endif
655 enddo
656 JoiY(IJK,maxFluxS) = 0d0
657 do s = 1, smax
658 if(s /= maxFluxS) JoiY(IJK,maxFluxS) = JoiY(IJK,maxFluxS) - JoiY(IJK,s)
659 enddo
660 endif
661
662
663
664
665
666
667 IF(.NOT.NO_K .AND. (.NOT.IP_AT_T(IJK) .OR. .NOT.SIP_AT_T(IJK))) THEN
668 IF(RO_g0 == ZERO) THEN
669 do s = 1, smax
670 rosN(s) = AVG_Z(ROP_S(IJK,s),ROP_S(IJKT,s),K)/ RO_S(IJK,s)
671 if(rosN(s) > zero_ep_s) then
672 w_s(ijk,s) = w_s(ijk,mmax) + JoiZ(IJK,s)/(rosN(s)*ro_s(IJK,s))
673 else
674 w_s(ijk,s) = w_s(ijk,mmax)
675 endif
676 enddo
677 ELSE
678 Vg = W_g(ijk) - W_s(ijk,mmax)
679 epgN = AVG_Z(EP_g(IJK),EP_g(IJKT),K)
680 rogN = AVG_Z(ROP_g(IJK),ROP_g(IJKT),K)
681 mugN = AVG_Z(MU_g(IJK),MU_g(IJKT),K)
682
683 do s = 1, smax
684 Ur(s) = w_g(ijk)-w_s(ijk,s)
685
686 (s) = (u_g(ijk)-u_s(ijk,s))**2 + (v_g(ijk)-v_s(ijk,s))**2
687 rosN(s) = AVG_Z(ROP_S(IJK,s),ROP_S(IJKT,s),K)
688 velup(s) = 0.d0
689 beta_cell(s) = beta_cell_Z(IJK,s)
690 dp(s) = D_P(IJK,s)
691
692 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
693 JoiM(s) = DELTAW(IJK,s)
694 ELSE
695 JoiM(s) = JoiMinusDragZ(ijk,s)
696 ENDIF
697
698 do kk = 1, smax
699 DijN_H(s,kk) = AVG_Z_S(DijF(IJK,s,kk),DijF(IJKT,s,kk),K)
700 DijN_A(s,kk) = AVG_Z(DijF(IJK,s,kk),DijF(IJKT,s,kk),K)
701 if(DijF_HarmT(IJK,kk))THEN
702 DijN(s,kk) = DijN_H(s,kk)
703 ELSE
704 DijN(s,kk) = DijN_A(s,kk)
705 ENDIF
706 if(s .eq. kk)then
707 beta_ij_cell(s,kk)=0.d0
708 else
709 beta_ij_cell(s,kk)=beta_ij_cell_Z(IJK,s,kk)
710 endif
711 enddo
712 enddo
713
714 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
715 vel=W_S(IJK,MMAX)
716 CALL VELOCITY_UPDATE(velup, smax, rosN, DijN, &
717 JoiM, beta_cell, beta_ij_cell,vel)
718 ELSE
719 CALL UrNEWT(ntrial, Ur, smax, ijk, tolx, tolf, &
720 epgN, rogN, mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
721 ENDIF
722
723
724 do s = 1, smax
725 IF(DRAG_TYPE_ENUM .eq. HYS)THEN
726 W_S(IJK,s)=velup(s)
727 JoiZ(IJK,s) = rosN(s) * (w_s(ijk,s)-w_s(ijk,mmax))
728 ELSE
729 w_s(ijk,s) = w_g(ijk) - Ur(s)
730 JoiZ(IJK,s) = rosN(s) * (w_s(ijk,s)-w_s(ijk,mmax))
731 ENDIF
732 enddo
733 ENDIF
734 ENDIF
735
736 if(smax==2) then
737 (IJK,2)=-JoiZ(IJK,1)
738 elseif(smax > 2) then
739 maxFlux = JoiZ(IJK,1)
740 maxFluxS = 1
741 do s = 2, smax
742 if( abs(JoiZ(IJK,s)) > abs(maxFlux) ) then
743 maxFlux = JoiZ(IJK,s)
744 maxFluxS = s
745 endif
746 enddo
747 JoiZ(IJK,maxFluxS) = 0d0
748 do s = 1, smax
749 if(s /= maxFluxS) JoiZ(IJK,maxFluxS) = JoiZ(IJK,maxFluxS) - JoiZ(IJK,s)
750 enddo
751 endif
752
753
754
755
756
757 ENDIF
758 CONTINUE
759
760 RETURN
761 END SUBROUTINE UpdateSpeciesVelocities
762
763
764
765
766
767
768
769
770
771 subroutine UrNEWT(ntrial, x, s, ijk, tolx, tolf, epgN, rogN, &
772 mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
773
774 Implicit NONE
775
776
777
778
779 integer, intent(in) :: s
780 integer, intent(in) :: ijk
781 integer, intent(in) :: ntrial
782 double precision, intent(in) :: tolx, tolf
783
784 DOUBLE PRECISION, intent(inout) :: X(s)
785 double precision, intent(in) :: epgN, rogN, mugN, Vg
786 double precision, intent(in) :: vrelSq(s), rosN(s), dp(s)
787 double precision, intent(in) :: DijN(s,s), JoiM(s)
788
789
790
791 INTEGER :: NP
792 PARAMETER (NP=15)
793
794
795
796 INTEGER :: I, K, INDX(s)
797 DOUBLE PRECISION :: D, ERRF, ERRX, FJAC(s,s), FVEC(s), P(s)
798
799
800 DO K = 1, NTRIAL
801 CALL Ur_JACOBI_EVAL(X, s, ijk, FVEC, FJAC, epgN, rogN, mugN, Vg, &
802 vrelSq, rosN, dp, DijN, JoiM)
803 ERRF = 0d0
804 DO I = 1, s
805 ERRF = ERRF + DABS(FVEC(I))
806 ENDDO
807
808 IF(ERRF <= TOLF) RETURN
809 DO I = 1, s
810 P(I) = -FVEC(I)
811 IF(rosN(I) .eq. 0.d0)THEN
812 P(I) = 0.d0
813 ENDIF
814 ENDDO
815
816 CALL LUDCMP(fjac, s, NP, indx, d, 'UrNewt')
817 CALL LUBKSB(fjac, s, NP, indx, p)
818
819 = 0d0
820 DO I = 1, s
821 ERRX = ERRX + DABS(P(I))
822 X(I) = X(I) + P(I)
823 ENDDO
824 IF(ERRX <= TOLX) RETURN
825 ENDDO
826
827 RETURN
828 END SUBROUTINE UrNEWT
829
830
831
832
833
834
835
836
837
838 SUBROUTINE Ur_JACOBI_EVAL(X, s, ijk, FVEC, FJAC, epgN, rogN, &
839 mugN, Vg, vrelSq, rosN, dp, DijN, JoiM)
840 Implicit NONE
841
842
843
844
845 integer, intent(in) :: s
846 integer, intent(in) :: ijk
847 double precision, intent(in) :: X(s)
848
849 DOUBLE PRECISION, INTENT(OUT) :: FVEC(s), FJAC(s,s)
850 DOUBLE PRECISION, INTENT(IN) :: epgN, rogN, mugN, Vg
851 DOUBLE PRECISION, INTENT(IN) :: vrelSq(s), rosN(s), dp(s)
852 DOUBLE PRECISION, INTENT(IN) :: DijN(s,s), JoiM(s)
853
854
855
856 double precision :: pi
857 parameter (pi=3.14159265458979323846d0)
858 double precision :: one
859 parameter (one=1.d0)
860 double precision :: zero
861 parameter (zero=0.d0)
862
863
864
865
866 INTEGER :: I, J
867
868 DOUBLE PRECISION :: RE_G, C_d, DgA, Vi, vrel
869 DOUBLE PRECISION :: FgsOni(s), dFgsdVi(s), sum(s)
870
871
872 DO i = 1, s
873 vrel = dsqrt(vrelSq(i) + x(i)**2)
874 Vi = pi * dp(i)**3 / 6d0
875 RE_G = dp(i)*vrel*rogN/mugN
876 IF(Re_G <= 1000D0 .and. Re_G> zero)THEN
877 C_d = (24.D0/Re_G) * (ONE + 0.15D0 * Re_G**0.687D0)
878 ELSE
879 C_d = 0.44D0
880 ENDIF
881 DgA = 0.75D0 * C_d * vrel * rogN / (epgN**2.65D0 * dp(i))
882 FgsOni(i) = DgA * Vi
883 IF(vrel == ZERO)THEN
884 dFgsdVi(i) = ZERO
885 ELSEIF(Re_G <= 1000D0)THEN
886 dFgsdVi(i) = 1.8549d0*mugN*Re_G**0.687d0*Vi*dabs(x(i))/ &
887 (dp(i)**2*epgN**2.65d0*vrel**2)
888 ELSE
889 dFgsdVi(i) = 0.33d0*rogN*Vi*dabs(x(i))/ &
890 (dp(i)*epgN**2.65d0*vrel)
891 ENDIF
892 ENDDO
893
894
895 DO i = 1, s
896 sum(i) = zero
897 DO j = 1, s
898 sum(i) = sum(i) + DijN(i,j) * FgsOni(j) * x(j)
899 ENDDO
900 ENDDO
901 DO i = 1, s
902 FVEC(i) = sum(i) - rosN(i)*x(i) + rosN(i)*Vg + JoiM(i)
903 ENDDO
904
905
906
907
908 DO i = 1, s
909 DO j = 1, s
910 if(j == i) THEN
911 FJAC(i,j) = DijN(i,i) * (FgsOni(i) + dFgsdVi(i)*x(i)) - rosN(i)
912 IF(rosN(i) .eq. 0.d0)THEN
913 FJAC(i,j) = 1.d0
914 ENDIF
915 else
916 FJAC(i,j) = DijN(i,j) * (FgsOni(j) + dFgsdVi(j)*x(j))
917 endif
918 ENDDO
919 ENDDO
920
921
922
923 RETURN
924 END SUBROUTINE Ur_JACOBI_EVAL
925
926
927
928
929
930
931
932
933
934 SUBROUTINE VELOCITY_UPDATE(FVEC, s, rosi, Diji, Joii, &
935 beta_celli, beta_ij_celli, &
936 velocity)
937 USE toleranc
938 USE fldvar
939 Implicit NONE
940
941
942
943
944 integer, intent (in) :: s
945
946 DOUBLE PRECISION, INTENT(IN) :: rosi(s), Diji(s,s), &
947 Joii(s), beta_celli(s), &
948 beta_ij_celli(s,s), velocity
949
950 DOUBLE PRECISION, INTENT(OUT) :: FVEC(s)
951
952
953
954 integer :: NP
955 PARAMETER (NP=15)
956
957
958
959
960 INTEGER :: i, k, l
961 INTEGER :: indx(s)
962 DOUBLE PRECISION :: D
963
964 DOUBLE PRECISION :: FJAC(s,s)
965
966
967 DO i=1,s
968 FVEC(i)= 0.d0
969 DO k=1,s
970 FJAC(i,k) = 0.d0
971 ENDDO
972 ENDDO
973
974 DO i=1,s
975
976 IF(rosi(i) .ne. 0.d0)THEN
977 FVEC(i) = rosi(i)*velocity+Joii(i)
978 ELSE
979 FVEC(i) = 0.d0
980 ENDIF
981
982 DO l=1,s
983
984 if(i .eq. l)then
985 FJAC(i,l)=FJAC(i,l)+rosi(i)
986 IF(rosi(i) .eq. 0.d0)THEN
987 (i,l) = 1.d0
988 ENDIF
989 endif
990
991 FJAC(i,l)=FJAC(i,l)-Diji(i,l)*beta_celli(l)
992
993 DO k=1,s
994
995 if(l .ne. k)then
996 FJAC(i,k)=FJAC(i,k)+Diji(i,l)*beta_ij_celli(l,k)
997 endif
998
999 enddo
1000
1001 enddo
1002 enddo
1003
1004 CALL LUDCMP(fjac, s, NP, indx, d, 'velocity_update')
1005 CALL LUBKSB(fjac, s, NP, indx, FVEC)
1006
1007 RETURN
1008
1009 END SUBROUTINE VELOCITY_UPDATE
1010
1011