c------------------------------------------------------------------------------- c ... set chem ref quantities subroutine chem_ref(crefl,wsprefl,wrefl,temprefl,r0l,wtl, 1 patml,boltzkl,avogadrol,cprefl,hrefl, 2 prefl,rgrefl,rhorefl,capdrefl, 3 htcndrefl,caplrefl,velrefl,vdiffrefl,amurefl, 4 tymrefl,cgrefl,grefl,damkohlrl) IMPLICIT NONE include 'parameter.par' include 'chem_ref.com' include 'chem.com' double precision r0l,patml,boltzkl,temprefl,prefl, 1 wrefl,rgrefl,rhorefl,wsprefl,capdrefl, 2 caplrefl,vdiffrefl,amurefl,tymrefl, 3 crefl,cgrefl,grefl,crefcgms,wsprefcg, 4 wrefcgms,avogadrol,cprefl,hrefl, 5 htcndrefl,velrefl,c1,c2,ctcon,eps0, 6 sig0,qradref,one,damkohlrl integer k,kthb,kfl,krev real*8 wtl(nsmx) c========================================================================== c Fill in chem_ref.com with reference values generated from c ref_gen.f r0 = r0l patm = patml boltzk = boltzkl avogadro = avogadrol tempref = temprefl cpref = cprefl href = hrefl pref = prefl wref = wrefl rgref = rgrefl rhoref = rhorefl wspref = wsprefl capdref = capdrefl htcndref = htcndrefl caplref = caplrefl velref = velrefl vdiffref = vdiffrefl amuref = amurefl tymref = tymrefl cref = crefl cgref = cgrefl gref = grefl damkohlr = damkohlrl c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... note. Ideally, I would like to keep all ref quantities in SI units. c ... however, capbref (or B_ref) is complicated. c ... chem.bin generated by chemkin gives me the dimensional B for each c ... reaction in cm-gmole-sec units c ... rather the m-kgmole-sec units I use above. c ... further, the exact units of B depend on the stoichiometry and number c ... of reactants of each reaction. Therefore, instead of converting B c ... fore each reaction to SI upon reading it in gckwk and then dividing c ... it by the SI B_ref here, it is easier to simply evaluate B_ref for c ... each reaction in cm-gmole-sec units. This is not too obtrusive c ... since B_ref is not used elsewhere c ... this is implemented by converting the cref and wspref used below c ... to cm-gmole-sec before using them crefcgms = cref*0.001 ! convert from kgmole/m^3 to gmole/cm^3 wsprefcg = wspref*0.001 ! convert fr kg/m^3.s to gm/cm^3.s wrefcgms = wref ! kg/kgmole = gm/gmole do k=1,numreac c1 = tempref**par(2,k) ! T_ref^{alpha_k} c2 = crefcgms**nursum(k) capbref(k)=wsprefcg/(wrefcgms*c1*c2) ! Ref Pre-exp Constant, B_ref(k) ! units depend on nursum(k) ! and par(2,k) end do c ... scale capbref(k) accordingly, if k is a 3rd body reaction c ... nb. if k is also a kfal=0 falloff reaction, its capbref c ... should NOT be scaled. It will be reset to the unscaled c ... valued in the following code block below for fall off reactions if(nthb.gt.0)then do kthb = 1,nthb k = ithb(kthb) capbref(k) = capbref(k)/crefcgms end do endif c ... fall-off reactions c ... nb, two possibilities c ... (1) (+M) reaction dependent on concentrations of all species c ... reaction is a formal third body reaction c ... i.e. it belongs to the set of ithb(kthb) reactions c ... but, the dep. on concentrations will be absorbed in pr c ... and klow, and will not multiply rop, hence reset capbref c ... by capbref = capbref *crefcgms c ... this case has kfal=0 c ... (2) (+particular_species) reaction dependent on concentration of c ... particular_species only c ... reaction is NOT a formal third body reaction, c ... it does NOT belong to the set of ithb(kthb) reactions c ... so do not reset its capbref, because it has not been changed c ... anyway c ... this case has kfal=1 c ... c ... and, for either case, the Bref for the low regime c ... capblref, is scaled using crefcgms and Tref as below. if(nfal.gt.0)then do kfl = 1, nfal k = ifal(kfl) ! reaction index if(kfal(kfl).eq.0)then capbref(k) = capbref(k)*crefcgms ! 3rd body with caplref endif ctcon = tempref**(par(2,k)-fpar(2,kfl)) capblref(k) = (capbref(k)/crefcgms)*ctcon end do endif c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... non-dimensional reaction rate data c ... par(1,k) is the dimensional B_k c ... par(3,k) is the dimensional E_k/R0 (units: K) do k=1,numreac capb(k) = par(1,k)/capbref(k) cape(k) = par(3,k)/tempref end do if(nfal.gt.0)then do kfl = 1, nfal k = ifal(kfl) ! reaction index capbl(k) = fpar(1,kfl)/capblref(k) capel(k) = fpar(3,kfl)/tempref end do endif c------------------------------------------------------------------------------- c ... define B_ref for reactions with specified reverse Arrhenius c ... parameters, and non-dimensionalize their B_k and E_k/R0 c ... kflagrev(k) is a flag, when 1->reverse parameters specified do k=1,numreac kflagrev(k) = 0 rexpo(k) = 0.0 end do do krev=1,nrev k = irev(krev) kflagrev(k) = 1 c1 = tempref**(par(2,k)-rpar(2,krev)) c2 = crefcgms**(nursum(k)-nupsum(k)) rcapbref(k) = capbref(k)*c1*c2 rcapb(k) = rpar(1,krev)/rcapbref(k) rexpo(k) = rpar(2,krev) rcape(k) = rpar(3,krev)/tempref end do c------------------------------------------------------------------------------- c ... for convenience/checking do k=1,numreac rr_f_ref(k) = capbref(k)*tempref**par(2,k) end do c------------------------------------------------------------------------------- do k=1,ns rg(k) = r0/wt(k) ! gas constant for each species (J/kg.K) end do c------------------------------------------------------------------------------- c------------------------------------------------------------------------------- c ... for radiation eps0 = 1.0d0 sig0 = 5.6696d-8 ! Stefan-Boltzmann constant (W/m^2.K^4) qradref = eps0 * sig0 * tempref**4 ! ref radiation heat flux (W/m^2) c------------------------------------------------------------------------------- do k=1,numreac ithb_type(k) = 0 ifal_type(k) = 0 end do if(nthb.gt.0)then do kthb=1,nthb k =ithb(kthb) ithb_type(k) = kthb enddo endif if(nfal.gt.0)then do kfl = 1, nfal k = ifal(kfl) ifal_type(k) = kfl enddo endif c------------------------------------------------------------------------------- c ... evaluate non-dimensional molar weights one = 1.0 do k=1,ns wtl(k) = wt(k) spwt(k) = wt(k)/wref ospwt(k) = one/spwt(k) end do c------------------------------------------------------------------------------- return end