File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/shear_viscosity.f
1
2
3
4
5
6
7
8
9
10
11
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)
32 double precision pi
33 parameter (pi=3.14159265458979323846d0)
34
35 do i=1,s
36 sum1(i) = 0.d0
37 (i) = 0.d0
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
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
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)
77 endif
78
79 enddo
80 bmat(i) = ni(i)*Ti(i)+sum2(i)
81 enddo
82
83 CALL LUDCMP(Amat, s, NP, indx, d, 'shear_viscosity')
84 CALL LUBKSB(Amat, s, NP, indx, bmat)
85
86 = 0.d0
87 do i=1,s
88 etajk(i) = bmat(i)
89 etakin = etakin + etajk(i)
90 enddo
91
92 etacol = 0
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
102
103 return
104 end subroutine shear_viscosity
105
106