subroutine drfdfk(ns,xmol,wmix,t,olmass,drfdm,drtdv,xml2, & sig,epsm,delm,chim,vstm,rstm,error) C ** routine to return mixture diffusion coefficinets ** C ** & thermal diffusion ratio per Kihara ** C ** C SNL / phpaul all rights reserved 1997 ** C last modified Oct/20/97 ** C **************************************************************** C inputs: C NS - number of species C XMOL - NS long vector containing molecular weights C WMIX - ns long vector containing mixture fractions C OLMASS - mixture molecular weight C T - temperature (K) C *** - set of molecular parameter arrays from DRINIT C returns for primary species: C DRFDM - NS long vector of mixture diffusion coefficients C (cm*cm/sec) as reduced to one atmosphere pressure C DRTDV - NS long thermal diffusion ratio (dimensionless) C per Kihara model and sum rule of Chapman and Cowling C **************************************************************** IMPLICIT NONE #include "defs.h" include 'drfmcm.for' C passed arrays logical error real*8 sig(KD),epsm(KD,KD),delm(KD,KD),chim(KD,KD) real*8 vstm(KD,KD),rstm(KD,KD),xmol(KD) real*8 wmix(KD),drfdm(KD),drtdv(KD) real*8 xml2(kd,kd) C internal working arrays real*8 sgx(KD),og2(KD),xs11(KD,KD),fmix(KD) real*8 dcor(KD,KD),alp(KD,KD) integer j, ns, k, jj, kk double precision delta, zz, t, zz01, zz02 double precision TR, VR, RR, SGR, sq2, cndif, o6, og1j double precision og2j, obj, ocj, oej, og1jk, og2jk, objk double precision ocjk, oejk, cncon, olmass, sum1, sum2 C square root of 2 data sq2/1.4142135623730951d0/ c delta DATA DELTA/1.0d-14/ C DELTA = a small number which is added to the mixture fractions C to prevent divide by zero under pure species conditions DATA CNDIF /2.6280D-3/ C CNDIF = scaling constant for Dij at 1 atm. in cm^2/s C for single species c ... o6 = 1./6. data o6/0.166666666666666667d0/ do 20 j = 1,ns fmix(j) = wmix(j)+delta ZZ = 1.D0 + CHIM(j,j) + DELM(j,j)/T zz01 = dexp(o6*dlog(zz)) zz02 = zz*zz TR = T/(EPSM(j,j)*zz02) VR = VSTM(j,j)*zz02 RR = RSTM(j,j)*zz01 SGR = sig(j)/zz01 CALL DRXIOM12(TR,RR,VR,OG1J,OG2J,OBJ,OCJ,OEJ,ERROR) sgx(j) = sgr * sgr xs11(j,j) = og1j * sgx(j) #ifndef FIRST_ORDER_DRFM og2(j) = og2j hn01 = 6.d0*ocj-5.d0 dcor(j,j) = 1.3d0*(hn01*hn01)/58.d0 alp(j,j) = 0.d0 #endif 20 continue C for the binary pairs do 40 jj = 2,ns do 30 kk = 1,jj-1 C generate primatives ZZ = 1.0D0 + CHIM(kk,jj) + DELM(kk,jj)/T zz01 = dexp(o6*dlog(zz)) zz02 = zz*zz TR = T/(EPSM(kk,jj)*zz02) VR = VSTM(kk,jj)*zz02 RR = RSTM(kk,jj)*zz01 SGR = 0.5d0*(sig(jj)+sig(kk))/zz01 CALL DRXIOM12(TR,RR,VR,OG1JK,OG2JK,OBJK,OCJK,OEJK,ERROR) xs11(kk,jj) = og1jk*sgr*sgr xs11(jj,kk) = xs11(kk,jj) #ifndef FIRST_ORDER_DRFM astar = og2jk/og1jk bstar = 5.d0-2.4d0*objk cstar = 6.d0*ocjk - 5.d0 C order according to moelcular weights Mj > Mk if(xmol(jj).gt.xmol(kk))then j = jj k = kk ff = 1.0d0 else j = kk k = jj ff = -1.0d0 endif C Kihara correction to the binary diffusion coefficient y = xmol(k)/xmol(j) cjk = fmix(j)/(fmix(j)+fmix(k)) an01 = 1.0d0+1.8d0*y aw = sq2*og1jk/(8.0d0*og2(k)*(an01*an01)) bw = 10.0d0*(1.0d0+1.8d0*y+3.0d0*y*y) dcor(jj,kk) = 1.3d0*(cstar*cstar)* & aw*cjk/(1.d0+(aw*bw-1.d0)*cjk) dcor(kk,jj) = dcor(jj,kk) #endif C thermal diffusion factor per Kihara C z = 1.d0/y C zp1 = z + 1.d0 C zm1 = z - 1.d0 C zr = zm1/zp1 C ckj =fmix(k)/(fmix(j)+fmix(k)) C hpj = sgx(j)*og2(j)/xs11(jj,kk) C hpk = sgx(k)*og2(k)/xs11(jj,kk) C s1 = z*DSQRT(2.d0/zp1)*hpj + C & (7.5d0*zm1 - 4.d0*z*astar) / (zp1*zp1) C q1 = (2.d0/zp1)*DSQRT(2.d0/zp1)*hpj* C & (0.5d0*bstar*z*z+3.d0+1.6d0*z*astar) C q12=15.d0*(zr*zr)*(0.5d0*bstar) C & + 4.d0*z*astar*(6.d0+bstar)/(zp1*zp1) C & + 1.6d0*zp1*hpj*hpk/DSQRT(z) C s2 = DSQRT(2.d0/(z+z*z))*hpk + C & (7.5d0*(z-z*z) - 4.d0*z*astar) / (zp1*zp1) C q2 = (2.d0/zp1)*DSQRT(2.d0/(z+z*z))*hpk* C & ((0.5d0*bstar)+3.d0*z*z+1.6d0*z*astar) C alp(jj,kk) = ff*cstar*(cjk*s1 - ckj*s2) / C & (cjk*cjk*q1 + ckj*ckj*q2 + cjk*ckj*q12) C alp(kk,jj) = (-1.0d0)*alp(jj,kk) 30 continue 40 continue C compute mixture diffusion coefficients & thermal diffusion ratio cncon = cndif*t*dsqrt(t)/olmass do 80 j = 1,ns sum1 = 0.d0 sum2 = 0.d0 C sum3 = 0.d0 if(j.gt.1)then do k = 1,j-1 sum1 = sum1 + fmix(k)*xmol(k) sum2 = sum2 + fmix(k)*xml2(k,j)*xs11(k,j) #ifndef FIRST_ORDER_DRFM & /(1.0d0+dcor(j,k)) #endif C sum3 = sum3 + wmix(k)*alp(j,k) end do endif if(j.lt.ns)then do k = j+1,ns sum1 = sum1 + fmix(k)*xmol(k) sum2 = sum2 + fmix(k)*xml2(k,j)*xs11(k,j) #ifndef FIRST_ORDER_DRFM & /(1.0d0+dcor(j,k)) #endif C sum3 = sum3 + wmix(k)*alp(j,k) end do endif drfdm(j) = cncon * sum1 / sum2 C drtdv(j) = wmix(j)*sum3 80 continue return end c******************************************************************************** SUBROUTINE drxiom12(TR,RS,VS,OG11,OG22,OBX,OCX,OEX,ERROR) C** Collision integrals using LJ for TR<10 and BM for TR>10 ** C** coefs. from Bzowski et al. JPCRD v19 p1179 (1990). ** C** Range for 0.2 10.0 A10 = DLOG(0.1D0*VS) AA = DLOG(VS/TR) C type I AT = 1.D0/(TR*TR) gn01 = g1(3)/a10 gn02 = g2(3)/a10 gn03 = g3(3)/a10 B1 = G1(1)+G1(2)/A10+gn01*gn01 B2 = (G2(1)+G2(2)/A10+gn02*gn02)*(-1.D3) B3 = (G3(1)+G3(2)/A10+gn03*gn03)*(1.D5) SUM1 = RS*RS*(BT(1)+AT*(BT(2)+AT*(BT(3)+AT*BT(4)))) SUM2 = AT*(B1+AT*(B2+AT*B3))/(A10*A10) OG11 = (SUM1 + SUM2) * (AA*AA) #ifndef FIRST_ORDER_DRFM SUM3 = RS*RS*AT*(BT(2)+AT*(2.D0*BT(3)+3.D0*AT*BT(4))) SUM4=AT*(B1+AT*(2.D0*B2+AT*3.D0*B3))/(A10*A10) HP1 = (SUM3+SUM4)/(SUM1+SUM2) OCX = 1.D0-((1.D0/AA)+ HP1 )*2.D0/3.D0 SUM5=RS*RS*AT*(2.D0*BT(2)+AT*(8.D0*BT(3)+18.D0*AT*BT(4))) SUM6=AT*(2.D0*B1+AT*(8.D0*B2+AT*18.D0*B3))/(A10*A10) HP2 = (SUM5+SUM6)/(SUM1+SUM2) +2.D0*(HP1*HP1) OBX = 4.D0*OCX-3.D0*(OCX*OCX)+((1.D0/(AA*AA))-HP2)*2.D0/3.D0 C type II AT = 1.D0/DLOG(TR) fn01 = f1(3)/a10 fn02 = f2(3)/a10 fn03 = f3(3)/a10 B1 = F1(1)+F1(2)/A10+fn01*fn01 B2 = (F2(1)+F2(2)/A10+fn02*fn02)*(-1.D0) B3 = (F3(1)+F3(2)/A10+fn03*fn03) SUM1=RS*RS*(DT(1)+AT*AT*(DT(3)+AT*(DT(4)+AT*DT(5)))) SUM2 = AT*AT*(B1+AT*(B2+AT*B3))/(A10*A10) OG22 = (SUM1 + SUM2) * (AA*AA) SUM3 = RS*RS*AT*AT*AT*(2.D0*DT(3)+AT*(3.D0*DT(4)+ & 4.D0*AT*DT(5))) SUM4 = AT*AT*AT*(2.D0*B1+AT*(3.D0*B2+AT*4.D0*B3))/(A10*A10) OEX = 1.D0-( 2.D0/AA + (SUM3+SUM4)/(SUM1+SUM2) )/4.D0 #endif RETURN END