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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: cooling_rate_tc(s,mi,sigmai,alpha,ni,n,v0,mu,sigma,chi,T,
4     !                                   zeta0,theta,Ti,zetau)
5     !
6     !  author:  C. Hrenya, Feb 2009
7     !
8     !  Purpose: find transport coefficient for cooling rate according to
9     !           GHD polydisperse KT
10     !
11     !  Literature/References:
12     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
13     !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
15     
16           subroutine cooling_rate_tc(s,mi,sigmai,alpha,ni,n,v0,mu,sigma, &
17                                     chi,T,zeta0,theta,Ti,zetau)
18           Implicit NONE
19     
20           integer s, indx(s)
21     
22           double precision mi(s),sigmai(s),alpha(s,s),ni(s),n,v0,mu(s,s), &
23                           sigma(s,s),chi(s,s),T,zeta0,theta(s), &
24                           Ti(s),zetau
25     
26           integer i,j
27           double precision kronecker(s,s),sum1(s),psi(s,s),eid_bar(s), &
28                           Amat(s,s),bmat(s),ejd(s),zeta10,zeta11
29           double precision d
30     
31           integer NP
32           parameter (NP=15)     !max no. of linear equations to solve
33     
34           double precision pi
35           parameter (pi=3.14159265458979323846d0)
36     
37           do i=1,s
38              do j=1,s
39                 if (i.eq.j) then
40                    kronecker(i,j) = 1.d0
41                 else
42                    kronecker(i,j) = 0.d0
43                 endif
44              enddo
45           enddo
46     
47           do i=1,s
48              sum1(i) = 0.d0     !calculate summation used in psi(i,i) - p 5 CMH notes
49           enddo
50           do i=1,s
51              do j=1,s
52                 if (i.ne.j) then
53                    sum1(i) = sum1(i)+dsqrt(pi)/2.d0*ni(j)*chi(i,j)* &
54                     sigma(i,j)**2*v0*mu(j,i)*(1.d0+alpha(i,j))/ &
55                     dsqrt(theta(i)*theta(j))/(theta(i)+theta(j))**2.5d0  &
56                     *(-2.d0*(90.d0*theta(j)**3+231.d0*theta(i)* &
57                     theta(j)**2+184.d0*theta(i)**2*theta(j)+40.d0* &
58                     theta(i)**3)+3.d0*mu(j,i)*(1.d0+alpha(i,j))* &
59                     (theta(i)+theta(j))*(70.d0*theta(j)**2+117.d0* &
60                     theta(i)*theta(j)+44.d0*theta(i)**2)-24.d0* &
61                     mu(j,i)**2*(1.d0+alpha(i,j))**2* &
62                     (theta(i)+theta(j))**2*(5.d0*theta(j)+4.d0*theta(i))  &
63                     +30.d0*mu(j,i)**3*(1.d0+alpha(i,j))**3*(theta(i)+ &
64                     theta(j))**3+5.d0*theta(j)*(theta(i)+theta(j))* &
65                     (2.d0*(4.d0*theta(i)+3.d0*theta(j))-3.d0*mu(j,i)* &
66                     (1.d0+alpha(i,j))*(theta(i)+theta(j))))
67                 endif
68              enddo
69           enddo
70     
71     ! find psi (p 5 CMH notes) and eid_bar (p 4) values
72           do i=1,s
73              eid_bar(i) = 0.d0
74           enddo
75           do i=1,s
76              do j=1,s
77                 if (i.eq.j) then
78                    psi(i,i) = -2.d0/15.d0*(sum1(i)+dsqrt(2.d0*pi) &
79                      /32.d0*ni(i)*chi(i,i)*sigmai(i)**2*(1.d0+ &
80                      alpha(i,i))/dsqrt(theta(i))*v0*(30.d0*alpha(i,i)**3 &
81                      -126.d0*alpha(i,i)**2+177.d0*alpha(i,i)+48.d0* &
82                      (3.d0*alpha(i,i)-7.d0)-137.d0)+3.d0/32.d0* &
83                      dsqrt(2.d0*pi)*ni(i)*chi(i,i)*sigmai(i)**2* &
84                      (1.d0+alpha(i,i))/dsqrt(theta(i))*v0  &
85                      *(10.d0*alpha(i,i)**3+22.d0*alpha(i,i)**2+ &
86                      11.d0*alpha(i,i)-3.d0))
87                  else
88                     psi(i,j) = -2.d0/15.d0*(dsqrt(pi)/2.d0* &
89                       ni(j)*chi(i,j)*sigma(i,j)**2*v0*mu(j,i)*(1.d0+ &
90                       alpha(i,j))*(theta(i)/theta(j))**1.5d0/(theta(i)+ &
91                       theta(j))**2.5d0*((2.d0*theta(j)+5.d0*theta(i))* &
92                       (2.d0*theta(j)+3.d0*mu(j,i)*(1.d0+alpha(i,j))  &
93                       *(theta(i)+theta(j)))-24.d0*mu(j,i)**2*(1.d0+ &
94                       alpha(i,j))**2*(theta(i)+theta(j))**2+30.d0* &
95                       mu(j,i)**3*(1.d0+alpha(i,j))**3*(theta(i)+ &
96                       theta(j))**3/theta(j)-5.d0*(theta(i)+theta(j))* &
97                       (2.d0*theta(j)+3.d0*mu(j,i)*(1.d0+alpha(i,j))* &
98                       (theta(i)+theta(j)))))
99                  endif
100                  eid_bar(i) = eid_bar(i) - pi/6.d0*ni(j)*chi(i,j)* &
101                       sigma(i,j)**3*mu(j,i)*(1.d0+alpha(i,j))   &
102                       *(40.d0*(mu(i,j)-1.d0)+4.d0*(19.d0+9.d0* &
103                       alpha(i,j))*mu(j,i)-48.d0*mu(j,i)**2*(theta(i)+ &
104                       theta(j))/theta(j)*(1.d0+alpha(i,j))**2   &
105                       +15.d0*mu(j,i)**3*(theta(i)+theta(j))**2/ &
106                       theta(j)**2*(1.d0+alpha(i,j))**3)
107              enddo
108           enddo
109     
110           do i=1,s
111              do j=1,s
112                 Amat(i,j) = psi(i,j) - 1.5d0*zeta0*kronecker(i,j)     !A matrix for solution of ejd (p 4 CMH notes)
113              enddo
114              bmat(i) = 2.d0*eid_bar(i)/15.d0                  !b matrix for solution of ejd (p 4 CMH notes)
115           enddo
116     
117           CALL LUDCMP(Amat, s, NP, indx, d, 'cooling_rate_tc')     ! solve system of s linear equations using
118           CALL LUBKSB(Amat, s, NP, indx, bmat)  ! LU decomposition
119     
120           do i=1,s
121              ejd(i) = bmat(i)
122           enddo
123     
124     ! determine summations for cooling rate transport coefficients (p 4 CMH notes)
125           zeta10 = 0.d0
126           zeta11 = 0.d0
127           do i=1,s
128              do j=1,s
129                 zeta10 = zeta10 - 2.d0*pi/(3.d0*n*T)*ni(i)*ni(j)*mu(j,i) &
130                   *sigma(i,j)**3*chi(i,j)*(1.d0-alpha(i,j)**2)*Ti(i)
131                 zeta11 = zeta11 + dsqrt(pi)*v0**3/(2.d0*n*T)*ejd(i)* &
132                   ni(i)*ni(j)*sigma(i,j)**2*chi(i,j)*mi(j)  &
133                   *mu(i,j)*(1.d0-alpha(i,j)**2)*dsqrt(theta(j)/ &
134                   (theta(i)**3*(theta(i)+theta(j))))
135              enddo
136           enddo
137     
138           zetau = zeta10 + zeta11         !transport coefficient for cooling rate (p 4 CMH notes)
139     
140           return
141           end subroutine cooling_rate_tc
142     
143