Counter diffusion and reaction in a slab - Multiphase Reactor Modeling using PyMRM

Reaction: $A+B\rightarrow *$

$\mathrm{div}(-D \, \mathrm{grad} (c_A) ) = -k \, c_A c_B$ , $\mathrm{div}(-D \, \mathrm{grad} (c_B)) = -k \, c_A c_B$ ,

with boundary conditions: $c_A(0)=1$ , $\frac{\partial c_A(L)}{\partial n} = 0$ , $\frac{\partial c_B(0)}{\partial n} = 0$ , $c_B(L)=1$

Imports¶

import numpy as np
from scipy.sparse import linalg as sla
import matplotlib.pyplot as plt
from pymrm import construct_grad, construct_div, NumJac

Parameters and boundary conditions¶

Physical and numerical parameters, and the boundary dictionaries in the a*dc/dn + b*c = d convention (outward normal): a fixed concentration $c_A=1$ at $x=0$ (with $c_B=0$ there) and the mirror image at $x=L$ .

# Physical parameters
num_c = 2
L = 1.0
D = 1
k = 500
bc_L = {'a': 0, 'b': 1, 'd': [[1, 0]]}
bc_R = {'a': 0, 'b': 1, 'd': [[0, 1]]}

# Numerical parameters
num_x = 100
shape = (num_x, num_c)

Reaction term¶

$A + B \rightarrow *$ with rate $k \, c_A c_B$ , consuming $A$ and $B$ in equal amounts.

def reaction(c, k):
    f = np.empty_like(c)
    r = k * c[:, 0] * c[:, 1]
    f[:, 0] = -r
    f[:, 1] = -r
    return f

Grid and operators¶

Uniform grid, the constant diffusion operator, and a NumJac object for the nonlinear reaction Jacobian.

x_f = np.linspace(0, L, num_x + 1)
x_c = 0.5 * (x_f[:-1] + x_f[1:])

numjac = NumJac(shape)

grad_mat, grad_bc = construct_grad(shape, x_f, x_c, bc=(bc_L, bc_R), axis=0)
div_mat = construct_div(shape, x_f, nu=0, axis=0)
jac_disp = div_mat @ (D * grad_mat)
jac_disp_bc = div_mat @ (D * grad_bc)

Solve (Newton iteration)¶

A fixed number of Newton iterations on the steady residual jac_disp @ c + jac_disp_bc + f_react = 0, with the reaction Jacobian rebuilt from NumJac at the current state each iteration.

c = np.zeros(shape)
for i in range(10):
    f_react, jac_react = numjac(lambda c: reaction(c, k), c)
    g = jac_disp @ c.reshape(-1, 1) + jac_disp_bc + f_react.reshape(-1, 1)
    lu = sla.splu(jac_disp + jac_react)
    c = c - lu.solve(g).reshape(shape)

Plot the concentration profiles¶

plt.plot(x_c, c[:, 0], label='Component 1')
plt.plot(x_c, c[:, 1], label='Component 2')
plt.xlabel('Position')
plt.ylabel('Concentration')
plt.legend()
plt.show()

Checking the result¶

The notebook stops after exactly 10 Newton iterations without checking convergence. pymrm.checks.residual_check judges the final state by a scale-free backward error on the same residual used above; this would catch a case where 10 iterations are not enough, for example with a much larger $k$ or with a Jacobian that does not match the residual and slows convergence down.

from pymrm.checks import residual_check


def residual(c):
    f_react, jac_react = numjac(lambda c: reaction(c, k), c.reshape(shape))
    g = jac_disp @ c.reshape(-1, 1) + jac_disp_bc + f_react.reshape(-1, 1)
    jac = jac_disp + jac_react
    return g, jac


check = residual_check(residual, c)
print(f"backward error {check['backward_error']:.3e}, ok={check['ok']}")
print("A tiny backward error confirms the 10 fixed iterations did converge for this k." if check['ok']
      else "The fixed iteration count did not converge: refine or add iterations.")

backward error 1.914e-16, ok=True
A tiny backward error confirms the 10 fixed iterations did converge for this k.