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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_IA_MOMSOURCE_V_S                                   C
4     !  Purpose: Determine source terms for V_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_V_S(M)
19     
20     !-----------------------------------------------
21     ! Modules
22     !-----------------------------------------------
23           USE param1, only: zero
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 :: I, J, K, IJK, IJKN, JP, IM,  KM, IJPK, IJMK,&
62                      IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
63                      IJKTN, IJKB, IJKBN, IJKM, IJPKM,&
64                      IPJK, IJKP
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
70     ! Bulk viscosity values
71           DOUBLE PRECISION :: XI_sL_pN, XI_sL_pS, LAMBDA_sL_pN, LAMBDA_sL_pS
72     ! Variables for drag calculations
73           DOUBLE PRECISION :: M_PM, M_PL, D_PM, D_PL, NU_PM_pN, NU_PM_pS, NU_PM_p, &
74                               NU_PL_pN, NU_PL_pS, NU_PL_p, T_PM_pN, T_PM_pS, &
75                               T_PL_pN, T_PL_pS, Fnu_s_p, FT_sM_p, FT_sL_p
76     ! Average volume fraction
77           DOUBLE PRECISION :: EPSA
78     ! Source terms (Surface)
79           DOUBLE PRECISION :: ssvx, ssvy, ssvz, ssx, ssy, ssz, ssbv
80     ! Source terms (Volumetric)
81           DOUBLE PRECISION :: DS1, DS2, DS3, DS4, DS1plusDS2
82     !-----------------------------------------------
83     
84     ! section largely based on tau_v_g:
85     
86           DO IJK = IJKSTART3, IJKEND3
87     
88     ! Skip walls where some values are undefined.
89              IF(WALL_AT(IJK)) cycle
90     
91               D_PM = D_P(IJK,M)
92               M_PM = (Pi/6d0) * D_PM**3 *RO_S(IJK,M)
93     
94               J = J_OF(IJK)
95               IJKN = NORTH_OF(IJK)
96               EPSA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
97               IF ( .NOT.SIP_AT_N(IJK) .AND. EPSA>DIL_EP_S) THEN
98     
99                    JP = JP1(J)
100                    I = I_OF(IJK)
101                    IM = IM1(I)
102                    K = K_OF(IJK)
103                    KM = KM1(K)
104                    IJPK = JP_OF(IJK)
105                    IJMK = JM_OF(IJK)
106                    IMJK = IM_OF(IJK)
107                    IJKM = KM_OF(IJK)
108                    IMJPK = IM_OF(IJPK)
109                    IJPKM = JP_OF(IJKM)
110     
111                    IJKW = WEST_OF(IJK)
112                    IJKE = EAST_OF(IJK)
113                    IJKNE = EAST_OF(IJKN)
114                    IJKNW = NORTH_OF(IJKW)
115     
116                    IJKB = BOTTOM_OF(IJK)
117                    IJKT = TOP_OF(IJK)
118                    IJKTN = NORTH_OF(IJKT)
119                    IJKBN = NORTH_OF(IJKB)
120     
121     ! additional required quantities:
122                    IPJK = IP_OF(IJK)
123                    IJKP = KP_OF(IJK)
124     
125     ! initialize variable
126                    STRESS_TERMS = ZERO
127                    DRAG_TERMS = ZERO
128     
129                    DO L = 1, MMAX
130                         IF (M .ne. L) THEN
131     
132     !--------------------- Sources from Stress Terms ---------------------
133     ! Surface Forces
134     ! standard shear stress terms (i.e. ~diffusion)
135                              MU_sL_pE = AVG_Y_H(AVG_X_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKE,M,L),I),&
136                                   AVG_X_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKNE,M,L),I),J)
137                              MU_sL_pW = AVG_Y_H(AVG_X_H(MU_sL_ip(IJKW,M,L),MU_sL_ip(IJK,M,L),IM),&
138                                   AVG_X_H(MU_sL_ip(IJKNW,M,L),MU_sL_ip(IJKN,M,L),IM),J)
139                              SSVX = MU_sL_pE*(V_S(IPJK,L)-V_S(IJK,L))*AYZ_V(IJK)*ODX_E(I)&
140                                   -MU_sL_pW*(V_S(IJK,L)-V_S(IMJK,L))*AYZ_V(IMJK)*ODX_E(IM)
141     
142                              MU_sL_pN = MU_sL_ip(IJKN,M,L)
143                              MU_sL_pS = MU_sL_ip(IJK,M,L)
144                              SSVY = MU_sL_pN*(V_S(IJPK,L)-V_S(IJK,L))*ODY(JP)*AXZ_V(IJK)&
145                                   -MU_sL_pS*(V_S(IJK,L)-V_S(IJMK,L))*ODY(J)*AXZ_V(IJMK)
146     
147                              MU_sL_pT = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJK,M,L),MU_sL_ip(IJKT,M,L),K),&
148                                   AVG_Z_H(MU_sL_ip(IJKN,M,L),MU_sL_ip(IJKTN,M,L),K),J)
149                              MU_sL_pB = AVG_Y_H(AVG_Z_H(MU_sL_ip(IJKB,M,L),MU_sL_ip(IJK,M,L),KM),&
150                                   AVG_Z_H(MU_sL_ip(IJKBN,M,L),MU_sL_ip(IJKN,M,L),KM),J)
151                              SSVZ = MU_sL_pT*(V_S(IJKP,L)-V_S(IJK,L))*AXY_V(IJK)*ODZ_T(K)*OX(I)&
152                                   -MU_sL_pB*(V_S(IJK,L)-V_S(IJKM,L))*AXY_V(IJKM)*ODZ_T(KM)*OX(I)
153     
154     ! bulk viscosity term
155                              XI_sL_pN = XI_sL_ip(IJKN,M,L)
156                              XI_sL_pS = XI_sL_ip(IJK,M,L)
157                              LAMBDA_sL_pN = -(2.d0/3.d0)*MU_sL_pN + XI_sL_pN
158                              LAMBDA_sL_pS = -(2.d0/3.d0)*MU_sL_pS + XI_sL_pS
159                              SSBV = (LAMBDA_sL_pN*TRD_S(IJKN,L)-LAMBDA_sL_pS*TRD_S(IJK,L))*AXZ(IJK)
160     
161     ! off diagonal shear stress terms
162                              SSX = MU_sL_pE*(U_S(IJPK,L)-U_S(IJK,L))*ODY_N(J)*AYZ_V(IJK)&
163                                   -MU_sL_pW*(U_S(IMJPK,L)-U_S(IMJK,L))*ODY_N(J)*AYZ_V(IMJK)
164                              SSY = SSVY
165                              SSZ = MU_sL_pT*(W_S(IJPK,L)-W_S(IJK,L))*AXY_V(IJK)*ODY_N(J)&
166                                   -MU_sL_pB*(W_S(IJPKM,L)-W_S(IJKM,L))*AXY_V(IJKM)*ODY_N(J)
167     !--------------------- End Sources from Stress Term ---------------------
168     
169     
170     !--------------------- Sources from Momentum Source Term ---------------------
171     ! Momentum source associated with the difference in the gradients in
172     ! number density of solids phase m and all other solids phases
173                              D_PL = D_P(IJK,L)
174                              M_PL = (Pi/6d0)* D_PL**3 *RO_S(IJK,L)
175     
176                              NU_PM_pN = ROP_S(IJKN,M)/M_PM
177                              NU_PM_pS = ROP_S(IJK,M)/M_PM
178                              NU_PM_p = AVG_Y(NU_PM_pS,NU_PM_pN,J)
179     
180                              NU_PL_pN = ROP_S(IJKN,L)/M_PL
181                              NU_PL_pS = ROP_S(IJK,L)/M_PL
182                              NU_PL_p = AVG_Y(NU_PL_pS,NU_PL_pN,J)
183     
184                              Fnu_s_p = AVG_Y(Fnu_s_ip(IJK,M,L),Fnu_s_ip(IJKN,M,L),J)
185                              DS1 = Fnu_s_p*NU_PL_p*(NU_PM_pN-NU_PM_pS)*ODY_N(J)
186                              DS2 = -Fnu_s_p*NU_PM_p*(NU_PL_pN-NU_PL_pS)*ODY_N(J)
187                              DS1plusDS2 = DS1 + DS2
188     
189     ! Momentum source associated with the gradient in granular
190     ! temperature of species M
191                              T_PM_pN = Theta_M(IJKN,M)
192                              T_PM_pS = Theta_M(IJK,M)
193     
194                              FT_sM_p = AVG_Y(FT_sM_ip(IJK,M,L),FT_sM_ip(IJKN,M,L),J)
195                              DS3 = FT_sM_p*(T_PM_pN-T_PM_pS)*ODY_N(J)
196     
197     ! Momentum source associated with the gradient in granular
198     ! temperature of species L
199                              T_PL_pN = Theta_M(IJKN,L)
200                              T_PL_pS = Theta_M(IJK,L)
201     
202                              FT_sL_p = AVG_Y(FT_sL_ip(IJK,M,L),FT_sL_ip(IJKN,M,L),J)
203                              DS4 = FT_sL_p*(T_PL_pN-T_PL_pS)*ODY_N(J)
204     !--------------------- End Sources from Momentum Source Term ---------------------
205     
206     
207     ! Add the terms
208                              STRESS_TERMS = STRESS_TERMS + SSVX + SSVY + SSVZ + &
209                                  SSBV + SSX + SSY + SSZ
210                              DRAG_TERMS = DRAG_TERMS + (DS1plusDS2+DS3+DS4)*VOL_V(IJK)
211     
212                         ELSE ! if m .ne. L
213     ! for m = l all stress terms should already be handled in existing routines
214     ! for m = l all drag terms should become zero
215                              STRESS_TERMS = STRESS_TERMS + ZERO
216                              DRAG_TERMS = DRAG_TERMS + ZERO
217     
218                         ENDIF ! if m .ne. L
219                    ENDDO     ! over L
220     
221                    KTMOM_V_S(IJK,M) = STRESS_TERMS + DRAG_TERMS
222               ELSE
223                    KTMOM_V_S(IJK,M) = ZERO
224               ENDIF     ! dilute
225           ENDDO        ! over ijk
226     
227           RETURN
228           END SUBROUTINE CALC_IA_MOMSOURCE_V_S
229