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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: bulk_viscosity(s,mi,sigmai,alpha,ni,v0,mu,sigma,chi,
4     !                                  beta,zeta0,theta,Ti,kappa,eta)
5     !
6     !  author:  C. Hrenya, Jan 2009
7     !
8     !  Purpose: find shear viscosity 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 shear_viscosity(s,mi,sigmai,alpha,ni,v0,mu,sigma,chi, &
16                                      beta,zeta0,theta,Ti,kappa,eta)
17           Implicit NONE
18     
19           integer s, indx(s)
20     
21           double precision mi(s),sigmai(s),alpha(s,s),ni(s),v0,mu(s,s), &
22                            sigma(s,s),chi(s,s),beta(s,s),zeta0,theta(s), &
23                            Ti(s),kappa,eta
24     
25           integer i,j
26           double precision sum1(s),sum2(s),tau(s,s),Amat(s,s),bmat(s), &
27                            etajk(s),etakin,etacol
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 tau(i,i) - p 10 CMH notes
37              sum2(i) = 0.d0     !calculate summation used in b vector - p 10 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)*chi(i,j)*sigma(i,j)**2* &
44                      mu(j,i)*(1.d0+alpha(i,j))*theta(i)**1.5d0 &
45                      /dsqrt(theta(j)) * (6.d0*beta(i,j)/ &
46                      (theta(i)**2*dsqrt(theta(i)+theta(j)))  &
47                      + (9.d0-3.d0*alpha(i,j))/2.d0*mu(j,i)/theta(i)**2* &
48                      dsqrt(theta(i)+theta(j))   &
49                      + 5.d0/(theta(i)*dsqrt(theta(i)+theta(j))))
50                endif
51                sum2(i) = sum2(i) + 2.d0*pi/15.d0*mi(i)*ni(i)*ni(j)*mu(j,i)* &
52                  sigma(i,j)**3*chi(i,j)*(1d0+alpha(i,j)) *(mu(j,i)* &
53                 (3.d0*alpha(i,j)-1.d0)*(Ti(i)/mi(i)+Ti(j)/mi(j)) &
54                 -4.d0*(Ti(i)-Ti(j))/(mi(i)+mi(j)))
55              enddo
56           enddo
57     
58     ! calculate tau (p 10-11 CMH notes)
59           do i=1,s
60              do j=1,s
61                 if (i .eq. j) then
62                    tau(i,i) = 4.d0*dsqrt(pi)/15.d0*v0 *(ni(i)*sigmai(i)**2* &
63                      chi(i,i)/dsqrt(2.d0*theta(i))*(9.d0-3.d0*alpha(i,i))* &
64                      (1.d0+alpha(i,i)) + 2.d0*sum1(i))
65     
66                    Amat(i,j) = tau(i,j) - 0.5d0*zeta0     !A matrix for solution of etajk (p 10 CMH notes)
67                 else
68                    tau(i,j) = 8.d0*dsqrt(pi)/15.d0*v0   &
69                      * ni(i)*chi(i,j)*sigma(i,j)**2*mu(i,j)*theta(j)**1.5d0 &
70                      /dsqrt(theta(i))*(1.d0+alpha(i,j))  * (6.d0*beta(i,j)/ &
71                      (theta(j)**2*dsqrt(theta(i)+theta(j)))   &
72                      + (9.d0-3.d0*alpha(i,j))/2.d0*mu(j,i)/theta(j)**2* &
73                      dsqrt(theta(i)+theta(j))   &
74                      - 5.d0/(theta(j)*dsqrt(theta(i)+theta(j))))
75     
76                    Amat(i,j) = tau(i,j)                 !A matrix for solution of etajk (p 10 CMH notes)
77                 endif
78     !
79              enddo
80              bmat(i) = ni(i)*Ti(i)+sum2(i)              !b matrix for solution of etajk (p 10 CMH notes)
81           enddo
82     
83           CALL LUDCMP(Amat, s, NP, indx, d, 'shear_viscosity')     ! solve system of s linear equations using
84           CALL LUBKSB(Amat, s, NP, indx, bmat)  ! LU decomposition
85     
86           etakin = 0.d0
87           do i=1,s
88              etajk(i) = bmat(i)
89              etakin = etakin + etajk(i)     !kinetic shear viscosity (p 11 CMH notes)
90           enddo
91     
92           etacol = 0                        !collisional shear viscosity (p 11 CMH notes)
93           do i=1,s
94              do j=1,s
95                 etacol = etacol + 4.d0*(pi)/15.d0*ni(j)*sigma(i,j)**3 &
96                             *chi(i,j)*mu(j,i)*(1d0+alpha(i,j))*etajk(i)
97              enddo
98           enddo
99           etacol = etacol + 0.6d0*kappa
100     
101           eta = etakin + etacol           !total shear viscosity (p 11 CMH notes)
102     
103           return
104           end subroutine shear_viscosity
105     
106