From c62542c24a6808136eb4a66d3ecada2ed1072ddc Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:03:05 -0400 Subject: [PATCH 01/16] Add binary diffusion --- src/Calculators/Diffusion.jl | 39 ++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/Calculators/Diffusion.jl b/src/Calculators/Diffusion.jl index 548fe87a..00d9419b 100644 --- a/src/Calculators/Diffusion.jl +++ b/src/Calculators/Diffusion.jl @@ -18,3 +18,42 @@ export ConstantDiffusivity end (sd::StokesDiffusivity)(;T::N,mu::Q,P::R=0.0) where {N,R,Q<:Number} = @fastmath kB*T/(6*Base.pi*mu*sd.r) export StokesDiffusivity + +""" +ChapmanEnskogBinaryDiffusivity: Calculate binary diffusion coefficient for species i and species j using Chapman-Enskog theory +and Lorentz-Berthelot combining rules +See Eq. 11.4 in Kee et al. 2018 +""" + +function getChapmanEnskogBinaryDiffusivity(transporti::TransportModel, transportj::TransportModel; T::N1, P::N2) where {N1<:Number,N2<:Number} + mi = transporti.m + mj = transportj.m + + mij = (mi * mj) / (mi + mj) + + epsiloni = transporti.potential.epsilon + epsilonj = transportj.potential.epsilon + sigmai = transporti.potential.sigma + sigmaj = transportj.potential.sigma + + reduced_polari = transporti.polarizability / sigmai^3 # Eq. 11.41 in Kee et al. 2018 + reduced_dipolej = transportj.dipolemoment / sqrt(epsilonj * sigmaj^3) # Eq. 11.42 in Kee et al. 2018 + xi = 1 + 0.25 * reduced_polari * reduced_dipolej * sqrt(epsiloni / epsilonj) # induction energy term Eq. 11.40 in Kee et al. 2018 + + epsilonij = xi^2 * sqrt(epsiloni * epsilonj) # Eq. 11.38 in Kee et al. 2018 + sigmaij = 0.5 * (sigmai + sigmaj) * xi^(-1 / 6) # Eq. 11.39 in Kee et al. 2018 + + omegaij = getcollisionintegral11(transporti.potential, transportj.potential, epsilonij; T=T) + + return 3 / 16 * sqrt(2 * pi * kB^3 * T^3 / mij) / (P * pi * sigmaij^2 * omegaij) +end + +function getChapmanEnskogSelfDiffusivity(transport::TransportModel; T::N1, P::N2) where {N1<:Number,N2<:Number} + m = transport.m + epsilon = transport.potential.epsilon + sigma = transport.potential.sigma + + omegaii = getcollisionintegral11(transport.potential, transport.potential, epsilon; T=T) + + return 3 / 8 * sqrt(pi * kB^3 * T^3 / m) / (P * pi * sigma^2 * omegaii) +end From 342ffc5a97177470984ccda80b69a003a1dfdb81 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:03:20 -0400 Subject: [PATCH 02/16] Add LennardJonesPotential --- src/Calculators/Potential.jl | 42 ++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) create mode 100644 src/Calculators/Potential.jl diff --git a/src/Calculators/Potential.jl b/src/Calculators/Potential.jl new file mode 100644 index 00000000..ef4459dc --- /dev/null +++ b/src/Calculators/Potential.jl @@ -0,0 +1,42 @@ +abstract type AbstractPotential end +export AbstractPotential + +struct EmptyPotential <: AbstractPotential end +export EmptyPotential + +""" +Lennard-Jones potential parameters +Used later in binary diffusion coefficient calculations +See Eq. 11.6 in Kee et al. 2018 +""" +@with_kw struct LennardJonesPotential <: AbstractPotential + epsilon::Float64 #L-J depth + sigma::Float64 #L-J diameter +end +export LennardJonesPotential + +function getcollisionintegral11(potentiali::LennardJonesPotential, potentialj::LennardJonesPotential, epsilonij::Float64; T::Number,) + + a1 = 1.0548 + a2 = 0.15504 + a3 = 0.55909 + a4 = 2.1705 + + reduced_T = T * kB / epsilonij + omegaij = a1 * reduced_T^(-a2) + (reduced_T + a3)^(-a4) # Eq. 11.6 in Kee et al. 2018 + + return omegaij +end + +function getcollisionintegral22(potentiali::LennardJonesPotential, potentialj::LennardJonesPotential, epsilonij::Float64; T::Number, ) + + b1 = 1.0413 + b2 = 0.11930 + b3 = 0.43628 + b4 = 1.6041 + + reduced_T = T * kB / epsilonij + omegaij = b1 * reduced_T^(-b2) + (reduced_T + b3)^(-b4) # Eq. 11.6 in Kee et al. 2018 + + return omegaij +end \ No newline at end of file From 70334f42c8ad8d34be252ca67573cc76036e8f16 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:03:32 -0400 Subject: [PATCH 03/16] Add ThermalConductivity --- src/Calculators/ThermalConductivity.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) create mode 100644 src/Calculators/ThermalConductivity.jl diff --git a/src/Calculators/ThermalConductivity.jl b/src/Calculators/ThermalConductivity.jl new file mode 100644 index 00000000..d9c19cd6 --- /dev/null +++ b/src/Calculators/ThermalConductivity.jl @@ -0,0 +1,13 @@ + +function getChapmanEnskogThermalConductivity(transport::TransportModel, thermo::NASA; T::N1) where {N1<:Number} + m = transport.m + potential = transport.potential + epsilon = potential.epsilon + sigma = potential.sigma + cp = getHeatCapacity(thermo, T) + cv = cp - R + + omega = getcollisionintegral22(potential, potential, epsilon; T=T) + return 25 / (32 * pi^0.5) * (kB * T / m)^0.5 * cv / (sigma^2 * Na * omega) +end +export ChapmanEnskogThermalConductivity From aa603ba6b9d11b2fef23a6bcd2f6f41da1bf2f41 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:03:41 -0400 Subject: [PATCH 04/16] Add transport model --- src/Calculators/Transport.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 src/Calculators/Transport.jl diff --git a/src/Calculators/Transport.jl b/src/Calculators/Transport.jl new file mode 100644 index 00000000..e130ceca --- /dev/null +++ b/src/Calculators/Transport.jl @@ -0,0 +1,12 @@ +abstract type AbstractTransportModel end +export AbstractTransportModel + +struct EmptyTransportModel <: AbstractTransportModel end + +@with_kw struct TransportModel{N1<:AbstractPotential} <: AbstractTransportModel + m::Float64 # molecular mass + potential::N1 + dipolemoment::Float64 + polarizability::Float64 +end +export TransportModel From bcdceb443c2bbe93db17c89e16b7f12d3ab39520 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:03:52 -0400 Subject: [PATCH 05/16] Add mixture averaged diffusion --- src/Flame1D/Diffusion1D.jl | 41 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 src/Flame1D/Diffusion1D.jl diff --git a/src/Flame1D/Diffusion1D.jl b/src/Flame1D/Diffusion1D.jl new file mode 100644 index 00000000..8f0ed4d0 --- /dev/null +++ b/src/Flame1D/Diffusion1D.jl @@ -0,0 +1,41 @@ +abstract type AbstractDiffusiveFlux end + +struct MixtureAveragedDiffusiveFlux <: AbstractDiffusiveFlux end + +function mixtureaverageddiffusivities!(mixdiffs, phase::IdealGas, C, cs, T, P) + mixdiffs .= 0.0 + + sum_cs = sum(cs) + for spcind in 1:length(phase.species) + if sum_cs == cs[spcind] # pure species + mixdiffs[spcind] = getChapmanEnskogSelfDiffusivity(phase.species[spcind].transport; T=T, P=P) + return mixdiffs + end + end + + for spcind in 1:length(phase.species) + cdivD = 0.0 + for otherspcind in 1:length(phase.species) + if spcind != otherspcind + bindarydiff = getChapmanEnskogBinaryDiffusivity(phase.species[spcind].transport, phase.species[otherspcind].transport; T=T, P=P) + cdivD += cs[otherspcind] / bindarydiff + end + end + mixdiffs[spcind] = (C - cs[spcind]) / cdivD + end + return mixdiffs +end + +function mixtureaveragedthermalconductivity(phase::IdealGas, cs::N1, T::N2, P::N3) where {N1<:AbstractArray, N2<:Number, N3<:Number} + mixtureaveragedlambdamult = 0.0 + mixtureaveragedlambdadiv = 0.0 + C = P / (R * T) + for spcind in 1:length(phase.species) + spc = phase.species[spcind] + lambda = getChapmanEnskogThermalConductivity(spc.transport, spc.thermo; T=T) + molfrac = cs[spcind] / C + mixtureaveragedlambdamult += molfrac * lambda + mixtureaveragedlambdadiv += molfrac / lambda + end + return 0.5 * (mixtureaveragedlambdamult + 1.0 / mixtureaveragedlambdadiv) +end \ No newline at end of file From cf8c0bc48f3224fbbab59753aba399632c8e9440 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:02 -0400 Subject: [PATCH 06/16] Add flame objects --- src/Flame1D/Flame1D.jl | 197 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 src/Flame1D/Flame1D.jl diff --git a/src/Flame1D/Flame1D.jl b/src/Flame1D/Flame1D.jl new file mode 100644 index 00000000..c18dcbf4 --- /dev/null +++ b/src/Flame1D/Flame1D.jl @@ -0,0 +1,197 @@ +abstract type AbstractFlame end +abstract type AxisymmetricFlame <: AbstractFlame end +abstract type FlatFlame <: AxisymmetricFlame end + +function getrho(cs, molecularweights::Array{Float64,1}) + return dot(cs, molecularweights) +end + +""" +`Burner Flame`: struct for 1D pre-mixed burner stabilized flame +- `phase`: IdealGas phase +- `inlet`: Inlet1D struct contains the burner condition + - `T`: burner temperature (K) + - `cs`: Array{Float64,1} contains the inlet molar concentration of species (mol/m^3) +- `indexes`: Array{Int64,1} contains [start index of species, end index of species, index of temperature, index of axial velocity] +- `P`: pressure (Pa) +- `mdot`: net mass flux (kg/m^2/s) +- `u0`: initial axial velocity (m/s) +- `rxnarray`: Array{Int64, 2} contains the reaction array of shape (number of reactions, 8) with indexes of participating species +- `p`: Array{Float64,1} contains the parameters +- `variable_index_dict`: Dict{String,Int64} contains the index of the variable in the solution vector +""" + +struct BurnerFlame{N1<:AbstractDiffusiveFlux} <: FlatFlame + phase::IdealGas + inlet::Inlet1D + molecularweights::Array{Float64,1} + indexes::Array{Int64,1} + P::Float64 + mdot::Float64 + u0::Float64 + rxnarray::Array{Int64, 2} + p::Array{Float64,1} + variable_index_dict::Dict{String,Int64} + diffusive_flux_model::N1 +end + +""" +`Burner Flame`: constructor for 1D pre-mixed burner stabilized flame +- `phase`: IdealGas phase +- `initialcond`: Dict{String,Float64} contains the initial conditions for the flame. + It takes in the following keys and values: + - species name: premixed gas molar concentration (mol/m^3) + - "P": pressure (Pa) + - "T": the burner temperature (K) + - "mdot": the net mass flux (kg/m^2/s) +""" + +function BurnerFlame(;phase::IdealGas, initialcond::Dict, diffusive_flux_model::String="mix") + @assert isa(phase, IdealGas) "phase must be an IdealGas phase" + @assert "T" in keys(initialcond) "Initial condition must contain the burner temperature, T (K)" + @assert "P" in keys(initialcond) "Initial condition must contain the pressure, P (Pa)" + @assert "mdot" in keys(initialcond) "Initial condition must contain the net mass flux, mdot (kg/m^2/s)" + @assert diffusive_flux_model in ["mix"] "diffusive flux model must be one of the following: mix" # currently only supporting mixture averaged diffusive flux + + # set inlet conditions + T = 0.0 + P = 0.0 + mdot = 0.0 + + cs = zeros(length(phase.species)) # inlet species molar concentration (mol/m^3) + spcnames = getfield.(phase.species, :name) + molecularweights = getfield.(phase.species, :molecularweight) + + for (key, value) in initialcond + if key == "T" + T = value + elseif key == "P" + P = value + elseif key == "mdot" + mdot = value + else + ind = findfirst(isequal(key), spcnames) + @assert !isnothing(ind) "Species $key not found in the mechanism" + cs[ind] = value + end + end + + rho = getrho(cs, molecularweights) # density (kg/m^3) + u = mdot / rho # axial velocity (m/s) + + inlet = Inlet1D(T, cs) + + variable_index_dict = phase.spcdict + variable_index_dict["T"] = length(spcnames) + 1 + variable_index_dict["u"] = length(spcnames) + 2 + + p = vcat(zeros(length(phase.species)), ones(length(phase.reactions))) + + if diffusive_flux_model == "mix" + diffusive_flux_model = MixtureAveragedDiffusiveFlux() + end + return BurnerFlame(phase, inlet, molecularweights, [1, length(phase.species), length(phase.species) + 1, length(phase.species) + 2], P, mdot, u, phase.rxnarray, p, variable_index_dict, diffusive_flux_model), p +end + +export BurnerFlame + +""" +`FreeFlame`: struct for 1D pre-mixed freely propagating adiabatic flame +- `phase`: IdealGas phase +- `inlet`: Inlet1D struct contains the inlet condition + - `T`: inlet temperature (K) + - `cs`: Array{Float64,1} contains the inlet molar concentration of species (mol/m^3) +- `indexes`: Array{Int64,1} contains [start index of species, end index of species, index of temperature, index of axial velocity] +- `P`: pressure (Pa) +- `rxnarray`: Array{Int64, 2} contains the reaction array of shape (number of reactions, 8) with indexes of participating species +- `p`: Array{Float64,1} contains the parameters +- `variable_index_dict`: Dict{String,Int64} contains the index of the variable in the solution vector +""" + +struct FreeFlame{N1<:AbstractDiffusiveFlux} <: FlatFlame + phase::IdealGas + inlet::Inlet1D + molecularweights::Array{Float64,1} + indexes::Array{Int64,1} + P::Float64 + rxnarray::Array{Int64, 2} + p::Array{Float64,1} + variable_index_dict::Dict{String,Int64} + diffusive_flux_model::N1 +end + +""" +`FreeFlame`: constructor for 1D pre-mixed freely propagating adiabatic flame +- `phase`: IdealGas phase +- `initialcond`: Dict{String,Float64} contains the initial conditions for the flame. + It takes in the following keys and values: + - species name: premixed gas molar concentration (mol/m^3) + - "P": pressure (Pa) + - "T": the inlet temperature (K) +""" + +function FreeFlame(;phase::IdealGas, initialcond::Dict, diffusive_flux_model::String="mix") + @assert isa(phase, IdealGas) "phase must be an IdealGas phase" + @assert "T" in keys(initialcond) "Initial condition must contain the burner temperature, T (K)" + @assert "P" in keys(initialcond) "Initial condition must contain the pressure, P (Pa)" + @assert diffusive_flux_model in ["mix"] "diffusive flux model must be one of the following: mix" # currently only supporting mixture averaged diffusive flux + + # set inlet conditions + T = 0.0 + P = 0.0 + + cs = zeros(length(phase.species)) # inlet species molar concentration (mol/m^3) + spcnames = getfield.(phase.species, :name) + molecularweights = getfield.(phase.species, :molecularweight) + + for (key, value) in initialcond + if key == "T" + T = value + elseif key == "P" + P = value + else + ind = findfirst(isequal(key), spcnames) + @assert !isnothing(ind) "Species $key not found in the mechanism" + cs[ind] = value + end + end + + variable_index_dict = phase.spcdict + variable_index_dict["T"] = length(spcnames) + 1 + variable_index_dict["u"] = length(spcnames) + 2 + + p = vcat(zeros(length(phase.species)), ones(length(phase.reactions))) + + inlet = Inlet1D(T, cs) + + if diffusive_flux_model == "mix" + diffusive_flux_model = MixtureAveragedDiffusiveFlux() + end + + return FreeFlame(phase, inlet, molecularweights, [1, length(phase.species), length(phase.species) + 1, length(phase.species) + 2], P, phase.rxnarray, p, variable_index_dict, diffusive_flux_model), p +end + +export FreeFlame + +function calcthermo(flame::N1, y, t, p, cell_size) where {N1<:AxisymmetricFlame} + @views cs = y[flame.indexes[1]:flame.indexes[2]] + T = y[flame.variable_index_dict["T"]] + P = flame.P + V = cell_size + C = P / (R * T) + N = P * V / (R * T) + ns = cs .* V + + cpdivR, hdivRT, sdivR = calcHSCpdless(flame.phase.vecthermo, T) + Gs = (hdivRT .- sdivR) .* (R * T) + Hs = hdivRT .* (R * T) + Cvave = dot(cpdivR, ns) * R / N - R + + kfs, krevs = getkfkrevs(flame.phase, T, P, C, N, ns, Gs, Array{Float64,1}(), V, 0.0) + return ns, cs, T, P, V, C, N, 0.0, kfs, krevs, Hs, Array{Float64,1}(), Gs, Array{Float64,1}(), Cvave, cpdivR, 0.0 +end + +function calcdomainderivatives!(flame::N1, rhs; t, T, P, Us, Hs, V, C, ns, N, Cvave) where {N1<:AxisymmetricFlame} + Cpave = Cvave + R + @views rhs[flame.variable_index_dict["T"]] = - dot(Hs, rhs[flame.indexes[1]:flame.indexes[2]])/(N * Cpave) +end \ No newline at end of file From 5f23b4c65bb9847957dd22ed9e3e87eaa101646c Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:10 -0400 Subject: [PATCH 07/16] Add 1D inlet interface --- src/Flame1D/Interface1D.jl | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 src/Flame1D/Interface1D.jl diff --git a/src/Flame1D/Interface1D.jl b/src/Flame1D/Interface1D.jl new file mode 100644 index 00000000..d559b466 --- /dev/null +++ b/src/Flame1D/Interface1D.jl @@ -0,0 +1,8 @@ +""" +`Inlet1D` is a struct that contains the inlet boundary condition for a 1D flame. + - `T`: inlet temperature (K) +""" +struct Inlet1D + T::Float64 + cs::Array{Float64,1} +end \ No newline at end of file From 26cab493b9754edef46cc2f5ef9805b03cef9455 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:17 -0400 Subject: [PATCH 08/16] Add Flame simulation --- src/Flame1D/Reactor1D.jl | 214 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 src/Flame1D/Reactor1D.jl diff --git a/src/Flame1D/Reactor1D.jl b/src/Flame1D/Reactor1D.jl new file mode 100644 index 00000000..43990277 --- /dev/null +++ b/src/Flame1D/Reactor1D.jl @@ -0,0 +1,214 @@ +using PreallocationTools +using NonlinearSolve +using SciMLNLSolve + +""" +`FlameReactor`: function to construct the ODE problem +- `flame`: Type of flame to solve for. Determines the boundary conditions. +- `p`: Array{Float64,1} contains the parameters +- `L`: Axial length (m) +""" + +struct FlameSimulationCache{X1,X2,X3,X4} + mixdiffusivities::X1 + diffusivefluxes_jphalf::X2 + diffusivefluxes_jmhalf::X3 + diffusivefluxes_j::X4 +end + +struct FlameReactor{F <: AbstractFlame, R} + flame::F + p::Array{Float64,1} + L::Float64 + chunk_size::Int + cell_centers::Array{Float64,1} + cell_sizes::Array{Float64,1} + nonlinearprob::NonlinearProblem + recommendedsolver::R +end + +function FlameReactor(flame::F, p::Array{Float64,1}, L::Float64; chunk_size::Int=12) where {F<:AbstractFlame,} + # function f!(residual, y, p) + # f_flame!(unflatten(residual), unflatten(y), p, 0.0, flame, cell_centers, cell_sizes, cache) + # end + # function jacobianyforwarddiff!(J, y, p) + # ForwardDiff.jacobian!(J,(residual, y) -> f!(residual, y, p), residualcache, y) + # end + function f!(residual, y, p, t) + f_flame!(unflatten(residual), unflatten(y), p, 0.0, flame, cell_centers, cell_sizes, cache) + end + function jacobianyforwarddiff!(J, y, p, t) + ForwardDiff.jacobian!(J,(residual, y) -> f!(residual, y, p, t), residualcache, y) + end + function unflatten(y::AbstractVector) + return reshape(y, num_variables, :) + end + # Create one dimensional grid + num_cells = 5 + num_variables = length(flame.phase.species) + 2 + zs = range(0.0, stop=L, length=num_cells+1) + cell_sizes = diff(zs) + cell_centers = zs[1:end-1] .+ cell_sizes ./ 2 # Cell centers + + mixdiffusivities = dualcache(zeros(length(flame.phase.species)), chunk_size) + diffusivefluxes_jphalf = dualcache(zeros(length(flame.phase.species)), chunk_size) + diffusivefluxes_jmhalf = dualcache(zeros(length(flame.phase.species)), chunk_size) + diffusivefluxes_j = dualcache(zeros(length(flame.phase.species)), chunk_size) + cache = FlameSimulationCache(mixdiffusivities, diffusivefluxes_jphalf, diffusivefluxes_jmhalf, diffusivefluxes_j) + residualcache = zeros(num_variables * num_cells) + + if isa(flame, BurnerFlame) + y0 = zeros(num_variables, num_cells) + @views y0[flame.indexes[1]:flame.indexes[2], :] .= flame.inlet.cs + @views y0[flame.variable_index_dict["T"], :] .= flame.inlet.T + @views y0[flame.variable_index_dict["u"], :] .= flame.u0 + y0 = reshape(y0, :) + end + + # residuals = zeros(num_variables * num_cells) + # f!(residuals, y0, p) + # println("y0") + # display(unflatten(y0)) + # println("residual") + # display(unflatten(residuals)) + + # nonlinearfcn = NonlinearFunction(f!; jac = jacobianyforwarddiff!) + # nonlinearprob = NonlinearProblem(nonlinearfcn, y0, p) + # recommendedsolver = NewtonRaphson() + + # sol = solve(nonlinearprob, recommendedsolver, abstol=1e-18, reltol=1e-6) + + # f!(residuals, sol.u, p) + # println("sol") + # display(unflatten(sol.u)) + # println("residual") + # display(unflatten(residuals)) + + residuals = zeros(num_variables * num_cells) + f!(residuals, y0, p, 0.0) + println("y0") + display(unflatten(y0)) + println("residual") + display(unflatten(residuals)) + + nonlinearfcn = ODEFunction(f!) + nonlinearprob = NonlinearProblem(nonlinearfcn, y0, p) + recommendedsolver = NLSolveJL() + + sol = solve(nonlinearprob, recommendedsolver, abstol=1e-18, reltol=1e-6) + + f!(residuals, sol.u, p, 0.0) + println("sol") + display(unflatten(sol.u)) + println("residual") + display(unflatten(residuals)) + + return FlameReactor(flame, p, L, chunk_size, cell_centers, cell_sizes, nonlinearprob, recommendedsolver) +end + +export FlameReactor + +function reaction!(residual, flame::F; t, ns, cs, T, P, V, C, N, kfs, krevs, Hs, Us, Cvave) where F <: AxisymmetricFlame + addreactionratecontributions!(residual, flame.rxnarray, cs, kfs, krevs) + @views residual[flame.indexes[1]:flame.indexes[2]] .*= V + calcdomainderivatives!(flame, residual; t=t, T=T, P=P, Us=Us, Hs=Hs, V=V, C=C, ns=ns, N=N, Cvave=Cvave) +end + +function convection!(residual, y_j, y_jm1, z_j, z_jm1, flame::F) where F <: AxisymmetricFlame + u_j = y_j[flame.variable_index_dict["u"]] + @views residual[flame.indexes[1]:flame.indexes[2]] .+= - u_j * (y_j[flame.indexes[1]:flame.indexes[2]] .- y_jm1[flame.indexes[1]:flame.indexes[2]]) ./ (z_j - z_jm1) + residual[flame.variable_index_dict["T"]] += - u_j * (y_j[flame.variable_index_dict["T"]] - y_jm1[flame.variable_index_dict["T"]]) / (z_j - z_jm1) +end + +function diffusion!(residual, y_jm1, y_jp1, z_j, z_jm1, z_jp1, dz_j, dz_jm1, dz_jp1, flame::F, cache::FlameSimulationCache; cs, T, P, N, Cvave, cpdivR) where F <: AxisymmetricFlame + @views cs_j = cs + @views cs_jm1 = y_jm1[flame.indexes[1]:flame.indexes[2]] + @views cs_jp1 = y_jp1[flame.indexes[1]:flame.indexes[2]] + T_j = T + T_jm1 = y_jm1[flame.variable_index_dict["T"]] + T_jp1 = y_jp1[flame.variable_index_dict["T"]] + + J_jphalf = get_tmp(cache.diffusivefluxes_jphalf, first(y_jp1)) + diffusive_flux!(J_jphalf, cs_jp1, cs_j, z_jp1, z_j, T_jp1, T_j, P, P, dz_jp1, dz_j, flame.diffusive_flux_model, flame, cache) + J_jmhalf = get_tmp(cache.diffusivefluxes_jmhalf, first(y_jm1)) + diffusive_flux!(J_jmhalf, cs_j, cs_jm1, z_j, z_jm1, T_j, T_jm1, P, P, dz_j, dz_jm1, flame.diffusive_flux_model, flame, cache) + @views residual[flame.indexes[1]:flame.indexes[2]] .+= (J_jphalf .- J_jmhalf) ./ dz_j + + J_j = get_tmp(cache.diffusivefluxes_j, first(cs)) + J_j .= (J_jphalf .+ J_jmhalf) ./ 2 + Cpave = Cvave + R + residual[flame.variable_index_dict["T"]] += -dot(J_j, cpdivR) * R * (T_jp1 - T_jm1) / (z_jp1 - z_jm1) / (N * Cpave) + + lambda_jphalf = mixtureaveragedthermalconductivity(flame.phase, cs_jp1, T_jp1, P) + lambda_jmhalf = mixtureaveragedthermalconductivity(flame.phase, cs_jm1, T_jm1, P) + residual[flame.variable_index_dict["T"]] += 2 / (z_jp1 - z_jm1) * (lambda_jphalf * (T_jp1 - T_j) / (z_jp1 - z_j) - lambda_jmhalf * (T_j - T_jm1) / (z_j - z_jm1)) / (N * Cpave) +end + +function diffusive_flux!(J_jphalf, cs_jp1, cs_j, z_jp1, z_j, T_jp1, T_j, P_jp1, P_j, dz_jp1, dz_j, diffusive_flux_model::MixtureAveragedDiffusiveFlux, flame::N1, cache::FlameSimulationCache) where N1 <: AxisymmetricFlame + C_jp1 = P_jp1 / (R * T_jp1) + C_j = P_j / (R * T_j) + C = (C_jp1 * dz_jp1 + C_j * dz_j) / (dz_jp1 + dz_j) + cs = (cs_jp1 * dz_jp1 + cs_j * dz_j) / (dz_jp1 + dz_j) + mixdiffs = get_tmp(cache.mixdiffusivities, first(cs)) + T = (T_jp1 * dz_jp1 + T_j * dz_j) / (dz_jp1 + dz_j) + P = (P_jp1 * dz_jp1 + P_j * dz_j) / (dz_jp1 + dz_j) + mixtureaverageddiffusivities!(mixdiffs, flame.phase, C, cs, T, P) + + J_jphalf .= - mixdiffs .* (cs_jp1 - cs_j) / (z_jp1 - z_j) + J_jphalf .-= cs ./ C .* sum(J_jphalf) # correction terms to ensure species conservation +end + +function masscontinuity!(residual, y_j, y_jm1, y_jp1, flame::BurnerFlame) + u_j = y_j[flame.variable_index_dict["u"]] + u_jm1 = y_jm1[flame.variable_index_dict["u"]] + @views cs_j = y_j[flame.indexes[1]:flame.indexes[2]] + @views cs_jm1 = y_jm1[flame.indexes[1]:flame.indexes[2]] + rho_j = getrho(cs_j, flame.molecularweights) + rho_jm1 = getrho(cs_jm1, flame.molecularweights) + residual[flame.variable_index_dict["u"]] = rho_j * u_j - rho_jm1 * u_jm1 +end + +function inletboundary!(residual, y_1, z_1, dz_1, flame::BurnerFlame, cache::FlameSimulationCache) + cs_0 = flame.inlet.cs + u_0 = flame.u0 + T_0 = flame.inlet.T + @views cs_1 = y_1[flame.indexes[1]:flame.indexes[2]] + u_1 = y_1[flame.variable_index_dict["u"]] + rho_1 = getrho(cs_1, flame.molecularweights) + T_1 = y_1[flame.variable_index_dict["T"]] + P = flame.P + + J_1mhalf = get_tmp(cache.diffusivefluxes_jmhalf, first(y_1)) + diffusive_flux!(J_1mhalf, cs_1, cs_0, z_1, z_1-dz_1, T_1, T_0, P, P, dz_1, dz_1, flame.diffusive_flux_model, flame, cache) + @views residual[flame.indexes[1]:flame.indexes[2]] .= u_0 .* cs_0 .- J_1mhalf .- u_1 .* cs_1 + residual[flame.variable_index_dict["T"]] = T_1 - T_0 + residual[flame.variable_index_dict["u"]] = rho_1 * u_1 - flame.mdot +end + +function outletboundary!(residual, y_N, z_N, y_Nm1, z_Nm1, flame::F) where F <: AxisymmetricFlame + @views cs_N = y_N[flame.indexes[1]:flame.indexes[2]] + T_N = y_N[flame.variable_index_dict["T"]] + u_N = y_N[flame.variable_index_dict["u"]] + rho_N = getrho(cs_N, flame.molecularweights) + @views cs_Nm1 = y_Nm1[flame.indexes[1]:flame.indexes[2]] + T_Nm1 = y_Nm1[flame.variable_index_dict["T"]] + u_Nm1 = y_Nm1[flame.variable_index_dict["u"]] + rho_Nm1 = getrho(cs_Nm1, flame.molecularweights) + + residual[flame.indexes[1]:flame.indexes[2]] .= (cs_N .- cs_Nm1) ./ (z_N - z_Nm1) + residual[flame.variable_index_dict["T"]] = (T_N .- T_Nm1) ./ (z_N - z_Nm1) + residual[flame.variable_index_dict["u"]] = rho_N * u_N - rho_Nm1 * u_Nm1 +end + +function f_flame!(residual, y, p, t, flame::F, cell_centers::Array{Float64,1}, cell_sizes::Array{Float64,1}, cache::FlameSimulationCache) where F <: AxisymmetricFlame + num_cells = length(cell_centers) + for j in 2:num_cells-1 + ns, cs, T, P, V, C, N, G, kfs, krevs, Hs, Us, Gs, diffs, Cvave, cpdivR, phi = calcthermo(flame, y[:, j], t, p, cell_sizes[j]) + @views reaction!(residual[:, j], flame; t=t, ns=ns, cs=cs, T=T, P=P, V=V, C=C, N=N, kfs=kfs, krevs=krevs, Hs=Hs, Us=Us, Cvave=Cvave) + @views diffusion!(residual[:, j], y[:, j-1], y[:, j+1], cell_centers[j], cell_centers[j-1], cell_centers[j+1], cell_sizes[j], cell_sizes[j-1], cell_sizes[j+1], flame, cache; cs, T, P, N, Cvave, cpdivR) + # @views convection!(residual[:, j], y[:, j], y[:, j-1], cell_centers[j], cell_centers[j-1], flame) + @views masscontinuity!(residual[:, j], y[:, j], y[:, j-1], y[:, j+1], flame) + end + @views inletboundary!(residual[:, 1], y[:, 1], cell_centers[1], cell_sizes[1], flame, cache) + @views outletboundary!(residual[:, end], y[:, end], cell_centers[end], y[:, end-1], cell_centers[end-1], flame) +end From 81441b6dec10c81a8fbbc53e835bae1c569cdaed Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:32 -0400 Subject: [PATCH 09/16] Update authors --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 7de0771a..3eda6e5f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ReactionMechanismSimulator" uuid = "c2d78dd2-25c4-5b79-bebc-be6c69dd440f" -authors = ["Matt Johnson "] +authors = ["Matt Johnson ", "Hao-Wei Pang "] version = "0.4.0" [deps] From bafa40a59f7dc7dab97e9e0182026015e3da2b52 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:39 -0400 Subject: [PATCH 10/16] Add packages --- Project.toml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Project.toml b/Project.toml index 3eda6e5f..41dfcaf3 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" @@ -26,11 +27,13 @@ QuartzImageIO = "dca85d43-d64c-5e67-8c65-017450d5d020" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710" SmoothingSplines = "102930c3-cf33-599f-b3b1-9a29a5acab30" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -50,6 +53,7 @@ IncompleteLU = "0.2.0" IterTools = "1.3.0" LsqFit = "0.12" ModelingToolkit = "8" +NonlinearSolve = "1" OrdinaryDiffEq = "6" Parameters = "0.12" PreallocationTools = "0" @@ -59,9 +63,11 @@ QuartzImageIO = "0.7" RecursiveArrayTools = "2.17" ReverseDiff = "1.9" SciMLBase = "1" +SciMLNLSolve = "0" SmoothingSplines = "0.3" SpecialFunctions = "1" StaticArrays = "1" +SteadyStateDiffEq = "1" Sundials = "4" Symbolics = "4" Tracker = "0.2" From 42208c899b677b32fbb964843621ee0fe91a994d Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:04:50 -0400 Subject: [PATCH 11/16] Add new modules --- src/Parse.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Parse.jl b/src/Parse.jl index 95e444b9..f7ae7f5a 100644 --- a/src/Parse.jl +++ b/src/Parse.jl @@ -7,7 +7,10 @@ module Calc include("Calculators/RateUncertainty.jl") include("Calculators/ThermoUncertainty.jl") include("Calculators/Thermo.jl") + include("Calculators/Potential.jl") + include("Calculators/Transport.jl") include("Calculators/Diffusion.jl") + include("Calculators/ThermalConductivity.jl") include("Calculators/Rate.jl") include("Calculators/Viscosity.jl") include("Calculators/Thermovec.jl") @@ -17,7 +20,10 @@ end module Spc include("Calculators/ThermoUncertainty.jl") include("Calculators/Thermo.jl") + include("Calculators/Potential.jl") + include("Calculators/Transport.jl") include("Calculators/Diffusion.jl") + include("Calculators/ThermalConductivity.jl") include("Calculators/kLAkH.jl") include("Species.jl") end @@ -25,6 +31,8 @@ module Rxn include("Calculators/RateUncertainty.jl") include("Calculators/ThermoUncertainty.jl") include("Calculators/Thermo.jl") + include("Calculators/Potential.jl") + include("Calculators/Transport.jl") include("Calculators/Diffusion.jl") include("Calculators/Rate.jl") include("Calculators/kLAkH.jl") From 8900761c7450a199adcc7048ddee15880d2440d3 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:05:00 -0400 Subject: [PATCH 12/16] Add molecular weight --- src/Parse.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/Parse.jl b/src/Parse.jl index f7ae7f5a..2c8ec04e 100644 --- a/src/Parse.jl +++ b/src/Parse.jl @@ -330,6 +330,10 @@ function readinputyml(fname::String) end end + if "transport" in keys(d) + d["transport"]["m"] = d["molecularweight"]/Na + end + spc = fcndict2obj(d,ymlunitsdict) push!(spclist,spc) if haskey(spcdict,spc.name) From a197f7cc57b25b2f5caf2c9976dac8bbed863ca2 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:05:21 -0400 Subject: [PATCH 13/16] Add molar weights --- src/Phase.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Phase.jl b/src/Phase.jl index 784ed5c8..5ff180d9 100644 --- a/src/Phase.jl +++ b/src/Phase.jl @@ -35,6 +35,7 @@ include("Reaction.jl") reversibility::Array{Bool,1} forwardability::Array{Bool,1} diffusionlimited::Bool = false + molecularweights::Array{Float64, 1} end function IdealGas(species,reactions; name="",diffusionlimited=false) @@ -50,10 +51,12 @@ function IdealGas(species,reactions; name="",diffusionlimited=false) electronchange = convert(Array{Float64,1},echangevec) reversibility = getfield.(rxns,:reversible) forwardability = getfield.(rxns,:forwardable) + molecularweights = getfield.(species,:molecularweight) return IdealGas(species=species,reactions=rxns,name=name, spcdict=Dict([sp.name=>i for (i,sp) in enumerate(species)]),stoichmatrix=M,Nrp=Nrp,rxnarray=rxnarray,veckinetics=vectuple, veckineticsinds=posinds, vecthermo=therm, otherreactions=otherrxns, electronchange=electronchange, - reversibility=reversibility,forwardability=forwardability,diffusionlimited=diffusionlimited,) + reversibility=reversibility,forwardability=forwardability,diffusionlimited=diffusionlimited, + molecularweights=molecularweights) end export IdealGas From d9d789526c5a91339f803cfa89a1950f3dabb430 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:05:27 -0400 Subject: [PATCH 14/16] Add new modules --- src/ReactionMechanismSimulator.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/ReactionMechanismSimulator.jl b/src/ReactionMechanismSimulator.jl index 3cecfa24..58a1c1f7 100644 --- a/src/ReactionMechanismSimulator.jl +++ b/src/ReactionMechanismSimulator.jl @@ -38,7 +38,10 @@ module ReactionMechanismSimulator include("Calculators/RateUncertainty.jl") include("Calculators/ThermoUncertainty.jl") include("Calculators/Thermo.jl") + include("Calculators/Potential.jl") + include("Calculators/Transport.jl") include("Calculators/Diffusion.jl") + include("Calculators/ThermalConductivity.jl") include("Calculators/Rate.jl") include("Calculators/Viscosity.jl") include("Calculators/Thermovec.jl") @@ -63,4 +66,8 @@ module ReactionMechanismSimulator include("Debugging.jl") include("Plotting.jl") include("fluxdiagrams.jl") + include("Flame1D/Diffusion1D.jl") + include("Flame1D/Interface1D.jl") + include("Flame1D/Flame1D.jl") + include("Flame1D/Reactor1D.jl") end From c1b90c555d595292e8e2cc1d34ba7db85f33e388 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:05:38 -0400 Subject: [PATCH 15/16] Add transport model to species --- src/Species.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Species.jl b/src/Species.jl index f8e086fb..12b5d3ee 100644 --- a/src/Species.jl +++ b/src/Species.jl @@ -5,7 +5,7 @@ export AbstractSpecies struct EmptySpecies <: AbstractSpecies end -struct Species{T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawConstant,N2<:AbstractLiquidVolumetricMassTransferCoefficient} <: AbstractSpecies +struct Species{T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawConstant,N2<:AbstractLiquidVolumetricMassTransferCoefficient,N3<:AbstractTransportModel} <: AbstractSpecies name::String index::Integer inchi::String @@ -23,6 +23,7 @@ struct Species{T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawCons comment::String isfragment::Bool isfragmentintermediate::Bool + transport::N3 end function Species(; name::String, index::Integer, inchi::String="", smiles::String="", adjlist::String="", @@ -31,8 +32,9 @@ function Species(; name::String, index::Integer, inchi::String="", smiles::Strin henrylawconstant::N1=EmptyHenryLawConstant(), liquidvolumetricmasstransfercoefficient::N2=EmptyLiquidVolumetricMassTransferCoefficient(), comment::String="",isfragment::Bool=false,isfragmentintermediate::Bool=false, - ) where {T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawConstant,N2<:AbstractLiquidVolumetricMassTransferCoefficient} - return Species(name, index, inchi, smiles, adjlist, thermo, atomnums, bondnum, diffusion, radius, radicalelectrons, molecularweight, henrylawconstant, liquidvolumetricmasstransfercoefficient, comment, isfragment, isfragmentintermediate) + transport::N3=EmptyTransportModel() + ) where {T<:AbstractThermo,N<:AbstractDiffusivity,N1<:AbstractHenryLawConstant,N2<:AbstractLiquidVolumetricMassTransferCoefficient, N3 <: AbstractTransportModel} + return Species(name, index, inchi, smiles, adjlist, thermo, atomnums, bondnum, diffusion, radius, radicalelectrons, molecularweight, henrylawconstant, liquidvolumetricmasstransfercoefficient, comment, isfragment, isfragmentintermediate, transport) end export Species From 4c25db9199c5f108dcf77b6e51dad5aa5a052963 Mon Sep 17 00:00:00 2001 From: hwpang Date: Fri, 9 Jun 2023 15:05:47 -0400 Subject: [PATCH 16/16] Add superminimal with transport model --- src/testing/superminimal_transport.rms | 663 +++++++++++++++++++++++++ 1 file changed, 663 insertions(+) create mode 100644 src/testing/superminimal_transport.rms diff --git a/src/testing/superminimal_transport.rms b/src/testing/superminimal_transport.rms new file mode 100644 index 00000000..a44d89b1 --- /dev/null +++ b/src/testing/superminimal_transport.rms @@ -0,0 +1,663 @@ +Phases: +- Species: + - name: Ar + radicalelectrons: 0 + smiles: '[Ar]' + thermo: + polys: + - Tmax: 1000.0 + Tmin: 200.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 4.37967 + type: NASApolynomial + - Tmax: 6000.0 + Tmin: 1000.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 4.37967 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.8845952811748727e-21 + sigma: 3.3300000000000004e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: He + radicalelectrons: 0 + smiles: '[He]' + thermo: + polys: + - Tmax: 1000.0 + Tmin: 200.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 0.928724 + type: NASApolynomial + - Tmax: 6000.0 + Tmin: 1000.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 0.928724 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.4082633281870966e-22 + sigma: 2.5760000000000007e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: Ne + radicalelectrons: 0 + smiles: '[Ne]' + thermo: + polys: + - Tmax: 1000.0 + Tmin: 200.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 3.35532 + type: NASApolynomial + - Tmax: 6000.0 + Tmin: 1000.0 + coefs: + - 2.5 + - 0.0 + - 0.0 + - 0.0 + - 0.0 + - -745.375 + - 3.35532 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 2.051646378123555e-21 + sigma: 3.758e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: N2 + radicalelectrons: 0 + smiles: N#N + thermo: + polys: + - Tmax: 1000.0 + Tmin: 200.0 + coefs: + - 3.53101 + - -0.000123661 + - -5.02999e-07 + - 2.43531e-09 + - -1.40881e-12 + - -1046.98 + - 2.96747 + type: NASApolynomial + - Tmax: 6000.0 + Tmin: 1000.0 + coefs: + - 2.95258 + - 0.0013969 + - -4.92632e-07 + - 7.8601e-11 + - -4.60755e-15 + - -923.949 + - 5.87189 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 1.7600000000000008e-30 + potential: + epsilon: 1.3465524862708356e-21 + sigma: 3.6210000000000007e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: H2 + radicalelectrons: 0 + smiles: '[H][H]' + thermo: + polys: + - Tmax: 1959.07 + Tmin: 100.0 + coefs: + - 3.43536393 + - 0.000212711953 + - -2.78628671e-07 + - 3.40270013e-10 + - -7.76039045e-14 + - -1031.35983 + - -3.90841661 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 1959.07 + coefs: + - 2.78818509 + - 0.000587615921 + - 1.5902213e-07 + - -5.52762536e-11 + - 4.3432812e-15 + - -596.155632 + - 0.112618494 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 7.900000000000003e-31 + potential: + epsilon: 5.246488890790464e-22 + sigma: 2.9199999999999997e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: O2 + radicalelectrons: 0 + smiles: '[O][O]' + thermo: + polys: + - Tmax: 1074.55 + Tmin: 100.0 + coefs: + - 3.53732243 + - -0.00121571647 + - 5.31620254e-06 + - -4.89446434e-09 + - 1.45846258e-12 + - -1038.58849 + - 4.68368183 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 1074.55 + coefs: + - 3.15382081 + - 0.00167804371 + - -7.69974236e-07 + - 1.51275462e-10 + - -1.08782414e-14 + - -1040.81728 + - 6.16755832 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 1.600000000000001e-30 + potential: + epsilon: 1.4828229409723015e-21 + sigma: 3.458000000000001e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: HO2 + radicalelectrons: 1 + smiles: '[O]O' + thermo: + polys: + - Tmax: 932.15 + Tmin: 100.0 + coefs: + - 4.04594488 + - -0.00173464779 + - 1.03766518e-05 + - -1.02202522e-08 + - 3.34908581e-12 + - -986.754245 + - 4.63581294 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 932.15 + coefs: + - 3.21023857 + - 0.00367941991 + - -1.27701572e-06 + - 2.18045259e-10 + - -1.46337935e-14 + - -910.368497 + - 8.1829188 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.4828229409723015e-21 + sigma: 3.458000000000001e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: OH + radicalelectrons: 1 + smiles: '[OH]' + thermo: + polys: + - Tmax: 1145.76 + Tmin: 100.0 + coefs: + - 3.51456841 + - 2.92732126e-05 + - -5.32149895e-07 + - 1.01947438e-09 + - -3.85939087e-13 + - 3414.25418 + - 2.10434749 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 1145.76 + coefs: + - 3.0719371 + - 0.000604020067 + - -1.39807196e-08 + - -2.13440522e-11 + - 2.4806113e-15 + - 3579.38798 + - 4.57801549 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.1045239770085187e-21 + sigma: 2.7500000000000003e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: H + radicalelectrons: 1 + smiles: '[H]' + thermo: + polys: + - Tmax: 4879.8 + Tmin: 100.0 + coefs: + - 2.5 + - -3.01680531e-12 + - 3.74582141e-15 + - -1.50856878e-18 + - 1.86626471e-22 + - 25474.2178 + - -0.444972899 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 4879.8 + coefs: + - 4.28461071 + - -0.00145494649 + - 4.44804306e-07 + - -6.04359642e-11 + - 3.07921551e-15 + - 23723.0923 + - -11.8931307 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 2.001945556980982e-21 + sigma: 2.0500000000000002e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: H2O + radicalelectrons: 0 + smiles: O + thermo: + polys: + - Tmax: 1130.24 + Tmin: 100.0 + coefs: + - 4.05763557 + - -0.000787932878 + - 2.90876518e-06 + - -1.47517698e-09 + - 2.12838441e-13 + - -30281.5866 + - -0.311363105 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 1130.24 + coefs: + - 2.84325231 + - 0.00275108244 + - -7.81029811e-07 + - 1.07243254e-10 + - -5.79389106e-15 + - -29958.6136 + - 5.91040933 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 6.1509219154539245e-30 + polarizability: 0.0 + potential: + epsilon: 7.90286938760371e-21 + sigma: 2.6050000000000003e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: H2O2 + radicalelectrons: 0 + smiles: OO + thermo: + polys: + - Tmax: 908.87 + Tmin: 100.0 + coefs: + - 3.73136061 + - 0.00335067714 + - 9.35045149e-06 + - -1.52101308e-08 + - 6.41593098e-12 + - -17721.1709 + - 5.4590992 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 908.87 + coefs: + - 5.41578065 + - 0.00261009268 + - -4.39898683e-07 + - 4.91103613e-11 + - -3.35202076e-15 + - -18302.9497 + - -4.02244574 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.4828229409723015e-21 + sigma: 3.458000000000001e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + - name: O(T) + radicalelectrons: 2 + smiles: '[O]' + thermo: + polys: + - Tmax: 4879.8 + Tmin: 100.0 + coefs: + - 2.5 + - -3.01680531e-12 + - 3.74582141e-15 + - -1.50856878e-18 + - 1.86626471e-22 + - 29230.2441 + - 4.09104286 + type: NASApolynomial + - Tmax: 5000.0 + Tmin: 4879.8 + coefs: + - 4.28461071 + - -0.00145494649 + - 4.44804306e-07 + - -6.04359642e-11 + - 3.07921551e-15 + - 27479.1187 + - -7.35711496 + type: NASApolynomial + type: NASA + transport: + dipolemoment: 0.0 + polarizability: 0.0 + potential: + epsilon: 1.1045239770085187e-21 + sigma: 2.7500000000000003e-10 + type: LennardJonesPotential + type: TransportModel + type: Species + name: phase +Reactions: +- kinetics: + A: 54500.00000000001 + Ea: 6276.0 + n: 0.0 + type: Arrhenius + products: + - H2 + radicalchange: -2 + reactants: + - H + - H + type: ElementaryReaction +- kinetics: + A: 290000000.00000006 + Ea: 236981.76 + n: 0.0 + type: Arrhenius + products: + - HO2 + - H + radicalchange: 2 + reactants: + - H2 + - O2 + type: ElementaryReaction +- kinetics: + A: 87900.00000000001 + Ea: 1882.8 + n: 1.0 + type: Arrhenius + products: + - HO2 + radicalchange: 0 + reactants: + - O2 + - H + type: ElementaryReaction +- kinetics: + A: 1.2350200000000002 + Ea: 25463.824 + n: 1.634 + type: Arrhenius + products: + - H2 + - HO2 + radicalchange: 0 + reactants: + - H + - H2O2 + type: ElementaryReaction +- kinetics: + A: 5005182.000000001 + Ea: 7096.063999999999 + n: 0.282 + type: Arrhenius + products: + - H2O2 + radicalchange: -2 + reactants: + - HO2 + - H + type: ElementaryReaction +- kinetics: + A: 17500.000000000004 + Ea: -13702.6 + n: 0.0 + type: Arrhenius + products: + - O2 + - H2O2 + radicalchange: -2 + reactants: + - HO2 + - HO2 + type: ElementaryReaction +- kinetics: + A: 7850000.000000001 + Ea: 0.0 + n: 0.0 + type: Arrhenius + products: + - H2O2 + radicalchange: -2 + reactants: + - OH + - OH + type: ElementaryReaction +- kinetics: + A: 1820.0000000000002 + Ea: 83972.88 + n: 1.21 + type: Arrhenius + products: + - H + - H2O + radicalchange: 0 + reactants: + - H2 + - OH + type: ElementaryReaction +- kinetics: + A: 162000000.00000003 + Ea: 627.6 + n: 0.0 + type: Arrhenius + products: + - H2O + radicalchange: -2 + reactants: + - OH + - H + type: ElementaryReaction +- kinetics: + A: 9300000.000000002 + Ea: 310118.08 + n: 0.0 + type: Arrhenius + products: + - HO2 + - OH + radicalchange: 2 + reactants: + - O2 + - H2O + type: ElementaryReaction +- kinetics: + A: 0.4994995000000001 + Ea: 26463.8 + n: 1.927 + type: Arrhenius + products: + - HO2 + - H2O + radicalchange: 0 + reactants: + - OH + - H2O2 + type: ElementaryReaction +- kinetics: + A: 43.58386000000001 + Ea: 5167.240000000001 + n: 1.88 + type: Arrhenius + products: + - HO2 + radicalchange: -2 + reactants: + - OH + - O(T) + type: ElementaryReaction +- kinetics: + A: 10000000.000000002 + Ea: 0.0 + n: 0.0 + type: Arrhenius + products: + - OH + radicalchange: -2 + reactants: + - H + - O(T) + type: ElementaryReaction +- kinetics: + A: 3680.3920000000007 + Ea: 24681.416 + n: 0.678 + type: Arrhenius + products: + - O2 + - OH + radicalchange: -2 + reactants: + - HO2 + - O(T) + type: ElementaryReaction +- kinetics: + A: 340.00000000000006 + Ea: 96022.8 + n: 1.5 + type: Arrhenius + products: + - OH + - H + radicalchange: 0 + reactants: + - H2 + - O(T) + type: ElementaryReaction +- kinetics: + A: 17400000.000000004 + Ea: 19874.0 + n: 0.0 + type: Arrhenius + products: + - HO2 + - OH + radicalchange: 0 + reactants: + - H2O2 + - O(T) + type: ElementaryReaction +- kinetics: + A: 5260.000000000001 + Ea: 74600.71999999999 + n: 1.2 + type: Arrhenius + products: + - OH + - OH + radicalchange: 0 + reactants: + - H2O + - O(T) + type: ElementaryReaction +Units: {}