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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name:  dufour_coeff(s,mi,alpha,T,ni,rho,v0,
4     !            mu,sigma,chi,beta,zeta0,theta,Ti,Dij,lambdai,gammaij,
5     !            omega,I_ilj,dTl_dnj,dzeta0_dnj,dchi0il_dnj,Dq)
6     !
7     !  author:  C. Hrenya, Feb 2009
8     !
9     !  Purpose: find dufour coefficient according to GHD polydisperse KT
10     !
11     !  Literature/References:
12     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
13     !
14     !  Modifications:
15     !     Sof: removed divisions by ni.
16     !
17     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
18     
19           subroutine dufour_coeff(s,mi,alpha,T,ni,rho,v0, &
20                  mu,sigma,chi,beta,zeta0,theta,Ti,Dij,lambdai,gammaij, &
21                  omega,I_ilj,dTl_dnj,dzeta0_dnj,dchi0il_dnj,Dq)
22           Implicit NONE
23     
24           integer s, indx(s)
25     
26           double precision mi(s),alpha(s,s),T,ni(s),rho,v0,mu(s,s), &
27                           sigma(s,s),chi(s,s),beta(s,s),zeta0,theta(s), &
28                           Ti(s),Dij(s,s),lambdai(s),gammaij(s,s), &
29                           omega(s,s),I_ilj(s,s,s),Dtl_dnj(s,s), &
30                           dzeta0_dnj(s),dchi0il_dnj(s,s,s),Dq(s,s)
31     
32           integer i,j,l,p,kk
33           double precision kronecker(s,s),integ1(s,s,s),integ2(s,s,s), &
34                           integ(s,s,s),sum1(s,s),sum2(s,s), &
35                           dq_bar(s,s),dqlj(s,s), &
36                           Dqkin(s,s),CipjT(s,s,s),Dqcol(s,s), &
37                           Amat(s,s),bmat(s,s), &
38                           Amat0(s,s),bmat0(s)
39           double precision d
40     
41           integer NP
42           parameter (NP=15)     !max no. of linear equations to solve
43     
44           double precision pi
45           parameter (pi=3.14159265458979323846d0)
46     
47           do i=1,s
48              do j=1,s
49                 if (i.eq.j) then
50                    kronecker(i,j) = 1.d0
51                 else
52                    kronecker(i,j) = 0.d0
53                 endif
54              enddo
55           enddo
56     
57     !determination of integral - p 18 of CMH notes (eq. 3) - via summation of two smaller integrals
58     ! Note that the new I_ilj(i,l,j) is now defined as old I_ilj(i,l,j)*ni(l)/ni(j)
59           do i=1,s
60              do l=1,s
61                 do j=1,s
62                    !integral 1 - p 18.1 of CMH notes
63                    integ1(i,l,j) = (kronecker(l,j)*ni(l) + 0.5d0*(ni(l)/chi(i,l) &
64                      *dchi0il_dnj(i,l,j)*ni(l)+I_ilj(i,l,j)*ni(j)))  &
65                      *pi*mi(i)*ni(i)*chi(i,l)*sigma(i,l)**3*mu(l,i) &
66                      *(1.d0+alpha(i,l))  &
67                      *((11.d0*mu(i,l)**2+(13.d0-9.d0*alpha(i,l))*mu(i,l) &
68                      *mu(l,i)+(5.d0+3.d0*alpha(i,l)**2-3.d0*alpha(i,l)) &
69                      *mu(l,i)**2)*Ti(i)**2/mi(i)**2   &
70                      +3.d0*mu(l,i)**2*(1.d0+alpha(i,l))**2*Ti(l)**2 &
71                      /mi(l)**2+(5.d0*mu(i,l)**2+(1.d0-9.d0*alpha(i,l)) &
72                      *mu(i,l)*mu(l,i)+(2.d0+3.d0*alpha(i,l)+6.d0 &
73                      *alpha(i,l)**2)*mu(l,i)**2)*Ti(i)*Ti(l)/mi(i)/mi(l)  &
74                      -5.d0*(Ti(i)/mi(i)+Ti(l)/mi(l))*Ti(i)/mi(i))
75                    !integral 2 - p 18.2 of CMH notes
76                    integ2(i,l,j) = 0.5d0*ni(j)/Ti(l)*dTl_dnj(l,j)  &
77                       *2.d0*pi*ni(i)*ni(l)*mu(i,l)*chi(i,l)*sigma(i,l)**3 &
78                       *Ti(l)*(1.d0+alpha(i,l))   &
79                       *(Ti(i)/mi(i)*(5.d0*(mu(i,l)**2-1.d0)  &
80                       +(1.d0-9.d0*alpha(i,l))*mu(i,l)*mu(l,i)  &
81                       +(2.d0+3.d0*alpha(i,l)+6.d0*alpha(i,l)**2) &
82                       *mu(l,i)**2)  &
83                       +6.d0*Ti(l)/mi(l)*mu(l,i)**2*(1.d0+alpha(i,l))**2)
84                    !summation
85                    integ(i,l,j) = integ1(i,l,j) + integ2(i,l,j)
86                 enddo
87              enddo
88           enddo
89     
90           sum1(1:s,1:s) = 0.d0           !calculate 1st summation used in dq_bar - p 18 CMH notes
91           sum2(1:s,1:s) = 0.d0           !calculate 2nd summation used in dq_bar - p 18 CMH notes
92           do i=1,s
93              do j=1,s
94                 do l=1,s  ! modification to ni(l) > 0 same as in thermal_mobility.f
95                    if(ni(l) > 0d0) sum1(i,j) = sum1(i,j) + mi(l)*(omega(i,l)-zeta0 &
96                               *kronecker(i,l))/ni(l)/Ti(l)*Dij(l,j)
97                    sum2(i,j) = sum2(i,j) + integ(i,l,j)
98                 enddo
99                 dq_bar(i,j) = -2.5d0*ni(j)*Ti(i)**3/mi(i)/T**2   &
100                             *(mi(j)/rho/Ti(i)*sum1(i,j)*ni(i)  &
101                             -0.4d0*mi(i)*T/Ti(i)**3*dzeta0_dnj(j) &           ! ni(i)/ni(i) cancels
102                             *lambdai(i) - dTl_dnj(i,j)/3.d0/Ti(i)**2*ni(i))  &
103                             +sum2(i,j)/3.d0/T**2
104              enddo
105           enddo
106     
107     !find dq,lj (p 18 CMH notes), which will later be used to find the kinetic contribution to Dq
108           do i=1,s
109              do l=1,s
110                 Amat(i,l) = gammaij(i,l)-1.5d0*zeta0*kronecker(i,l)      !A matrix for solution of dq,lj (p 18 CMH notes)
111                 bmat(i,l) = dq_bar(i,l)                                  !b matrix for solution of dq,lj (p 18 CMH notes)
112              enddo
113           enddo
114     
115     ! this extra kk loop and addition of Amat0 and bmat0 is necessary
116     ! since x & b in Ax=b are s by s matrices rather than vectors of length s,
117     ! whereas LUBSKB is specific to x & b vectors of length s
118           do kk=1,s
119              do i=1,s
120                 do j=1,s
121                     Amat0(i,j) = Amat(i,j)
122                 enddo
123              enddo
124              do i=1,s
125                 bmat0(i) = bmat(i,kk)
126              enddo
127     
128              CALL LUDCMP(Amat0, s, NP, indx, d, 'dufour_coeff') ! solve system of s linear equations using
129              CALL LUBKSB(Amat0, s, NP, indx, bmat0) ! LU decomposition
130     
131              do i=1,s
132                 dqlj(i,kk) = bmat0(i)
133              enddo
134           enddo
135     
136     !calculate kinetic contribution to the Dufour coefficient - p18.2 CMH notes
137           do l=1,s
138              do j=1,s
139                 Dqkin(l,j) = dqlj(l,j)   &
140                   +2.5d0/T**2*mi(j)*ni(j)*Ti(l)/rho*Dij(l,j)
141              enddo
142           enddo
143     
144     !evaluate coefficient CipjT needed to evaluate collisional component - p 18.3 CMH notes
145           do i=1,s
146              do p=1,s
147                 do j=1,s
148                    CipjT(i,p,j) = 8.d0*dsqrt(pi)/3.d0*ni(i)*ni(p)*v0**3   &
149                      /dsqrt(theta(i)+theta(p))/(theta(i)*theta(p))**1.5d0  &
150                      *(kronecker(j,p)*beta(i,p)*(theta(i)+theta(p))  &
151                      -0.5d0*theta(i)*theta(p)*(1.d0+(mu(p,i)* &
152                      (theta(i)+theta(p))  &
153                      -2.d0*beta(i,p))/theta(p))*ni(j)/Ti(p)*dTl_dnj(p,j))  &
154                      +2.d0*dsqrt(pi)/3.d0*ni(i)*ni(p)*v0**3* &
155                      (1.d0-alpha(i,p))  &
156                      *(mu(p,i)-mu(i,p))*((theta(i)+theta(p))/theta(i) &
157                      /theta(p))**1.5d0   &
158                      *(kronecker(j,p)+1.5d0*theta(i)/(theta(i)+theta(p))  &
159                      *ni(j)/Ti(p)*dTl_dnj(p,j))
160                 enddo
161              enddo
162           enddo
163     
164     !find collisional contribution - p 18.3 CMH notes
165           Dqcol(1:s,1:s) = 0.d0
166           do i=1,s
167              do j=1,s
168                 do p=1,s
169                    Dqcol(i,j) = Dqcol(i,j)  &
170                      + (1.d0+alpha(i,p))/8.d0*mi(p)*mu(i,p)*sigma(i,p)**3 &
171                      *chi(i,p)*(4.d0*pi/5.d0*(1.d0-alpha(i,p))* &
172                      (mu(i,p)-mu(p,i))  &
173                      *ni(i)*(2.d0/mi(p)*Dqkin(p,j)  &
174                      +5.d0*Ti(i)/T**2*mi(j)*ni(j)/rho/mi(i)*Dij(p,j))  &
175                      +16.d0*pi/5.d0*ni(i)*(2.d0*mu(p,i)/mi(p)*Dqkin(p,j)  &
176                      -5.d0*(2.d0*mu(i,p)-mu(p,i))*Ti(i)/T**2*ni(j)*mi(j) &
177                      /rho/mi(i)*Dij(p,j)) &
178                      -sigma(i,j)/T**2*CipjT(i,p,j))
179                 enddo
180                 !final total Dufour coefficient - p 18 CMH notes
181                 Dq(i,j) = Dqkin(i,j) + Dqcol(i,j)
182              enddo
183           enddo
184     
185           return
186           end subroutine dufour_coeff
187     
188