File: /nfs/home/0/users/jenkins/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(IER)
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     ! Dummy arguments
32     !-----------------------------------------------
33     ! Error index
34           INTEGER, INTENT(INOUT) :: IER
35     !-----------------------------------------------
36     ! Local variables
37     !-----------------------------------------------
38     ! Solids phase index
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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
54     !                                                                      C
55     !  Subroutine: CALC_IA_MOMSOURCE_U_S                                   C
56     !  Purpose: Determine source terms for U_S momentum equation arising   C
57     !           from IA kinetic theory constitutive relations for stress   C
58     !           and solid-solid drag                                       C
59     !                                                                      C
60     !  Literature/Document References:                                     C
61     !    Idir, Y.H., "Modeling of the multiphase mixture of particles      C
62     !      using the kinetic theory approach," PhD Thesis, Illinois        C
63     !      Institute of Technology, Chicago, Illinois, 2004                C
64     !    Iddir, Y.H., & H. Arastoopour, "Modeling of multitype particle    C
65     !      flow using the kinetic theory approach," AIChE J., Vol 51,      C
66     !      No 6, June 2005                                                 C
67     !                                                                      C
68     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
69     
70           SUBROUTINE CALC_IA_MOMSOURCE_U_S(M)
71     
72     !-----------------------------------------------
73     ! Modules
74     !-----------------------------------------------
75           USE param1, only: half, zero
76           USE constant, only: pi
77     
78     ! trace of D_s at i, j, k
79           USE visc_s, only: trD_s
80     
81     ! number of solids phases
82           USE physprop, only: mmax
83     
84     ! x,y,z-components of solids velocity
85           USE fldvar, only: u_s, v_s, w_s
86     ! particle diameter, bulk density, material density
87           USE fldvar, only: d_p, rop_s, ro_s
88     ! granular temperature
89           USE fldvar, only: theta_m, ep_s
90     ! dilute threshold
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     ! Dummy arguments
104     !-----------------------------------------------
105     ! solids phase index
106           INTEGER, INTENT(IN) :: M
107     !-----------------------------------------------
108     ! Local variables
109     !-----------------------------------------------
110     ! Temporary variable
111           DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
112     ! Indices
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     ! solids phase index
118           INTEGER :: L
119     ! Viscosity values for stress calculations
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     ! Bulk viscosity values for calculations
123           DOUBLE PRECISION :: XI_sL_pE, XI_sL_pW, LAMBDA_sL_pE, LAMBDA_sL_pW
124     ! Variables for drag calculations
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     ! Average volume fraction
129           DOUBLE PRECISION :: EPSA
130     ! Source terms (Surface)
131           DOUBLE PRECISION :: ssux, ssuy, ssuz, ssx, ssy, ssz, ssbv
132     ! Source terms (Volumetric)
133           DOUBLE PRECISION :: tauzz, DS1, DS2, DS3, DS4, DS1plusDS2
134     !-----------------------------------------------
135     
136     ! section largely based on tau_u_g:
137     
138           DO IJK = IJKSTART3, IJKEND3
139     
140     ! Skip walls where some values are undefined.
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     ! additional required quantities:
176                 IJPK = JP_OF(IJK)
177                 IJKP = KP_OF(IJK)
178     
179     ! initialize variable
180                 STRESS_TERMS = ZERO
181                 DRAG_TERMS = ZERO
182     
183                 DO L = 1,MMAX
184                    IF (M .ne. L) THEN
185     
186     !--------------------- Sources from Stress Terms ---------------------
187     ! Surface Forces
188     ! standard shear stress terms (i.e. ~diffusion)
189                       MU_sL_pE = 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     ! bulk viscosity term
209                       XI_sL_pE = 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     ! off diagonal shear stress terms
216                       SSX = 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     ! special terms for cylindrical coordinates
223                       IF (CYLINDRICAL) THEN
224     
225     ! modify Ssz: (1/x) (d/dz) (tau_zz)
226     !                    integral of (1/x) (d/dz) (mu*(-w/x))
227     !                    (normally part of tau_u_s) - explicit
228                          SSZ = 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     ! tau_zz/x terms: (tau_zz/x)
232     !                    integral of -(2mu/x)*(1/x)*dw/dz
233     !                    (normally part of tau_u_s) - explicit
234                          MU_sL_p = 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     !                    integral of -(2mu/x)*(1/x)*u
241     !                    (normally part of source_u_s)
242                          tauzz = 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     !--------------------- End Sources from Stress Term ---------------------
248     
249     
250     !--------------------- Sources from Momentum Source Term ---------------------
251     ! Momentum source associated with the difference in the gradients in
252     ! number density of solids phase m and all other solids phases
253                       D_PL = 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     ! Momentum source associated with the gradient in granular
270     ! temperature of species M
271                       T_PM_pE = 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     ! Momentum source associated with the gradient in granular
278     ! temperature of species L
279                       T_PL_pE = 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     !--------------------- End Sources from Momentum Source Term ---------------------
285     
286     
287     ! Add the terms
288                       STRESS_TERMS = 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 ! if m .ne. L
293     ! for m = l all stress terms should already be handled in existing routines
294     ! for m = l all drag terms should become zero
295                       STRESS_TERMS = STRESS_TERMS + ZERO
296                       DRAG_TERMS = DRAG_TERMS + ZERO
297     
298                    ENDIF ! if m .ne. L
299                 ENDDO     ! over L
300     
301                 KTMOM_U_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
302              ELSE
303                 KTMOM_U_S(IJK,M) = ZERO
304              ENDIF     ! dilute
305           ENDDO        ! over ijk
306     
307           RETURN
308           END SUBROUTINE CALC_IA_MOMSOURCE_U_S
309