File: /nfs/home/0/users/jenkins/mfix.git/model/kintheory_w_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_KTMOMSOURCE_W_S                                    C
4     !  Purpose: Determine source terms for W_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_W_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_w_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_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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
49     !                                                                      C
50     !  Subroutine: CALC_IA_MOMSOURCE_W_S                                   C
51     !  Purpose: Determine source terms for W_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_W_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 :: 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     ! Phase index
113           INTEGER :: L
114     ! Viscosity values
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
118           DOUBLE PRECISION :: XI_sL_pT, XI_sL_pB, LAMBDA_sL_pT, LAMBDA_sL_pB
119     ! Variables for drag calculations
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     ! Average velocity gradients
124           DOUBLE PRECISION :: duodz
125     ! Average volume fraction
126           DOUBLE PRECISION :: EPSA
127     ! Source terms (Surface)
128           DOUBLE PRECISION :: sswx, sswy, sswz, ssx, ssy, ssz, ssbv, tauxz_x
129     ! Source terms (Volumetric)
130           DOUBLE PRECISION :: tauxz_ox, DS1, DS2, DS3, DS4, DS1plusDS2
131     !-----------------------------------------------
132     
133     ! section largely based on tau_w_g:
134     
135           DO IJK = IJKSTART3, IJKEND3
136     
137     ! Skip walls where some values are undefined.
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     ! additional required quantities:
171                    IPJK = IP_OF(IJK)
172                    IJPK = JP_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 = 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     ! bulk viscosity term
204                              XI_sL_pT = 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     ! off diagonal shear stress terms
211                              SSX = 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                              !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
215     
216                              SSY = 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     ! special terms for cylindrical coordinates
222                              IF (CYLINDRICAL) THEN
223     ! tau_zz term: modify Ssz: (1/x) (d/dz) (tau_zz)
224     !                             integral of (1/x) (d/dz) (2mu*(u/x))
225     !                             (normally part of tau_w_s) - explicit
226                                   SSZ = 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     ! tau_xz/x terms: (tau_xz/x)
232     !                             integral of (1/x)*(mu/x)*du/dz
233     !                             (normally part of tau_w_s) - explicit
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     !                             integral of (1/x)*mu*dw/dx
245     !                             (normally part of source_w_s)
246                                   tauxz_ox = 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     !                             integral of (1/x)*(-mu/x)*w
251     !                             (normally part of source_w_s)
252                                   tauxz_ox = tauxz_ox - (OX(I)*OX(I)*MU_sL_p*W_S(IJK,L))
253     
254     !                             multiply all tau_xz/x terms by volume:
255                                   tauxz_ox = tauxz_ox*VOL_W(IJK)
256     
257     ! x*tau_xz term: (1/x) (d/dz) (x*tau_xz)
258     !                             integral of (1/x)*d( (-x)*(mu/x)*w )/dx
259     !                             (normally part of source_w_s)
260                                   tauxz_x = -( ( 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                                        !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
265                              ELSE
266                                   tauxz_ox = ZERO
267                                   tauxz_x = ZERO
268                              ENDIF
269     !--------------------- End Sources from Stress Term ---------------------
270     
271     
272     !--------------------- Sources from Momentum Source Term ---------------------
273     ! Momentum source associated with the difference in the gradients in
274     ! number density of solids phase m and all other solids phases
275                              D_PL = 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     ! Momentum source associated with the gradient in granular
292     ! temperature of species M
293                              T_PM_pT = 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     ! Momentum source associated with the gradient in granular
300     ! temperature of species L
301                              T_PL_pT = 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     !--------------------- End Sources from Momentum Source Term ---------------------
307     
308     
309     ! Add the terms
310                         STRESS_TERMS = 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 ! if m .ne. L
315     ! for m = l all stress terms should already be handled in existing routines
316     ! for m = l all drag terms should become zero
317                              STRESS_TERMS = STRESS_TERMS + ZERO
318                              DRAG_TERMS = DRAG_TERMS + ZERO
319                         ENDIF ! if m .ne. L
320                    ENDDO     ! over L
321     
322                    KTMOM_W_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
323               ELSE
324                    KTMOM_W_S(IJK,M) = ZERO
325               ENDIF     ! dilute
326           ENDDO     ! over ijk
327     
328           RETURN
329           END SUBROUTINE CALC_IA_MOMSOURCE_W_S
330