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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: thermal_conductivity(s,mi,T,sigmai,alpha,ni,rho,v0,mu,
4     !                     sigma,chi,beta,zeta0,theta,Ti,DT,lambda,
5     !                     omega,gamma,lambdai)
6     !
7     !  author:  C. Hrenya, Jan 2009
8     !
9     !  Purpose: find thermal conductivity according to GHD polydisperse KT
10     !
11     !  Literature/References:
12     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
13     !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
15     
16           subroutine thermal_conductivity(s,mi,T,sigmai,alpha,ni,rho,v0,mu, &
17                                     sigma,chi,beta,zeta0,theta,Ti,DT,lambda, &
18                                    omega,gamma,lambdai)
19           Implicit NONE
20     
21           integer s, indx(s)
22     
23           double precision mi(s),sigmai(s),T,alpha(s,s),ni(s),rho,v0, &
24                            mu(s,s),sigma(s,s),chi(s,s),beta(s,s),zeta0, &
25                            theta(s),Ti(s),DT(s),lambda
26     
27           integer i,j
28           double precision E(s,s),A(s,s),C(s,s),F(s,s),sum1(s),sum2(s), &
29                            sum3(s),sum4(s),omega(s,s),lambdai_bar(s), &
30                            gamma(s,s),Amat(s,s),bmat(s),lambdai(s), &
31                            lambdakini(s),lambdakin,lambdacol,CT(s,s)
32           double precision d
33     
34           integer NP
35           parameter (NP=15)     !max no. of linear equations to solve
36           double precision pi
37           parameter (pi=3.14159265458979323846d0)
38     
39     !calculate quantities on p.14-15 of CMH notes
40           do i=1,s
41              do j=1,s
42                 ! CMH edit (9/18/09):  re-introduced mu(j,i)**2 instead of mu(j,i)
43                 ! in first line to correct error
44                 E(i,j) = 2.d0*mu(j,i)**2/theta(i)**2*(theta(i) &
45                         +theta(j))**2  &
46                         *(2.d0*alpha(i,j)**2-3.d0*alpha(i,j)+4.d0) &
47                         *(5.d0*theta(i)+8.d0*theta(j)) &
48                         -mu(j,i)*(theta(i)+theta(j))  &
49                         *(beta(i,j)/theta(i)**2*(5.d0*theta(i)+8.d0*theta(j)) &
50                         *(14.d0*alpha(i,j)-22.d0)-theta(j)/theta(i) &
51                         *(20.d0+3.d0*(15.d0-7.d0*alpha(i,j))+9.d0 &
52                         *(1.d0-alpha(i,j))-28.d0*alpha(i,j))  &
53                         -25.d0*(1d0-alpha(i,j))) +18.d0*(beta(i,j)/theta(i))**2 &
54                         *(5.d0*theta(i)+8.d0*theta(j))+2.d0*beta(i,j)  &
55                         /theta(i)*(25.d0*theta(i)+66.d0*theta(j))+5.d0*theta(j)  &
56                         /theta(i)*(11.d0*theta(i)+6.d0*theta(j))-5.d0*(theta(i) &
57                         +theta(j))*theta(j)/theta(i)**2*(5.d0*theta(i) &
58                         +6.d0*theta(j))
59     
60                 ! CMH edit (9/18/09): corrected error in final line,
61                 ! by making theta(j)/theta(i) instead of its reciprocal
62                 A(i,j) = 5.d0*(2.d0*beta(i,j)+theta(j))+mu(j,i) &
63                         *(theta(i)+theta(j))*(5.d0*(1.d0-alpha(i,j)) &
64                         -(14.d0*alpha(i,j)-22.d0)*beta(i,j)/theta(i)) &
65                         +18.d0*beta(i,j)**2/theta(i)+2.d0*mu(j,i)**2 &
66                         *(2.d0*alpha(i,j)**2-3.d0*alpha(i,j)+4.d0) &
67                         *(theta(i)+theta(j))**2/theta(i) &
68                         -5.d0*theta(j)/theta(i)*(theta(i)+theta(j))
69     
70                 C(i,j) = 5.d0*(2.d0*beta(i,j)-theta(i))+mu(j,i) &
71                          *(theta(i)+theta(j))*(5.d0*(1.d0-alpha(i,j)) &
72                          +(14.d0*alpha(i,j)-22.d0)*beta(i,j)/theta(j)) &
73                          -18.d0*beta(i,j)**2/theta(j)-2.d0*mu(j,i)**2 &
74                          *(2.d0*alpha(i,j)**2-3.d0*alpha(i,j)+4.d0) &
75                          *(theta(i)+theta(j))**2/theta(j) &
76                          +5.d0*(theta(i)+theta(j))
77     
78                 F(i,j) = 2.d0*(mu(j,i)*(theta(i)+theta(j)) &
79                         /theta(j))**2*(2.d0*alpha(i,j)**2-3.d0 &
80                         *alpha(i,j)+4.d0)*(8.d0*theta(i)+5.d0*theta(j)) &
81                         -mu(j,i)*(theta(i)+theta(j))*(beta(i,j) &
82                         /theta(j)**2*(8.d0*theta(i)+5.d0*theta(j)) &
83                         *(14.d0*alpha(i,j)-22.d0)+theta(i)/theta(j) &
84                         *(20.d0+3.d0*(15.d0-7.d0*alpha(i,j))+9.d0 &
85                         *(1.d0-alpha(i,j))-28.d0*alpha(i,j))+25.d0 &
86                         *(1.d0-alpha(i,j)))+18.d0 &
87                         *(beta(i,j)/theta(j))**2*(6.d0*theta(i) &
88                         +5.d0*theta(j)) -2.d0*beta(i,j)/theta(j)*(66.d0 &
89                         *theta(i)+25.d0*theta(j))+5.d0*theta(i)/theta(j) &
90                         *(6.d0*theta(i)+11.d0*theta(j))-5.d0*(theta(i) &
91                         +theta(j))/theta(j)*(6.d0*theta(i)+5.d0*theta(j))
92              enddo
93           enddo
94     
95     ! calculate summations used in gamma(i,i) and omega(i,i) - p. 13 and 16 of CMH notes, respectively
96           do i=1,s
97              sum1(i) = 0.d0     !calculate summation used in gamma(i,i)
98              sum2(i) = 0.d0     !calculate summation used in omega(i,i)
99           enddo
100     
101           do i=1,s
102              do j=1,s
103                 if (j .ne. i) then
104                    sum1(i) = sum1(i)+ni(j)*chi(i,j)*sigma(i,j)**2 &
105                       *v0*mu(j,i)*(1.d0+alpha(i,j))*(theta(i) &
106                       /(theta(j)*(theta(i)+theta(j))))**1.5d0 &
107                       *(E(i,j)-5.d0*(theta(i)+theta(j))/theta(i) &
108                       *A(i,j))
109                    sum2(i) = sum2(i) + ni(j)*chi(i,j)*sigma(i,j)**2 &
110                       *v0*mu(j,i)*(1.d0+alpha(i,j))*dsqrt(theta(i) &
111                       /((theta(i)+theta(j))*theta(j)**3))*A(i,j)
112                 endif
113              enddo
114           enddo
115     
116     ! calculate summations used in lambda_bar - p. 15 of CMH notes
117           do i=1,s
118              sum3(i) = 0.d0     !calculate 1st summation used in b vector (lambda_bar)
119              sum4(i) = 0.d0     !calculate 2nd summation used in b vector (lambda_bar)
120           enddo
121     
122           do i=1,s
123              do j=1,s
124                 !first find omega's (p 16 CMH notes)
125                 if (j .eq. i) then
126                    omega(i,i) = 4.d0/3.d0*dsqrt(pi/2.d0)*sigmai(i)**2 &
127                      *ni(i)*chi(i,i)*v0*(1.d0-alpha(i,i)**2) &
128                      /dsqrt(theta(i)) + 4.d0*dsqrt(pi)/15d0*sum2(i)
129     
130                    if(ni(j)>0d0) sum3(i) = sum3(i)+  &
131                                  (omega(i,j)-zeta0)/(ni(j)*Ti(j))*DT(j)
132                 else
133                    omega(i,j) = 4.d0*dsqrt(pi)/15.d0*ni(j)*chi(i,j) &
134                      *sigma(i,j)**2*v0*mu(j,i)*(1.d0+alpha(i,j))   &
135                      *dsqrt(theta(i)/((theta(i)+theta(j)) &
136                      *theta(j)**3))*C(i,j)
137     
138                    if(ni(j)>0d0) sum3(i) = sum3(i)+  &
139                                  omega(i,j)/(ni(j)*Ti(j))*DT(j)
140                 endif
141                 sum4(i) = sum4(i) + 2.d0*pi*ni(i)*ni(j)*mu(i,j) &
142                      *chi(i,j)*sigma(i,j)**3*Ti(j)*(1.d0+alpha(i,j))   &
143                      *(Ti(i)/mi(i)*(5.d0*(mu(i,j)**2-1.d0)  &
144                      +(1.d0-9.d0*alpha(i,j))*mu(i,j)*mu(j,i)  &
145                      +(2.d0+3.d0*alpha(i,j)+6.d0*alpha(i,j)**2) &
146                      *mu(j,i)**2) +6.d0*Ti(j)/mi(j)*mu(j,i)**2 &
147                      *(1.d0+alpha(i,j))**2)
148              enddo
149           enddo
150     
151     !calculate gamma (p13 of CMH notes) an lambda_bar (p 15 of CMH notes)
152           do i=1,s
153              lambdai_bar(i) = -2.5d0*rho*ni(i)*Ti(i)**3 &
154                 /(Ti(i)*T*mi(i))*sum3(i)+2.5d0*ni(i)*Ti(i)**2 &
155                 /(mi(i)*T)+sum4(i)/(6.d0*T)
156              do j=1,s
157                 if (i .eq. j) then
158                    gamma(i,i) = 16.d0/15.d0*dsqrt(pi)*sigmai(i)**2 &
159                       *ni(i)*chi(i,i)*v0/dsqrt(2.d0*theta(i)) &
160                       *(1.d0+alpha(i,i))*(1.d0+33.d0/16.d0 &
161                       *(1.d0-alpha(i,i)))  &
162                        +2.d0/15.d0*dsqrt(pi)*sum1(i)
163     
164                    Amat(i,j) = gamma(i,j) - 2.d0*zeta0 !A matrix for solution of lambdai (p 13 CMH notes)
165                 else
166                    gamma(i,j) = -2.d0/15.d0*dsqrt(pi)*ni(i)*chi(i,j) &
167                       *sigma(i,j)**2*v0*mu(i,j)*(1.d0+alpha(i,j))   &
168                       *(theta(j)/(theta(i)*(theta(i)+theta(j))))**1.5d0   &
169                       *(F(i,j)+5.d0*((theta(i)+theta(j))/theta(j)) &
170                       *C(i,j))
171     
172                    Amat(i,j) = gamma(i,j)             !A matrix for solution of lambdai (p 13 CMH notes)
173                 endif
174              enddo
175              bmat(i) = lambdai_bar(i)                 !B vector for solution of lambdai (p 13 CMH notes)
176           enddo
177     
178           CALL LUDCMP(Amat, s, NP, indx, d, 'thermal_conductivity')     ! solve system of s linear equations using
179           CALL LUBKSB(Amat, s, NP, indx, bmat)  ! LU decomposition
180     
181           do i=1,s
182              lambdai(i) = bmat(i)
183           enddo
184     
185           lambdakin = 0             !kinetic contribution to conductivity (p 16 CMH notes)
186           do i = 1,s
187              lambdakini(i) = lambdai(i) + 2.5d0/T*rho*Ti(i)*DT(i)/mi(i)
188              lambdakin = lambdakin + lambdakini(i)
189           enddo
190     
191           lambdacol = 0d0            !collisional contribution to conductivity (p 17 CMH notes)
192           do i=1,s
193              do j=1,s
194                 CT(i,j) = -4.d0/3.d0*dsqrt(pi)*ni(i)*ni(j)*v0**3 &
195                    /dsqrt((theta(i)+theta(j))*(theta(i)*theta(j))**3)  &
196                    *(2.d0*beta(i,j)**2+theta(i)*theta(j)  &
197                    + (theta(i)+theta(j)) *((theta(i)+theta(j))*mu(i,j) &
198                    *mu(j,i)+beta(i,j)*(1.d0+mu(j,i))))-dsqrt(pi)*ni(i) &
199                    *ni(j)*v0**3*(1.d0-alpha(i,j))*(mu(j,i)-mu(i,j))  &
200                    *((theta(i)+theta(j))/theta(i)/theta(j))**1.5d0 &
201                    *(mu(j,i)+beta(i,j)/(theta(i)+theta(j)))
202                  lambdacol = lambdacol + (1.d0+alpha(i,j))/8.d0*mi(j) &
203                    *mu(i,j)*sigma(i,j)**3*chi(i,j)  *(4.d0*pi/5.d0 &
204                    *(1.d0-alpha(i,j))*(mu(i,j)-mu(j,i))*ni(i)  &
205                    *(2.d0/mi(j)*lambdakini(j)+5.d0*Ti(i)/(mi(i)*mi(j) &
206                    *T)*rho*DT(j)) +48.d0*pi/15.d0*ni(i)*(2.d0*mu(i,j) &
207                    /mi(j)*lambdakini(j)-5.d0*Ti(i)/(mi(i)*mi(j)*T) &
208                    *(2.d0*mu(i,j)-mu(j,i))*rho*DT(j))  &
209                    -sigma(i,j)*CT(i,j)/T)
210              enddo
211           enddo
212     
213           lambda = lambdakin + lambdacol                  !conductivity (p 13 CMH notes)
214     
215           return
216           end subroutine thermal_conductivity
217