Reactor with surface chemistry and several phases

[Mikael Larsen]

Hi Cantera folks

I am studying calcination of lime in a hydrogen atmosphere

CaCO3 → CaO + CO2

I have produced my own data (attached)

I have defined the interphase of CaO(S) with CaCO3_s and CaO_s plus a gas phase as adjacent phases.

CaCO3_s and CaO_s are bulk phases with fixed stoichiometry.

Everything is placed in a IdealGasConstantPressureReactor.

When modelling I get the right gas composition approaching the equilibrium.

However, the CaCO3_s and CaO_s seem to be treated as reservoirs. I can not enter  amounts of them and r.get_state() only gives mass fraction of gas-species and coverages of interface species. I need to get the amounts of these two bulk phases.

I have tried to introduce a mixture in the reactor without success and also tried to use the insert function.

Here is the code

import numpy as np
import cantera as ct

p = ct.one_atm # pressure
T = 1080.0 # surface temperature

gas = ct.Solution("CaCO3-calcination-surface.yaml", "gas")
CaCO3_phase = ct.Solution("CaCO3-calcination-surface.yaml","CaCO3_s")
CaO_phase = ct.Solution("CaCO3-calcination-surface.yaml", "CaO_s")
surf_phase = ct.Interface("CaCO3-calcination-surface.yaml", "CaO_surf", adjacent=[gas, CaCO3_phase, CaO_phase])

surf_phase.TPX = T, p, 'CaO(S): 1, CaCO3(S):1, CO2(S):0 '
gas.TPX = T, p, 'H2:1, CO2:0'
CaCO3_phase.TP = T, p
CaO_phase.TP = T, p

r = ct.IdealGasConstPressureReactor(gas, energy='off')
surf = ct.ReactorSurface(surf_phase, r = r, A=1)
net = ct.ReactorNet([r])

states = ct.SolutionArray(gas, extra=['t','surf_coverages'])

print('{:10s} {:10s} {:10s} {:14s}'.format('t [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))

dt_max = 1
t_end = 1 * dt_max
t = 0
while t < t_end:
t = net.step()
states.append(r.thermo.state, t=t, surf_coverages=surf.coverages)
print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(net.time, r.T, r.thermo.P, r.thermo.u))

print(gas.species_names)
print(r.get_state())

# plot major and minor gas species separately
minor_idx = []
major_idx = []
for i, name in enumerate(gas.species_names):
mean = np.mean(states(name).X)
if mean >= 0.0001:
major_idx.append(i)
elif mean >= 1e-20:
minor_idx.append(i)

import matplotlib.pyplot as plt
plt.clf()
plt.subplot(2, 2, 1)
plt.plot(states.t, states.T)
plt.xlabel('Time (s)')
plt.ylabel('Temperature (K)')
plt.subplot(2, 2, 2)
for j in major_idx:
plt.plot(states.t, states.X[:,j], label=gas.species_name(j))
plt.legend(fontsize=8, loc = 'best')
plt.xlabel('Time (s)')
plt.ylabel('Mole Fraction')
plt.subplot(2, 2, 3)
for k, i in enumerate(minor_idx):
plt.plot(states.t, states.X[:,i], label=gas.species_name(i))
plt.legend(fontsize=8, loc = 'best')
plt.xlabel('Time (s)')
plt.ylabel('Mole Fraction')
plt.subplot(2, 2, 4)
for i, name in enumerate(surf_phase.species_names):
plt.plot(states.t, states.surf_coverages[:, i], label=name)
plt.legend(fontsize=8, loc = 'best')
plt.xlabel('Time (s)')
plt.ylabel('Site Fraction')

plt.tight_layout()
plt.show()

[Ray Speth]

Hi,

This is an interesting question. The way to add these additional bulk phases is to introduce new reactor objects for each phase, and add these reactors to the network. The one slightly tricky requirement is that you also need to define a new ReactorSurface connecting those reactors to the surface phase:

r = ct.IdealGasConstPressureReactor(gas, energy='off') rCaCO3 = ct.ConstPressureReactor(CaCO3_phase, energy='off') rCaO = ct.ConstPressureReactor(CaO_phase, energy='off') rCaCO3.volume = 0.0005 rCaO.volume = 0.0005 surf = ct.ReactorSurface(surf_phase, r=r, A=1) surf = ct.ReactorSurface(surf_phase, r=rCaCO3, A=1) surf = ct.ReactorSurface(surf_phase, r=rCaO, A=1) net = ct.ReactorNet([r, rCaCO3, rCaO])

You can then modify the integration loop to track the mass of each phase:

masses = [] while t < t_end: t = net.step() states.append(r.thermo.state, t=t, surf_coverages=surf.coverages) masses.append((r.mass, rCaCO3.mass, rCaO.mass))

And plot the result at the end:

masses = np.array(masses) fig, ax = plt.subplots() ax.plot(states.t * 1000, masses[:,0], label='gas') ax.plot(states.t * 1000, masses[:,1], label='CaCO3') ax.plot(states.t * 1000, masses[:,2], label='CaO') ax.set(xlim=(0.0, 0.2), xlabel='time [ms]', ylabel='mass [kg]') ax.legend()

I’ll note that at present, I think this only works correctly as long as you use distinct Solution objects for each Reactor in the ReactorNet, due to the issue discussed in Enhancement #201. Also, there’s nothing that stops the reactions involving the solid phase species if the mass of one of those phases goes to zero, so you can end up with negative mass in those reactors. I’ll try to write up another Enhancement proposal to track the improvements that could be made here.

Regards,
Ray

​

[Mikael Larsen]

Hi Ray

Thanks at lot - that works beautifully.

However, I don't understand when trying to assign masses to the reactor for the bulk phases instead of volume (rCaCO3.mass = 1) I get the following error message "AttributeError: attribute 'mass' of 'cantera.reactor.ReactorBase' objects is not writable"

Regards

Mikael

[Mikael Larsen]

and thanks for   Enhancement #201. 

" This is the main problem: While the thermo property in Python does this synchronization, there's nothing that prevents the user from inadvertently accessing the Solution object directly without using the thermo property."

I learned it the hard way!

[Mikael Larsen]

I encountered a problem with this solution. Using a H2 atmosphere everything work fine, however switching to a N2 atmosphere I get the following error message when time reaches 1.5e-4 s:

Traceback (most recent call last):
  File "C:\Users\rml\PycharmProjects\ThermoChemistry\calc-surface1.py", line 81, in <module>
    surf_states.append(r.thermo.state, t=t, surf_coverages=surf.coverages, CaOS_prod=surf_prod[0], CaCO3S_prod=surf_prod[1], CO2S_prod=surf_prod[2])
  File "C:\Users\rml\AppData\Local\miniconda3\Lib\site-packages\cantera\composite.py", line 806, in append
    self._phase.state = state
    ^^^^^^^^^^^^^^^^^
  File "build\\python\\cantera\\thermo.pyx", line 1355, in cantera.thermo.ThermoPhase.state.__set__
cantera._utils.CanteraError:
*******************************************************************************
CanteraError thrown by Phase::assignDensity:
density must be positive. density = -20519257829809.63
*******************************************************************************

Until then the simulations are very similar in both cases - the first two graphs are simulations in H2 atmosphere and below is the simulations in N2.

Anyway, looking at the net_production_rates it is clear that the simulations struggles to go to a new steady state, which especially is seen for the net_production_rate for CO2(S) 

[Mikael Larsen]

I encountered a problem with this solution. Using a H2 atmosphere everything work fine, however .........

and here is the code

import numpy as np
import cantera as ct

p = ct.one_atm # pressure

T = 1073.0 # surface temperature

gas = ct.Solution("CaCO3-calcination-surface.yaml", "gas")
CaCO3_phase = ct.Solution("CaCO3-calcination-surface.yaml","CaCO3_s")
CaO_phase = ct.Solution("CaCO3-calcination-surface.yaml", "CaO_s")
surf_phase = ct.Interface("CaCO3-calcination-surface.yaml", "CaO_surf", adjacent=[gas, CaCO3_phase, CaO_phase])

surf_phase.TPX = T, p, 'CaO(S): 1, CaCO3(S):1, CO2(S):0 '

gas.TPX = T, p, 'N2:1, CO2:0'

CaCO3_phase.TP = T, p
CaO_phase.TP = T, p

r = ct.IdealGasConstPressureReactor(gas, energy='off')
rCaCO3 = ct.ConstPressureReactor(CaCO3_phase, energy='off')
rCaO = ct.ConstPressureReactor(CaO_phase, energy='off')

rCaCO3.volume = 0.0002
rCaO.volume = 0.0000000001

surf = ct.ReactorSurface(surf_phase, r=r, A=1)
surf = ct.ReactorSurface(surf_phase, r=rCaCO3, A=1)
surf = ct.ReactorSurface(surf_phase, r=rCaO, A=1)
net = ct.ReactorNet([r, rCaCO3, rCaO])

masses = []
gas_states = ct.SolutionArray(gas, extra=['t', 'r_prod', 'CO2_prod'])
surf_states = ct.SolutionArray(surf_phase, extra=['t', 'surf_coverages', 'CaOS_prod', 'CaCO3S_prod', 'CO2S_prod'])

print('{:10s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s} {:10s}'.format('t [s]', 'Mass [kg]', 'CaO(S)', 'CaCO3(S)', 'CO2(S)', 'CO2', 'CaCO3', 'CaO'))

t_end = 1
t = 0

while t < t_end:
t = net.step()

surf_fw_prog = surf_phase.forward_rates_of_progress # Array
surf_phase_prod = surf_phase.net_production_rates
surf_phase_conc = surf_phase.concentrations
gas_prod = r.kinetics.net_production_rates
surf_prod = surf.kinetics.net_production_rates
print(f' {net.time:10.3e} {r.mass:10.3f} {surf_prod[0]:10.3e} {surf_prod[1]:10.3e} {surf_prod[2]:10.3e} {surf_prod[18]:10.3e} {surf_prod[35]:10.3e} {surf_prod[36]:10.3e} ')
gas_states.append(r.thermo.state, t=t, r_prod=gas_prod[1], CO2_prod=gas_prod[15])

surf_states.append(r.thermo.state, t=t, surf_coverages=surf.coverages, CaOS_prod=surf_prod[0], CaCO3S_prod=surf_prod[1], CO2S_prod=surf_prod[2])

masses.append((r.mass, rCaCO3.mass, rCaO.mass))

# plot major and minor gas species separately

CO2 = gas.species_index('CO2')

minor_idx = []
major_idx = []
for i, name in enumerate(gas.species_names):

mean = np.mean(gas_states(name).X)

if mean >= 0.0001:
major_idx.append(i)
elif mean >= 1e-20:
minor_idx.append(i)

# minor_idx.append(CO2)

masses = np.array(masses)

import matplotlib.pyplot as plt
# plt.clf()
fig, (ax1, ax2, ax3, ax4, ax5, ax6) = plt.subplots(6, 1)

ax1.loglog(gas_states.t, masses[:, 0], label='gas')
ax1.loglog(gas_states.t, masses[:, 1], label='CaCO3')
ax1.loglog(gas_states.t, masses[:, 2], label='CaO')
ax1.set(xlabel='time [ms]', ylabel='mass [kg]')
ax1.legend()

for j in major_idx:
ax2.loglog(gas_states.t, gas_states.X[:, j], label=gas.species_name(j))
ax2.set_ylim([1e-10, 1])
ax2.legend(fontsize=8, loc ='best')
ax2.set_xlabel('Time (s)')
ax2.set_ylabel('Mole Fraction')

for k, i in enumerate(minor_idx):

ax3.loglog(gas_states.t, gas_states.X[:, i], label=gas.species_name(i))
ax3.set(ylim=(1e-20, 1e-10), xlabel='Time (s)', ylabel='Mole Fraction')
ax3.legend(fontsize=8, loc ='best')

ax4.loglog(gas_states.t, gas_states.r_prod, label=gas.species_name(1))
ax4.loglog(gas_states.t, gas_states.CO2_prod, label=gas.species_name(15))
ax4.legend(fontsize=8, loc = 'best')
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('Net production rates (kmol/m² s)')

ax5.loglog(surf_states.t, surf_states.CaOS_prod, label=surf_phase.species_name(0))
ax5.loglog(surf_states.t, surf_states.CaCO3S_prod, label=surf_phase.species_name(1))
ax5.loglog(surf_states.t, surf_states.CO2S_prod, label=surf_phase.species_name(2))
ax5.legend(fontsize=8, loc ='best')
ax5.set_xlabel('Time (s)')
ax5.set_ylabel('Net production rates (kmol/m² s)')

for i, name in enumerate(surf_phase.species_names):

ax6.loglog(surf_states.t, surf_states.surf_coverages[:, i], label=name)
ax6.legend(fontsize=8)
ax6.set_xlabel('Time (s)')
ax6.set_ylabel('Site Fraction')

# fig.tight_layout()
plt.show()