File: /nfs/home/0/users/jenkins/mfix.git/model/kintheory_v_s.f
1
2
3
4
5
6
7
8
9
10
11 SUBROUTINE CALC_KTMOMSOURCE_V_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_v_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_V_s(:,M) = ZERO
44 IF (KT_TYPE_ENUM == IA_2005) THEN
45 CALL CALC_IA_MOMSOURCE_V_S (M)
46 ENDIF
47 ENDDO
48
49 RETURN
50 END SUBROUTINE CALC_KTMOMSOURCE_V_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_V_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, K, IJK, IJKN, JP, IM, KM, IJPK, IJMK,&
114 IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
115 IJKTN, IJKB, IJKBN, IJKM, IJPKM,&
116 IPJK, IJKP
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
122
123 DOUBLE PRECISION :: XI_sL_pN, XI_sL_pS, LAMBDA_sL_pN, LAMBDA_sL_pS
124
125 DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pN, NU_PM_pS, NU_PM_p, &
126 NU_PL_pN, NU_PL_pS, NU_PL_p, T_PM_pN, T_PM_pS, &
127 T_PL_pN, T_PL_pS, Fnu_s_p, FT_sM_p, FT_sL_p
128
129 DOUBLE PRECISION :: EPSA
130
131 DOUBLE PRECISION :: ssvx, ssvy, ssvz, ssx, ssy, ssz, ssbv
132
133 DOUBLE PRECISION :: 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 J = J_OF(IJK)
147 IJKN = NORTH_OF(IJK)
148 EPSA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
149 IF ( .NOT.SIP_AT_N(IJK) .AND. EPSA>DIL_EP_S) THEN
150
151 JP = JP1(J)
152 I = I_OF(IJK)
153 IM = IM1(I)
154 K = K_OF(IJK)
155 KM = KM1(K)
156 IJPK = JP_OF(IJK)
157 IJMK = JM_OF(IJK)
158 IMJK = IM_OF(IJK)
159 IJKM = KM_OF(IJK)
160 IMJPK = IM_OF(IJPK)
161 IJPKM = JP_OF(IJKM)
162
163 IJKW = WEST_OF(IJK)
164 IJKE = EAST_OF(IJK)
165 IJKNE = EAST_OF(IJKN)
166 IJKNW = NORTH_OF(IJKW)
167
168 IJKB = BOTTOM_OF(IJK)
169 IJKT = TOP_OF(IJK)
170 IJKTN = NORTH_OF(IJKT)
171 IJKBN = NORTH_OF(IJKB)
172
173
174 = IP_OF(IJK)
175 IJKP = KP_OF(IJK)
176
177
178 = ZERO
179 DRAG_TERMS = ZERO
180
181 DO L = 1, MMAX
182 IF (M .ne. L) THEN
183
184
185
186
187 = AVG_Y_H(AVG_X_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I),&
188 AVG_X_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKNE,M,L),I),J)
189 MU_sL_pW = AVG_Y_H(AVG_X_H(MU_sL_ip(IJKW,M,L),MU_sL_ip(IJK,M,L),IM),&
190 AVG_X_H(MU_sL_ip(IJKNW,M,L),MU_sL_ip(IJKN,M,L),IM),J)
191 SSVX = MU_sL_pE*(V_S(IPJK,L)-V_S(IJK,L))*AYZ_V(IJK)*ODX_E(I)&
192 -MU_sL_pW*(V_S(IJK,L)-V_S(IMJK,L))*AYZ_V(IMJK)*ODX_E(IM)
193
194 MU_sL_pN = MU_sL_ip(IJKN,M,L)
195 MU_sL_pS = MU_sL_ip(IJK,M,L)
196 SSVY = MU_sL_pN*(V_S(IJPK,L)-V_S(IJK,L))*ODY(JP)*AXZ_V(IJK)&
197 -MU_sL_pS*(V_S(IJK,L)-V_S(IJMK,L))*ODY(J)*AXZ_V(IJMK)
198
199 MU_sL_pT = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
200 AVG_Z_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKTN,M,L),K),J)
201 MU_sL_pB = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
202 AVG_Z_H(MU_sL_ip(IJKBN,M,L),MU_sL_ip(IJKN,M,L),KM),J)
203 SSVZ = MU_sL_pT*(V_S(IJKP,L)-V_S(IJK,L))*AXY_V(IJK)*ODZ_T(K)*OX(I)&
204 -MU_sL_pB*(V_S(IJK,L)-V_S(IJKM,L))*AXY_V(IJKM)*ODZ_T(KM)*OX(I)
205
206
207 = XI_sL_ip(IJKN,M,L)
208 XI_sL_pS = XI_sL_ip(IJK,M,L)
209 LAMBDA_sL_pN = -(2.d0/3.d0)*MU_sL_pN + XI_sL_pN
210 LAMBDA_sL_pS = -(2.d0/3.d0)*MU_sL_pS + XI_sL_pS
211 SSBV = (LAMBDA_sL_pN*TRD_S(IJKN,L)-LAMBDA_sL_pS*TRD_S(IJK,L))*AXZ(IJK)
212
213
214 = MU_sL_pE*(U_S(IJPK,L)-U_S(IJK,L))*ODY_N(J)*AYZ_V(IJK)&
215 -MU_sL_pW*(U_S(IMJPK,L)-U_S(IMJK,L))*ODY_N(J)*AYZ_V(IMJK)
216 SSY = SSVY
217 SSZ = MU_sL_pT*(W_S(IJPK,L)-W_S(IJK,L))*AXY_V(IJK)*ODY_N(J)&
218 -MU_sL_pB*(W_S(IJPKM,L)-W_S(IJKM,L))*AXY_V(IJKM)*ODY_N(J)
219
220
221
222
223
224
225 = D_P(IJK,L)
226 M_PL = (Pi/6d0)* D_PL**3 *RO_S(IJK,L)
227
228 NU_PM_pN = ROP_S(IJKN,M)/M_PM
229 NU_PM_pS = ROP_S(IJK,M)/M_PM
230 NU_PM_p = AVG_Y(NU_PM_pS,NU_PM_pN,J)
231
232 NU_PL_pN = ROP_S(IJKN,L)/M_PL
233 NU_PL_pS = ROP_S(IJK,L)/M_PL
234 NU_PL_p = AVG_Y(NU_PL_pS,NU_PL_pN,J)
235
236 Fnu_s_p = AVG_Y(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKN,M,L),J)
237 DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pN-NU_PM_pS)*ODY_N(J)
238 DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pN-NU_PL_pS)*ODY_N(J)
239 DS1plusDS2 = DS1 + DS2
240
241
242
243 = Theta_M(IJKN,M)
244 T_PM_pS = Theta_M(IJK,M)
245
246 FT_sM_p = AVG_Y(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKN,M,L),J)
247 DS3 = FT_sM_p*(T_PM_pN-T_PM_pS)*ODY_N(J)
248
249
250
251 = Theta_M(IJKN,L)
252 T_PL_pS = Theta_M(IJK,L)
253
254 FT_sL_p = AVG_Y(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKN,M,L),J)
255 DS4 = FT_sL_p*(T_PL_pN-T_PL_pS)*ODY_N(J)
256
257
258
259
260 = STRESS_TERMS + SSVX + SSVY + SSVZ + &
261 SSBV + SSX + SSY + SSZ
262 DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_V(IJK)
263
264 ELSE
265
266
267 = STRESS_TERMS + ZERO
268 DRAG_TERMS = DRAG_TERMS + ZERO
269
270 ENDIF
271 ENDDO
272
273 (IJK,M) = STRESS_TERMS + DRAG_TERMS
274 ELSE
275 KTMOM_V_S(IJK,M) = ZERO
276 ENDIF
277 ENDDO
278
279 RETURN
280 END SUBROUTINE CALC_IA_MOMSOURCE_V_S
281