MFIX  2016-1
ghd.f
Go to the documentation of this file.
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 
dufour_coeff
subroutine dufour_coeff(s, mi, alpha, T, ni, rho, v0, mu, sigma, chi, beta, zeta0, theta, Ti, Dij, lambdai, gammaij, omega, I_ilj, dTl_dnj, dzeta0_dnj, dchi0il_dnj, Dq)
Definition: dufour_coeff.f:22
ghd_model
subroutine ghd_model(S, SIGMAI, IJK, alpha, MI, phii, T, Zeta0, zetau, Ti, P, Kappa, Eta, DT, DF, Lambda, Lij, Dij, Dq)
Definition: ghd.f:16
cooling_rate_tc
subroutine cooling_rate_tc(s, mi, sigmai, alpha, ni, n, v0, mu, sigma, chi, T, zeta0, theta, Ti, zetau)
Definition: cooling_rate_tc.f:18
drag
Definition: drag_mod.f:11
cooling_rate
subroutine cooling_rate(s, mi, ni, n, m, Ti, T, chi, sigmai, alpha, rhoi, xve
Definition: cooling_rate.f:15
ordinary_diff
subroutine ordinary_diff(s, mi, sigmai, alpha, phii, T, phi, ni, n, rhoi, rho, m, mu, sigma, chi, zeta0, theta, Ti, DT, nu, Dij, I_ilj, dTl_dnj, dzeta0_dnj, dchi0il_dnj)
Definition: ordinary_diff.f:19
chi_ij_ghd
subroutine chi_ij_ghd(s, i, j, sigmai, phi, ni, chi_ij)
Definition: chi_ij_GHD.f:2
thermal_conductivity
subroutine thermal_conductivity(s, mi, T, sigmai, alpha, ni, rho, v0, mu,
Definition: thermal_conductivity.f:17
shear_viscosity
subroutine shear_viscosity(s, mi, sigmai, alpha, ni, v0, mu, sigma, chi, beta, zeta0, theta, Ti, kappa, eta)
Definition: shear_viscosity.f:17
thermal_diffusivity
subroutine thermal_diffusivity(s, alpha, ni, mi, rho, v0, mu, sigma, chi, zeta0, theta, Ti, p, DT, nu)
Definition: thermal_diffusivity.f:17
thermal_mobility
subroutine thermal_mobility(s, mi, alpha, ni, mu, sigma, chi, zeta0, theta, Ti, DF, gammaij, omega, Lij)
Definition: thermal_mobility.f:20
pressure
subroutine pressure(s, alpha, ni, n, mu, sigma, chi, T, Ti, p)
Definition: pressure.f:15
drag::f_gs
double precision, dimension(:,:), allocatable f_gs
Definition: drag_mod.f:14
bulk_viscosity
subroutine bulk_viscosity(s, mi, alpha, ni, v0, mu, sigma, chi, theta, kappa)
Definition: bulk_viscosity.f:16
mass_mobility
subroutine mass_mobility(s, mi, ni, rho, zeta0, theta, nu, DF)
Definition: mass_mobility.f:15