Question on Multicomponent Reactive Transport Problem

[nirxg]

Hello everyone,

I’d like to start by saying thanks to the OpenPNM team for developing and maintaining such a helpful tool.

I am currently working on a simulation of methane reforming using OpenPNM. Below, I have attached the case setup and the flowchart of my solution procedure. I am using a segregated solver approach for this problem.

The issue I encountered is that the mole fraction field of hydrogen obtained from the simulation appears to be lower than expected compared to carbon monoxide. According to many references, the hydrogen mole fraction at the outlet is usually about three times higher than that of carbon dioxide, but in my results the hydrogen concentration seems significantly lower.

I would greatly appreciate any insights or suggestions regarding possible causes of this discrepancy or directions I might explore to resolve it. I would also be happy to provide more details or share files if that would help with troubleshooting.

Thank you very much for your time and support!

Best regards,
nirxg

import openpnm as op
import numpy as np
import matplotlib.pyplot as plt
from tqdm import tqdm

k1 = 1.0  # reaction rate constant
iters = 100  # number of iterations for the simulation

pn = op.network.Cubic(shape=[10, 10, 1], spacing=1.5e-1)
pn['pore.diameter'] = 1e-1
pn['throat.diameter'] = 1e-1
pn.add_model(propname='throat.diffusive_size_factors',
             model=op.models.geometry.diffusive_size_factors.spheres_and_cylinders)
pn.add_model(propname='throat.hydraulic_size_factors',
             model=op.models.geometry.hydraulic_size_factors.spheres_and_cylinders)

# CH4 + H2O -> CO + 3H2
CH4 = op.phase.StandardGas(network=pn, species='74-82-8', name='methane')
H2O = op.phase.StandardGas(network=pn, species='7732-18-5', name='water')
CO = op.phase.StandardGas(network=pn, species='630-08-0', name='carbon monoxide')
H2 = op.phase.StandardGas(network=pn, species='1333-74-0', name='hydrogen')

CH4['pore.mass_fraction'] = 0.4
H2O['pore.mass_fraction'] = 0.6
for phase in [CH4, H2O, CO, H2]:
    phase['pore.temperature'] = 543.0
    phase['pore.diffusivity'] = 5.0e-5
    phase.regenerate_models()
    phase.add_model(propname='throat.diffusive_conductance',
                    model=op.models.physics.diffusive_conductance.generic_diffusive)
    phase.add_model(propname='throat.hydraulic_conductance',
                    model=op.models.physics.hydraulic_conductance.generic_hydraulic) 
    
Mixture = op.phase.StandardGasMixture(network=pn, components=[CH4, H2O, CO, H2])
Mixture['pore.temperature'] = 543.0
for phase, init in zip(['methane', 'water', 'carbon monoxide', 'hydrogen'],
                       [0.4, 0.6, 0.0, 0.0]):
    Mixture.y(phase, init)
Mixture.regenerate_models()
Mixture.add_model(propname='throat.hydraulic_conductance',
                  model=op.models.physics.hydraulic_conductance.generic_hydraulic)

sf_Mixture = op.algorithms.StokesFlow(network=pn, phase=Mixture)
ad_CH4 = op.algorithms.AdvectionDiffusion(network=pn, phase=CH4)
ad_H2O = op.algorithms.AdvectionDiffusion(network=pn, phase=H2O)
ad_CO = op.algorithms.AdvectionDiffusion(network=pn, phase=CO)
ad_H2 = op.algorithms.AdvectionDiffusion(network=pn, phase=H2)

for phase in [ad_CH4, ad_H2O, ad_CO, ad_H2]:
    phase.settings['quantity'] = 'pore.mass_fraction'

inlet = pn.pores('left')
outlet = pn.pores('right')
reaction = list(set(pn.Ps) - set(inlet) - set(outlet))

for i in tqdm(range(iters)):

    Mixture.regenerate_models()
    sf_Mixture.set_value_BC(pores=inlet, values=1.0 + 101325.0, mode='overwrite')
    sf_Mixture.set_value_BC(pores=outlet, values=101325, mode='overwrite')
    sf_Mixture.run()

    for phase in [CH4, H2O, CO, H2]:
        # update ad_dif_conductance
        phase['pore.pressure'] = sf_Mixture['pore.pressure']
        phase.add_model(propname='throat.ad_dif_conductance',
                        model=op.models.physics.ad_dif_conductance.ad_dif,
                        s_scheme='exponential')
    
    # calculate the reaction rate (R = k1 * w_CH4/Mw_CH4 * w_H2O/Mw_H2O)
    R = k1 * CH4['pore.mass_fraction']/CH4.params['molecular_weight'] * \
        H2O['pore.mass_fraction']/H2O.params['molecular_weight']

    for alg, phase, bc, val, name in zip(
        [ad_CH4, ad_H2O, ad_CO, ad_H2],
        [CH4, H2O, CO, H2],
        [0.4, 0.6, 0.0, 0.0],
        [-1, -1, 1, 3],
        ['methane', 'water', 'carbon monoxide', 'hydrogen']
    ):

        model = op.models.physics.source_terms.linear
        phase['pore.A1'] = 0
        phase['pore.A2'] = val * R * phase.params['molecular_weight']
        phase.add_model(propname='pore.reaction',
                        model=model, A1='pore.A1', A2='pore.A2',
                        X='pore.mass_fraction', regen_mode='deferred')
    
        alg.set_value_BC(pores=inlet, values=bc, mode='overwrite')
        alg.set_outflow_BC(pores=outlet, mode='overwrite')
        alg.set_source(pores=reaction, propname='pore.reaction', mode='add')
        alg.run()
        Mixture[f'pore.mole_fraction.{name}'] = alg['pore.mass_fraction']
    
fluxin = ad_CH4.rate(pores=inlet)[0]
fluxout = -ad_CH4.rate(pores=outlet)[0]
print(f'conversion rate:{(fluxin - fluxout) / fluxin}')

phases = [CH4, H2O, CO, H2]
MW = np.array([p.params['molecular_weight'] for p in phases])   # kg/mol
W  = np.vstack([p['pore.mass_fraction'] for p in phases]).T     # (Np, ncomp)
n_tilde = W / MW                                               
den = n_tilde.sum(axis=1, keepdims=True) + 1e-30                
X = n_tilde / den                                              

for i, p in enumerate(phases):
    p['pore.mole_fraction_from_w'] = X[:, i]

nx, ny, nz = 10, 10, 1
def to_img(vec):
    if nz == 1:
        return vec.reshape(nx, ny)
    else:
        return vec.reshape(nx, ny, nz)[:, :, nz//2]

fig, axes = plt.subplots(1, len(phases) + 1, figsize=(4*(len(phases)+1), 4))
total_img = np.zeros((nx, ny))

for i, p in enumerate(phases):
    img = to_img(X[:, i])
    total_img += img
    im = axes[i].imshow(img.T, origin='lower', cmap='viridis')
    axes[i].set_title(f"{p.name} (mole fraction)")
    plt.colorbar(im, ax=axes[i], shrink=0.7)

im = axes[-1].imshow(total_img.T, origin='lower', cmap='viridis')
axes[-1].set_title("Sum of mole fractions")
plt.colorbar(im, ax=axes[-1], shrink=0.7)

plt.tight_layout()
plt.show()
print("max |sum(x)-1|:", float(np.max(np.abs(X.sum(axis=1) - 1))))

[jgostick]

Thank you very much for your clear and detailed question. Without running your code I can take a few guesses. Despite our efforts, we often find that 'reusing' an algorithm inside a for-loop can lead to unexpected results. Basically, not all the old values get overwritten. This is clearly something that we should focus on fixing, but in the meantime I suggest you try creating new algorithm objects inside the loop. You can remove the old ones from the Project to avoid memory leaks (hopefully?). Try this and see if it helps.

[nirxg]

Thank you very much for your suggestion. After trying it, I did observe some effect, but the results are still quite different from what I expected. I found that this is due to the influence of convection: if only diffusion is considered, the outlet hydrogen concentration would be 3 times that of carbon monoxide; with convection included, this ratio drops to around 1.6. Could you please let me know if this result aligns with reality?