File: RELATIVE:/../../../mfix.git/model/GhdTheory/thermal_diffusivity.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: thermal_diffusivity(s,alpha,ni,mi,rho,v0,mu,sigma,chi,
4     !                                       zeta0,theta,Ti,p,DT,nu)
5     !
6     !  author:  C. Hrenya, Jan 2009
7     !
8     !  Purpose: find thermal diffusivity according to GHD polydisperse KT
9     !
10     !  Literature/References:
11     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
12     !
13     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
14     
15           subroutine thermal_diffusivity(s,alpha,ni,mi,rho,v0,mu,sigma, &
16                                          chi,zeta0,theta,Ti,p,DT,nu)
17     
18           Implicit NONE
19     
20           integer s, indx(s)
21     
22           double precision alpha(s,s),ni(s),mi(s),rho,v0,mu(s,s), &
23                            sigma(s,s),chi(s,s),zeta0,theta(s), &
24                            Ti(s),p,DT(s)
25     
26           integer i,j
27           double precision sum1(s),sum2(s),nu(s,s),Amat(s,s),bmat(s)
28           double precision d
29     
30           integer NP
31           parameter (NP=15)     !max no. of linear equations to solve
32           double precision pi
33           parameter (pi=3.14159265458979323846d0)
34     
35           do i=1,s
36              sum1(i) = 0.d0     !calculate summation used in nu(i,i) - p 7 CMH notes
37              sum2(i) = 0.d0     !calculate summation used in b vector - p 7 CMH notes
38           enddo
39     
40           do i=1,s
41              do j=1,s
42                 if (i .ne. j) then
43                    sum1(i) = sum1(i) + ni(j)*sigma(i,j)**2*chi(i,j) &
44                         *mu(j,i)*v0*(1.d0+alpha(i,j)) &
45                         *dsqrt((theta(i)+theta(j))/(theta(i)*theta(j)))
46                 endif
47                 sum2(i) = sum2(i) + ni(j)*mu(i,j)*chi(i,j)* &
48                           sigma(i,j)**3*Ti(j)*(1.d0+alpha(i,j))
49              enddo
50           enddo
51     
52           do i=1,s
53              do j=1,s
54                 if (i .eq. j) then
55                    nu(i,i) = 4.d0*dsqrt(pi)/3.d0*sum1(i)
56     
57                    Amat(i,j) = nu(i,j)-zeta0                    !A matrix for solution of DT (p 7 CMH notes)
58                 else
59                    nu(i,j) = -4.d0*dsqrt(pi)/3.d0*ni(i)*sigma(i,j)**2 &
60                       *chi(i,j)*mu(i,j)*v0*(1.d0+alpha(i,j))   &
61                       *dsqrt((theta(i)+theta(j))/(theta(i)*theta(j)))
62     
63                    Amat(i,j) = nu(i,j)                          !A matrix for solution of DT (p 7 CMH notes)
64                 endif
65              enddo
66              bmat(i) = -p*ni(i)*mi(i)/rho**2*(1.d0-rho*Ti(i) &  !b matrix for solution of DT (p 7 CMH notes)
67                           /(mi(i)*p))+2d0*pi/3d0*ni(i)/rho*sum2(i)
68           enddo
69     
70           CALL LUDCMP(Amat, s, NP, indx, d, 'thermal_diffusivity')     ! solve system of s linear equations using
71           CALL LUBKSB(Amat, s, NP, indx, bmat)  ! LU decomposition
72     
73           do i=1,s
74              DT(i) = bmat(i)
75           enddo
76     
77           return
78           end subroutine thermal_diffusivity
79     
80