Skip to content

Commit 2f76997

Browse files
mjohnson541mjohnson541
authored
Merge pull request #194 from ReactionMechanismGenerator/two_phase_pyrolysis
Interfaces for two phase pyrolysis and bug fix
2 parents c5a836f + 956c03d commit 2f76997

8 files changed

Lines changed: 514571 additions & 24 deletions

File tree

‎iJulia/Vapor-Liquid Two-Phase Low-Temperature Pyrolysis.ipynb‎

Lines changed: 1659 additions & 0 deletions
Large diffs are not rendered by default.

‎src/Calculators/kLAkH.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ abstract type AbstractLiquidVolumetricMassTransferCoefficient end
2828
export AbstractLiquidVolumetricMassTransferCoefficient
2929

3030
struct EmptyLiquidVolumetricMassTransferCoefficient <: AbstractLiquidVolumetricMassTransferCoefficient end
31-
(eh::EmptyLiquidVolumetricMassTransferCoefficient)(;T::N) where {N<:Number} = 0
31+
(eh::EmptyLiquidVolumetricMassTransferCoefficient)(;T::N) where {N<:Number} = 0.0
3232
export EmptyLiquidVolumetricMassTransferCoefficient
3333

3434
@with_kw struct ConstantLiquidVolumetricMassTransferCoefficient{N<:Number} <: AbstractLiquidVolumetricMassTransferCoefficient

‎src/Domain.jl‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1503,6 +1503,14 @@ end
15031503
dydt[d.indexes[3]] -= inter.Vout(t)
15041504
end
15051505
end
1506+
for inter in interfaces
1507+
if isa(inter,VolumeMaintainingOutlet) && d == inter.domain #VolumeMaintainingOutlet has to be evaluated after dVdt has been modified by everything else
1508+
@inbounds dVdt = dydt[d.indexes[3]]
1509+
@inbounds flow = P*dVdt/(R*T)
1510+
@views @inbounds dydt[d.indexes[1]:d.indexes[2]] .-= flow * ns/N
1511+
@inbounds dydt[d.indexes[3]] -= dVdt
1512+
end
1513+
end
15061514
end
15071515

15081516
@inline function calcdomainderivatives!(d::ConstantVDomain{W,Y},dydt::K,interfaces::Z12;t::Z10,T::Z4,P::Z9,Us::Z,Hs::Z11,V::Z2,C::Z3,ns::Z5,N::Z6,Cvave::Z7) where {Z12,Z11,Z10,Z9,W<:IdealGas,Z7,K,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5}
@@ -1586,6 +1594,14 @@ end
15861594
dydt[d.indexes[4]] -= inter.Vout(t)
15871595
end
15881596
end
1597+
for inter in interfaces
1598+
if isa(inter,VolumeMaintainingOutlet) && d == inter.domain #VolumeMaintainingOutlet has to be evaluated after dVdt has been modified by everything else
1599+
@inbounds dVdt = dydt[d.indexes[4]]
1600+
@inbounds flow = P*dVdt/(R*T)
1601+
@views @inbounds dydt[d.indexes[1]:d.indexes[2]] .-= flow * ns/N
1602+
@inbounds dydt[d.indexes[4]] -= dVdt
1603+
end
1604+
end
15891605
end
15901606

15911607
@inline function calcdomainderivatives!(d::ParametrizedTPDomain{W,Y},dydt::K,interfaces::Z12;t::Z10,T::Z4,P::Z9,Us::Z,Hs::Z11,V::Z2,C::Z3,ns::Z5,N::Z6,Cvave::Z7) where {Z11,Z10,Z9,W<:IdealGas,Z7,K,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5,Z12}
@@ -1619,6 +1635,14 @@ end
16191635
dydt[d.indexes[3]] -= inter.Vout(t)
16201636
end
16211637
end
1638+
for inter in interfaces
1639+
if isa(inter,VolumeMaintainingOutlet) && d == inter.domain #VolumeMaintainingOutlet has to be evaluated after dVdt has been modified by everything else
1640+
@inbounds dVdt = dydt[d.indexes[3]]
1641+
@inbounds flow = P*dVdt/(R*T)
1642+
@views @inbounds dydt[d.indexes[1]:d.indexes[2]] .-= flow * ns/N
1643+
@inbounds dydt[d.indexes[3]] -= dVdt
1644+
end
1645+
end
16221646
end
16231647

16241648
@inline function calcdomainderivatives!(d::ParametrizedVDomain{W,Y},dydt::K,interfaces::Z12;t::Z10,T::Z4,P::Z9,Us::Z,Hs::Z11,V::Z2,C::Z3,ns::Z5,N::Z6,Cvave::Z7) where {Z11,Z10,Z9,W<:IdealGas,Z7,K,Y<:Integer,Z6,Z,Z2,Z3,Z4,Z5,Z12}
@@ -1705,6 +1729,14 @@ end
17051729
dydt[d.indexes[4]] -= inter.Vout(t)
17061730
end
17071731
end
1732+
for inter in interfaces
1733+
if isa(inter,VolumeMaintainingOutlet) && d == inter.domain #VolumeMaintainingOutlet has to be evaluated after dVdt has been modified by everything else
1734+
@inbounds dVdt = dydt[d.indexes[4]]
1735+
@inbounds flow = P*dVdt/(R*T)
1736+
@views @inbounds dydt[d.indexes[1]:d.indexes[2]] .-= flow * ns/N
1737+
@inbounds dydt[d.indexes[4]] -= dVdt
1738+
end
1739+
end
17081740
end
17091741
export calcdomainderivatives!
17101742

‎src/Interface.jl‎

Lines changed: 65 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -164,7 +164,7 @@ function getstoichmatrix(domain1,domain2,rxns)
164164
end
165165

166166
function getinterfacereactioninds(domain1,domain2,reactions)
167-
indices = zeros(Int64,(6,length(reactions)))
167+
indices = zeros(Int64,(8,length(reactions)))
168168
N1 = length(domain1.phase.species)
169169
for (i,rxn) in enumerate(reactions)
170170
for (j,r) in enumerate(rxn.reactants)
@@ -183,7 +183,7 @@ function getinterfacereactioninds(domain1,domain2,reactions)
183183
isfirst = false
184184
ind = findfirst(isequal(r),domain2.phase.species)
185185
end
186-
indices[j+3,i] = isfirst ? ind : ind+N1
186+
indices[j+4,i] = isfirst ? ind : ind+N1
187187
end
188188
end
189189
return indices
@@ -317,3 +317,66 @@ struct VolumetricFlowRateOutlet{V,F1<:Function} <: AbstractBoundaryInterface
317317
Vout::F1
318318
end
319319
export VolumetricFlowRateOutlet
320+
321+
"""
322+
VolumeMaintainingOutlet is designed for gas phase domain such that the flow rate of this outlet will adjust to maintain the volume of the
323+
domain to be constant. This is particularly useful to simulate any vapor-liquid phase system where the gas phase outlet
324+
is determined by the amount of evaporation.
325+
"""
326+
struct VolumeMaintainingOutlet{V} <: AbstractBoundaryInterface
327+
domain::V
328+
end
329+
330+
export VolumeMaintainingOutlet
331+
332+
struct VaporLiquidMassTransferInternalInterfaceConstantT{D1,D2,B} <: AbstractInternalInterface
333+
domaingas::D1
334+
domainliq::D2
335+
ignoremasstransferspcnames::Array{String,1}
336+
ignoremasstransferspcinds::B
337+
kLAs::Array{Float64,1}
338+
kHs::Array{Float64,1}
339+
parameterindexes::Array{Int64,1}
340+
domaininds::Array{Int64,1}
341+
p::Array{Float64,1}
342+
end
343+
344+
function VaporLiquidMassTransferInternalInterfaceConstantT(domaingas,domainliq,ignoremasstransferspcnames)
345+
@assert isa(domaingas.phase,IdealGas)
346+
@assert isa(domainliq.phase,IdealDiluteSolution)
347+
@assert getfield.(domaingas.phase.species,:name) == getfield.(domainliq.phase.species,:name)
348+
T = domainliq.T
349+
phase = domainliq.phase
350+
ignoremasstransferspcinds = getinterfaceignoremasstransferspcinds(domaingas,domainliq,ignoremasstransferspcnames)
351+
kLAs = [kLA(T=T) for kLA in getfield.(phase.species,:liquidvolumetricmasstransfercoefficient)]
352+
kHs = [kH(T=T) for kH in getfield.(phase.species,:henrylawconstant)]
353+
return VaporLiquidMassTransferInternalInterfaceConstantT(domaingas,domainliq,ignoremasstransferspcnames,ignoremasstransferspcinds,kLAs,kHs,[1,length(domainliq.phase.species)],[0,0],ones(length(domainliq.phase.species))),ones(length(domainliq.phase.species))
354+
end
355+
export VaporLiquidMassTransferInternalInterfaceConstantT
356+
357+
function getkLAkHs(vl::VaporLiquidMassTransferInternalInterfaceConstantT,Tgas,Tliq)
358+
return vl.kLAs, vl.kHs
359+
end
360+
361+
function evaluate(vl::VaporLiquidMassTransferInternalInterfaceConstantT,dydt,Vgas,Vliq,Tgas,Tliq,N1,N2,P1,P2,Cvave1,Cvave2,ns1,ns2,Us1,Us2,cstot,p)
362+
kLAs, kHs = getkLAkHs(vl,Tgas,Tliq)
363+
@views @inbounds @fastmath evap = kLAs.*cstot[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]]*Vliq
364+
@views @inbounds @fastmath cond = kLAs./kHs.*cstot[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]]*Vliq
365+
netevap = (evap .- cond)
366+
@simd for ind in vl.ignoremasstransferspcinds
367+
@inbounds netevap[ind] = 0.0
368+
end
369+
@views @inbounds @fastmath dydt[vl.domaingas.indexes[1]:vl.domaingas.indexes[2]] .+= netevap
370+
@views @inbounds @fastmath dydt[vl.domainliq.indexes[1]:vl.domainliq.indexes[2]] .-= netevap
371+
end
372+
export evaluate
373+
374+
function getinterfaceignoremasstransferspcinds(domaingas,domainliq,ignoremasstransferspcnames)
375+
indices = zeros(Int64,length(ignoremasstransferspcnames))
376+
spcnamesliq = getfield.(domainliq.phase.species,:name)
377+
for (i,name) in enumerate(ignoremasstransferspcnames)
378+
indliq = findfirst(isequal(name),spcnamesliq)
379+
indices[i] = domainliq.indexes[1]-1+indliq
380+
end
381+
return indices
382+
end

‎src/Reactor.jl‎

Lines changed: 47 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -103,14 +103,16 @@ function Reactor(domains::T,y0s::W1,tspan::W2,interfaces::Z=Tuple(),ps::X=SciMLB
103103
for i in 1:length(domain.constantspeciesinds)
104104
domain.constantspeciesinds[i] += k-1
105105
end
106-
for (thermovar,ind) in domain.thermovariabledict
107-
domain.thermovariabledict[thermovar] += k-1
108-
end
109106
domain.indexes[1] = k
110107
k += Nspcs
111108
domain.indexes[2] = k-1
112109
y0[domain.indexes[1]:domain.indexes[2]] = y0s[j][1:Nspcs]
113110
for m = 3:2+Ntherm
111+
for (therm,ind) in domain.thermovariabledict
112+
if ind == domain.indexes[m]
113+
domain.thermovariabledict[therm] = n
114+
end
115+
end
114116
domain.indexes[m] = n
115117
y0[n] = y0s[j][Nspcs+m-2]
116118
n -= 1
@@ -161,6 +163,15 @@ function Reactor(domains::T,y0s::W1,tspan::W2,interfaces::Z=Tuple(),ps::X=SciMLB
161163
inter.parameterindexes[2] = length(p)+length(ps[i+length(domains)])
162164
inter.rxnarray .= getinterfacereactioninds(inter.domain1,inter.domain2,inter.reactions)
163165
p = vcat(p,ps[i+length(domains)])
166+
elseif isa(inter, VaporLiquidMassTransferInternalInterfaceConstantT)
167+
indgas = findfirst(isequal(inter.domaingas),domains)
168+
indliq = findfirst(isequal(inter.domainliq),domains)
169+
inter.domaininds[1] = indgas
170+
inter.domaininds[2] = indliq
171+
inter.parameterindexes[1] = length(p)+1
172+
inter.parameterindexes[2] = length(p)+length(ps[i+length(domains)])
173+
inter.ignoremasstransferspcinds .= getinterfaceignoremasstransferspcinds(inter.domaingas,inter.domainliq,inter.ignoremasstransferspcnames)
174+
p = vcat(p,ps[i+length(domains)])
164175
end
165176
end
166177

@@ -394,15 +405,19 @@ export getrate
394405
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
395406
elseif @inbounds rarray[3,i] == 0
396407
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
397-
else
408+
elseif @inbounds rarray[4,i] == 0
398409
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
410+
else
411+
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]]
399412
end
400-
if @inbounds rarray[5,i] == 0
401-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
402-
elseif @inbounds rarray[6,i] == 0
403-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
413+
if @inbounds rarray[6,i] == 0
414+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]
415+
elseif @inbounds rarray[7,i] == 0
416+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]
417+
elseif @inbounds rarray[8,i] == 0
418+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]
404419
else
405-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
420+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]]
406421
end
407422
@fastmath R = fR - rR
408423

@@ -462,29 +477,39 @@ end
462477
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]
463478
elseif @inbounds rarray[3,i] == 0
464479
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]
465-
else
480+
elseif @inbounds rarray[4,i] == 0
466481
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]
482+
else
483+
@inbounds @fastmath fR = kfs[i]*cs[rarray[1,i]]*cs[rarray[2,i]]*cs[rarray[3,i]]*cs[rarray[4,i]]
467484
end
468-
if @inbounds rarray[5,i] == 0
469-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]
470-
elseif @inbounds rarray[6,i] == 0
471-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]
485+
if @inbounds rarray[6,i] == 0
486+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]
487+
elseif @inbounds rarray[7,i] == 0
488+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]
489+
elseif @inbounds rarray[8,i] == 0
490+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]
472491
else
473-
@inbounds @fastmath rR = krevs[i]*cs[rarray[4,i]]*cs[rarray[5,i]]*cs[rarray[6,i]]
492+
@inbounds @fastmath rR = krevs[i]*cs[rarray[5,i]]*cs[rarray[6,i]]*cs[rarray[7,i]]*cs[rarray[8,i]]
474493
end
475494
@fastmath R = (fR - rR)*V
476495
@inbounds @fastmath dydt[rarray[1,i]] -= R
477496
if @inbounds rarray[2,i] != 0
478497
@inbounds @fastmath dydt[rarray[2,i]] -= R
479498
if @inbounds rarray[3,i] != 0
480499
@inbounds @fastmath dydt[rarray[3,i]] -= R
500+
if @inbounds rarray[4,i] != 0
501+
@inbounds @fastmath dydt[rarray[4,i]] -= R
502+
end
481503
end
482504
end
483-
@inbounds @fastmath dydt[rarray[4,i]] += R
484-
if @inbounds rarray[5,i] != 0
485-
@inbounds @fastmath dydt[rarray[5,i]] += R
486-
if @inbounds rarray[6,i] != 0
487-
@inbounds @fastmath dydt[rarray[6,i]] += R
505+
@inbounds @fastmath dydt[rarray[5,i]] += R
506+
if @inbounds rarray[6,i] != 0
507+
@inbounds @fastmath dydt[rarray[6,i]] += R
508+
if @inbounds rarray[7,i] != 0
509+
@inbounds @fastmath dydt[rarray[7,i]] += R
510+
if @inbounds rarray[8,i] != 0
511+
@inbounds @fastmath dydt[rarray[8,i]] += R
512+
end
488513
end
489514
end
490515
end
@@ -604,6 +629,8 @@ end
604629
for (i,inter) in enumerate(interfaces)
605630
if isa(inter,AbstractReactiveInternalInterface)
606631
evaluate(inter,dydt,domains,vT[inter.domaininds[1]],vT[inter.domaininds[2]],vphi[inter.domaininds[1]],vphi[inter.domaininds[2]],vGs[inter.domaininds[1]],vGs[inter.domaininds[2]],cstot,p)
632+
elseif isa(inter,VaporLiquidMassTransferInternalInterfaceConstantT)
633+
evaluate(inter,dydt,vV[inter.domaininds[1]],vV[inter.domaininds[2]],vT[inter.domaininds[1]],vT[inter.domaininds[2]],vN[inter.domaininds[1]],vN[inter.domaininds[2]],vP[inter.domaininds[1]],vP[inter.domaininds[2]],vCvave[inter.domaininds[1]],vCvave[inter.domaininds[2]],vns[inter.domaininds[1]],vns[inter.domaininds[2]],vUs[inter.domaininds[1]],vUs[inter.domaininds[2]],cstot,p)
607634
end
608635
end
609636
for (i,domain) in enumerate(domains)

‎src/Simulation.jl‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,8 @@ function interofdomain(inter,domain)
105105
return domain==inter.domain
106106
elseif hasfield(typeof(inter),:domain1)
107107
return domain==inter.domain1 || domain==inter.domain2
108+
elseif hasfield(typeof(inter),:domaingas)
109+
return domain==inter.domaingas || domain==inter.domainliq
108110
end
109111
end
110112

@@ -684,7 +686,7 @@ calculate the rates of all reactions at time t
684686
function rates(ssys::Q,t::X) where {Q<:SystemSimulation,X<:Real}
685687
rts = zeros(length(ssys.reactions))
686688
domains = getfield.(ssys.sims,:domain)
687-
Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([length(inter.reactions) for inter in ssys.interfaces if hasproperty(inter,:reactions)])
689+
Nrxns = sum([length(sim.domain.phase.reactions) for sim in ssys.sims])+sum([hasproperty(inter,:reactions) ? length(inter.reactions) : 0 for inter in ssys.interfaces])
688690
Nspcs = sum([length(sim.domain.phase.species) for sim in ssys.sims])
689691
cstot = zeros(Nspcs)
690692
vns = Array{Any,1}(undef,length(domains))

‎src/TestReactors.jl‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,52 @@ end;
7272

7373
end;
7474

75+
@testset "Test vapor-liquid phase multi-domain reactor simulation with VaporLiquidMassTransferInternalInterfaceConstantT and VolumeMaintainingOutlet interface" begin
76+
input_file = "../src/testing/TdependentkLAkH.rms"
77+
phaseDict = readinput(input_file)
78+
liqspcs = phaseDict["phase"]["Species"]
79+
liqrxns = phaseDict["phase"]["Reactions"]
80+
solvent = phaseDict["Solvents"][1]
81+
gasspcs = liqspcs
82+
gasrxns = []
83+
liqspcnames = getfield.(liqspcs,:name)
84+
gasspcnames = getfield.(gasspcs,:name)
85+
86+
gas = IdealGas(gasspcs,gasrxns;name="gas");
87+
liq = IdealDiluteSolution(liqspcs,liqrxns,solvent;name="liq",diffusionlimited=true)
88+
89+
Vliq = 1.0
90+
Vgas = 1.0
91+
T = 25 + 273.15
92+
octaneconc = 6478
93+
tf = 3600*24
94+
95+
initialconds = Dict("octane"=>octaneconc*Vliq,"T"=>T,"V"=>Vliq)
96+
domainliq,y0liq,pliq = ConstantTVDomain(phase=liq,initialconds=initialconds);
97+
98+
P = 101300.0
99+
N2N = P*Vgas/R/T*0.8
100+
oxygenN = P*Vgas/R/T*0.2
101+
initialconds = Dict("N2"=>N2N,"oxygen"=>oxygenN,"T"=>T,"P"=>P)
102+
domaingas,y0gas,pgas = ConstantTPDomain(phase=gas,initialconds=initialconds);
103+
initialconds = Dict("N2"=>N2N,"oxygen"=>oxygenN,"T"=>T,"P"=>P)
104+
inletgas = Inlet(domaingas,initialconds,x->42)
105+
outletgas = VolumeMaintainingOutlet(domaingas)
106+
107+
ignoremasstransferspcnames = ["octane"] #only allowing octane to be consumed via reactions, to avoid the liquid phase to dry out completely
108+
vl,pinter = VaporLiquidMassTransferInternalInterfaceConstantT(domaingas,domainliq,ignoremasstransferspcnames);
109+
110+
domains = (domainliq,domaingas)
111+
interfaces = [vl,inletgas,outletgas]
112+
react,y0,p = Reactor(domains,(y0liq,y0gas),(0.0,tf),interfaces,(pliq,pgas,pinter));
113+
sol = solve(react.ode,react.recommendedsolver,abstol=1e-18,reltol=1e-6);
114+
115+
name = "oxygen"
116+
ind = findfirst(x->x==name,liqspcnames)
117+
@test sol(tf)[ind] ≈ 1.5533140456432624e-7 rtol=1e-5 #test there are oxygen dissolved into the liquid
118+
119+
end
120+
75121
@testset "Test liquid phase Parametrized T Constant V reactor jacobian" begin
76122
#Parametrized T constant V Ideal Dilute Liquid
77123
initialconds = Dict(["ts"=>[0.,600.,1200.],"T"=>[450.0,490.,500.],"V"=>1.0e-6*1e6,"octane"=>6.154e-3*1e6,"oxygen"=>4.953e-6*1e6])

0 commit comments

Comments
 (0)