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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_KTMOMSOURCE_V_S                                    C
4     !  Purpose: Determine source terms for V_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_V_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_v_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_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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
54     !                                                                      C
55     !  Subroutine: CALC_IA_MOMSOURCE_V_S                                   C
56     !  Purpose: Determine source terms for V_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_V_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, 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     ! Phase index
118           INTEGER :: L
119     ! Viscosity values
120           DOUBLE PRECISION :: MU_sL_pE, MU_sL_pW, MU_sL_pN, MU_sL_pS, MU_sL_pT,&
121                               MU_sL_pB
122     ! Bulk viscosity values
123           DOUBLE PRECISION :: XI_sL_pN, XI_sL_pS, LAMBDA_sL_pN, LAMBDA_sL_pS
124     ! Variables for drag calculations
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     ! Average volume fraction
129           DOUBLE PRECISION :: EPSA
130     ! Source terms (Surface)
131           DOUBLE PRECISION :: ssvx, ssvy, ssvz, ssx, ssy, ssz, ssbv
132     ! Source terms (Volumetric)
133           DOUBLE PRECISION :: DS1, DS2, DS3, DS4, DS1plusDS2
134     !-----------------------------------------------
135     
136     ! section largely based on tau_v_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               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     ! additional required quantities:
174                    IPJK = IP_OF(IJK)
175                    IJKP = KP_OF(IJK)
176     
177     ! initialize variable
178                    STRESS_TERMS = ZERO
179                    DRAG_TERMS = ZERO
180     
181                    DO L = 1, MMAX
182                         IF (M .ne. L) THEN
183     
184     !--------------------- Sources from Stress Terms ---------------------
185     ! Surface Forces
186     ! standard shear stress terms (i.e. ~diffusion)
187                              MU_sL_pE = 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     ! bulk viscosity term
207                              XI_sL_pN = 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     ! off diagonal shear stress terms
214                              SSX = 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     !--------------------- End Sources from Stress Term ---------------------
220     
221     
222     !--------------------- Sources from Momentum Source Term ---------------------
223     ! Momentum source associated with the difference in the gradients in
224     ! number density of solids phase m and all other solids phases
225                              D_PL = 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     ! Momentum source associated with the gradient in granular
242     ! temperature of species M
243                              T_PM_pN = 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     ! Momentum source associated with the gradient in granular
250     ! temperature of species L
251                              T_PL_pN = 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     !--------------------- End Sources from Momentum Source Term ---------------------
257     
258     
259     ! Add the terms
260                              STRESS_TERMS = STRESS_TERMS + SSVX + SSVY + SSVZ + &
261                                  SSBV + SSX + SSY + SSZ
262                              DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_V(IJK)
263     
264                         ELSE ! if m .ne. L
265     ! for m = l all stress terms should already be handled in existing routines
266     ! for m = l all drag terms should become zero
267                              STRESS_TERMS = STRESS_TERMS + ZERO
268                              DRAG_TERMS = DRAG_TERMS + ZERO
269     
270                         ENDIF ! if m .ne. L
271                    ENDDO     ! over L
272     
273                    KTMOM_V_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
274               ELSE
275                    KTMOM_V_S(IJK,M) = ZERO
276               ENDIF     ! dilute
277           ENDDO        ! over ijk
278     
279           RETURN
280           END SUBROUTINE CALC_IA_MOMSOURCE_V_S
281