File: N:\mfix\model\kintheory_w_s.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_IA_MOMSOURCE_W_S                                   C
4     !  Purpose: Determine source terms for W_S momentum equation arising   C
5     !           from IA kinetic theory constitutive relations for stress   C
6     !           and solid-solid drag                                       C
7     !                                                                      C
8     !  Literature/Document References:                                     C
9     !    Idir, Y.H., "Modeling of the multiphase mixture of particles      C
10     !      using the kinetic theory approach," PhD Thesis, Illinois        C
11     !      Institute of Technology, Chicago, Illinois, 2004                C
12     !    Iddir, Y.H., & H. Arastoopour, "Modeling of multitype particle    C
13     !      flow using the kinetic theory approach," AIChE J., Vol 51,      C
14     !      No 6, June 2005                                                 C
15     !                                                                      C
16     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
17     
18           SUBROUTINE CALC_IA_MOMSOURCE_W_S(M)
19     
20     !-----------------------------------------------
21     ! Modules
22     !-----------------------------------------------
23           USE param1, only: half, zero, undefined
24           USE constant, only: pi
25     
26     ! trace of D_s at i, j, k
27           USE visc_s, only: trD_s
28     
29     ! number of solids phases
30           USE physprop, only: mmax
31     
32     ! x,y,z-components of solids velocity
33           USE fldvar, only: u_s, v_s, w_s
34     ! particle diameter, bulk density, material density
35           USE fldvar, only: d_p, rop_s, ro_s
36     ! granular temperature
37           USE fldvar, only: theta_m, ep_s
38     ! dilute threshold
39           USE toleranc, only: dil_ep_s
40     
41           Use kintheory
42     
43           USE geometry
44           USE indices
45           USE compar
46           USE fun_avg
47           USE functions
48     
49           IMPLICIT NONE
50     !-----------------------------------------------
51     ! Dummy arguments
52     !-----------------------------------------------
53     ! solids phase index
54           INTEGER, INTENT(IN) :: M
55     !-----------------------------------------------
56     ! Local variables
57     !-----------------------------------------------
58     ! Temporary variable
59           DOUBLE PRECISION :: STRESS_TERMS, DRAG_TERMS
60     ! Indices
61           INTEGER :: IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
62                      IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
63                      IMJKP, K, IJKT, JM, KP, IJKM, IPJK,&
64                      IJPK
65     ! Phase index
66           INTEGER :: L
67     ! Viscosity values
68           DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
69                               MU_sL_pB, MU_sL_p
70     ! Bulk viscosity values
71           DOUBLE PRECISION :: XI_sL_pT, XI_sL_pB, LAMBDA_sL_pT, LAMBDA_sL_pB
72     ! Variables for drag calculations
73           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pT, NU_PM_pB, NU_PM_p, &
74                               NU_PL_pT, NU_PL_pB, NU_PL_p, T_PM_pT, T_PM_pB,&
75                               T_PL_pT, T_PL_pB, Fnu_s_p, FT_sM_p, FT_sL_p
76     ! Average velocity gradients
77           DOUBLE PRECISION :: duodz
78     ! Average volume fraction
79           DOUBLE PRECISION :: EPSA
80     ! Source terms (Surface)
81           DOUBLE PRECISION :: sswx, sswy, sswz, ssx, ssy, ssz, ssbv, tauxz_x
82     ! Source terms (Volumetric)
83           DOUBLE PRECISION :: tauxz_ox, DS1, DS2, DS3, DS4, DS1plusDS2
84     !-----------------------------------------------
85     
86     ! section largely based on tau_w_g:
87     
88           DO IJK = IJKSTART3, IJKEND3
89     
90     ! Skip walls where some values are undefined.
91              IF(WALL_AT(IJK)) cycle
92     
93               D_PM = D_P(IJK,M)
94               M_PM = (Pi/6d0)*(D_PM**3.)*RO_S(IJK,M)
95     
96               K = K_OF(IJK)
97               IJKT = TOP_OF(IJK)
98               EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
99               IF ( .NOT.SIP_AT_T(IJK) .AND. EPSA>DIL_EP_S) THEN
100     
101                    J = J_OF(IJK)
102                    I = I_OF(IJK)
103                    KP = KP1(K)
104                    IM = IM1(I)
105                    JM = JM1(J)
106                    IMJK = IM_OF(IJK)
107                    IJMK = JM_OF(IJK)
108                    IJKM = KM_OF(IJK)
109                    IJKP = KP_OF(IJK)
110                    IJMKP = JM_OF(IJKP)
111                    IMJKP = KP_OF(IMJK)
112     
113                    IJKN = NORTH_OF(IJK)
114                    IJKS = SOUTH_OF(IJK)
115                    IJKNT = TOP_OF(IJKN)
116                    IJKST = TOP_OF(IJKS)
117     
118                    IJKE = EAST_OF(IJK)
119                    IJKW = WEST_OF(IJK)
120                    IJKTE = EAST_OF(IJKT)
121                    IJKTW = WEST_OF(IJKT)
122     
123     ! additional required quantities:
124                    IPJK = IP_OF(IJK)
125                    IJPK = JP_OF(IJK)
126     
127     ! initialize variable
128                    STRESS_TERMS = ZERO
129                    DRAG_TERMS = ZERO
130     
131                    DO L = 1, MMAX
132                         IF (M .ne. L) THEN
133     
134     !--------------------- Sources from Stress Terms ---------------------
135     ! Surface Forces
136     ! standard shear stress terms (i.e. ~diffusion)
137                              MU_sL_pE = AVG_Z_H(AVG_X_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I),&
138                                   AVG_X_H(MU_sL_ip(IJKT,M,L),MU_sL_ip(IJKTE,M,L),I),K)
139                              MU_sL_pW = AVG_Z_H(AVG_X_H(MU_sL_ip(IJKW,M,L),MU_sL_ip(IJK,M,L),IM),&
140                                   AVG_X_H(MU_sL_ip(IJKTW,M,L),MU_sL_ip(IJKT,M,L),IM),K)
141                              SSWX = MU_sL_pE*(W_S(IPJK,L)-W_S(IJK,L))*AYZ_W(IJK)*ODX_E(I)&
142                                   -MU_sL_PW*(W_S(IJK,L)-W_S(IMJK,L))*AYZ_W(IMJK)*ODX_E(IM)
143     
144                              MU_sL_pN = AVG_Z_H(AVG_Y_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKN,M,L), J),&
145                                   AVG_Y_H(MU_sL_ip(IJKT,M,L),MU_sL_ip(IJKNT,M,L), J), K)
146                              MU_sL_pS = AVG_Z_H(AVG_Y_H(MU_sL_ip(IJKS,M,L),MU_sL_ip(IJK,M,L), JM),&
147                                   AVG_Y_H(MU_sL_ip(IJKST,M,L),MU_sL_ip(IJKT,M,L), JM), K)
148                              SSWY = MU_sL_pN*(W_S(IJPK,L)-W_S(IJK,L))*AXZ_W(IJK)*ODY_N(J)&
149                                   -MU_sL_pS*(W_S(IJK,L)-W_S(IJMK,L))*AXZ_W(IJKM)*ODY_N(JM)
150     
151                              MU_sL_pT = MU_sL_ip(IJKT,M,L)
152                              MU_sL_pB = MU_sL_ip(IJK,M,L)
153                              SSWZ = MU_sL_pT*(W_S(IJKP,L)-W_S(IJK,L))*AXY_W(IJK)*ODZ(KP)*OX(I)&
154                                   -MU_sL_pB*(W_S(IJK,L)-W_S(IJKM,L))*AXY_W(IJKM)*ODZ(K)*OX(I)
155     
156     ! bulk viscosity term
157                              XI_sL_pT = XI_sL_ip(IJKT,M,L)
158                              XI_sL_pB = XI_sL_ip(IJK,M,L)
159                              LAMBDA_sL_pT = -(2d0/3d0)*MU_sL_pT + XI_sL_pT
160                              LAMBDA_sL_pB = -(2d0/3d0)*MU_sL_pB + XI_sL_pB
161                              SSBV = (LAMBDA_sL_pT*TRD_S(IJKT,L)-LAMBDA_sL_pB*TRD_S(IJK,L))*AXY(IJK)
162     
163     ! off diagonal shear stress terms
164                              SSX = MU_sL_pE*(U_S(IJKP,L)-U_S(IJK,L))*OX_E(I)*AYZ_W(IJK)*ODZ_T(K)&
165                                   -MU_sL_pW*(U_S(IMJKP,L)-U_S(IMJK,L))*(DY(J)*HALF*(DZ(K)+&
166                                   DZ(KP)))*ODZ_T(K)
167                              !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
168     
169                              SSY = MU_sL_pN*(V_S(IJKP,L)-V_S(IJK,L))*AXZ_W(IJK)*ODZ_T(K)*OX(I)&
170                                   -MU_sL_pS*(V_S(IJMKP,L)-V_S(IJMK,L))*AXZ_W(IJMK)*ODZ_T(K)*OX(I)
171     
172                              SSZ = SSWZ
173     
174     ! special terms for cylindrical coordinates
175                              IF (CYLINDRICAL) THEN
176     ! tau_zz term: modify Ssz: (1/x) (d/dz) (tau_zz)
177     !                             integral of (1/x) (d/dz) (2mu*(u/x))
178     !                             (normally part of tau_w_s) - explicit
179                                   SSZ = SSZ +(&
180                                        (MU_sL_pT*(U_S(IJKP,L)+U_S(IMJKP,L))*OX(I)*AXY_W(IJK))-&
181                                        (MU_sL_pB*(U_S(IJK,L)+U_S(IMJK,L))*OX(I)*AXY_W(IJKM)))
182     
183                                   MU_sL_p = AVG_Z(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K)
184     ! tau_xz/x terms: (tau_xz/x)
185     !                             integral of (1/x)*(mu/x)*du/dz
186     !                             (normally part of tau_w_s) - explicit
187                                   IF (OX_E(IM) /= UNDEFINED) THEN
188                                        DUODZ = (U_S(IMJKP,L)-U_S(IMJK,L))*OX_E(IM)*ODZ_T(K)
189                                   ELSE
190                                        DUODZ = ZERO
191                                   ENDIF
192     
193                                   tauxz_ox = MU_sL_p*OX(I)*HALF*(&
194                                        (OX_E(I)*(U_S(IJKP,L)-U_S(IJK,L))*ODZ_T(K))+&
195                                        DUODZ)
196     
197     !                             integral of (1/x)*mu*dw/dx
198     !                             (normally part of source_w_s)
199                                   tauxz_ox = tauxz_ox + (MU_sL_p*OX(I)*HALF*&
200                                        ( (W_S(IPJK,L)-W_S(IJK,L))*ODX_E(I) + &
201                                        (W_S(IJK,L)-W_S(IMJK,L))*ODX_E(IM) ) )
202     
203     !                             integral of (1/x)*(-mu/x)*w
204     !                             (normally part of source_w_s)
205                                   tauxz_ox = tauxz_ox - (OX(I)*OX(I)*MU_sL_p*W_S(IJK,L))
206     
207     !                             multiply all tau_xz/x terms by volume:
208                                   tauxz_ox = tauxz_ox*VOL_W(IJK)
209     
210     ! x*tau_xz term: (1/x) (d/dz) (x*tau_xz)
211     !                             integral of (1/x)*d( (-x)*(mu/x)*w )/dx
212     !                             (normally part of source_w_s)
213                                   tauxz_x = -( ( MU_sL_pE*OX_E(I)*HALF*(W_S(IPJK,L)+&
214                                        W_S(IJK,L))*AYZ_W(IJK) )-( MU_sL_pW*HALF*&
215                                        (W_S(IJK,L)+W_S(IMJK,L))*&
216                                        (DY(J)*HALF*(DZ(K)+DZ(KP)) ) ) )
217                                        !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
218                              ELSE
219                                   tauxz_ox = ZERO
220                                   tauxz_x = ZERO
221                              ENDIF
222     !--------------------- End Sources from Stress Term ---------------------
223     
224     
225     !--------------------- Sources from Momentum Source Term ---------------------
226     ! Momentum source associated with the difference in the gradients in
227     ! number density of solids phase m and all other solids phases
228                              D_PL = D_P(IJK,L)
229                              M_PL = (Pi/6d0)*(D_PL**3.)*RO_S(IJK,L)
230     
231                              NU_PM_pT = ROP_S(IJKT,M)/M_PM
232                              NU_PM_pB = ROP_S(IJK,M)/M_PM
233                              NU_PM_p = AVG_Z(NU_PM_pB,NU_PM_pT,K)
234     
235                              NU_PL_pT = ROP_S(IJKT,L)/M_PL
236                              NU_PL_pB = ROP_S(IJK,L)/M_PL
237                              NU_PL_p = AVG_Z(NU_PL_pB,NU_PL_pT,K)
238     
239                              Fnu_s_p = AVG_Z(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKT,M,L),K)
240                              DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pT-NU_PM_pB)*OX(I)*ODZ_T(K)
241                              DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pT-NU_PL_pB)*OX(I)*ODZ_T(K)
242                              DS1plusDS2 = DS1 + DS2
243     
244     ! Momentum source associated with the gradient in granular
245     ! temperature of species M
246                              T_PM_pT = Theta_M(IJKT,M)
247                              T_PM_pB = Theta_M(IJK,M)
248     
249                              FT_sM_p = AVG_Z(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKT,M,L),K)
250                              DS3 = FT_sM_p*(T_PM_pT-T_PM_pB)*OX(I)*ODZ_T(K)
251     
252     ! Momentum source associated with the gradient in granular
253     ! temperature of species L
254                              T_PL_pT = Theta_M(IJKT,L)
255                              T_PL_pB = Theta_M(IJK,L)
256     
257                              FT_sL_p = AVG_Z(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKT,M,L),K)
258                              DS4 = FT_sL_p*(T_PL_pT-T_PL_pB)*OX(I)*ODZ_T(K)
259     !--------------------- End Sources from Momentum Source Term ---------------------
260     
261     
262     ! Add the terms
263                         STRESS_TERMS = STRESS_TERMS + SSWX + SSWY + SSWZ + &
264                             SSBV + SSX + SSY + SSZ + tauxz_ox + tauxz_x
265                         DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_W(IJK)
266     
267                         ELSE ! if m .ne. L
268     ! for m = l all stress terms should already be handled in existing routines
269     ! for m = l all drag terms should become zero
270                              STRESS_TERMS = STRESS_TERMS + ZERO
271                              DRAG_TERMS = DRAG_TERMS + ZERO
272                         ENDIF ! if m .ne. L
273                    ENDDO     ! over L
274     
275                    KTMOM_W_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
276               ELSE
277                    KTMOM_W_S(IJK,M) = ZERO
278               ENDIF     ! dilute
279           ENDDO     ! over ijk
280     
281           RETURN
282           END SUBROUTINE CALC_IA_MOMSOURCE_W_S
283