File: RELATIVE:/../../../mfix.git/model/kintheory_u_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_KTMOMSOURCE_U_S                                    C
4     !  Purpose: Determine source terms for U_S momentum equation arising   C
5     !           from kinetic theory constitutive relations for stress      C
6     !           and solid-solid drag                                       C
7     !                                                                      C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10     
11           SUBROUTINE CALC_KTMOMSOURCE_U_S()
12     
13     !-----------------------------------------------
14     ! Modules
15     !-----------------------------------------------
16     
17           USE param1, only: zero
18     
19     ! kinetic theories
20           USE run, only: kt_type_enum
21           USE run, only: ia_2005
22     
23     ! number of solids phases
24           USE physprop, only: smax
25     
26     ! solids source term
27           USE kintheory, only: ktmom_u_s
28     
29           IMPLICIT NONE
30     !-----------------------------------------------
31     ! Local variables
32     !-----------------------------------------------
33     ! Solids phase index
34           INTEGER :: M
35     !-----------------------------------------------
36     
37           DO M = 1, SMAX
38              KTMOM_U_s(:,M) = ZERO
39              IF (KT_TYPE_ENUM == IA_2005) THEN
40                 CALL CALC_IA_MOMSOURCE_U_S (M)
41              ENDIF
42           ENDDO
43     
44           RETURN
45           END SUBROUTINE CALC_KTMOMSOURCE_U_S
46     
47     
48     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
49     !                                                                      C
50     !  Subroutine: CALC_IA_MOMSOURCE_U_S                                   C
51     !  Purpose: Determine source terms for U_S momentum equation arising   C
52     !           from IA kinetic theory constitutive relations for stress   C
53     !           and solid-solid drag                                       C
54     !                                                                      C
55     !  Literature/Document References:                                     C
56     !    Idir, Y.H., "Modeling of the multiphase mixture of particles      C
57     !      using the kinetic theory approach," PhD Thesis, Illinois        C
58     !      Institute of Technology, Chicago, Illinois, 2004                C
59     !    Iddir, Y.H., & H. Arastoopour, "Modeling of multitype particle    C
60     !      flow using the kinetic theory approach," AIChE J., Vol 51,      C
61     !      No 6, June 2005                                                 C
62     !                                                                      C
63     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
64     
65           SUBROUTINE CALC_IA_MOMSOURCE_U_S(M)
66     
67     !-----------------------------------------------
68     ! Modules
69     !-----------------------------------------------
70           USE param1, only: half, zero
71           USE constant, only: pi
72     
73     ! trace of D_s at i, j, k
74           USE visc_s, only: trD_s
75     
76     ! number of solids phases
77           USE physprop, only: mmax
78     
79     ! x,y,z-components of solids velocity
80           USE fldvar, only: u_s, v_s, w_s
81     ! particle diameter, bulk density, material density
82           USE fldvar, only: d_p, rop_s, ro_s
83     ! granular temperature
84           USE fldvar, only: theta_m, ep_s
85     ! dilute threshold
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     ! Dummy arguments
99     !-----------------------------------------------
100     ! solids phase index
101           INTEGER, INTENT(IN) :: M
102     !-----------------------------------------------
103     ! Local variables
104     !-----------------------------------------------
105     ! Temporary variable
106           DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
107     ! Indices
108           INTEGER :: I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
109                      IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
110                      IJKB, IJKBE, IJKM, IPJKM, &
111                      IJPK, IJKP, IM, IJKW
112     ! solids phase index
113           INTEGER :: L
114     ! Viscosity values for stress calculations
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     ! Bulk viscosity values for calculations
118           DOUBLE PRECISION :: XI_sL_pE, XI_sL_pW, LAMBDA_sL_pE, LAMBDA_sL_pW
119     ! Variables for drag calculations
120           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pE, NU_PM_pW, NU_PM_p, &
121                               NU_PL_pE, NU_PL_pW, NU_PL_p, T_PM_pE, T_PM_pW, &
122                               T_PL_pE, T_PL_pW, Fnu_s_p, FT_sM_p, FT_sL_p
123     ! Average volume fraction
124           DOUBLE PRECISION :: EPSA
125     ! Source terms (Surface)
126           DOUBLE PRECISION :: ssux, ssuy, ssuz, ssx, ssy, ssz, ssbv
127     ! Source terms (Volumetric)
128           DOUBLE PRECISION :: tauzz, DS1, DS2, DS3, DS4, DS1plusDS2
129     !-----------------------------------------------
130     
131     ! section largely based on tau_u_g:
132     
133           DO IJK = IJKSTART3, IJKEND3
134     
135     ! Skip walls where some values are undefined.
136              IF(WALL_AT(IJK)) cycle
137     
138              D_PM = D_P(IJK,M)
139              M_PM = (Pi/6d0)*D_PM**3 * RO_S(IJK,M)
140     
141              I = I_OF(IJK)
142              IJKE = EAST_OF(IJK)
143              EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
144              IF ( .NOT.SIP_AT_E(IJK) .AND. EPSA>DIL_EP_S) THEN
145     
146                 IP = IP1(I)
147                 IM = Im1(I)
148                 IJKW  = WEST_OF(IJK)
149                 J = J_OF(IJK)
150                 JM = JM1(J)
151                 K = K_OF(IJK)
152                 KM = KM1(K)
153                 IPJK = IP_OF(IJK)
154                 IMJK = IM_OF(IJK)
155                 IJMK = JM_OF(IJK)
156                 IJKM = KM_OF(IJK)
157                 IPJMK = JM_OF(IPJK)
158                 IPJKM = IP_OF(IJKM)
159     
160                 IJKN = NORTH_OF(IJK)
161                 IJKS = SOUTH_OF(IJK)
162                 IJKNE = EAST_OF(IJKN)
163                 IJKSE = EAST_OF(IJKS)
164     
165                 IJKT = TOP_OF(IJK)
166                 IJKB = BOTTOM_OF(IJK)
167                 IJKTE = EAST_OF(IJKT)
168                 IJKBE = EAST_OF(IJKB)
169     
170     ! additional required quantities:
171                 IJPK = JP_OF(IJK)
172                 IJKP = KP_OF(IJK)
173     
174     ! initialize variable
175                 STRESS_TERMS = ZERO
176                 DRAG_TERMS = ZERO
177     
178                 DO L = 1,MMAX
179                    IF (M .ne. L) THEN
180     
181     !--------------------- Sources from Stress Terms ---------------------
182     ! Surface Forces
183     ! standard shear stress terms (i.e. ~diffusion)
184                       MU_sL_pE = MU_sL_ip(IJKE,M,L)
185                       MU_sL_pW = MU_sL_ip(IJK,M,L)
186                       SSUX = MU_sL_pE*(U_S(IPJK,L)-U_S(IJK,L))*AYZ_U(IJK)*ODX(IP)&
187                            -MU_sL_PW*(U_S(IJK,L)-U_S(IMJK,L))*AYZ_U(IMJK)*ODX(I)
188     
189                       MU_sL_pN = AVG_X_H(AVG_Y_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKN,M,L), J),&
190                            AVG_Y_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKNE,M,L), J), I)
191                       MU_sL_pS = AVG_X_H(AVG_Y_H(MU_sL_ip(IJKS,M,L),MU_sL_ip(IJK,M,L), JM),&
192                            AVG_Y_H(MU_sL_ip(IJKSE,M,L),MU_sL_ip(IJKE,M,L), JM), I)
193                       SSUY = MU_sL_pN*(U_S(IJPK,L)-U_S(IJK,L))*AXZ_U(IJK)*ODY_N(J)&
194                            -MU_sL_pS*(U_S(IJK,L)-U_S(IJMK,L))*AXZ_U(IJMK)*ODY_N(JM)
195     
196                       MU_sL_pT = AVG_X_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
197                            AVG_Z_H(MU_sL_ip(IJKE,M,L),MU_sL_ip(IJKTE,M,L),K),I)
198                       MU_sL_pB = AVG_X_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
199                            AVG_Z_H(MU_sL_ip(IJKBE,M,L),MU_sL_ip(IJKE,M,L),KM),I)
200                       SSUZ = MU_sL_pT*(U_S(IJKP,L)-U_S(IJK,L))*AXY_U(IJK)*ODZ_T(K)*OX_E(I)&
201                            -MU_sL_pB*(U_S(IJK,L)-U_S(IJKM,L))*AXY_U(IJKM)*ODZ_T(KM)*OX_E(I)
202     
203     ! bulk viscosity term
204                       XI_sL_pE = XI_sL_ip(IJKE,M,L)
205                       XI_sL_pW = XI_sL_ip(IJK,M,L)
206                       LAMBDA_sL_pE = -(2.d0/3.d0)*Mu_sL_pE + XI_sL_pE
207                       LAMBDA_sL_pW = -(2.d0/3.d0)*Mu_sL_pW + XI_sL_pW
208                       SSBV = (LAMBDA_sL_pE*TRD_S(IJKE,L)-LAMBDA_sL_pW*TRD_S(IJK,L))*AYZ(IJK)
209     
210     ! off diagonal shear stress terms
211                       SSX = SSUX
212                       SSY = MU_sL_pN*(V_S(IPJK,L)-V_S(IJK,L))*AXZ_U(IJK)*ODX_E(I)&
213                            -MU_sL_pS*(V_S(IPJMK,L)-V_S(IJMK,L))*AXZ_U(IJMK)*ODX_E(I)
214                       SSZ = MU_sL_pT*(W_S(IPJK,L)-W_S(IJK,L))*AXY_U(IJK)*ODX_E(I)&
215                            -MU_sL_pB*(W_S(IPJKM,L)-W_S(IJKM,L))*AXY_U(IJKM)*ODX_E(I)
216     
217     ! special terms for cylindrical coordinates
218                       IF (CYLINDRICAL) THEN
219     
220     ! modify Ssz: (1/x) (d/dz) (tau_zz)
221     !                    integral of (1/x) (d/dz) (mu*(-w/x))
222     !                    (normally part of tau_u_s) - explicit
223                          SSZ = SSZ - (MU_sL_pT*(W_S(IPJK,L)+W_S(IJK,L))*HALF*OX_E(I)*AXY_U(IJK)&
224                               -MU_sL_pB*(W_S(IPJKM,L)+W_S(IJKM,L))*HALF*OX_E(I)*AXY_U(IJKM))
225     
226     ! tau_zz/x terms: (tau_zz/x)
227     !                    integral of -(2mu/x)*(1/x)*dw/dz
228     !                    (normally part of tau_u_s) - explicit
229                          MU_sL_p = AVG_X(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I)
230                          tauzz = -2.d0*MU_sL_p*OX_E(I)*HALF*(&
231                               ((W_S(IPJK,L)-W_S(IPJKM,L))*OX(IP)*ODZ(K))+&
232                               ((W_S(IJK,L)-W_S(IJKM,L))*OX(I)*ODZ(K)) &
233                               ) * VOL_U(IJK)
234     
235     !                    integral of -(2mu/x)*(1/x)*u
236     !                    (normally part of source_u_s)
237                          tauzz = tauzz + (-2.d0*MU_sL_p*OX_E(I)*OX_E(I)*&
238                               U_S(IJK,L)*VOL_U(IJK))
239                       ELSE
240                          tauzz = ZERO
241                       ENDIF
242     !--------------------- End Sources from Stress Term ---------------------
243     
244     
245     !--------------------- Sources from Momentum Source Term ---------------------
246     ! Momentum source associated with the difference in the gradients in
247     ! number density of solids phase m and all other solids phases
248                       D_PL = D_P(IJK,L)
249                       M_PL = (PI/6.d0)*D_PL**3 * RO_S(IJK,L)
250     
251                       NU_PM_pE = ROP_S(IJKE,M)/M_PM
252                       NU_PM_pW = ROP_S(IJK,M)/M_PM
253                       NU_PM_p = AVG_X(NU_PM_pW,NU_PM_pE,I)
254     
255                       NU_PL_pE = ROP_S(IJKE,L)/M_PL
256                       NU_PL_pW = ROP_S(IJK,L)/M_PL
257                       NU_PL_p = AVG_X(NU_PL_pW,NU_PL_pE,I)
258     
259                       Fnu_s_p = AVG_X(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKE,M,L),I)
260                       DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pE-NU_PM_pW)*ODX_E(I)
261                       DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pE-NU_PL_pW)*ODX_E(I)
262                       DS1plusDS2 = DS1 + DS2
263     
264     ! Momentum source associated with the gradient in granular
265     ! temperature of species M
266                       T_PM_pE = Theta_M(IJKE,M)
267                       T_PM_pW = Theta_M(IJK,M)
268     
269                       FT_sM_p = AVG_X(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKE,M,L),I)
270                       DS3 = FT_sM_p*(T_PM_pE-T_PM_pW)*ODX_E(I)
271     
272     ! Momentum source associated with the gradient in granular
273     ! temperature of species L
274                       T_PL_pE = Theta_M(IJKE,L)
275                       T_PL_pW = Theta_M(IJK,L)
276     
277                       FT_sL_p = AVG_X(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKE,M,L),I)
278                       DS4 = FT_sL_p*(T_PL_pE-T_PL_pW)*ODX_E(I)
279     !--------------------- End Sources from Momentum Source Term ---------------------
280     
281     
282     ! Add the terms
283                       STRESS_TERMS = STRESS_TERMS + SSUX + SSUY + SSUZ + &
284                            SSBV + SSX + SSY + SSZ + tauzz
285                       DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_U(IJK)
286     
287                    ELSE ! if m .ne. L
288     ! for m = l all stress terms should already be handled in existing routines
289     ! for m = l all drag terms should become zero
290                       STRESS_TERMS = STRESS_TERMS + ZERO
291                       DRAG_TERMS = DRAG_TERMS + ZERO
292     
293                    ENDIF ! if m .ne. L
294                 ENDDO     ! over L
295     
296                 KTMOM_U_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
297              ELSE
298                 KTMOM_U_S(IJK,M) = ZERO
299              ENDIF     ! dilute
300           ENDDO        ! over ijk
301     
302           RETURN
303           END SUBROUTINE CALC_IA_MOMSOURCE_U_S
304