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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  Subroutine name: ghd_model
4     !
5     !  Purpose: find transport coefficients of GHD polydisperse KT
6     !           for known inputs.
7     !
8     !  Literature / References
9     !     C. Hrenya handnotes and Garzo, Hrenya, Dufty papers (PRE, 2007)
10     !
11     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
12     
13           SUBROUTINE  GHD_MODEL (S, SIGMAI,IJK, alpha, MI, phii, T, Zeta0, &
14                                  zetau, Ti, P, Kappa, Eta, DT, DF, Lambda, &
15                                  Lij, Dij, Dq)
16     
17     !-----------------------------------------------
18     ! Modules
19     !-----------------------------------------------
20           USE drag
21           Implicit NONE
22     !-----------------------------------------------
23     ! Local parameters
24     !-----------------------------------------------
25           double precision :: pi
26           parameter (pi=3.14159265458979323846d0)
27     !-----------------------------------------------
28     ! Dummy arguments
29     !-----------------------------------------------
30           integer :: s                           !number of species
31           double precision sigmai(s)             !diameter of each species
32           integer :: ijk                         !fluid cell index
33           double precision :: alpha(s,s)         !restitution coeff of each pair
34           double precision :: mi(s)              !mass of each species
35           double precision :: phii(s)            !solids volume fraction of each species
36           double precision :: T                  !granular (mixture) temperature
37           double precision :: zeta0              !zeta0 is cooling rate
38           double precision :: zetau              !cooling rate transport coefficient (1st order)
39           double precision :: Ti(s)              !species temperature
40           double precision :: p                  !pressure
41           double precision :: kappa              !bulk viscosity
42           double precision :: eta                !shear viscosity
43           double precision :: DT(s)              !thermal diffusivity
44           double precision :: DF(s,s)            !mass mobility coefficient
45           double precision :: lambda             !thermal conductivity
46           double precision :: Lij(s,s)           !thermal mobility
47           double precision :: Dij(s,s)           !ordinary diffusion
48           double precision :: Dq(s,s)            !Dufour coefficient
49     !-----------------------------------------------
50     ! Local variables
51     !-----------------------------------------------
52           integer i, j, k
53           double precision ni(s)                 !number density of each species
54           double precision phi                   !total solids fraction
55           double precision n                     !total number density
56           double precision rhoi(s)               !mi(i)*ni(i) of each species
57           double precision rho                   !m*n
58           double precision m                     !average mass
59           double precision v0                    !commonly used quantities: v0, mu, sigma
60           double precision mu(s,s)
61           double precision sigma(s,s)
62           double precision chi_ij
63           double precision chi(s,s)              !radial distribution function at contact for each pair
64           double precision group1(s,s), group2(s,s)
65           double precision theta(s)
66           double precision beta(s,s)             !defined on p11 CMH notes
67           double precision Beta_tot, niTi, scale_fac
68           double precision nu(s,s)               !commonly used quantity: p7 cmh notes
69           double precision I_ilj(s,s,s)          !commonly used quantity: p6.1 CMH notes
70           double precision dTl_dnj(s,s)          !partial derivative of Tl wrt nj: p 6 CMH notes
71           double precision dzeta0_dnj(s)         !partial derivative of zeta 0 wrt nj: p 6 CMH notes
72           double precision dchi0il_dnj(s,s,s)    !partial derivative of chi_il wrt nj: p 6 CMH notes
73           double precision :: omega(s,s)         !commonly used quantity: p16 CMH notes
74           double precision :: gammaij(s,s)       !commonly used quantity: p13 CMH notes
75           double precision :: lambdai(s)         !commonly used quantity: p13 CMH notes
76     !-----------------------------------------------
77     
78     
79     ! Initializing
80           phi = 0.d0
81           n   = 0.d0
82           rho = 0.d0
83           m   = 0.d0
84           niTi = 0.d0
85           Beta_tot = 0.d0
86           do i=1,s
87              phi = phi + phii(i)
88              ni(i) = 6.d0*phii(i) / (pi*sigmai(i)**3)
89              n = n+ ni(i)
90              rhoi(i) = mi(i)*ni(i)
91              rho = rho + rhoi(i)
92              m = m + mi(i)
93              Beta_tot = Beta_tot + F_GS(IJK,s)
94           enddo
95           do i=1,s
96             if(n .eq. 0) then  ! do not do any calculation if total solids concentration is zero.
97               Ti(:)      = T
98               zeta0      = 0.d0
99               zetau      = 0.d0
100               p          = 0.d0
101               kappa      = 0.d0
102               eta        = 0.d0
103               DT(:)      = 0.d0
104               DF(:,:)    = 0.d0
105               lambda     = 0.d0
106               Lij(:,:)   = 0.d0
107               Dij(:,:)   = 0.d0
108               RETURN
109             endif
110           enddo
111           m = m/dble(s)
112           v0 = dsqrt(2.d0*T/m)
113     
114           do i=1,s
115              do j=1,s
116                 mu(i,j) = mi(i)/(mi(i)+mi(j))
117                 sigma(i,j) = 0.5d0*(sigmai(i)+sigmai(j))
118                 call chi_ij_GHD(s,i,j,sigmai,phi,ni,chi_ij)
119                 chi(i,j) = chi_ij
120              enddo
121           enddo
122     !---------------------------------------------------------------------
123     ! ZEROTH-ORDER COOLING RATE (zeta0)
124     !---------------------------------------------------------------------
125     ! This subroutine solves the nonlinear algebraic equations for theta
126           call cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi,theta)
127     
128     ! Ti and zeta0 are calculated from theta
129           do i=1,s
130              Ti(i) = mi(i)*T/(m*theta(i))
131              niTi = niTi + ni(i)*Ti(i)
132           enddo
133     
134           do j=1,s
135                group1(1,j)  = chi(1,j)*ni(j)*mu(j,1)* &
136                                 sigma(1,j)**2*(1d0+alpha(1,j))
137                group2(1,j)  = mu(j,1)/2.d0*(1d0+alpha(1,j))
138           enddo
139           zeta0 = 0d0;
140           do k=1,s
141              zeta0 = zeta0 + group1(1,k)*dsqrt((theta(1)+theta(k)) &
142                            /(theta(1)*theta(k))) &
143                     * ((1.d0-group2(1,k)*(theta(1)+theta(k))/theta(k)))
144           enddo
145           zeta0 = 8.d0/3.d0*dsqrt(2.d0*pi*T/m)*zeta0
146     
147     ! this quantity is needed for determination of several other transport coefficients
148           do i=1,s
149              do j=1,s
150                 beta(i,j) = mu(i,j)*theta(j)-mu(j,i)*theta(i);  !p 11 CMH notes
151              enddo
152           enddo
153     
154     !---------------------------------------------------------------------
155     ! COOLING RATE TRANSPORT COEFFICIENT (zetaU)
156     !---------------------------------------------------------------------
157           call cooling_rate_tc(s,mi,sigmai,alpha,ni,n,v0,mu,sigma,chi,T, &
158                               zeta0,theta,Ti,zetau)
159     !---------------------------------------------------------------------
160     ! PRESSURE (p)
161     !---------------------------------------------------------------------
162           call pressure (s,alpha,ni,n,mu,sigma,chi,T,Ti,p)
163     
164     !---------------------------------------------------------------------
165     ! BULK VISCOSITY (kappa)
166     !---------------------------------------------------------------------
167           call bulk_viscosity(s,mi,alpha,ni,v0,mu,sigma,chi,theta,kappa)
168     
169     !---------------------------------------------------------------------
170     ! SHEAR VISCOSITY (eta)
171     !---------------------------------------------------------------------
172           call shear_viscosity(s,mi,sigmai,alpha,ni,v0,mu,sigma,chi, &
173                                beta,zeta0,theta,Ti,kappa,eta)
174     
175     !---------------------------------------------------------------------
176     ! THERMAL DIFFUSION COEFFICIENT (DT) & nu
177     !---------------------------------------------------------------------
178           call thermal_diffusivity(s,alpha,ni,mi,rho,v0,mu,sigma,chi, &
179                                    zeta0,theta,Ti,p,DT,nu)
180     
181     !---------------------------------------------------------------------
182     ! MASS MOBILITY COEFFICIENT (DF)
183     !---------------------------------------------------------------------
184           call mass_mobility(s,mi,ni,rho,zeta0,theta,nu,DF)
185     
186     !---------------------------------------------------------------------
187     ! ORDINARY DIFFUSION (Dij)
188     !---------------------------------------------------------------------
189           call ordinary_diff(s,mi,sigmai,alpha,phii,T,phi,ni,n,rhoi,rho, &
190                  m,mu,sigma,chi,zeta0,theta,Ti,DT,nu,Dij,I_ilj,dTl_dnj, &
191                  dzeta0_dnj,dchi0il_dnj)
192     
193     !---------------------------------------------------------------------
194     ! CONDUCTIVITY (lambda)
195     !---------------------------------------------------------------------
196           call thermal_conductivity(s,mi,T,sigmai,alpha,ni,rho,v0,mu,&
197                              sigma,chi,beta,zeta0,theta,Ti,DT,&
198                              lambda,omega,gammaij,lambdai)
199     
200     !---------------------------------------------------------------------
201     ! THERMAL MOBILITY COEFFICIENT (Lij)
202     !---------------------------------------------------------------------
203           call thermal_mobility(s,mi,alpha,ni,mu,sigma,chi,zeta0,&
204                                 theta,Ti,DF,gammaij,omega,Lij)
205     
206     !---------------------------------------------------------------------
207     ! DUFOUR COEFFICIENT (Dq)
208     !---------------------------------------------------------------------
209           call dufour_coeff(s,mi,alpha,T,ni,rho,v0, &
210                  mu,sigma,chi,beta,zeta0,theta,Ti,Dij,lambdai,gammaij, &
211                  omega,I_ilj,dTl_dnj,dzeta0_dnj,dchi0il_dnj,Dq)
212     
213     
214     ! SCALE ALL TRANSPORT COEFFICIENTS FOR THE LOW NUMBER DENSITY LIMIT
215           scale_fac = 1d0+2d0*eta*Beta_tot/(rho*niTi)
216     
217     
218           lambda = lambda/scale_fac
219           eta = eta/scale_fac
220           kappa = kappa/scale_fac
221           do i=1,s
222               DT(i) = DT(i)/scale_fac
223             do j=1,s
224               Lij(i,j) = Lij(i,j)/scale_fac
225               Dq(i,j) = Dq(i,j)/scale_fac
226               Dij(i,j) = Dij(i,j)/scale_fac
227     
228     ! do not scale mass mobility as fluxes need be computed correctly in dilute limit.
229     !          DF(i,j) = DF(i,j)/scale_fac
230             enddo
231           enddo
232     
233           RETURN
234           END SUBROUTINE GHD_MODEL
235     
236