File: /nfs/home/0/users/jenkins/mfix.git/model/cartesian_grid/calc_vort_out.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14 SUBROUTINE CALC_VORTICITY
15
16 USE param
17 USE param1
18 USE parallel
19 USE constant
20 USE run
21 USE toleranc
22 USE geometry
23 USE indices
24 USE compar
25 USE sendrecv
26 USE fldvar
27 USE quadric
28 USE cutcell
29 USE functions
30
31 IMPLICIT NONE
32 INTEGER :: I,J,K,IJK,IP,IM,JP,JM,KP,KM,IJKE,IJKN,IJKT,IJKW,IJKS,IJKB
33 DOUBLE PRECISION :: DU_DX,DU_DY,DU_DZ,DV_DX,DV_DY,DV_DZ,DW_DX,DW_DY,DW_DZ
34 DOUBLE PRECISION :: OMEGA_X,OMEGA_Y,OMEGA_Z
35 DOUBLE PRECISION :: LAMBDA_1,LAMBDA_2,LAMBDA_3
36 DOUBLE PRECISION,DIMENSION(3,3) :: OMEGA,SS,AA
37 DOUBLE PRECISION,DIMENSION(4) :: POLY
38
39 DO IJK = IJKSTART3, IJKEND3
40
41 I = I_OF(IJK)
42 J = J_OF(IJK)
43 K = K_OF(IJK)
44
45 IM = I - 1
46 IP = I + 1
47
48 JM = J - 1
49 JP = J + 1
50
51 KM = K - 1
52 KP = K + 1
53
54
55 IJKE = FUNIJK(IP,J,K)
56 IJKN = FUNIJK(I,JP,K)
57 IJKT = FUNIJK(I,J,KP)
58
59 IJKW = FUNIJK(IM,J,K)
60 IJKS = FUNIJK(I,JM,K)
61 IJKB = FUNIJK(I,J,KM)
62
63
64
65 IF (FLUID_AT(IJK).AND.INTERIOR_CELL_AT(IJK)) THEN
66
67 DU_DX = ZERO
68 DU_DY = ZERO
69 DU_DZ = ZERO
70
71 DV_DX = ZERO
72 DV_DY = ZERO
73 DV_DZ = ZERO
74
75 DW_DX = ZERO
76 DW_DY = ZERO
77 DW_DZ = ZERO
78
79
80
81
82
83
84
85
86
87 IF ((FLUID_AT(IJKE)).AND.(FLUID_AT(IJKW))) THEN
88
89 DU_DX = (U_g(IJKE) - U_g(IJKW)) / (X_U(IJKE) - X_U(IJKW))
90
91 ELSE IF ((FLUID_AT(IJKE)).AND.(.NOT.FLUID_AT(IJKW))) THEN
92
93 DU_DX = (U_g(IJKE) - U_g(IJK)) / (X_U(IJKE) - X_U(IJK))
94
95 ELSE IF ((FLUID_AT(IJKW)).AND.(.NOT.FLUID_AT(IJKE))) THEN
96
97 DU_DX = (U_g(IJK) - U_g(IJKW)) / (X_U(IJK) - X_U(IJKW))
98
99 END IF
100
101
102
103
104
105 IF ((FLUID_AT(IJKN)).AND.(FLUID_AT(IJKS))) THEN
106
107 DU_DY = (U_g(IJKN) - U_g(IJKS)) / (Y_U(IJKN) - Y_U(IJKS))
108
109 ELSE IF ((FLUID_AT(IJKN)).AND.(.NOT.FLUID_AT(IJKS))) THEN
110
111 DU_DY = (U_g(IJKN) - U_g(IJK)) / (Y_U(IJKN) - Y_U(IJK))
112
113 ELSE IF ((FLUID_AT(IJKS)).AND.(.NOT.FLUID_AT(IJKN))) THEN
114
115 DU_DY = (U_g(IJK) - U_g(IJKS)) / (Y_U(IJK) - Y_U(IJKS))
116
117 END IF
118
119
120
121
122
123 IF(DO_K) THEN
124 IF ((FLUID_AT(IJKT)).AND.(FLUID_AT(IJKB))) THEN
125
126 DU_DZ = (U_g(IJKT) - U_g(IJKB)) / (Z_U(IJKT) - Z_U(IJKB))
127
128 ELSE IF ((FLUID_AT(IJKT)).AND.(.NOT.FLUID_AT(IJKB))) THEN
129
130 DU_DZ = (U_g(IJKT) - U_g(IJK)) / (Z_U(IJKT) - Z_U(IJK))
131
132 ELSE IF ((FLUID_AT(IJKB)).AND.(.NOT.FLUID_AT(IJKT))) THEN
133
134 DU_DZ = (U_g(IJK) - U_g(IJKB)) / (Z_U(IJK) - Z_U(IJKB))
135
136 END IF
137 ENDIF
138
139
140
141
142
143 IF ((FLUID_AT(IJKE)).AND.(FLUID_AT(IJKW))) THEN
144
145 DV_DX = (V_g(IJKE) - V_g(IJKW)) / (X_V(IJKE) - X_V(IJKW))
146
147 ELSE IF ((FLUID_AT(IJKE)).AND.(.NOT.FLUID_AT(IJKW))) THEN
148
149 DV_DX = (V_g(IJKE) - V_g(IJK)) / (X_V(IJKE) - X_V(IJK))
150
151 ELSE IF ((FLUID_AT(IJKW)).AND.(.NOT.FLUID_AT(IJKE))) THEN
152
153 DV_DX = (V_g(IJK) - V_g(IJKW)) / (X_V(IJK) - X_V(IJKW))
154
155 END IF
156
157
158
159
160
161 IF ((FLUID_AT(IJKN)).AND.(FLUID_AT(IJKS))) THEN
162
163 DV_DY = (V_g(IJKN) - V_g(IJKS)) / (Y_V(IJKN) - Y_V(IJKS))
164
165 ELSE IF ((FLUID_AT(IJKN)).AND.(.NOT.FLUID_AT(IJKS))) THEN
166
167 DV_DY = (V_g(IJKN) - V_g(IJK)) / (Y_V(IJKN) - Y_V(IJK))
168
169 ELSE IF ((FLUID_AT(IJKS)).AND.(.NOT.FLUID_AT(IJKN))) THEN
170
171 DV_DY = (V_g(IJK) - V_g(IJKS)) / (Y_V(IJK) - Y_V(IJKS))
172
173 END IF
174
175
176
177
178
179 IF(DO_K) THEN
180 IF ((FLUID_AT(IJKT)).AND.(FLUID_AT(IJKB))) THEN
181
182 DV_DZ = (V_g(IJKT) - V_g(IJKB)) / (Z_V(IJKT) - Z_V(IJKB))
183
184 ELSE IF ((FLUID_AT(IJKT)).AND.(.NOT.FLUID_AT(IJKB))) THEN
185
186 DV_DZ = (V_g(IJKT) - V_g(IJK)) / (Z_V(IJKT) - Z_V(IJK))
187
188 ELSE IF ((FLUID_AT(IJKB)).AND.(.NOT.FLUID_AT(IJKT))) THEN
189
190 DV_DZ = (V_g(IJK) - V_g(IJKB)) / (Z_V(IJK) - Z_V(IJKB))
191
192 END IF
193 ENDIF
194
195
196
197
198
199 IF(DO_K) THEN
200 IF ((FLUID_AT(IJKE)).AND.(FLUID_AT(IJKW))) THEN
201
202 DW_DX = (W_g(IJKE) - W_g(IJKW)) / (X_W(IJKE) - X_W(IJKW))
203
204 ELSE IF ((FLUID_AT(IJKE)).AND.(.NOT.FLUID_AT(IJKW))) THEN
205
206 DW_DX = (W_g(IJKE) - W_g(IJK)) / (X_W(IJKE) - X_W(IJK))
207
208 ELSE IF ((FLUID_AT(IJKW)).AND.(.NOT.FLUID_AT(IJKE))) THEN
209
210 DW_DX = (W_g(IJK) - W_g(IJKW)) / (X_W(IJK) - X_W(IJKW))
211
212 END IF
213
214
215
216
217
218 IF ((FLUID_AT(IJKN)).AND.(FLUID_AT(IJKS))) THEN
219
220 DW_DY = (W_g(IJKN) - W_g(IJKS)) / (Y_W(IJKN) - Y_W(IJKS))
221
222 ELSE IF ((FLUID_AT(IJKN)).AND.(.NOT.FLUID_AT(IJKS))) THEN
223
224 DW_DY = (W_g(IJKN) - W_g(IJK)) / (Y_W(IJKN) - Y_W(IJK))
225
226 ELSE IF ((FLUID_AT(IJKS)).AND.(.NOT.FLUID_AT(IJKN))) THEN
227
228 DW_DY = (W_g(IJK) - W_g(IJKS)) / (Y_W(IJK) - Y_W(IJKS))
229
230 END IF
231
232
233
234
235
236 IF ((FLUID_AT(IJKT)).AND.(FLUID_AT(IJKB))) THEN
237
238 DW_DZ = (W_g(IJKT) - W_g(IJKB)) / (Z_W(IJKT) - Z_W(IJKB))
239
240 ELSE IF ((FLUID_AT(IJKT)).AND.(.NOT.FLUID_AT(IJKB))) THEN
241
242 DW_DZ = (W_g(IJKT) - W_g(IJK)) / (Z_W(IJKT) - Z_W(IJK))
243
244 ELSE IF ((FLUID_AT(IJKB)).AND.(.NOT.FLUID_AT(IJKT))) THEN
245
246 DW_DZ = (W_g(IJK) - W_g(IJKB)) / (Z_W(IJK) - Z_W(IJKB))
247
248 END IF
249 ENDIF
250
251
252
253
254
255 = DW_DY - DV_DZ
256 OMEGA_Y = DU_DZ - DW_DX
257 OMEGA_Z = DV_DX - DU_DY
258
259
260 VORTICITY(IJK) = DSQRT(OMEGA_X**2 + OMEGA_Y**2 + OMEGA_Z**2)
261
262
263
264
265
266
267
268
269
270 (1,1) = DU_DX
271 OMEGA(1,2) = HALF * (DU_DY + DV_DX)
272 OMEGA(1,3) = HALF * (DU_DZ + DW_DX)
273
274 OMEGA(2,1) = OMEGA(1,2)
275 OMEGA(2,2) = DV_DY
276 OMEGA(2,3) = HALF * (DV_DZ + DW_DY)
277
278 OMEGA(3,1) = OMEGA(1,3)
279 OMEGA(3,2) = OMEGA(2,3)
280 OMEGA(3,3) = DW_DZ
281
282
283 SS(1,1) = ZERO
284 SS(1,2) = HALF * (DU_DY - DV_DX)
285 SS(1,3) = HALF * (DU_DZ - DW_DX)
286
287 SS(2,1) = - SS(1,2)
288 SS(2,2) = ZERO
289 SS(2,3) = HALF * (DV_DZ - DW_DY)
290
291 SS(3,1) = - SS(1,3)
292 SS(3,2) = - SS(2,3)
293 SS(3,3) = ZERO
294
295
296 AA = MATMUL(OMEGA,OMEGA) + MATMUL(SS,SS)
297
298
299
300
301
302 (1) = - ONE
303 POLY(2) = AA(1,1) + AA(2,2) + AA(3,3)
304 POLY(3) = AA(2,1)*AA(1,2) + AA(3,1)*AA(1,3) + AA(3,2)*AA(2,3) &
305 -( AA(1,1)*AA(2,2) + AA(1,1)*AA(3,3) + AA(2,2)*AA(3,3) )
306 POLY(4) = AA(1,1)*AA(2,2)*AA(3,3) + AA(1,2)*AA(2,3)*AA(3,1) + AA(2,1)*AA(3,2)*AA(1,3) &
307 -( AA(3,1)*AA(2,2)*AA(1,3) + AA(1,2)*AA(2,1)*AA(3,3) + AA(2,3)*AA(3,2)*AA(1,1) )
308
309
310
311
312
313
314
315
316 CALL BAIRSTOW(POLY,LAMBDA_1,LAMBDA_2,LAMBDA_3)
317
318 LAMBDA2(IJK) = LAMBDA_2
319
320 ENDIF
321
322 END DO
323
324 RETURN
325
326
327 END SUBROUTINE CALC_VORTICITY
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344 SUBROUTINE BAIRSTOW(A,X1,X2,X3)
345 USE param1
346
347 IMPLICIT NONE
348 INTEGER :: N,NMAX,ITERMAX,J,ITER
349
350 PARAMETER(NMAX = 4)
351
352 DOUBLE PRECISION, DIMENSION(NMAX):: A,B,C
353 DOUBLE PRECISION:: R,S,DELTA,DELTAR,DELTAS,DENOM
354 DOUBLE PRECISION:: EPS
355 DOUBLE PRECISION:: X1,X2,X3
356 DOUBLE PRECISION:: BUFFER1
357
358 LOGICAL :: TEST1, TEST2, TEST3
359
360 N = 3
361 = 200
362 = 1.0D-6
363
364
365
366
367
368
369
370 = -0.9D0
371 S = -1.0D0
372
373 ITER = 0
374
375
376
377
378
379 CONTINUE
380
381 B(1) = A(1)
382 B(2) = A(2) + R * B(1)
383 C(1) = B(1)
384 C(2) = B(2) + R * C(1)
385 B(3) = A(3) + R * B(2) + S * B(1)
386 C(3) = B(3) + R * C(2) + S * C(1)
387 B(4) = A(4) + R * B(3) + S * B(2)
388
389
390
391
392
393 = C(N-1) * C(N-1) - C(N)*C(N-2)
394
395 DELTAR = (-B(N)*C(N-1) + B(N+1)*C(N-2) ) / DENOM
396 DELTAS = (-C(N-1)*B(N+1) + B(N)*C(N) ) / DENOM
397
398 R = R + DELTAR
399 S = S + DELTAS
400
401 ITER = ITER +1
402
403
404
405
406
407 = (ABS( B(N) ) > EPS)
408 TEST2 = (ABS( B(N+1)) > EPS)
409 TEST3 = (ITER <= ITERMAX)
410
411 IF(TEST1.AND.TEST2.AND.TEST3) GOTO 20
412
413 IF (.NOT.TEST3) THEN
414
415
416
417 = UNDEFINED
418 X2 = UNDEFINED
419 X3 = UNDEFINED
420 ENDIF
421
422
423
424
425
426
427
428 = R*R + 4.0D0 * S
429
430 IF (DELTA.LT.0.0D0) THEN
431
432
433
434 = UNDEFINED
435 X2 = UNDEFINED
436
437 ELSE
438
439 X1 = ( R + DSQRT(DELTA)) / ( 2.0D0 )
440 X2 = ( R - DSQRT(DELTA)) / ( 2.0D0 )
441
442 ENDIF
443
444
445
446
447
448
449
450 = N - 2
451
452 DO J=1,N+1
453 A(J) = B(J)
454 END DO
455
456 IF(DABS(A(1))<1.0D-9) THEN
457 X3 = UNDEFINED
458 ELSE
459 X3 = -A(2)/A(1)
460 ENDIF
461
462
463
464
465
466 IF(X1>X2) THEN
467 BUFFER1 = X1
468 X1 = X2
469 X2 = BUFFER1
470 ENDIF
471
472 IF(X2>X3) THEN
473 BUFFER1 = X2
474 X2 = X3
475 X3 = BUFFER1
476 ENDIF
477
478 IF(X1>X2) THEN
479 BUFFER1 = X1
480 X1 = X2
481 X2 = BUFFER1
482 ENDIF
483
484 RETURN
485
486 END SUBROUTINE BAIRSTOW
487
488
489