File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/thermal_conductivity.f
1
2
3
4
5
6
7
8
9
10
11
12
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)
36 double precision pi
37 parameter (pi=3.14159265458979323846d0)
38
39
40 do i=1,s
41 do j=1,s
42
43
44 (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
61
62 (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
96 do i=1,s
97 sum1(i) = 0.d0
98 (i) = 0.d0
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
117 do i=1,s
118 sum3(i) = 0.d0
119 (i) = 0.d0
120 enddo
121
122 do i=1,s
123 do j=1,s
124
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
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
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)
173 endif
174 enddo
175 bmat(i) = lambdai_bar(i)
176 enddo
177
178 CALL LUDCMP(Amat, s, NP, indx, d, 'thermal_conductivity')
179 CALL LUBKSB(Amat, s, NP, indx, bmat)
180
181 do i=1,s
182 lambdai(i) = bmat(i)
183 enddo
184
185 lambdakin = 0
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
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
214
215 return
216 end subroutine thermal_conductivity
217