File: RELATIVE:/../../../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()
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     ! Local variables
32     !-----------------------------------------------
33     ! Solids phase index
34           INTEGER :: M
35     !-----------------------------------------------
36     
37           DO M = 1, SMAX
38              KTMOM_V_s(:,M) = ZERO
39              IF (KT_TYPE_ENUM == IA_2005) THEN
40                 CALL CALC_IA_MOMSOURCE_V_S (M)
41              ENDIF
42           ENDDO
43     
44           RETURN
45           END SUBROUTINE CALC_KTMOMSOURCE_V_S
46     
47     
48     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
49     !                                                                      C
50     !  Subroutine: CALC_IA_MOMSOURCE_V_S                                   C
51     !  Purpose: Determine source terms for V_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_V_S(M)
66     
67     !-----------------------------------------------
68     ! Modules
69     !-----------------------------------------------
70           USE param1, only: 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, K, IJK, IJKN, JP, IM,  KM, IJPK, IJMK,&
109                      IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
110                      IJKTN, IJKB, IJKBN, IJKM, IJPKM,&
111                      IPJK, IJKP
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
117     ! Bulk viscosity values
118           DOUBLE PRECISION :: XI_sL_pN, XI_sL_pS, LAMBDA_sL_pN, LAMBDA_sL_pS
119     ! Variables for drag calculations
120           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pN, NU_PM_pS, NU_PM_p, &
121                               NU_PL_pN, NU_PL_pS, NU_PL_p, T_PM_pN, T_PM_pS, &
122                               T_PL_pN, T_PL_pS, Fnu_s_p, FT_sM_p, FT_sL_p
123     ! Average volume fraction
124           DOUBLE PRECISION :: EPSA
125     ! Source terms (Surface)
126           DOUBLE PRECISION :: ssvx, ssvy, ssvz, ssx, ssy, ssz, ssbv
127     ! Source terms (Volumetric)
128           DOUBLE PRECISION :: DS1, DS2, DS3, DS4, DS1plusDS2
129     !-----------------------------------------------
130     
131     ! section largely based on tau_v_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               J = J_OF(IJK)
142               IJKN = NORTH_OF(IJK)
143               EPSA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
144               IF ( .NOT.SIP_AT_N(IJK) .AND. EPSA>DIL_EP_S) THEN
145     
146                    JP = JP1(J)
147                    I = I_OF(IJK)
148                    IM = IM1(I)
149                    K = K_OF(IJK)
150                    KM = KM1(K)
151                    IJPK = JP_OF(IJK)
152                    IJMK = JM_OF(IJK)
153                    IMJK = IM_OF(IJK)
154                    IJKM = KM_OF(IJK)
155                    IMJPK = IM_OF(IJPK)
156                    IJPKM = JP_OF(IJKM)
157     
158                    IJKW = WEST_OF(IJK)
159                    IJKE = EAST_OF(IJK)
160                    IJKNE = EAST_OF(IJKN)
161                    IJKNW = NORTH_OF(IJKW)
162     
163                    IJKB = BOTTOM_OF(IJK)
164                    IJKT = TOP_OF(IJK)
165                    IJKTN = NORTH_OF(IJKT)
166                    IJKBN = NORTH_OF(IJKB)
167     
168     ! additional required quantities:
169                    IPJK = IP_OF(IJK)
170                    IJKP = KP_OF(IJK)
171     
172     ! initialize variable
173                    STRESS_TERMS = ZERO
174                    DRAG_TERMS = ZERO
175     
176                    DO L = 1, MMAX
177                         IF (M .ne. L) THEN
178     
179     !--------------------- Sources from Stress Terms ---------------------
180     ! Surface Forces
181     ! standard shear stress terms (i.e. ~diffusion)
182                              MU_sL_pE = AVG_Y_H(AVG_X_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I),&
183                                   AVG_X_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKNE,M,L),I),J)
184                              MU_sL_pW = AVG_Y_H(AVG_X_H(MU_sL_ip(IJKW,M,L),MU_sL_ip(IJK,M,L),IM),&
185                                   AVG_X_H(MU_sL_ip(IJKNW,M,L),MU_sL_ip(IJKN,M,L),IM),J)
186                              SSVX = MU_sL_pE*(V_S(IPJK,L)-V_S(IJK,L))*AYZ_V(IJK)*ODX_E(I)&
187                                   -MU_sL_pW*(V_S(IJK,L)-V_S(IMJK,L))*AYZ_V(IMJK)*ODX_E(IM)
188     
189                              MU_sL_pN = MU_sL_ip(IJKN,M,L)
190                              MU_sL_pS = MU_sL_ip(IJK,M,L)
191                              SSVY = MU_sL_pN*(V_S(IJPK,L)-V_S(IJK,L))*ODY(JP)*AXZ_V(IJK)&
192                                   -MU_sL_pS*(V_S(IJK,L)-V_S(IJMK,L))*ODY(J)*AXZ_V(IJMK)
193     
194                              MU_sL_pT = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
195                                   AVG_Z_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKTN,M,L),K),J)
196                              MU_sL_pB = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
197                                   AVG_Z_H(MU_sL_ip(IJKBN,M,L),MU_sL_ip(IJKN,M,L),KM),J)
198                              SSVZ = MU_sL_pT*(V_S(IJKP,L)-V_S(IJK,L))*AXY_V(IJK)*ODZ_T(K)*OX(I)&
199                                   -MU_sL_pB*(V_S(IJK,L)-V_S(IJKM,L))*AXY_V(IJKM)*ODZ_T(KM)*OX(I)
200     
201     ! bulk viscosity term
202                              XI_sL_pN = XI_sL_ip(IJKN,M,L)
203                              XI_sL_pS = XI_sL_ip(IJK,M,L)
204                              LAMBDA_sL_pN = -(2.d0/3.d0)*MU_sL_pN + XI_sL_pN
205                              LAMBDA_sL_pS = -(2.d0/3.d0)*MU_sL_pS + XI_sL_pS
206                              SSBV = (LAMBDA_sL_pN*TRD_S(IJKN,L)-LAMBDA_sL_pS*TRD_S(IJK,L))*AXZ(IJK)
207     
208     ! off diagonal shear stress terms
209                              SSX = MU_sL_pE*(U_S(IJPK,L)-U_S(IJK,L))*ODY_N(J)*AYZ_V(IJK)&
210                                   -MU_sL_pW*(U_S(IMJPK,L)-U_S(IMJK,L))*ODY_N(J)*AYZ_V(IMJK)
211                              SSY = SSVY
212                              SSZ = MU_sL_pT*(W_S(IJPK,L)-W_S(IJK,L))*AXY_V(IJK)*ODY_N(J)&
213                                   -MU_sL_pB*(W_S(IJPKM,L)-W_S(IJKM,L))*AXY_V(IJKM)*ODY_N(J)
214     !--------------------- End Sources from Stress Term ---------------------
215     
216     
217     !--------------------- Sources from Momentum Source Term ---------------------
218     ! Momentum source associated with the difference in the gradients in
219     ! number density of solids phase m and all other solids phases
220                              D_PL = D_P(IJK,L)
221                              M_PL = (Pi/6d0)* D_PL**3 *RO_S(IJK,L)
222     
223                              NU_PM_pN = ROP_S(IJKN,M)/M_PM
224                              NU_PM_pS = ROP_S(IJK,M)/M_PM
225                              NU_PM_p = AVG_Y(NU_PM_pS,NU_PM_pN,J)
226     
227                              NU_PL_pN = ROP_S(IJKN,L)/M_PL
228                              NU_PL_pS = ROP_S(IJK,L)/M_PL
229                              NU_PL_p = AVG_Y(NU_PL_pS,NU_PL_pN,J)
230     
231                              Fnu_s_p = AVG_Y(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKN,M,L),J)
232                              DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pN-NU_PM_pS)*ODY_N(J)
233                              DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pN-NU_PL_pS)*ODY_N(J)
234                              DS1plusDS2 = DS1 + DS2
235     
236     ! Momentum source associated with the gradient in granular
237     ! temperature of species M
238                              T_PM_pN = Theta_M(IJKN,M)
239                              T_PM_pS = Theta_M(IJK,M)
240     
241                              FT_sM_p = AVG_Y(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKN,M,L),J)
242                              DS3 = FT_sM_p*(T_PM_pN-T_PM_pS)*ODY_N(J)
243     
244     ! Momentum source associated with the gradient in granular
245     ! temperature of species L
246                              T_PL_pN = Theta_M(IJKN,L)
247                              T_PL_pS = Theta_M(IJK,L)
248     
249                              FT_sL_p = AVG_Y(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKN,M,L),J)
250                              DS4 = FT_sL_p*(T_PL_pN-T_PL_pS)*ODY_N(J)
251     !--------------------- End Sources from Momentum Source Term ---------------------
252     
253     
254     ! Add the terms
255                              STRESS_TERMS = STRESS_TERMS + SSVX + SSVY + SSVZ + &
256                                  SSBV + SSX + SSY + SSZ
257                              DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_V(IJK)
258     
259                         ELSE ! if m .ne. L
260     ! for m = l all stress terms should already be handled in existing routines
261     ! for m = l all drag terms should become zero
262                              STRESS_TERMS = STRESS_TERMS + ZERO
263                              DRAG_TERMS = DRAG_TERMS + ZERO
264     
265                         ENDIF ! if m .ne. L
266                    ENDDO     ! over L
267     
268                    KTMOM_V_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
269               ELSE
270                    KTMOM_V_S(IJK,M) = ZERO
271               ENDIF     ! dilute
272           ENDDO        ! over ijk
273     
274           RETURN
275           END SUBROUTINE CALC_IA_MOMSOURCE_V_S
276