! ! ********* C1-C3 MECHANISM ********* ! Ranzi, E., Frassoldati, A., Stagni, A., Pelucchi, M., Cuoci, A., Faravelli, T., Reduced kinetic schemes of complex reaction systems: Fossil and biomass-derived transportation fuels (2014) International Journal of Chemical Kinetics, 46 (9), pp. 512-542, DOI: 10.1002/kin.20867 ! Ranzi, E., Cavallotti, C., Cuoci, A., Frassoldati, A., Pelucchi, M., Faravelli, T., New reaction classes in the kinetic modeling of low temperature oxidation of n-alkanes (2015) Combustion and Flame, 162 (5), pp. 1679-1691, DOI: 10.1016/j.combustflame.2014.11.030 ! Bagheri, G., Ranzi, E., Pelucchi, M., Parente, A., Frassoldati, A., Faravelli, T., Comprehensive kinetic study of combustion technologies for low environmental impact: MILD and OXY-fuel combustion of methane (2020) Combustion and Flame 212, pp. 142-155, DOI: 10.1016/j.combustflame.2019.10.014 ! ! ********* NOx SUB-MECHANISM ********* ! Song, Y., Marrodán, L., Vin, N., Herbinet, O., Assaf, E., Fittschen, C., Stagni, A., Faravelli, T. & Battin-Leclerc, F., The sensitizing effects of NO2 and NO on methane low temperature oxidation in a jet stirred reactor. (2019) Proceedings of the Combustion Institute, 37(1), 667-675. DOI: 10.1016/j.proci.2018.06.115 ! ! ********* REDUCED USING THE ARCANE SOFTWARE ********* ! Cazères, Q., Pepiot, P., Riber, E., & Cuenot, B., A fully automatic procedure for the analytical reduction of chemical kinetics mechanisms for Computational Fluid Dynamics applications. (2021) Fuel, 303, 121247. DOI: 10.1016/j.fuel.2021.121247 ! ! ********* REDUCED AT ETH ZURICH BY QUENTIN MALE ********* ! Malé, Q., Pandey, K., & Noiray, N., The LEAF concept operated with hydrogen: Flame topology and NOx formation. (2024) Proceedings of the Combustion Institute. ! !-------------------------------------------------------------------------------------------------- ! Copyright (c) CERFACS (all rights reserved) !-------------------------------------------------------------------------------------------------- ! FILE ARC_S12R38QSS4.f90 !> @file ARC_S12R38QSS4.f90 !! Module for calculating the ARC source terms. !! @details !! @authors Q. Malé !! @date 2023/10/18 !! @since !! @note !-------------------------------------------------------------------------------------------------- ! !-------------------------------------------------------------------------------------------------- ! MODULE mod_customkinetics !> @details Generated by ARCANE custom kinetics routine to compute the chemical source terms. !! @authors Q. Cazères, P. Pepiot !! @date 2019/01/24 !-------------------------------------------------------------------------------------------------- module mod_customkinetics implicit none integer, parameter :: pr = selected_real_kind(15,307) ! Ideal gas constant real(pr), parameter :: Rcst = 8.3144621_pr ! Number of elements in the chemical system integer, parameter :: ne = 3 ! Number of non-qss and qss species and reactions integer, parameter :: nspec = 12 integer, parameter :: nqss = 4 integer, parameter :: nreac = 38 integer, parameter :: nreac_reverse = 38 ! Actual expression of each reaction character(len=65), dimension(nreac + nreac_reverse) :: reacexp ! Number of thirdbodies integer, parameter :: nTB = 3 integer, parameter :: nFO = 3 integer, parameter :: nTB_reverse = 3 integer, parameter :: nFO_reverse = 3 ! Index of elements integer, parameter :: eN = 1 integer, parameter :: eH = 2 integer, parameter :: eO = 3 ! Index of species integer, parameter :: sN2 = 1 integer, parameter :: sH2 = 2 integer, parameter :: sH = 3 integer, parameter :: sO2 = 4 integer, parameter :: sO = 5 integer, parameter :: sH2O = 6 integer, parameter :: sOH = 7 integer, parameter :: sHO2 = 8 integer, parameter :: sNO = 9 integer, parameter :: sNO2 = 10 integer, parameter :: sN = 11 integer, parameter :: sNH = 12 integer, parameter :: sqssN2O = 1 integer, parameter :: sqssHONO = 2 integer, parameter :: sqssNNH = 3 integer, parameter :: sqssNH2 = 4 ! Index of reactions integer, parameter :: r1f = 1 integer, parameter :: r2f = 2 integer, parameter :: r3f = 3 integer, parameter :: r4f = 4 integer, parameter :: r5f = 5 integer, parameter :: r6f = 6 integer, parameter :: r7f = 7 integer, parameter :: r8f = 8 integer, parameter :: r9f = 9 integer, parameter :: r10f = 10 integer, parameter :: r11f = 11 integer, parameter :: r12f = 12 integer, parameter :: r13f = 13 integer, parameter :: r14f = 14 integer, parameter :: r15f = 15 integer, parameter :: r16f = 16 integer, parameter :: r17f = 17 integer, parameter :: r18f = 18 integer, parameter :: r19f = 19 integer, parameter :: r20f = 20 integer, parameter :: r21f = 21 integer, parameter :: r22f = 22 integer, parameter :: r23f = 23 integer, parameter :: r24f = 24 integer, parameter :: r25f = 25 integer, parameter :: r26f = 26 integer, parameter :: r27f = 27 integer, parameter :: r28f = 28 integer, parameter :: r29f = 29 integer, parameter :: r30f = 30 integer, parameter :: r31f = 31 integer, parameter :: r32f = 32 integer, parameter :: r33f = 33 integer, parameter :: r34f = 34 integer, parameter :: r35f = 35 integer, parameter :: r36f = 36 integer, parameter :: r37f = 37 integer, parameter :: r38f = 38 integer, parameter :: r1b = 39 integer, parameter :: r2b = 40 integer, parameter :: r3b = 41 integer, parameter :: r4b = 42 integer, parameter :: r5b = 43 integer, parameter :: r6b = 44 integer, parameter :: r7b = 45 integer, parameter :: r8b = 46 integer, parameter :: r9b = 47 integer, parameter :: r10b = 48 integer, parameter :: r11b = 49 integer, parameter :: r12b = 50 integer, parameter :: r13b = 51 integer, parameter :: r14b = 52 integer, parameter :: r15b = 53 integer, parameter :: r16b = 54 integer, parameter :: r17b = 55 integer, parameter :: r18b = 56 integer, parameter :: r19b = 57 integer, parameter :: r20b = 58 integer, parameter :: r21b = 59 integer, parameter :: r22b = 60 integer, parameter :: r23b = 61 integer, parameter :: r24b = 62 integer, parameter :: r25b = 63 integer, parameter :: r26b = 64 integer, parameter :: r27b = 65 integer, parameter :: r28b = 66 integer, parameter :: r29b = 67 integer, parameter :: r30b = 68 integer, parameter :: r31b = 69 integer, parameter :: r32b = 70 integer, parameter :: r33b = 71 integer, parameter :: r34b = 72 integer, parameter :: r35b = 73 integer, parameter :: r36b = 74 integer, parameter :: r37b = 75 integer, parameter :: r38b = 76 ! Index of third body species integer, parameter :: mM1 = 1 integer, parameter :: mM5 = 2 integer, parameter :: mM13 = 3 integer, parameter :: mM12 = 4 integer, parameter :: mM29 = 5 integer, parameter :: mM36 = 6 ! Index of third body reactions integer, parameter :: TBr1f = 1 integer, parameter :: TBr5f = 2 integer, parameter :: TBr13f = 3 integer, parameter :: TBr1b = 4 integer, parameter :: TBr5b = 5 integer, parameter :: TBr13b = 6 ! Index of fall off reactions integer, parameter :: FOr12f = 1 integer, parameter :: FOr29f = 2 integer, parameter :: FOr36f = 3 integer, parameter :: FOr12b = 4 integer, parameter :: FOr29b = 5 integer, parameter :: FOr36b = 6 ! Molar mass real(pr), parameter, dimension(nspec) :: W_sp =(/ & 2.801348e-02_pr, & ! N2 2.015880e-03_pr, & ! H2 1.007940e-03_pr, & ! H 3.199880e-02_pr, & ! O2 1.599940e-02_pr, & ! O 1.801528e-02_pr, & ! H2O 1.700734e-02_pr, & ! OH 3.300674e-02_pr, & ! HO2 3.000614e-02_pr, & ! NO 4.600554e-02_pr, & ! NO2 1.400674e-02_pr, & ! N 1.501468e-02_pr & ! NH !4.401288e-02_pr, & ! N2O !4.701348e-02_pr, & ! HONO !2.902142e-02_pr, & ! NNH !1.602262e-02_pr & ! NH2 /) contains ! ----------------------------------------------- ! ! Subroutine for pressure dependent coefficients ! ! ----------------------------------------------- ! real(pr) function getlindratecoeff(Tloc,k0,kinf,fc,concin,Ploc) implicit none real(pr) :: Tloc,k0,kinf,fc,redP,Ploc real(pr) :: ntmp,ccoeff,dcoeff,lgknull real(pr) :: f real(pr) :: conc, concin if (concin.gt.0.0_pr) then conc = concin else conc = Ploc / ( Rcst * Tloc ) end if redP = abs(k0) * conc / max(abs(kinf), tiny(1.0_pr)) + tiny(1.0_pr) ntmp = 0.75_pr - 1.27_pr * log10( fc ) ccoeff = - 0.4_pr - 0.67_pr * log10( fc ) dcoeff = 0.14_pr lgknull = log10(redP) f = (lgknull+ccoeff)/(ntmp-dcoeff*(lgknull+ccoeff)) f = fc**(1.0_pr / ( f * f + 1.0_pr )) getlindratecoeff = kinf * f * redP / ( 1.0_pr + redP ) end function getlindratecoeff ! ----------------------------------------------- ! ! Evaluate thirdbodies ! ! ----------------------------------------------- ! subroutine get_thirdbodies(M,c) implicit none real(pr), dimension(nspec) :: c real(pr), dimension(nTB + nFO) :: M M(mM1) = (1.500000e+00_pr)*c(sH2) & + (1.100000e+01_pr)*c(sH2O) & + sum(c) M(mM5) = (-2.700000e-01_pr)*c(sH2) & + (2.650000e+00_pr)*c(sH2O) & + sum(c) M(mM13) = sum(c) M(mM12) = (3.000000e-01_pr)*c(sH2) & + (9.000000e+00_pr)*c(sH2O) & + sum(c) M(mM29) = (9.000000e+00_pr)*c(sH2O) & + (7.000000e-01_pr)*c(sN2) & + (5.000000e-01_pr)*c(sO2) & + sum(c) M(mM36) = (1.100000e+01_pr)*c(sH2O) & + (7.000000e-01_pr)*c(sN2) & + (4.000000e-01_pr)*c(sO2) & + sum(c) end subroutine get_thirdbodies ! ----------------------------------------------- ! ! Evaluate rate coefficients ! ! ----------------------------------------------- ! subroutine get_rate_coefficients(k,M,Tloc,Ploc) implicit none real(pr), dimension(nreac + nreac_reverse) :: k real(pr), dimension(nFO + nFO_reverse) :: k_0 real(pr), dimension(nFO + nFO_reverse) :: k_inf real(pr), dimension(nFO + nFO_reverse) :: FC real(pr), dimension(nTB + nFO) :: M real(pr) :: Tloc,Ploc,R_T_inv,T_log,redP ! Rate coefficients R_T_inv = 1.0_pr/(Rcst*Tloc) T_log = log(Tloc) k(r1f) = (4.577000e+13_pr)*exp((-4.368096e+05_pr)*R_T_inv + (-1.400000e+00_pr)*T_log ) k(r2f) = (5.080000e-02_pr)*exp((-2.632573e+04_pr)*R_T_inv + (2.670000e+00_pr)*T_log ) k(r3f) = (4.380000e+07_pr)*exp((-2.924616e+04_pr)*R_T_inv ) k(r4f) = (1.140000e+08_pr)*exp((-6.395662e+04_pr)*R_T_inv ) k(r5f) = (3.500000e+10_pr)*exp((-2.000000e+00_pr)*T_log ) k(r6f) = (6.700000e+01_pr)*exp((-6.270477e+04_pr)*R_T_inv + (1.704000e+00_pr)*T_log ) k(r7f) = (7.079000e+07_pr)*exp((-1.234280e+03_pr)*R_T_inv ) k(r8f) = (1.140200e+04_pr)*exp((-2.317016e+03_pr)*R_T_inv + (1.083000e+00_pr)*T_log ) k(r9f) = (3.250000e+07_pr) k(r10f) = (7.000000e+06_pr)*exp((4.572945e+03_pr)*R_T_inv ) k(r11f) = (4.500000e+08_pr)*exp((-4.572945e+04_pr)*R_T_inv ) k_0(FOr12f) =(1.740000e+07_pr)*exp((-1.230000e+00_pr)*T_log ) k_inf(FOr12f) =(4.650000e+06_pr)*exp((4.400000e-01_pr)*T_log ) FC(FOr12f) = ((1.0_pr - 6.700000e-01_pr)*exp(-Tloc/(1.000000e-30_pr)) + (6.700000e-01_pr)*exp(-Tloc/(1.000000e+30_pr)))& + exp(-(1.000000e+30_pr)/Tloc) k(r12f) = getlindratecoeff(Tloc, k_0(FOr12f), k_inf(FOr12f),FC(FOr12f), M(mM12),Ploc) k(r13f) = (1.000000e+04_pr) k(r14f) = (4.000000e+07_pr)*exp((-1.527160e+04_pr)*R_T_inv ) k(r15f) = (6.200000e+09_pr)*exp((-1.053322e+04_pr)*R_T_inv + (-1.150000e+00_pr)*T_log ) k(r16f) = (3.500000e+07_pr)*exp((-7.232881e+03_pr)*R_T_inv ) k(r17f) = (7.000000e+07_pr) k(r18f) = (1.570000e+01_pr)*exp((2.426720e+03_pr)*R_T_inv + (1.740000e+00_pr)*T_log ) k(r19f) = (1.000000e+07_pr) k(r20f) = (2.010000e+09_pr)*exp((-2.372328e+04_pr)*R_T_inv + (-1.380000e+00_pr)*T_log ) k(r21f) = (1.800000e+08_pr)*exp((1.020896e+03_pr)*R_T_inv + (-3.510000e-01_pr)*T_log ) k(r22f) = (2.830000e+07_pr) k(r23f) = (9.027000e+03_pr)*exp((-2.719600e+04_pr)*R_T_inv + (1.000000e+00_pr)*T_log ) k(r24f) = (4.280000e+07_pr)*exp((-6.568880e+03_pr)*R_T_inv ) k(r25f) = (1.000000e+09_pr) k(r26f) = (1.900000e+08_pr)*exp((9.204800e+01_pr)*R_T_inv + (-2.740000e-01_pr)*T_log ) k(r27f) = (5.200000e+05_pr)*exp((1.711256e+03_pr)*R_T_inv + (3.880000e-01_pr)*T_log ) k(r28f) = (2.110000e+06_pr)*exp((2.008320e+03_pr)*R_T_inv ) k_0(FOr29f) =(4.720000e+12_pr)*exp((-6.485200e+03_pr)*R_T_inv + (-2.870000e+00_pr)*T_log ) k_inf(FOr29f) =(1.300000e+09_pr)*exp((-7.500000e-01_pr)*T_log ) FC(FOr29f) = ((1.0_pr - 7.500000e-01_pr)*exp(-Tloc/(1.000000e+03_pr)) + (7.500000e-01_pr)*exp(-Tloc/(1.000000e+05_pr)))& + exp(-(1.000000e+30_pr)/Tloc) k(r29f) = getlindratecoeff(Tloc, k_0(FOr29f), k_inf(FOr29f),FC(FOr29f), M(mM29),Ploc) k(r30f) = (3.090000e+17_pr)*exp((-6.782264e+03_pr)*R_T_inv + (-4.170000e+00_pr)*T_log ) k(r31f) = (1.890000e-03_pr)*exp((-5.952577e+03_pr)*R_T_inv + (2.830000e+00_pr)*T_log ) k(r32f) = (4.300000e+03_pr)*exp((-1.702888e+04_pr)*R_T_inv + (9.800000e-01_pr)*T_log ) k(r33f) = (1.700000e+06_pr)*exp((2.175680e+03_pr)*R_T_inv ) k(r34f) = (2.000000e+05_pr)*exp((4.426672e+03_pr)*R_T_inv + (8.400000e-01_pr)*T_log ) k(r35f) = (3.920000e+06_pr)*exp((9.957920e+02_pr)*R_T_inv ) k_0(FOr36f) =(6.020000e+08_pr)*exp((-2.403457e+05_pr)*R_T_inv ) k_inf(FOr36f) =(9.900000e+10_pr)*exp((-2.422578e+05_pr)*R_T_inv ) redP = abs(k_0(FOr36f)) * M(mM36) / max(abs(k_inf(FOr36f)), tiny(1.0_pr)) + tiny(1.0_pr) k(r36f) = k_inf(FOr36f) * redP / ( 1.0_pr + redP ) k(r37f) = (5.000000e+08_pr)*exp((-7.573040e+04_pr)*R_T_inv ) k(r38f) = (6.620000e+07_pr)*exp((-1.114199e+05_pr)*R_T_inv ) k(r1b) = (2.598477e+08_pr)*exp((-4.054129e+03_pr)*R_T_inv + (-1.784129e+00_pr)*T_log ) k(r2b) = (2.502565e-02_pr)*exp((-2.013635e+04_pr)*R_T_inv + (2.658110e+00_pr)*T_log ) k(r3b) = (1.527374e+08_pr)*exp((-9.105995e+04_pr)*R_T_inv + (3.945820e-02_pr)*T_log ) k(r4b) = (4.965349e+05_pr)*exp((6.110298e+03_pr)*R_T_inv + (3.373898e-01_pr)*T_log ) k(r5b) = (2.149816e+16_pr)*exp((-4.945693e+05_pr)*R_T_inv + (-1.576413e+00_pr)*T_log ) k(r6b) = (9.465100e+00_pr)*exp((5.298392e+03_pr)*R_T_inv + (1.652652e+00_pr)*T_log ) k(r7b) = (4.315121e+04_pr)*exp((-1.538968e+05_pr)*R_T_inv + (6.251890e-01_pr)*T_log ) k(r8b) = (3.239183e+03_pr)*exp((-2.312358e+05_pr)*R_T_inv + (1.382689e+00_pr)*T_log ) k(r9b) = (4.548409e+06_pr)*exp((-2.227294e+05_pr)*R_T_inv + (2.877992e-01_pr)*T_log ) k(r10b) = (6.934637e+06_pr)*exp((-2.861597e+05_pr)*R_T_inv + (3.391470e-01_pr)*T_log ) k(r11b) = (4.457981e+08_pr)*exp((-3.364621e+05_pr)*R_T_inv + (3.391470e-01_pr)*T_log ) k_0(FOr12b) =(1.078839e+13_pr)*exp((-2.038366e+05_pr)*R_T_inv + (-1.145560e+00_pr)*T_log ) k_inf(FOr12b) =(2.883105e+12_pr)*exp((-2.038366e+05_pr)*R_T_inv + (5.244405e-01_pr)*T_log ) FC(FOr12b) = ((1.0_pr - 6.700000e-01_pr)*exp(-Tloc/(1.000000e-30_pr)) + (6.700000e-01_pr)*exp(-Tloc/(1.000000e+30_pr)))& + exp(-(1.000000e+30_pr)/Tloc) k(r12b) = getlindratecoeff(Tloc, k_0(FOr12b), k_inf(FOr12b),FC(FOr12b), M(mM12),Ploc) k(r13b) = (1.423517e+12_pr)*exp((-2.739036e+05_pr)*R_T_inv + (-2.529494e-01_pr)*T_log ) k(r14b) = (1.853949e+07_pr)*exp((-5.997735e+04_pr)*R_T_inv + (4.196312e-02_pr)*T_log ) k(r15b) = (1.214733e+10_pr)*exp((-1.373899e+03_pr)*R_T_inv + (-1.298055e+00_pr)*T_log ) k(r16b) = (8.669342e+07_pr)*exp((-1.110774e+05_pr)*R_T_inv + (7.427599e-02_pr)*T_log ) k(r17b) = (2.400132e+09_pr)*exp((-2.995616e+05_pr)*R_T_inv + (-2.112206e-01_pr)*T_log ) k(r18b) = (1.356091e+02_pr)*exp((-1.632316e+05_pr)*R_T_inv + (1.853734e+00_pr)*T_log ) k(r19b) = (6.960099e+08_pr)*exp((-3.057510e+05_pr)*R_T_inv + (-1.993309e-01_pr)*T_log ) k(r20b) = (3.001774e+08_pr)*exp((-2.532180e+05_pr)*R_T_inv + (-1.253831e+00_pr)*T_log ) k(r21b) = (3.704061e+14_pr)*exp((-1.530077e+05_pr)*R_T_inv + (-1.448680e+00_pr)*T_log ) k(r22b) = (7.952135e+08_pr)*exp((-2.019065e+05_pr)*R_T_inv + (-2.736069e-01_pr)*T_log ) k(r23b) = (1.104805e+03_pr)*exp((-1.590356e+05_pr)*R_T_inv + (1.063783e+00_pr)*T_log ) k(r24b) = (8.017840e+07_pr)*exp((-3.208938e+05_pr)*R_T_inv + (1.034297e-01_pr)*T_log ) k(r25b) = (3.070011e+03_pr)*exp((-3.308966e+04_pr)*R_T_inv + (-2.829555e-02_pr)*T_log ) k(r26b) = (4.556481e+13_pr)*exp((-2.016122e+05_pr)*R_T_inv + (-1.193552e+00_pr)*T_log ) k(r27b) = (6.060014e+04_pr)*exp((-4.596444e+04_pr)*R_T_inv + (5.661281e-01_pr)*T_log ) k(r28b) = (6.566553e+07_pr)*exp((-3.035686e+04_pr)*R_T_inv + (-2.325011e-01_pr)*T_log ) k_0(FOr29b) =(2.091027e+22_pr)*exp((-3.127540e+05_pr)*R_T_inv + (-3.355450e+00_pr)*T_log ) k_inf(FOr29b) =(5.759183e+18_pr)*exp((-3.062688e+05_pr)*R_T_inv + (-1.235450e+00_pr)*T_log ) FC(FOr29b) = ((1.0_pr - 7.500000e-01_pr)*exp(-Tloc/(1.000000e+03_pr)) + (7.500000e-01_pr)*exp(-Tloc/(1.000000e+05_pr)))& + exp(-(1.000000e+30_pr)/Tloc) k(r29b) = getlindratecoeff(Tloc, k_0(FOr29b), k_inf(FOr29b),FC(FOr29b), M(mM29),Ploc) k(r30b) = (4.602619e+29_pr)*exp((-2.162417e+05_pr)*R_T_inv + (-5.342724e+00_pr)*T_log ) k(r31b) = (1.141067e-05_pr)*exp((-1.089513e+05_pr)*R_T_inv + (3.529163e+00_pr)*T_log ) k(r32b) = (1.773189e-03_pr)*exp((-3.021387e+05_pr)*R_T_inv + (2.576311e+00_pr)*T_log ) k(r33b) = (3.579064e+04_pr)*exp((-1.626368e+05_pr)*R_T_inv + (7.386209e-01_pr)*T_log ) k(r34b) = (3.917383e+00_pr)*exp((-1.158707e+05_pr)*R_T_inv + (1.697690e+00_pr)*T_log ) k(r35b) = (1.762817e+04_pr)*exp((-1.893685e+05_pr)*R_T_inv + (5.203003e-01_pr)*T_log ) k_0(FOr36b) =(7.706559e-03_pr)*exp((-7.173106e+04_pr)*R_T_inv + (8.912561e-01_pr)*T_log ) k_inf(FOr36b) =(1.267358e+00_pr)*exp((-7.364315e+04_pr)*R_T_inv + (8.912561e-01_pr)*T_log ) redP = abs(k_0(FOr36b)) * M(mM36) / max(abs(k_inf(FOr36b)), tiny(1.0_pr)) + tiny(1.0_pr) k(r36b) = k_inf(FOr36b) * redP / ( 1.0_pr + redP ) k(r37b) = (5.554150e+02_pr)*exp((-3.336819e+05_pr)*R_T_inv + (1.263496e+00_pr)*T_log ) k(r38b) = (1.103035e+03_pr)*exp((-2.569530e+05_pr)*R_T_inv + (8.864592e-01_pr)*T_log ) return end subroutine get_rate_coefficients ! ----------------------------------------------- ! ! Evaluate reaction rates ! ! ----------------------------------------------- ! subroutine get_reaction_rates(w,k,m,c,cqss) implicit none real(pr), dimension(nspec) :: c real(pr), dimension(nqss) :: cqss real(pr), dimension(nreac + nreac_reverse) :: w,k real(pr), dimension(nTB + nFO) :: m w(r1f) = k(r1f) * c(sH2) * m(mM1) w(r2f) = k(r2f) * c(sH2) * c(sO) w(r3f) = k(r3f) * c(sH2) * c(sOH) w(r4f) = k(r4f) * c(sH) * c(sO2) w(r5f) = k(r5f) * c(sH) * c(sOH) * m(mM5) w(r6f) = k(r6f) * c(sH2O) * c(sO) w(r7f) = k(r7f) * c(sH) * c(sHO2) w(r8f) = k(r8f) * c(sH) * c(sHO2) w(r9f) = k(r9f) * c(sHO2) * c(sO) w(r10f) = k(r10f) * c(sHO2) * c(sOH) w(r11f) = k(r11f) * c(sHO2) * c(sOH) w(r12f) = k(r12f) * c(sH) * c(sO2) w(r13f) = k(r13f) * c(sO) * c(sOH) * m(mM13) w(r14f) = k(r14f) * c(sH) * cqss(sqssNH2) w(r15f) = k(r15f) * cqss(sqssNH2) * c(sNO) w(r16f) = k(r16f) * c(sH) * c(sNH) w(r17f) = k(r17f) * c(sNH) * c(sO) w(r18f) = k(r18f) * c(sNH) * c(sOH) w(r19f) = k(r19f) * c(sNH) * c(sOH) w(r20f) = k(r20f) * c(sNH) * c(sO2) w(r21f) = k(r21f) * c(sNH) * c(sNO) w(r22f) = k(r22f) * c(sN) * c(sOH) w(r23f) = k(r23f) * c(sN) * c(sO2) w(r24f) = k(r24f) * c(sN) * c(sNO) w(r25f) = k(r25f) * cqss(sqssNNH) w(r26f) = k(r26f) * cqss(sqssNNH) * c(sO) w(r27f) = k(r27f) * cqss(sqssNNH) * c(sO) w(r28f) = k(r28f) * c(sHO2) * c(sNO) w(r29f) = k(r29f) * c(sNO) * c(sO) w(r30f) = k(r30f) * c(sNO) * c(sOH) w(r31f) = k(r31f) * c(sH) * cqss(sqssHONO) w(r32f) = k(r32f) * c(sH) * cqss(sqssHONO) w(r33f) = k(r33f) * cqss(sqssHONO) * c(sOH) w(r34f) = k(r34f) * c(sH) * c(sNO2) w(r35f) = k(r35f) * c(sNO2) * c(sO) w(r36f) = k(r36f) * cqss(sqssN2O) w(r37f) = k(r37f) * c(sH) * cqss(sqssN2O) w(r38f) = k(r38f) * cqss(sqssN2O) * c(sO) w(r1b) = k(r1b) * c(sH)**2.0_pr * m(mM1) w(r2b) = k(r2b) * c(sH) * c(sOH) w(r3b) = k(r3b) * c(sH) * c(sH2O) w(r4b) = k(r4b) * c(sO) * c(sOH) w(r5b) = k(r5b) * c(sH2O) * m(mM5) w(r6b) = k(r6b) * c(sOH)**2.0_pr w(r7b) = k(r7b) * c(sOH)**2.0_pr w(r8b) = k(r8b) * c(sH2) * c(sO2) w(r9b) = k(r9b) * c(sO2) * c(sOH) w(r10b) = k(r10b) * c(sH2O) * c(sO2) w(r11b) = k(r11b) * c(sH2O) * c(sO2) w(r12b) = k(r12b) * c(sHO2) w(r13b) = k(r13b) * c(sHO2) * m(mM13) w(r14b) = k(r14b) * c(sH2) * c(sNH) w(r15b) = k(r15b) * cqss(sqssNNH) * c(sOH) w(r16b) = k(r16b) * c(sH2) * c(sN) w(r17b) = k(r17b) * c(sH) * c(sNO) w(r18b) = k(r18b) * c(sH2O) * c(sN) w(r19b) = k(r19b) * c(sH2) * c(sNO) w(r20b) = k(r20b) * c(sNO) * c(sOH) w(r21b) = k(r21b) * c(sH) * cqss(sqssN2O) w(r22b) = k(r22b) * c(sH) * c(sNO) w(r23b) = k(r23b) * c(sNO) * c(sO) w(r24b) = k(r24b) * c(sN2) * c(sO) w(r25b) = k(r25b) * c(sH) * c(sN2) w(r26b) = k(r26b) * c(sH) * cqss(sqssN2O) w(r27b) = k(r27b) * c(sNH) * c(sNO) w(r28b) = k(r28b) * c(sNO2) * c(sOH) w(r29b) = k(r29b) * c(sNO2) w(r30b) = k(r30b) * cqss(sqssHONO) w(r31b) = k(r31b) * c(sH2) * c(sNO2) w(r32b) = k(r32b) * c(sH2O) * c(sNO) w(r33b) = k(r33b) * c(sH2O) * c(sNO2) w(r34b) = k(r34b) * c(sNO) * c(sOH) w(r35b) = k(r35b) * c(sNO) * c(sO2) w(r36b) = k(r36b) * c(sN2) * c(sO) w(r37b) = k(r37b) * c(sN2) * c(sOH) w(r38b) = k(r38b) * c(sNO)**2.0_pr return end subroutine get_reaction_rates ! ----------------------------------------------- ! ! Evaluate production rates ! ! ----------------------------------------------- ! subroutine get_production_rates(cdot,w) implicit none real(pr), dimension(nspec) :: cdot real(pr), dimension(nreac + nreac_reverse) :: w cdot(sN2) = 0.0_pr & + w(r24f) & - w(r24b) & + w(r25f) & - w(r25b) & + w(r36f) & - w(r36b) & + w(r37f) & - w(r37b) cdot(sH2) = 0.0_pr & - w(r1f) & + w(r1b) & - w(r2f) & + w(r2b) & - w(r3f) & + w(r3b) & + w(r8f) & - w(r8b) & + w(r14f) & - w(r14b) & + w(r16f) & - w(r16b) & + w(r19f) & - w(r19b) & + w(r31f) & - w(r31b) cdot(sH) = 0.0_pr & + 2.0_pr * w(r1f) & - 2.0_pr * w(r1b) & + w(r2f) & - w(r2b) & + w(r3f) & - w(r3b) & - w(r4f) & + w(r4b) & - w(r5f) & + w(r5b) & - w(r7f) & + w(r7b) & - w(r8f) & + w(r8b) & - w(r12f) & + w(r12b) & - w(r14f) & + w(r14b) & - w(r16f) & + w(r16b) & + w(r17f) & - w(r17b) & + w(r21f) & - w(r21b) & + w(r22f) & - w(r22b) & + w(r25f) & - w(r25b) & + w(r26f) & - w(r26b) & - w(r31f) & + w(r31b) & - w(r32f) & + w(r32b) & - w(r34f) & + w(r34b) & - w(r37f) & + w(r37b) cdot(sO2) = 0.0_pr & - w(r4f) & + w(r4b) & + w(r8f) & - w(r8b) & + w(r9f) & - w(r9b) & + w(r10f) & - w(r10b) & + w(r11f) & - w(r11b) & - w(r12f) & + w(r12b) & - w(r20f) & + w(r20b) & - w(r23f) & + w(r23b) & + w(r35f) & - w(r35b) cdot(sO) = 0.0_pr & - w(r2f) & + w(r2b) & + w(r4f) & - w(r4b) & - w(r6f) & + w(r6b) & - w(r9f) & + w(r9b) & - w(r13f) & + w(r13b) & - w(r17f) & + w(r17b) & + w(r23f) & - w(r23b) & + w(r24f) & - w(r24b) & - w(r26f) & + w(r26b) & - w(r27f) & + w(r27b) & - w(r29f) & + w(r29b) & - w(r35f) & + w(r35b) & + w(r36f) & - w(r36b) & - w(r38f) & + w(r38b) cdot(sH2O) = 0.0_pr & + w(r3f) & - w(r3b) & + w(r5f) & - w(r5b) & - w(r6f) & + w(r6b) & + w(r10f) & - w(r10b) & + w(r11f) & - w(r11b) & + w(r18f) & - w(r18b) & + w(r32f) & - w(r32b) & + w(r33f) & - w(r33b) cdot(sOH) = 0.0_pr & + w(r2f) & - w(r2b) & - w(r3f) & + w(r3b) & + w(r4f) & - w(r4b) & - w(r5f) & + w(r5b) & + 2.0_pr * w(r6f) & - 2.0_pr * w(r6b) & + 2.0_pr * w(r7f) & - 2.0_pr * w(r7b) & + w(r9f) & - w(r9b) & - w(r10f) & + w(r10b) & - w(r11f) & + w(r11b) & - w(r13f) & + w(r13b) & + w(r15f) & - w(r15b) & - w(r18f) & + w(r18b) & - w(r19f) & + w(r19b) & + w(r20f) & - w(r20b) & - w(r22f) & + w(r22b) & + w(r28f) & - w(r28b) & - w(r30f) & + w(r30b) & - w(r33f) & + w(r33b) & + w(r34f) & - w(r34b) & + w(r37f) & - w(r37b) cdot(sHO2) = 0.0_pr & - w(r7f) & + w(r7b) & - w(r8f) & + w(r8b) & - w(r9f) & + w(r9b) & - w(r10f) & + w(r10b) & - w(r11f) & + w(r11b) & + w(r12f) & - w(r12b) & + w(r13f) & - w(r13b) & - w(r28f) & + w(r28b) cdot(sNO) = 0.0_pr & - w(r15f) & + w(r15b) & + w(r17f) & - w(r17b) & + w(r19f) & - w(r19b) & + w(r20f) & - w(r20b) & - w(r21f) & + w(r21b) & + w(r22f) & - w(r22b) & + w(r23f) & - w(r23b) & - w(r24f) & + w(r24b) & + w(r27f) & - w(r27b) & - w(r28f) & + w(r28b) & - w(r29f) & + w(r29b) & - w(r30f) & + w(r30b) & + w(r32f) & - w(r32b) & + w(r34f) & - w(r34b) & + w(r35f) & - w(r35b) & + 2.0_pr * w(r38f) & - 2.0_pr * w(r38b) cdot(sNO2) = 0.0_pr & + w(r28f) & - w(r28b) & + w(r29f) & - w(r29b) & + w(r31f) & - w(r31b) & + w(r33f) & - w(r33b) & - w(r34f) & + w(r34b) & - w(r35f) & + w(r35b) cdot(sN) = 0.0_pr & + w(r16f) & - w(r16b) & + w(r18f) & - w(r18b) & - w(r22f) & + w(r22b) & - w(r23f) & + w(r23b) & - w(r24f) & + w(r24b) cdot(sNH) = 0.0_pr & + w(r14f) & - w(r14b) & - w(r16f) & + w(r16b) & - w(r17f) & + w(r17b) & - w(r18f) & + w(r18b) & - w(r19f) & + w(r19b) & - w(r20f) & + w(r20b) & - w(r21f) & + w(r21b) & + w(r27f) & - w(r27b) return end subroutine get_production_rates ! ----------------------------------------------- ! ! Evaluation of QSS concentrations ! ! ----------------------------------------------- ! subroutine get_qss(cqss,c,k,M) implicit none real(pr), dimension(nqss) :: cqss real(pr), dimension(nspec) :: c real(pr), dimension(nreac + nreac_reverse) :: k real(pr), dimension(nTB + nFO) :: M integer :: index real(pr) :: N2O_ct real(pr) :: N2O_num real(pr) :: N2O_denom real(pr) :: N2O_NNH real(pr) :: HONO_ct real(pr) :: HONO_num real(pr) :: HONO_denom real(pr) :: NNH_ct real(pr) :: NNH_num real(pr) :: NNH_denom real(pr) :: NNH_N2O real(pr) :: NNH_NH2 real(pr) :: NH2_ct real(pr) :: NH2_num real(pr) :: NH2_denom real(pr) :: NH2_NNH real(pr) :: a_d_a real(pr) :: a_g_a real(pr) :: a_h_b real(pr) :: A_A_A real(pr) :: A_B_A real(pr) :: A_B_B NH2_denom = tiny(1.0_pr) + ( 0.0_pr & + k(r14f)* c(sH) & + k(r15f)* c(sNO) ) NH2_num = ( 0.0_pr & + k(r14b)* c(sH2) * c(sNH) ) NH2_ct = NH2_num / NH2_denom NH2_NNH = - ( 0.0_pr & + k(r15b) * c(sOH) ) / NH2_denom NNH_denom = tiny(1.0_pr) + ( 0.0_pr & + k(r25f)& + k(r26f)* c(sO) & + k(r27f)* c(sO) & + k(r15b)* c(sOH) ) NNH_num = ( 0.0_pr & + k(r25b)* c(sH) * c(sN2) & + k(r27b)* c(sNH) * c(sNO) ) NNH_ct = NNH_num / NNH_denom NNH_N2O = - ( 0.0_pr & + k(r26b) * c(sH) ) / NNH_denom NNH_NH2 = - ( 0.0_pr & + k(r15f) * c(sNO) ) / NNH_denom N2O_denom = tiny(1.0_pr) + ( 0.0_pr & + k(r36f)& + k(r37f)* c(sH) & + k(r38f)* c(sO) & + k(r21b)* c(sH) & + k(r26b)* c(sH) ) N2O_num = ( 0.0_pr & + k(r21f)* c(sNH) * c(sNO) & + k(r36b)* c(sN2) * c(sO) & + k(r37b)* c(sN2) * c(sOH) & + k(r38b) *c(sNO)** 2.0_pr ) N2O_ct = N2O_num / N2O_denom N2O_NNH = - ( 0.0_pr & + k(r26f) * c(sO) ) / N2O_denom a_d_a = (1.0_pr) - (N2O_NNH) * (NNH_N2O) a_g_a = NH2_NNH a_h_b = (1.0_pr) - (NNH_NH2) * (a_g_a) / (a_d_a) A_A_A = (NNH_ct) - (N2O_ct) * (NNH_N2O) A_B_A = NH2_ct A_B_B = (A_B_A) - (A_A_A) * (a_g_a) / (a_d_a) cqss(sqssNH2) = ( A_B_B ) / ( a_h_b ) cqss(sqssNNH) = (A_A_A - (NNH_NH2) * cqss(sqssNH2)) / (a_d_a) cqss(sqssN2O) = N2O_ct - (N2O_NNH) * cqss(sqssNNH) HONO_denom = tiny(1.0_pr) + ( 0.0_pr & + k(r31f)* c(sH) & + k(r32f)* c(sH) & + k(r33f)* c(sOH) & + k(r30b) ) HONO_num = ( 0.0_pr & + k(r30f)* c(sNO) * c(sOH) & + k(r31b)* c(sH2) * c(sNO2) & + k(r32b)* c(sH2O) * c(sNO) & + k(r33b)* c(sH2O) * c(sNO2) ) HONO_ct = HONO_num / HONO_denom cqss(sqssHONO) = HONO_ct cqss = max(cqss, tiny(1.0_pr)) cqss = min(cqss, 1e03_pr) ! Necessary for very low concentration variations ! Needs to be integrated only for Cantera do index=1,nqss if (cqss(index) /= cqss(index)) then cqss(index) = tiny(1.0_pr) end if end do return end subroutine get_qss ! ----------------------------------------------- ! ! Mass fractions to concentrations ! ! ----------------------------------------------- ! subroutine y2c(y, W_sp, P, T, c) implicit none real(pr),dimension(nspec) :: W_sp real(pr),dimension(nspec) :: c, y real(pr) :: rho, P, T, inv_W_g integer :: k ! Gas molecular weight inverse inv_W_g = 0.0_pr do k =1, nspec inv_W_g = inv_W_g + y(k) / W_sp(k) end do ! Gas density rho = P / (Rcst * inv_W_g * T) ! Conversion c = y * rho / W_sp return end subroutine y2c end module mod_customkinetics ! ---------------------------------------------------- ! ! routine to get the species source terms from the ARC ! ! ! ! INPUTS ! ! P: Pressure ! ! T: Temperature ! ! y: species mass fractions ! ! OUTPUT ! ! wdot: species source terms ! ! ! ! ---------------------------------------------------- ! subroutine customkinetics(P, T, y, wdot) use mod_customkinetics implicit none real(pr), dimension(nspec) :: y, c, wdot, cdot real(pr), dimension(nqss) :: cqss real(pr), dimension(nreac + nreac_reverse) :: w,k real(pr), dimension(nTB + nFO) :: M real(pr) :: P, T, rho ! Convert to concentrations call y2c(y, W_sp, P, T, c) ! Evaluate QSS concentrations and reaction rates call get_thirdbodies(M,c) call get_rate_coefficients(k, M, T, P) call get_qss(cqss, c, k, M) call get_reaction_rates(w, k, M, c, cqss) ! Evaluate production rates call get_production_rates(wdot, w) return end subroutine customkinetics