c******************************************************************************* c This is a subroutine that caculates the forward and backward rates of reaction c given species(mass fraction), temperature and pressure (thermodynamic). The c array(data_out) is 2xNumReactions long. c c This subroutine originally came from Habib Najm's dflame code. c Programmers : c H.N. Najm (SNL, hnnajm@ca.sandia.gov). c S. Lefantzi (SNL, slefant@ca.sandia.gov) c J. Lee (SNL, jclee@ca.sandia.gov) c c This subroutine is called from Thermoanalysis.cpp. It is not called from any other c fortran subroutines in the ChemicalRates component distribution. c c October 30 : Jerry Lee found an oversight - while the reaction rates were being c calculated, the pressure falloffs were not being applied. About 200 c below (look for "added by jclee") Jerry added the pressure falloff. c Other changes are prefaced by "jclee". The code seems to be correct c and is being checked in. c c******************************************************************************* c ... evaluate prod rate of sp ksp subroutine getrfrb(pressure,temperature,nspecies, 1 ysp,data_out) IMPLICIT NONE include 'parameter.par' include 'chem.com' include 'chem_ref.com' integer k, nspecies integer ksp,isp,ip3 integer nyp, nnn double precision one,ctot,xik double precision crprod,cpprod,cvp real*8 ysp(nspecies) !Y(species) real*8 data_out(2*numreac) !Return forward and backward reaction rates real*8 cgspl(nsmx),qfac(idim),rr_f(idim),rr_r(idim) real*8 cspl(nsmx) ! jclee Oct.13.2004 double precision pressure, temperature double precision temp,otc,vlntemp double precision sum,V,rho_k,rho integer kr,kthb,ispthb,kfl,kfc,ifp double precision vt,tc,vot,cthb, 1 rklow,pr,prlog,fc,xp,fcent,fclog,xn, 2 cprlog,flog,clog,pcor,rop c press is dimensional pressure [N/m^2] and temp is non-dimensional temperature #ifdef NON_DIM temp = temperature #else temp = temperature / tempref #endif c------------------------------------------------------------------------------- c Evaluate the cell density of the mixture: c 1. rho(mix) = m(mix) / V(mix) c 2. m(mix) = 1 kgr c 3. V(mix) = V(1) + V(2) +...+ V(ns), for species 1..ns c 4. V(k) = m(k) / rho(k) c 5. m(k) = ys(k,i,j) * m(mix) = ys(k,i,j) c 6. rho(k) = pressure / ( rg(k)*temp(i,j)*tempref ) c 7. rg(k) = r0 / wt V = 0.0 sum = 0.0 do k = 1,nspecies rg(k) = r0 / (1000.0 * wt(k)) #ifdef NON_DIM rho_k= pressure*pref / (1000.0 * rg(k) * temperature * tempref) #else rho_k= pressure / (1000.0 * rg(k) * temperature ) #endif V = ysp(k) / rho_k sum = sum + V end do ! end of k=1,nspecies loop c non-dimensionalize rho rho = (1.d0 / sum) / rhoref c------------------------------------------------------------------------------- c find the cspl's call y2c(rho,ysp,cspl) c--------------------------------------------------------------------------- c assign the input data tc = temp otc = 1.0/temp vlntemp = log(temp) c------------------------------------------------------------------------------- c ... prelude one = 1.0 call cgsp(temp,cgspl) ! gibbs free energies ctot = 0.0 do ksp=1,ns ctot = ctot + cspl(ksp) end do do k=1,numreac ! loop on all reaction steps; jclee Oct.13.2004 kr = k ! jclee Oct.13.2004 c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... initialize qfac array, product of concentrations to be used in c ... third body and fall off reactions qfac(k) = 1.0d0 c------------------------------------------------------------------------------- c ... evaluate forward reaction rate, rr_f if(itfdep(k).ne.0)then rr_f(k) = capb(k)*exp(par(2,k)*vlntemp - & cape(k)*otc) else rr_f(k) = capb(k) endif c------------------------------------------------------------------------------- c ... nspec(k) > 0 iff reaction is reversible c ... xik = SUM_i=1,ns nu(i,k)*G_ik c ... evaluate (1/equilibrium_reaction_rate) = orr_e c ... evaluate reverse reaction rate, rr_r c ... kflagrev(k) = 1 -> reverse Arrhenius rate parameters are defined c ... use these rather than orr_e to find rr_r if(nspec(k).gt.0)then if(kflagrev(k).eq.1)then if(itrdep(k).ne.0)then rr_r(k) = rcapb(k)*exp(rexpo(k)*vlntemp-rcape(k)*otc) else rr_r(k) = rcapb(k) endif else xik = 0.0 do isp=1,nreac(k) xik = xik + anu(isp,k)*cgspl(nunk(isp,k)) end do do ip3=4,nprod(k)+3 xik = xik + anu(ip3,k)*cgspl(nunk(ip3,k)) end do rr_r(k) = rr_f(k)*exp(xik*otc+anusum(k)*vlntemp) endif else rr_r(k) = 0.0d0 ! reaction irreversible-> rr_r=0 endif c------------------------------------------------------------------------------- c ... find reactants concentrations product PROD(c_j**nu_jk^prime)_j=1,ns c ... it is assumed below that nyp is a positive integer (non-zero) c ... which should be true in general crprod = 1.0d0 do isp = 1,nreac(k) nyp = - nu(isp,k) cvp = cspl(nunk(isp,k)) do nnn=1,nyp crprod = crprod * cvp end do end do c------------------------------------------------------------------------------- c ... for reverse reaction, c ... find products concentrations product PROD(c_j**nu_jk^double-prime)_j=1,ns c ... nspec(k) > 0 if reaction is reversible c ... it is assumed below that nu(ip3,k) is a positive integer (non-zero) c ... which should be true in general if(nspec(k).gt.0)then cpprod = 1.d0 do ip3 = 4,nprod(k)+3 cvp = cspl(nunk(ip3,k)) do nnn=1,nu(ip3,k) cpprod = cpprod * cvp end do end do else cpprod = 0.0d0 endif c------------------------------------------------------------------------------- c ... evaluate reaction rate-of-progress rr_f(k) = rr_f(k)*crprod rr_r(k) = -rr_r(k)*cpprod c------------------------------------------------------------------------------- c---------------------------added by jclee for the thrid-body reaction and Fall off reactions c---------------------------copied from asp rop_pack.f-------------------------- c ... evaluate third-body-scaled reaction rates-of-progress c ... nb. qfac is the product of concentrations, with possible 3rd c ... body efficiencies. c ... if the reaction is also a fall-off reaction with kfal=0, then c ... qfac is just used to find pr below, and then reset to 1.0 c ... otherwise, it is actually used to scale the rate of progress c ... later on below if(nthb.gt.0)then kthb = jthb(kr) if(kthb.gt.0)then qfac(kr) = ctot do ispthb = 1,ntbs(kthb) ksp = nktb(ispthb,kthb) qfac(kr) = qfac(kr)+(aik(ispthb,kthb)-1.d0)*cspl(ksp) end do endif endif c------------------------------------------------------------------------------- c ... evaluate fall-off pressure corrections c ... vt is dimensional temperature (K), vot = 1/vt c ... these are used in the computation of the dimensionless F if(nfal.gt.0)then kfl = jfal(kr) if(kfl.gt.0)then vt = tc*tempref vot = 1.0d0/vt kfc = kfal(kfl) ifp = ifop(kfl) if(kfc.eq.0)then cthb = qfac(kr) qfac(kr)= 1.0d0 else cthb = cspl(kfc) endif rklow = capbl(kr)*exp(fpar(2,kfl)*vlntemp-capel(kr)*otc) pr = rklow * cthb / rr_f(kr)*crprod ! jclee Oct.13.2004 prlog = log10(pr) if(ifp.eq.1)then ! Lindemann form fc = 1.0d0 else ! ifop = 2, 3, or 4 if(ifp.eq.2)then ! SRI form xp = 1.0/(1.0+prlog*prlog) fc = ( + (fpar(4,kfl) * exp(-fpar(5,kfl)*vot) + + exp(-vt/fpar(6,kfl)) + )**xp + ) * fpar(7,kfl) * vt**fpar(8,kfl) else ! 6 or 7-param TROE form fcent = (1.0-fpar(4,kfl))*exp(-vt/fpar(5,kfl)) + + fpar(4,kfl)*exp(-vt/fpar(6,kfl)) if(ifp.eq.4)then ! 7-parameter TROE form fcent = fcent + exp(-fpar(7,kfl)*vot) endif fclog = log10(fcent) xn = 0.75 - 1.27*fclog cprlog = prlog - (0.4+0.67*fclog) clog = cprlog/(xn-0.14*cprlog) flog = fclog/(1.0+clog*clog) fc = 10.0**flog endif endif pcor = fc*pr/(1.0+pr) rr_f(kr) = rr_f(kr) * pcor rr_r(kr) = rr_r(kr) * pcor endif endif c------------------------------------------------------------------------------- c ... modify ropl_f and ropl_r due to 3rd bodies, given possible mod to qfac bec falloff if(ithb_type(k) .ne. 0)then rr_f(k) = qfac(k) * rr_f(k) rr_r(k) = qfac(k) * rr_r(k) endif c------------------END of modification by jerry for the third-body reaction and Fall off reactions end do !end of reactions loop on k c ... assign the output data do k=1,2*numreac if (k.le.numreac) then data_out(k) = rr_f(k) else data_out(k) = rr_r(k-numreac) end if end do c------------------------------------------------------------------------------- return end c*******************************************************************************