File: /nfs/home/0/users/jenkins/mfix.git/model/kintheory_u_s.f
1
2
3
4
5
6
7
8
9
10
11 SUBROUTINE CALC_KTMOMSOURCE_U_S(IER)
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_u_s
28
29 IMPLICIT NONE
30
31
32
33
34 INTEGER, INTENT(INOUT) :: IER
35
36
37
38
39 INTEGER :: M
40
41
42 DO M = 1, SMAX
43 KTMOM_U_s(:,M) = ZERO
44 IF (KT_TYPE_ENUM == IA_2005) THEN
45 CALL CALC_IA_MOMSOURCE_U_S (M)
46 ENDIF
47 ENDDO
48
49 RETURN
50 END SUBROUTINE CALC_KTMOMSOURCE_U_S
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70 SUBROUTINE CALC_IA_MOMSOURCE_U_S(M)
71
72
73
74
75 USE param1, only: half, zero
76 USE constant, only: pi
77
78
79 USE visc_s, only: trD_s
80
81
82 USE physprop, only: mmax
83
84
85 USE fldvar, only: u_s, v_s, w_s
86
87 USE fldvar, only: d_p, rop_s, ro_s
88
89 USE fldvar, only: theta_m, ep_s
90
91 USE toleranc, only: dil_ep_s
92
93 Use kintheory
94
95 USE geometry
96 USE indices
97 USE compar
98 USE fun_avg
99 USE functions
100
101 IMPLICIT NONE
102
103
104
105
106 INTEGER, INTENT(IN) :: M
107
108
109
110
111 DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
112
113 INTEGER :: I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
114 IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
115 IJKB, IJKBE, IJKM, IPJKM, &
116 IJPK, IJKP, IM, IJKW
117
118 INTEGER :: L
119
120 DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
121 MU_sL_pB, Mu_sL_p
122
123 DOUBLE PRECISION :: XI_sL_pE, XI_sL_pW, LAMBDA_sL_pE, LAMBDA_sL_pW
124
125 DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pE, NU_PM_pW, NU_PM_p, &
126 NU_PL_pE, NU_PL_pW, NU_PL_p, T_PM_pE, T_PM_pW, &
127 T_PL_pE, T_PL_pW, Fnu_s_p, FT_sM_p, FT_sL_p
128
129 DOUBLE PRECISION :: EPSA
130
131 DOUBLE PRECISION :: ssux, ssuy, ssuz, ssx, ssy, ssz, ssbv
132
133 DOUBLE PRECISION :: tauzz, DS1, DS2, DS3, DS4, DS1plusDS2
134
135
136
137
138 DO IJK = IJKSTART3, IJKEND3
139
140
141 IF(WALL_AT(IJK)) cycle
142
143 D_PM = D_P(IJK,M)
144 M_PM = (Pi/6d0)*D_PM**3 * RO_S(IJK,M)
145
146 I = I_OF(IJK)
147 IJKE = EAST_OF(IJK)
148 EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
149 IF ( .NOT.SIP_AT_E(IJK) .AND. EPSA>DIL_EP_S) THEN
150
151 IP = IP1(I)
152 IM = Im1(I)
153 IJKW = WEST_OF(IJK)
154 J = J_OF(IJK)
155 JM = JM1(J)
156 K = K_OF(IJK)
157 KM = KM1(K)
158 IPJK = IP_OF(IJK)
159 IMJK = IM_OF(IJK)
160 IJMK = JM_OF(IJK)
161 IJKM = KM_OF(IJK)
162 IPJMK = JM_OF(IPJK)
163 IPJKM = IP_OF(IJKM)
164
165 IJKN = NORTH_OF(IJK)
166 IJKS = SOUTH_OF(IJK)
167 IJKNE = EAST_OF(IJKN)
168 IJKSE = EAST_OF(IJKS)
169
170 IJKT = TOP_OF(IJK)
171 IJKB = BOTTOM_OF(IJK)
172 IJKTE = EAST_OF(IJKT)
173 IJKBE = EAST_OF(IJKB)
174
175
176 = JP_OF(IJK)
177 IJKP = KP_OF(IJK)
178
179
180 = ZERO
181 DRAG_TERMS = ZERO
182
183 DO L = 1,MMAX
184 IF (M .ne. L) THEN
185
186
187
188
189 = MU_sL_ip(IJKE,M,L)
190 MU_sL_pW = MU_sL_ip(IJK,M,L)
191 SSUX = MU_sL_pE*(U_S(IPJK,L)-U_S(IJK,L))*AYZ_U(IJK)*ODX(IP)&
192 -MU_sL_PW*(U_S(IJK,L)-U_S(IMJK,L))*AYZ_U(IMJK)*ODX(I)
193
194 MU_sL_pN = AVG_X_H(AVG_Y_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKN,M,L), J),&
195 AVG_Y_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKNE,M,L), J), I)
196 MU_sL_pS = AVG_X_H(AVG_Y_H(MU_sL_ip(IJKS,M,L),MU_sL_ip(IJK,M,L), JM),&
197 AVG_Y_H(MU_sL_ip(IJKSE,M,L),MU_sL_ip(IJKE,M,L), JM), I)
198 SSUY = MU_sL_pN*(U_S(IJPK,L)-U_S(IJK,L))*AXZ_U(IJK)*ODY_N(J)&
199 -MU_sL_pS*(U_S(IJK,L)-U_S(IJMK,L))*AXZ_U(IJMK)*ODY_N(JM)
200
201 MU_sL_pT = AVG_X_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
202 AVG_Z_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKTE,M,L),K),I)
203 MU_sL_pB = AVG_X_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
204 AVG_Z_H(MU_sL_ip(IJKBE,M,L),MU_sL_ip(IJKE,M,L),KM),I)
205 SSUZ = MU_sL_pT*(U_S(IJKP,L)-U_S(IJK,L))*AXY_U(IJK)*ODZ_T(K)*OX_E(I)&
206 -MU_sL_pB*(U_S(IJK,L)-U_S(IJKM,L))*AXY_U(IJKM)*ODZ_T(KM)*OX_E(I)
207
208
209 = XI_sL_ip(IJKE,M,L)
210 XI_sL_pW = XI_sL_ip(IJK,M,L)
211 LAMBDA_sL_pE = -(2.d0/3.d0)*Mu_sL_pE + XI_sL_pE
212 LAMBDA_sL_pW = -(2.d0/3.d0)*Mu_sL_pW + XI_sL_pW
213 SSBV = (LAMBDA_sL_pE*TRD_S(IJKE,L)-LAMBDA_sL_pW*TRD_S(IJK,L))*AYZ(IJK)
214
215
216 = SSUX
217 SSY = MU_sL_pN*(V_S(IPJK,L)-V_S(IJK,L))*AXZ_U(IJK)*ODX_E(I)&
218 -MU_sL_pS*(V_S(IPJMK,L)-V_S(IJMK,L))*AXZ_U(IJMK)*ODX_E(I)
219 SSZ = MU_sL_pT*(W_S(IPJK,L)-W_S(IJK,L))*AXY_U(IJK)*ODX_E(I)&
220 -MU_sL_pB*(W_S(IPJKM,L)-W_S(IJKM,L))*AXY_U(IJKM)*ODX_E(I)
221
222
223 IF (CYLINDRICAL) THEN
224
225
226
227
228 = SSZ - (MU_sL_pT*(W_S(IPJK,L)+W_S(IJK,L))*HALF*OX_E(I)*AXY_U(IJK)&
229 -MU_sL_pB*(W_S(IPJKM,L)+W_S(IJKM,L))*HALF*OX_E(I)*AXY_U(IJKM))
230
231
232
233
234 = AVG_X(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I)
235 tauzz = -2.d0*MU_sL_p*OX_E(I)*HALF*(&
236 ((W_S(IPJK,L)-W_S(IPJKM,L))*OX(IP)*ODZ(K))+&
237 ((W_S(IJK,L)-W_S(IJKM,L))*OX(I)*ODZ(K)) &
238 ) * VOL_U(IJK)
239
240
241
242 = tauzz + (-2.d0*MU_sL_p*OX_E(I)*OX_E(I)*&
243 U_S(IJK,L)*VOL_U(IJK))
244 ELSE
245 tauzz = ZERO
246 ENDIF
247
248
249
250
251
252
253 = D_P(IJK,L)
254 M_PL = (PI/6.d0)*D_PL**3 * RO_S(IJK,L)
255
256 NU_PM_pE = ROP_S(IJKE,M)/M_PM
257 NU_PM_pW = ROP_S(IJK,M)/M_PM
258 NU_PM_p = AVG_X(NU_PM_pW,NU_PM_pE,I)
259
260 NU_PL_pE = ROP_S(IJKE,L)/M_PL
261 NU_PL_pW = ROP_S(IJK,L)/M_PL
262 NU_PL_p = AVG_X(NU_PL_pW,NU_PL_pE,I)
263
264 Fnu_s_p = AVG_X(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKE,M,L),I)
265 DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pE-NU_PM_pW)*ODX_E(I)
266 DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pE-NU_PL_pW)*ODX_E(I)
267 DS1plusDS2 = DS1 + DS2
268
269
270
271 = Theta_M(IJKE,M)
272 T_PM_pW = Theta_M(IJK,M)
273
274 FT_sM_p = AVG_X(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKE,M,L),I)
275 DS3 = FT_sM_p*(T_PM_pE-T_PM_pW)*ODX_E(I)
276
277
278
279 = Theta_M(IJKE,L)
280 T_PL_pW = Theta_M(IJK,L)
281
282 FT_sL_p = AVG_X(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKE,M,L),I)
283 DS4 = FT_sL_p*(T_PL_pE-T_PL_pW)*ODX_E(I)
284
285
286
287
288 = STRESS_TERMS + SSUX + SSUY + SSUZ + &
289 SSBV + SSX + SSY + SSZ + tauzz
290 DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_U(IJK)
291
292 ELSE
293
294
295 = STRESS_TERMS + ZERO
296 DRAG_TERMS = DRAG_TERMS + ZERO
297
298 ENDIF
299 ENDDO
300
301 (IJK,M) = STRESS_TERMS + DRAG_TERMS
302 ELSE
303 KTMOM_U_S(IJK,M) = ZERO
304 ENDIF
305 ENDDO
306
307 RETURN
308 END SUBROUTINE CALC_IA_MOMSOURCE_U_S
309