File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/ghd.f
1
2
3
4
5
6
7
8
9
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
19
20 USE drag
21 Implicit NONE
22
23
24
25 double precision :: pi
26 parameter (pi=3.14159265458979323846d0)
27
28
29
30 integer :: s
31 double precision sigmai(s)
32 integer :: ijk
33 double precision :: alpha(s,s)
34 double precision :: mi(s)
35 double precision :: phii(s)
36 double precision :: T
37 double precision :: zeta0
38 double precision :: zetau
39 double precision :: Ti(s)
40 double precision :: p
41 double precision :: kappa
42 double precision :: eta
43 double precision :: DT(s)
44 double precision :: DF(s,s)
45 double precision :: lambda
46 double precision :: Lij(s,s)
47 double precision :: Dij(s,s)
48 double precision :: Dq(s,s)
49
50
51
52 integer i, j, k
53 double precision ni(s)
54 double precision phi
55 double precision n
56 double precision rhoi(s)
57 double precision rho
58 double precision m
59 double precision v0
60 double precision mu(s,s)
61 double precision sigma(s,s)
62 double precision chi_ij
63 double precision chi(s,s)
64 double precision group1(s,s), group2(s,s)
65 double precision theta(s)
66 double precision beta(s,s)
67 double precision Beta_tot, niTi, scale_fac
68 double precision nu(s,s)
69 double precision I_ilj(s,s,s)
70 double precision dTl_dnj(s,s)
71 double precision dzeta0_dnj(s)
72 double precision dchi0il_dnj(s,s,s)
73 double precision :: omega(s,s)
74 double precision :: gammaij(s,s)
75 double precision :: lambdai(s)
76
77
78
79
80 = 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
97 (:) = 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
124
125
126 call cooling_rate(s,mi,ni,n,m,Ti,T,chi,sigmai,alpha,rhoi,theta)
127
128
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
148 do i=1,s
149 do j=1,s
150 beta(i,j) = mu(i,j)*theta(j)-mu(j,i)*theta(i);
151 enddo
152 enddo
153
154
155
156
157 call cooling_rate_tc(s,mi,sigmai,alpha,ni,n,v0,mu,sigma,chi,T, &
158 zeta0,theta,Ti,zetau)
159
160
161
162 call pressure (s,alpha,ni,n,mu,sigma,chi,T,Ti,p)
163
164
165
166
167 call bulk_viscosity(s,mi,alpha,ni,v0,mu,sigma,chi,theta,kappa)
168
169
170
171
172 call shear_viscosity(s,mi,sigmai,alpha,ni,v0,mu,sigma,chi, &
173 beta,zeta0,theta,Ti,kappa,eta)
174
175
176
177
178 call thermal_diffusivity(s,alpha,ni,mi,rho,v0,mu,sigma,chi, &
179 zeta0,theta,Ti,p,DT,nu)
180
181
182
183
184 call mass_mobility(s,mi,ni,rho,zeta0,theta,nu,DF)
185
186
187
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
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
202
203 call thermal_mobility(s,mi,alpha,ni,mu,sigma,chi,zeta0,&
204 theta,Ti,DF,gammaij,omega,Lij)
205
206
207
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
215 = 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
229
230 enddo
231 enddo
232
233 RETURN
234 END SUBROUTINE GHD_MODEL
235
236