File: /nfs/home/0/users/jenkins/mfix.git/model/kintheory_w_s.f
1
2
3
4
5
6
7
8
9
10
11 SUBROUTINE CALC_KTMOMSOURCE_W_S()
12
13
14
15
16
17 USE param1, only: zero
18
19
20 USE run, only: kt_type_enum
21 USE run, only: ia_2005
22
23
24 USE physprop, only: smax
25
26
27 USE kintheory, only: ktmom_w_s
28
29 IMPLICIT NONE
30
31
32
33
34 INTEGER :: M
35
36
37 DO M = 1, SMAX
38 KTMOM_W_s(:,M) = ZERO
39 IF (KT_TYPE_ENUM == IA_2005) THEN
40 CALL CALC_IA_MOMSOURCE_W_S (M)
41 ENDIF
42 ENDDO
43
44 RETURN
45 END SUBROUTINE CALC_KTMOMSOURCE_W_S
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65 SUBROUTINE CALC_IA_MOMSOURCE_W_S(M)
66
67
68
69
70 USE param1, only: half, zero
71 USE constant, only: pi
72
73
74 USE visc_s, only: trD_s
75
76
77 USE physprop, only: mmax
78
79
80 USE fldvar, only: u_s, v_s, w_s
81
82 USE fldvar, only: d_p, rop_s, ro_s
83
84 USE fldvar, only: theta_m, ep_s
85
86 USE toleranc, only: dil_ep_s
87
88 Use kintheory
89
90 USE geometry
91 USE indices
92 USE compar
93 USE fun_avg
94 USE functions
95
96 IMPLICIT NONE
97
98
99
100
101 INTEGER, INTENT(IN) :: M
102
103
104
105
106 DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
107
108 INTEGER :: IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
109 IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
110 IMJKP, K, IJKT, JM, KP, IJKM, IPJK,&
111 IJPK
112
113 INTEGER :: L
114
115 DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
116 MU_sL_pB, MU_sL_p
117
118 DOUBLE PRECISION :: XI_sL_pT, XI_sL_pB, LAMBDA_sL_pT, LAMBDA_sL_pB
119
120 DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pT, NU_PM_pB, NU_PM_p, &
121 NU_PL_pT, NU_PL_pB, NU_PL_p, T_PM_pT, T_PM_pB,&
122 T_PL_pT, T_PL_pB, Fnu_s_p, FT_sM_p, FT_sL_p
123
124 DOUBLE PRECISION :: duodz
125
126 DOUBLE PRECISION :: EPSA
127
128 DOUBLE PRECISION :: sswx, sswy, sswz, ssx, ssy, ssz, ssbv, tauxz_x
129
130 DOUBLE PRECISION :: tauxz_ox, DS1, DS2, DS3, DS4, DS1plusDS2
131
132
133
134
135 DO IJK = IJKSTART3, IJKEND3
136
137
138 IF(WALL_AT(IJK)) cycle
139
140 D_PM = D_P(IJK,M)
141 M_PM = (Pi/6d0)*(D_PM**3.)*RO_S(IJK,M)
142
143 K = K_OF(IJK)
144 IJKT = TOP_OF(IJK)
145 EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
146 IF ( .NOT.SIP_AT_T(IJK) .AND. EPSA>DIL_EP_S) THEN
147
148 J = J_OF(IJK)
149 I = I_OF(IJK)
150 KP = KP1(K)
151 IM = IM1(I)
152 JM = JM1(J)
153 IMJK = IM_OF(IJK)
154 IJMK = JM_OF(IJK)
155 IJKM = KM_OF(IJK)
156 IJKP = KP_OF(IJK)
157 IJMKP = JM_OF(IJKP)
158 IMJKP = KP_OF(IMJK)
159
160 IJKN = NORTH_OF(IJK)
161 IJKS = SOUTH_OF(IJK)
162 IJKNT = TOP_OF(IJKN)
163 IJKST = TOP_OF(IJKS)
164
165 IJKE = EAST_OF(IJK)
166 IJKW = WEST_OF(IJK)
167 IJKTE = EAST_OF(IJKT)
168 IJKTW = WEST_OF(IJKT)
169
170
171 = IP_OF(IJK)
172 IJPK = JP_OF(IJK)
173
174
175 = ZERO
176 DRAG_TERMS = ZERO
177
178 DO L = 1, MMAX
179 IF (M .ne. L) THEN
180
181
182
183
184 = AVG_Z_H(AVG_X_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I),&
185 AVG_X_H(MU_sL_ip(IJKT,M,L),MU_sL_ip(IJKTE,M,L),I),K)
186 MU_sL_pW = AVG_Z_H(AVG_X_H(MU_sL_ip(IJKW,M,L),MU_sL_ip(IJK,M,L),IM),&
187 AVG_X_H(MU_sL_ip(IJKTW,M,L),MU_sL_ip(IJKT,M,L),IM),K)
188 SSWX = MU_sL_pE*(W_S(IPJK,L)-W_S(IJK,L))*AYZ_W(IJK)*ODX_E(I)&
189 -MU_sL_PW*(W_S(IJK,L)-W_S(IMJK,L))*AYZ_W(IMJK)*ODX_E(IM)
190
191 MU_sL_pN = AVG_Z_H(AVG_Y_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKN,M,L), J),&
192 AVG_Y_H(MU_sL_ip(IJKT,M,L),MU_sL_ip(IJKNT,M,L), J), K)
193 MU_sL_pS = AVG_Z_H(AVG_Y_H(MU_sL_ip(IJKS,M,L),MU_sL_ip(IJK,M,L), JM),&
194 AVG_Y_H(MU_sL_ip(IJKST,M,L),MU_sL_ip(IJKT,M,L), JM), K)
195 SSWY = MU_sL_pN*(W_S(IJPK,L)-W_S(IJK,L))*AXZ_W(IJK)*ODY_N(J)&
196 -MU_sL_pS*(W_S(IJK,L)-W_S(IJMK,L))*AXZ_W(IJKM)*ODY_N(JM)
197
198 MU_sL_pT = MU_sL_ip(IJKT,M,L)
199 MU_sL_pB = MU_sL_ip(IJK,M,L)
200 SSWZ = MU_sL_pT*(W_S(IJKP,L)-W_S(IJK,L))*AXY_W(IJK)*ODZ(KP)*OX(I)&
201 -MU_sL_pB*(W_S(IJK,L)-W_S(IJKM,L))*AXY_W(IJKM)*ODZ(K)*OX(I)
202
203
204 = XI_sL_ip(IJKT,M,L)
205 XI_sL_pB = XI_sL_ip(IJK,M,L)
206 LAMBDA_sL_pT = -(2d0/3d0)*MU_sL_pT + XI_sL_pT
207 LAMBDA_sL_pB = -(2d0/3d0)*MU_sL_pB + XI_sL_pB
208 SSBV = (LAMBDA_sL_pT*TRD_S(IJKT,L)-LAMBDA_sL_pB*TRD_S(IJK,L))*AXY(IJK)
209
210
211 = MU_sL_pE*(U_S(IJKP,L)-U_S(IJK,L))*OX_E(I)*AYZ_W(IJK)*ODZ_T(K)&
212 -MU_sL_pW*(U_S(IMJKP,L)-U_S(IMJK,L))*(DY(J)*HALF*(DZ(K)+&
213 DZ(KP)))*ODZ_T(K)
214
215
216 = MU_sL_pN*(V_S(IJKP,L)-V_S(IJK,L))*AXZ_W(IJK)*ODZ_T(K)*OX(I)&
217 -MU_sL_pS*(V_S(IJMKP,L)-V_S(IJMK,L))*AXZ_W(IJMK)*ODZ_T(K)*OX(I)
218
219 SSZ = SSWZ
220
221
222 IF (CYLINDRICAL) THEN
223
224
225
226 = SSZ +(&
227 (MU_sL_pT*(U_S(IJKP,L)+U_S(IMJKP,L))*OX(I)*AXY_W(IJK))-&
228 (MU_sL_pB*(U_S(IJK,L)+U_S(IMJK,L))*OX(I)*AXY_W(IJKM)))
229
230 MU_sL_p = AVG_Z(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K)
231
232
233
234 IF (OX_E(IM) /= UNDEFINED) THEN
235 DUODZ = (U_S(IMJKP,L)-U_S(IMJK,L))*OX_E(IM)*ODZ_T(K)
236 ELSE
237 DUODZ = ZERO
238 ENDIF
239
240 tauxz_ox = MU_sL_p*OX(I)*HALF*(&
241 (OX_E(I)*(U_S(IJKP,L)-U_S(IJK,L))*ODZ_T(K))+&
242 DUODZ)
243
244
245
246 = tauxz_ox + (MU_sL_p*OX(I)*HALF*&
247 ( (W_S(IPJK,L)-W_S(IJK,L))*ODX_E(I) + &
248 (W_S(IJK,L)-W_S(IMJK,L))*ODX_E(IM) ) )
249
250
251
252 = tauxz_ox - (OX(I)*OX(I)*MU_sL_p*W_S(IJK,L))
253
254
255 = tauxz_ox*VOL_W(IJK)
256
257
258
259
260 = -( ( MU_sL_pE*OX_E(I)*HALF*(W_S(IPJK,L)+&
261 W_S(IJK,L))*AYZ_W(IJK) )-( MU_sL_pW*HALF*&
262 (W_S(IJK,L)+W_S(IMJK,L))*&
263 (DY(J)*HALF*(DZ(K)+DZ(KP)) ) ) )
264
265 ELSE
266 tauxz_ox = ZERO
267 tauxz_x = ZERO
268 ENDIF
269
270
271
272
273
274
275 = D_P(IJK,L)
276 M_PL = (Pi/6d0)*(D_PL**3.)*RO_S(IJK,L)
277
278 NU_PM_pT = ROP_S(IJKT,M)/M_PM
279 NU_PM_pB = ROP_S(IJK,M)/M_PM
280 NU_PM_p = AVG_Z(NU_PM_pB,NU_PM_pT,K)
281
282 NU_PL_pT = ROP_S(IJKT,L)/M_PL
283 NU_PL_pB = ROP_S(IJK,L)/M_PL
284 NU_PL_p = AVG_Z(NU_PL_pB,NU_PL_pT,K)
285
286 Fnu_s_p = AVG_Z(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKT,M,L),K)
287 DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pT-NU_PM_pB)*OX(I)*ODZ_T(K)
288 DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pT-NU_PL_pB)*OX(I)*ODZ_T(K)
289 DS1plusDS2 = DS1 + DS2
290
291
292
293 = Theta_M(IJKT,M)
294 T_PM_pB = Theta_M(IJK,M)
295
296 FT_sM_p = AVG_Z(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKT,M,L),K)
297 DS3 = FT_sM_p*(T_PM_pT-T_PM_pB)*OX(I)*ODZ_T(K)
298
299
300
301 = Theta_M(IJKT,L)
302 T_PL_pB = Theta_M(IJK,L)
303
304 FT_sL_p = AVG_Z(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKT,M,L),K)
305 DS4 = FT_sL_p*(T_PL_pT-T_PL_pB)*OX(I)*ODZ_T(K)
306
307
308
309
310 = STRESS_TERMS + SSWX + SSWY + SSWZ + &
311 SSBV + SSX + SSY + SSZ + tauxz_ox + tauxz_x
312 DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_W(IJK)
313
314 ELSE
315
316
317 = STRESS_TERMS + ZERO
318 DRAG_TERMS = DRAG_TERMS + ZERO
319 ENDIF
320 ENDDO
321
322 (IJK,M) = STRESS_TERMS + DRAG_TERMS
323 ELSE
324 KTMOM_W_S(IJK,M) = ZERO
325 ENDIF
326 ENDDO
327
328 RETURN
329 END SUBROUTINE CALC_IA_MOMSOURCE_W_S
330