A new approach for simulating inhomogeneous chemical kinetics

In this paper, inhomogeneous chemical kinetics are simulated by describing the concentrations of interacting chemical species by a linear expansion of basis functions in such a manner that the coupled reaction and diffusion processes are propagated through time efficiently by tailor-made numerical methods. The approach is illustrated through modelling \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document}α- and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma$$\end{document}γ-radiolysis in thin layers of water and at their solid interfaces from the start of the chemical phase until equilibrium was established. The method’s efficiency is such that hundreds of such systems can be modelled in a few hours using a single core of a typical laptop, allowing the investigation of the effects of the underlying parameter space. Illustrative calculations showing the effects of changing dose-rate and water-layer thickness are presented. Other simulations are presented which show the approach’s capability to solve problems with spherical symmetry (an approximation to an isolated radiolytic spur), where the hollowing out of an initial Gaussian distribution is observed, in line with previous calculations. These illustrative simulations show the generality and the computational efficiency of this approach to solving reaction-diffusion problems. Furthermore, these example simulations illustrate the method’s suitability for simulating solid-fluid interfaces, which have received a lot of experimental attention in contrast to the lack of computational studies.

the system of PDEs is given as where ρ(x,t) is the concentration of the the chemical species given as the subscript, k is the reaction rate coefficient with reactants given as the subscript, and products given as the superscripts.The first term in each of the PDEs describes the diffusion of that species through the medium, and R(ρ(x,t)) can be written in the form We can now form a stoichiometric matrix, M, and reaction vector, v, where each column in M represents a reaction, and each row represents a single chemical species.Each element in M is denoted by ±n, where n is the number of the product (+) or reactant (−) present in the reaction.Each row in v represents the reaction rate constants and species for each reaction.Now the reaction vector R(ρ(x,t)) can be represented as M • v.
The Jacobian matrix, J R (ρ(x,t)), required for the calculation of a single time step using Kahan's method can now be rewritten as M • J v (ρ(x,t)), where J v (ρ(x,t)) i, j = ∂ v i ∂ ρ(x,t) j .For this example, J v (ρ(x,t)) is given by

Discrete Transforms
Discrete transforms are used to transform from coefficient space to value space, whilst inverse discrete transforms are used to transform from value space to coefficient space.Given the concentration distribution for an arbitrary basis for N distinct spatial values x 0 , x 1 , ..., x N , this can be rewritten in the form of a matrix multiplication where Given a basis {φ n } N n=0 , the transform can be made fast if we chose points {x k } N k=0 in which the inverse of Φ can be calculated without the need for linear algebra.Below is an example using the trigonometric basis with Dirichlet boundary conditions at x = 0 and Neumann boundary conditions at x = L, (SinDN basis), All spatial points in our software package are defined as N+1 .Choosing N = 1 for simplicity, (S7) becomes where sin( 9π 8 ) = − sin( π 8 ).The inverse, Φ −1 , can be found as This can be rewritten as showing that for this example Φ −1 = 1 Det(Φ) Φ.
Readily available algorithms can perform DCT and DST in O(N log N) operations instead of the O(N 2 ) operations required by a naive approach 1 .
For implementation in our software, we choose the discrete transform best suited to our basis functions.Given the concentration distribution using the SinDN basis and the spatial points 2 ) N , we can rewrite (S13) as where subscript k refers to the k th point in the spatial domain.
We can now compare (S14) to the library of discrete sine transforms 2 , and see that it takes the form of the Type-IV Discrete Sine Transform (DST-IV) and so we can express the transform from coefficient space to value space as To transform from value space to coefficient space, the inverse of the function given in (S16) must be found.As DST-IV is the inverse of itself 2 , this is simply given by Analogous processes need to be applied for each basis.S2.Diffusion coefficients and reaction rate coefficients used for all simulations in the plutonium stewardship section.These simulations all use step size ∆t = 10 − 3, N = 100 spectral terms for water thicknesses, L, of 1-20 monolayers where 1 monolayer is assumed to have thickness of 0.25nm.Initial distributions for all species were set to ρ(x) = 0.

Hollowing out Effect
Table S3 provides all information required to recreate the hollowing out effect simulations.S3.This table contains all parameters gathered from Burns et al 5 to run the hollowing out effect simulations.The first two panels give information on the initial conditions used, where b 2 = 2r 2 0 at t = 0.The bottom two panels provide the diffusion coefficients and reaction rate coefficients required to simulate the system.Our simulations used an adaptive time stepping algorithm with 4606 timesteps, N = 40 spectral terms and a spatial domain of r = [0,100] nm.

Table S1 .
Table collecting all data used to simulate incident radiation as zeroth order reactions.Note that the H 3 O + G-value missing from Spinks and Woods is assumed to be the same as the e − aq G-value due to comparison with Kreipl et al.TableS2provides all other information required to recreate these simulations .