WIP: flame

Project.toml

@@ -1,6 +1,6 @@
 name = "ReactionMechanismSimulator"
 uuid = "c2d78dd2-25c4-5b79-bebc-be6c69dd440f"
-authors = ["Matt Johnson <mjohnson541@gmail.com>"]
+authors = ["Matt Johnson <mjohnson541@gmail.com>", "Hao-Wei Pang <hao4wei3pang2@gmail.com>"]
 version = "0.4.0"
 
 [deps]
@@ -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"

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

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

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

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

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