# A physics-constrained deep learning surrogate model of the runaway electron
avalanche growth rate

Jonathan Arnaud<EMAIL_ADDRESS>Nuclear Engineering Program, Department of
Materials Science and Engineering, University of Florida, Gainesville, FL
32611, United States of America  Tyler Mark  Nuclear Engineering Program,
Department of Materials Science and Engineering, University of Florida,
Gainesville, FL 32611, United States of America  Christopher J. McDevitt
<EMAIL_ADDRESS>Nuclear Engineering Program, Department of Materials Science
and Engineering, University of Florida, Gainesville, FL 32611, United States
of America

(August 28, 2024)

###### Abstract

A surrogate model of the runaway electron avalanche growth rate in a magnetic
fusion plasma is developed. This is accomplished by employing a physics-
informed neural network (PINN) to learn the parametric solution of the adjoint
to the relativistic Fokker-Planck equation. The resulting PINN is able to
evaluate the runaway probability function across a broad range of parameters
in the absence of any synthetic or experimental data. This surrogate of the
adjoint relativistic Fokker-Planck equation is then used to infer the
avalanche growth rate as a function of the electric field, synchrotron
radiation and effective charge. Predictions of the avalanche PINN are compared
against first principle calculations of the avalanche growth rate with
excellent agreement observed across a broad range of parameters.

##  I Introduction

The unintentional generation of a large relativistic electron population
continues to pose a substantial obstacle to the success of the tokamak reactor
concept. Such runaway electrons (RE) may be inadvertently generated by the
strong electric fields coinciding with a tokamak disruption and obtain
energies of several MeV Hender et al. (2007). Due to their high energy and
often localized impact, REs have the potential to induce substantial damage to
plasma facing components Matthews et al. (2016). Obtaining a robust
description of their formation processes has thus emerged as a topic of
immediate importance to tokamak devices in addition to its intrinsic interest
to plasma physics.

One of the primary challenges in describing REs is the multi-physics nature of
a fusion plasma. In particular, alongside a description of RE kinetics, an
accurate treatment of RE formation requires the self-consistent evolution of
the background magnetohydrodynamic (MHD) equilibrium, impurity transport, and
radiative losses. As a result, a promising approach toward developing an
integrated description of RE formation and evolution involves the
identification of reduced models of RE kinetics that can be coupled to a
broader plasma physics framework Hoppe et al. (2021); Bandaru et al. (2021);
Liu et al. (2021); Sainterme and Sovinec (2024). Developing such a reduced RE
module has, however, posed a challenge to the plasma physics community. While
several reduced models of runaway generation processes such as Dreicer
generation Dreicer (1959); Kruskal and Bernstein (1962); Connor and Hastie
(1975); Hesslow et al. (2019a), hot tail generation Helander et al. (2004);
Smith and Verwichte (2008); Yang et al. (2023), or avalanche generation
Rosenbluth and Putvinski (1997); Martín-Solís et al. (2015); Aleynikov and
Breizman (2015); Hesslow et al. (2019b) have been developed, these reduced
models often struggle to quantitatively describe the inherently kinetic
physics that characterize RE generation. Furthermore, many of these models
were derived under highly simplified assumptions, where their generalization
to more realistic plasma conditions is non-trivial.

The rapidly evolving field of deep learning suggests a new pathway to
developing reduced RE models. While often computationally intensive to train,
the online deployment of deep learning based models is typically orders of
magnitude faster than traditional plasma physics codes thus providing an
efficient surrogate that may be called by a broader plasma physics framework.
In the present paper, our aim will be to develop a physics-informed neural
network (PINN) to provide an accurate reduced model of the RE avalanche. In
contrast to purely data driven paradigms, this deep learning approach embeds
the physics model into the training of a deep neural network, enabling the
PINN to make predictions in the absence of synthetic or experimental data. In
so doing, the trained PINN encodes the underlying kinetic solution, allowing
for greater interpretability, along with predicting quantities of interest
(QoI) such as the rate of RE generation.

This approach was recently used to predict the number of hot tail seed
electrons in an axisymmetric model of the thermal quench McDevitt (2023). Our
aim in the current paper is to develop a PINN to describe the avalanche
mechanism of REs SOKOLOV (1979); Jayakumar et al. (1993). In contrast to the
hot tail seed mechanism, the avalanche mechanism of RE formation requires a
pre-existing ‘seed’ RE population to be present. Once established, large-angle
collisions of this RE seed with the cold background plasma will result in the
initially cold electrons being scattered to energies up to half of the initial
energy of the seed electron. For a sufficiently large electric field, these
‘secondary’ electrons will be accelerated to relativistic energies resulting
in the exponential growth of the original seed. Such a process is of
particular importance to reactor scale tokamak plasmas due to the ability of
this mechanism to convert the Ohmic plasma current into RE current, even in
the limit where only a minuscule seed population of REs is present Martín-
Solís et al. (2017); Vallhagen et al. (2020).

The remainder of this paper is organized as follows: Section II provides a
brief overview of the physics-informed deep learning framework employed. The
adjoint of the relativistic Fokker-Planck equation, together with its
solution, is described in Sec. III. Section IV evaluates the avalanche growth
rate and threshold across a broad range of parameters and verifies the
predictions against those from a traditional RE solver. Conclusions along with
a brief discussion are given in Sec. V.

##  II Physics-Constrained Deep Learning

Physics-constrained deep learning methods have emerged as a powerful means to
efficiently describe complex physical processes. In addition to exploiting
available data, such methods seek to embed physical constraints into the
training of a neural network Lagaris et al. (1998); Karpatne et al. (2017);
Karniadakis et al. (2021); Lusch et al. (2018); Wang et al. (2020), thus
providing a natural means of avoiding overfitting, along with allowing for
greater generalizability to unseen parameter regimes. Physics-informed neural
networks Karniadakis et al. (2021) have emerged as a particularly prominent
example. A PINN, in its simplest form, incorporates the partial differential
equation (PDE), boundary, and initial conditions into the loss function,
yielding

Loss
$\displaystyle=\frac{1}{N_{PDE}}\sum^{N_{PDE}}_{i}\mathcal{R}^{2}\left(\mathbf{p}_{i},t_{i};\bm{\lambda}_{i}\right)+\frac{1}{N_{bdy}}\sum^{N_{bdy}}_{i}\left[P_{i}-P\left(\mathbf{p}_{i},t_{i};\bm{\lambda}_{i}\right)\right]^{2}$
$\displaystyle+\frac{1}{N_{init}}\sum^{N_{init}}_{i}\left[P_{i}-P\left(\mathbf{p}_{i},t=0;\bm{\lambda}_{i}\right)\right]^{2},$
(1)

where $P\left(\mathbf{p},t;\bm{\lambda}\right)$ is the dependent variable (the
runaway probability function for the present work),
$\mathcal{R}\left(\mathbf{p},t;\bm{\lambda}\right)$ is the residual of the
PDE, $\mathbf{p}$ and $t$ are the independent variables (momentum and time),
and $\bm{\lambda}$ represents parameters of the physical system. Here, the
first term represents the loss against the PDE, the second term represents the
loss against the boundary conditions, and the third term represents the loss
against the initial conditions for time-dependent problems. Noting that
derivatives of the neural network output $P$ with respect to its inputs
($\mathbf{p},t;\bm{\lambda}$) can be evaluated by automatic differentiation, a
feature provided by standard machine learning libraries Abadi et al. (2016);
Paszke et al. (2017), no discretization of the PDE is required. As a result,
PINNs are inherently mesh-free and only require specification of a
distribution of training points. After minimization of the loss function, the
dependent variable $P$ will satisfy the PDE, boundary, and initial conditions
up to the loss defined by Eq. (1). Hence, a PINN provides a means of solving
PDEs, where the value of the loss achieved after training provides an estimate
of the accuracy of the solution.

A powerful property of PINNs is that they may be used to learn the parametric
solution to a PDE Sun et al. (2020); McDevitt and Tang (2024). In particular,
since the parameters of the physical problem $\bm{\lambda}$ are inputs into
the neural network, after minimization of the loss, the PINN can predict
$P\left(\mathbf{p},t;\bm{\lambda}\right)$ across a broad range of parameters
$\bm{\lambda}$. While obtaining a parametric solution of a PDE often requires
extensive offline training, the online execution of a PINN is rapid, where an
individual prediction typically requires a few microseconds. Physics-informed
neural networks thus provide a framework for developing efficient surrogate
models of a PDE. A significant limitation of the above approach, however, is
that PINNs often fail to train when treating the challenging PDEs that
characterize many scientific and engineering applications Wang et al. (2022).
A primary aim of this paper will therefore be to develop custom output layers
to the PINN that enable it to robustly evaluate the adjoint to the
relativistic Fokker-Planck equation across a broad range of plasma conditions
in the absence of synthetic or experimental data.

When carrying out the training of the PINN, a single Nvidia A100 GPU will be
used. A fully connected feedforward neural network is employed, containing six
hidden layers, each having a width of 64 neurons. Roughly a million training
points are used and distributed across the input space
$\left(\mathbf{p},t,\bm{\lambda}\right)$ according to a Hammersley
distribution. During training a small number of training points are added to
regions where the residual is maximal. An independent set of test points are
also applied to verify accuracy of predictions away from training points. The
python script used for training the PINN is written using the DeepXDE library
Lu et al. (2021) with TensorFlow Abadi et al. (2016) as the backend.

##  III Steady State Runaway Probability Function

###  III.1 Adjoint Relativistic Fokker-Planck Equation

The adjoint to the relativistic Fokker-Planck equation has been treated in a
variety of contexts, including wave-driven currents in magnetized fusion
plasmas Antonsen Jr and Chu (1982); Taguchi (1983); Fisch (1986); Karney and
Fisch (1986) and runaway electron formation Karney and Fisch (1986); Liu et
al. (2016, 2017); Zhang and del Castillo-Negrete (2017); McDevitt (2023). In
the absence of synchrotron radiation, the adjoint to the steady state
relativistic Fokker-Planck equation can be written as:

$\displaystyle\left[-E_{\|}\xi-C_{F}\right]\frac{\partial P}{\partial
p}-\left(1-\xi^{2}\right)\frac{E_{\|}}{p}\frac{\partial
P}{\partial\xi}=-\frac{\nu_{D}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
P}{\partial\xi}\right],$ (2a) where the collisional coefficients are taken to
be: $\nu_{D}=\left(1+Z_{eff}\right)\frac{\gamma}{p^{3}},\quad
C_{F}=\frac{1+p^{2}}{p^{2}}.$ (2b)

Here, the relativistic momentum $p$ is normalized as $p\to
p/\left(m_{e}c\right)$, the electron’s pitch is defined by $\xi\equiv
p_{\|}/p$, time is normalized as $t\to t/\tau_{c}$, where $\tau_{c}\equiv
4\pi\epsilon^{2}_{0}m^{2}_{e}c^{3}/\left(e^{4}n_{e}\ln\Lambda\right)$ is the
collision time of a relativistic electron, the collisional coefficients
$\nu_{D}$ and $C_{F}$ are normalized to $\tau_{c}$, and the parallel electric
field is normalized to the Connor-Hastie electric field $E_{\|}\to
E_{\|}/E_{c}$, where $E_{c}\equiv m_{e}c/\left(e\tau_{c}\right)$ Connor and
Hastie (1975). Energy diffusion has been neglected due to this term being
exceptionally small for conditions typical of a tokamak disruption. In
particular, noting that the energy diffusivity scales with
$T_{e}/\left(m_{e}c^{2}\right)$ for the high energies characteristic of REs,
we expect this approximation to be well satisfied for a low temperature post
thermal quench plasma. Furthermore, the collision frequencies used in Eq. (2)
assume the limit $v>v_{Te}$, where $v_{Te}$ is the electron thermal velocity.
For the parameters of interest, the critical speed for an electron to run away
will be much larger than $v_{Te}$, hence we anticipate that this approximation
will not substantially impact our results.

Concerning the boundary conditions, we will enforce $P=0$ at $p=p_{min}$ and
$P=1$ along the upper boundary, where the energy flux $U_{p}\equiv-E_{\|}\xi-
C_{F}$ is positive. Specifically, the high energy boundary condition enforces
that $P$ is unity for values of the pitch $\xi$ where electric field
acceleration exceeds collisional drag such that the electron will be
accelerated out of the simulation domain. Further noting that for values of
the pitch near $\left|\xi\right|\approx 1$, Eq. (2a) is nearly hyperbolic, we
anticipate that electrons accelerated through the high energy boundary will be
accelerated to arbitrarily high energies, and thus have a zero probability of
returning to the simulation domain. Hence, any electron located at the high
energy boundary $p_{max}$ with $U_{p}>0$, is treated as a RE. With these
boundary conditions, the quantity $P\left(p,\xi\right)$ indicates the
probability that an electron initially located at $\left(p,\xi\right)$ will
run away at a later time, and is often referred to as the runaway probability
function (RPF) Karney and Fisch (1986).

When near the threshold electric field for avalanche generation, synchrotron
radiation substantially impacts the RPF. While adding synchrotron radiation
leads to a modest modification of Eq. (2), it does complicate the physical
interpretation of the RPF. In particular, with the inclusion of synchrotron
radiation, the adjoint equation takes the form:

$\displaystyle\left[-E_{\|}\xi-C_{F}-\alpha\gamma
p\left(1-\xi^{2}\right)\right]\frac{\partial P}{\partial
p}+\left(1-\xi^{2}\right)\left[-\frac{E_{\|}}{p}+\alpha\frac{\xi}{\gamma}\right]\frac{\partial
P}{\partial\xi}=-\frac{\nu_{D}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
P}{\partial\xi}\right],$ (3)

where the strength of synchrotron radiation is set by the parameter
$\alpha\equiv\tau_{c}/\tau_{s}$, where $\tau_{s}\equiv
6\pi\epsilon_{0}m^{3}_{e}c^{3}/\left(e^{4}B^{2}\right)$. Here, Eq. (3) is
solved similarly as before with a boundary condition of $P=0$ at $p=p_{min}$
along with the condition that $P=1$ on the high energy boundary $p=p_{max}$
when the energy flux, now defined by $U^{(\alpha)}_{p}\equiv-E_{\|}\xi-
C_{F}-\alpha\gamma p\left(1-\xi^{2}\right)$, is positive. While this problem
formulation is directly analogous to the case that neglects synchrotron
radiation, the physical interpretation of the RPF is slightly modified.
Specifically, as shown in Refs. Andersson et al. (2001); Decker et al. (2016);
Guo et al. (2017); McDevitt et al. (2018), synchrotron radiation damping and
pitch-angle scattering results in electrons obtaining a saturated energy,
achieved via the formation of a circulation pattern in momentum space centered
about an O-point. Thus, an electron accelerated through the high energy
boundary will have a non-zero probability of returning to the simulation
domain after a finite time, rather than being accelerated to arbitrarily high
energy as was the case when synchrotron radiation was neglected. For the case
with synchrotron radiation, the RPF should thus be given the narrower
interpretation as the probability that an electron reaches the high energy
boundary of the simulation domain before slowing down to the low energy bulk,
rather than the probability of an electron being accelerated to arbitrarily
high energy.

###  III.2 Embedding Physical Constraints into the PINN

Our aim in this section will be to develop a PINN framework capable of
robustly representing solutions to Eq. (3) across the three-dimensional
parameter space $\left(E_{\|},Z_{eff},\alpha\right)$. A key component of our
description is enforcing a subset of physical properties of the RPF as hard
constraints. In particular, we will (i) enforce the low energy boundary
condition $P=0$ at $p=p_{min}$, (ii) constrain the RPF to have a range between
zero and one, and (iii) recover the limit that the RPF vanishes when
$\left|E_{\|}\right|<1$. These three constraints are enforced by introducing a
customized output layer to the neural network of the form:

$P^{\prime}\left(p,\xi\right)\equiv\frac{1}{2}\left[1+\tanh\left(\frac{E_{\|}-1}{\Delta
E}\right)\right]\left(\frac{p-p_{min}}{p_{max}-p_{min}}\right)P_{NN}\left(p,\xi\right),$
(4a)
$P\left(p,\xi\right)\equiv\tanh\left({P^{\prime}}^{2}\left(p,\xi\right)\right).$
(4b)

Here, $P_{NN}$ is the output of the hidden layers of the neural network, and
$\Delta E$ is a hyperparameter whose value should satisfy $\Delta E<1$, where
for all cases in this paper it is taken to be $\Delta E=0.1$. From Eq. (4) it
can be verified that regardless of the value of $P_{NN}$, the predicted RPF
$P\left(p,\xi\right)$: (i) vanishes at the low energy boundary, (ii) has a
range of zero to one, and (iii) vanishes for $E_{\|}<1$ (only positive
electric fields are considered when training the RPF PINN). We note in passing
that the output layer defined by Eq. (4) results in both the value, and the
first derivative of $P$ vanishing at $p=p_{min}$. While this latter condition
is not required when defining the RPF, it will nevertheless be satisfied as
long as the value of $p_{min}$ chosen is below the momentum $p_{crit}\equiv
1/\sqrt{E_{\|}-1}$ where collisional drag exceeds electric field acceleration.
Specifically, since energy diffusion is neglected in the present analysis,
electrons located below $p_{crit}$ will have zero probability of running away.
Thus, for $p<p_{crit}$ the solution of the RPF will be a constant with a
magnitude of zero, and thus is consistent with pure Dirichlet and Neumann
boundary conditions at $p=p_{min}$

The loss function employed is taken to have the form:

$\text{Loss}=\frac{1}{N_{PDE}}\sum^{N_{PDE}}_{i}\left[\left(\frac{p^{2}_{i}}{1+p^{2}_{i}}\right)\mathcal{R}\left(p_{i},\xi_{i};\bm{\lambda}_{i}\right)\right]^{2}+\frac{1}{N_{bdy}}\sum^{N_{bdy}}_{i}\left[1-P\left(p_{i},\xi_{i};\bm{\lambda}_{i}\right)\right]^{2},$
(5)

where $\mathcal{R}\left(p_{i},\xi_{i};\bm{\lambda}_{i}\right)$ is the residual
of Eq. (3) and $\bm{\lambda}$ represents the physics parameters
$\left(E_{\|},Z_{eff},\alpha\right)$. Here, the first term in Eq. (5)
penalizes deviations from the governing PDE given by Eq. (3), where the
prefactor $p^{2}_{i}/\left(1+p^{2}_{i}\right)$ removes the low energy
divergence of the pitch-angle scattering operator. The second term in the loss
function defined by Eq. (5) enforces the high energy boundary condition, i.e.
$P=1$ at $p_{max}$ when $U^{(\alpha)}_{p}>0$. In particular, the boundary
points $N_{bdy}$, will only be applied at locations that satisfy both
$U^{(\alpha)}_{p}>0$ and $p=p_{max}$.

###  III.3 Parametric Dependence of the Runaway Probability Function

Figure 1: Loss history for a feedforward neural network with six hidden layers
each with a width of 64 neurons, along with roughly 1,000,000 training points.
15,000 epochs were performed with the ADAM optimizer, with the remaining steps
performed using L-BFGS. The model was trained across $E_{\|}\in(1,10)$,
$Z_{eff}\in(1,10)$, and $\alpha\in(0,0.2)$. The range of $p$ was chosen such
that the low energy boundary was $10\;\text{keV}$ and the high energy boundary
was $5\;\text{MeV}$.

In this section, we will seek to obtain solutions to the PINN in the 5D space
defined by the two independent coordinates $\left(p,\xi\right)$ and the three
physics parameters $\left(E_{\|},Z_{eff},\alpha\right)$. The loss history of
the PINN is shown in Fig. 1. Here, after roughly 100,000 epochs, the loss
associated with the PDE drops below $10^{-7}$, along with the boundary loss
dropping to $\approx 10^{-9}$, with the loss dropping more slowly for the
remaining $\sim 400,000$ epochs. The test loss of the PDE reaches $\approx
10^{-9}$ by the end of the training indicating that an accurate solution was
found. The periodic spikes in the training data are due to additional training
points sampled after every 100,000 epochs at locations where the residual is
maximal. The test loss is only updated after every 100,000 epochs when using
the L-BFGS optimizer, leading to the sharp variations evident in the dashed
curves.

Figure 2: Runaway probability functions for different values of the physics
parameters $\left(E_{\|},Z_{eff},\alpha\right)$. Panel (a) is for
$\left(E_{\|}=3,Z_{eff}=1,\alpha=0\right)$, panel (b) is for
$\left(E_{\|}=3,Z_{eff}=10,\alpha=0\right)$, panel (c) is for
$\left(E_{\|}=3,Z_{eff}=1,\alpha=0.2\right)$, and panel (d) is for
$\left(E_{\|}=10,Z_{eff}=1,\alpha=0\right)$.

Four example predictions of the RPF are shown in Fig. 2. Here, the RPF
vanishes at low energies, where drag exceeds electric field acceleration, but
increases at higher energy due to the $\left(1+p^{2}\right)/p^{2}$ drop in the
collisional drag. In particular, for the parameters indicated in Fig. 2(a),
the $P=0.5$ contour is located at an approximate energy of $200\;\text{keV}$
for $\xi=-1$. Considering a case with a large $Z_{eff}$ [see Fig. 2(b) with
$Z_{eff}=10$], the $P=0.5$ contour shifts to higher energy, with a more
gradual transition between the $P\approx 0$ and $P\approx 1$ regions. The
impact of synchrotron radiation on the RPF is shown in Fig. 2(c), where the
magnitude of synchrotron radiation was taken to be $\alpha=0.2$. Compared to
an otherwise identical case, but without synchrotron radiation [Fig. 2(a)], it
is evident that the location of the $P=0.5$ contour at $\xi=-1$ has only
shifted slightly. This is due to synchrotron radiation vanishing for $\xi=-1$
and having a modest magnitude at low energies. In contrast, the RPF at high
energies is more substantially impacted, where the region with $P\approx 1$ is
now largely localized to negative values of pitch for the energy range
considered. Finally, increasing the electric field to $E_{\|}=10$, with
$Z_{eff}$ = 1 and $\alpha=0$, results in a drop in the location of the $P=0.5$
contour [compare Figs. 2(a) and (d)], due to the electric field being able to
overcome collisional drag for a larger range of energies.

###  III.4 Critical Energy to Run Away

Figure 3: $\log_{10}E_{crit}$, where $E_{crit}$ is the critical energy to run
away in eV. Panel (a) is for $\alpha=0$, panel (b) is for $\alpha=0.1$, and
panel (c) is for $\alpha=0.2$.

The location of the $P=0.5$ contour $\left(p_{crit},\xi_{crit}\right)$
provides a useful reference point for identifying the critical energy and
pitch above which electrons are likely to run away. The parametric solution of
the PDE given by the PINN thus provides $\left(p_{crit},\xi_{crit}\right)$
across the entire domain by simply extracting the $P=0.5$ surface. To show the
dependence of ($E_{\|}$, $Z_{eff}$, $\alpha$) on the critical energy to run
away, we evaluate the energy at which $P=0.5$ and $\xi=-1$ across the entire
training region, which is shown in Fig. 3. Here, the $\log_{10}$ of the
critical energy in units of eV is plotted, where the empty regions represent
scenarios where $E_{\|}$ is below the avalanche threshold $E_{av}$ (see Sec.
IV.1 below). The dependence of $E_{\|}$ and $Z_{eff}$ on the critical energy
can be seen in Fig. 3 (a) for $\alpha=0$, where increasing the electric field
from $E_{\|}=2$ to $E_{\|}=10$ decreases the critical energy from $\sim$ 1 MeV
to $\sim$ 36 keV. At a given electric field (e.g. $E_{\|}=4$), the critical
energy increases by almost an order of magnitude as $Z_{eff}$ is increased
from one to ten, indicating that higher $Z$ elements in the plasma increase
the critical energy for REs. Moreover, the addition of synchrotron radiation
[compare Figs. 3 (b) and (c)], shows that the region where $E_{\|}<E_{av}$ is
increased (larger white region). As $\alpha$ is increased further [compare
Figs. 3 (a) and (c)], however, we also see that the shape and range of the
critical energy varies modestly, indicating the critical energy is modestly
impacted by $\alpha$ when $E_{\|}\gg 1$.

##  IV Surrogate Model of the Avalanche Growth Rate

###  IV.1 Secondary Source of Runaway Electrons

Given the rate and distribution of secondary electrons generated by large-
angle collisions, the RPF described in Sec. III can be used to evaluate the
avalanche growth rate. In particular, denoting the source of secondary
electrons as $S\left(p,\xi\right)$, the rate that REs form due to the
avalanche mechanism can be expressed as Liu et al. (2017):

$\left.\frac{dn_{RE}}{dt}\right|_{av}=\int
d^{3}pS\left(p,\xi\right)P\left(p,\xi\right),$ (6)

where $n_{RE}$ is the density of REs. Here, $S\left(p,\xi\right)$ indicates
the rate and momentum space distribution of REs, whereas $P\left(p,\xi\right)$
indicates the probability that an electron at a given momentum space location
$\left(p,\xi\right)$ runs away. By integrating over momentum space,
$S\left(p,\xi\right)P\left(p,\xi\right)$ will thus indicate the expected rate
that REs form due to the avalanche mechanism. The primary challenge with
evaluating Eq. (6) is due to $S\left(p,\xi\right)$ depending on the primary RE
distribution $f_{e}$, i.e.

$S\left(p,\xi\right)=\int
d^{3}p^{\prime}S_{0}\left(p^{\prime},\xi^{\prime},p,\xi\right)f_{e}\left(p^{\prime},\xi^{\prime}\right),$
(7)

where $S_{0}\left(p^{\prime},\xi^{\prime},p,\xi\right)$ is defined by:

$S_{0}\left(p^{\prime},\xi^{\prime},p,\xi\right)=n_{e}cr^{2}_{e}\frac{v^{\prime}}{2\pi
p^{2}}\frac{d\sigma_{M}\left(p^{\prime},p\right)}{dp}\Pi\left(p^{\prime},\xi^{\prime},p;\xi\right).$
(8)

where, $r_{e}=e^{2}/\left(4\pi\epsilon_{0}m_{e}c^{2}\right)$ is the classical
electron radius, $d\sigma_{M}/dp$ is the Møller cross section Møller (1932);
Ashkin et al. (1954), $\Pi\left(p^{\prime},\xi^{\prime},p;\xi\right)$
describes the pitch-angle dependence of secondary electron generation (see
Ref. Boozer (2015) for an explicit expression), $d^{3}p=2\pi p^{2}dpd\xi$, and
all variables have been dedimensionalized $p\to p/m_{e}c$, $v^{\prime}\to
v^{\prime}/c$, $\sigma_{M}\to\sigma_{M}/r^{2}_{e}$. Noting that the solution
of the adjoint relativistic Fokker-Planck equation does not directly yield the
RE distribution $f_{e}\left(p^{\prime},\xi^{\prime}\right)$, a closure
relation will need to be introduced to evaluate the rate of RE generation via
avalanching. The simplest closure, introduced in Ref. Rosenbluth and Putvinski
(1997), involves taking the limit where REs are assumed to have asymptotically
high energies and a pitch of $\xi=-1$. While idealized, this closure has been
shown to provide a good approximation to the full Møller source evaluated
using a self-consistently computed RE distribution McDevitt et al. (2018). In
the limit of $p^{\prime}\to\infty$ and $\xi^{\prime}=-1$, Eq. (7) asymptotes
to

$S\left(p,\xi\right)=n_{e}n_{RE}cr^{2}_{e}\frac{v}{\gamma^{2}-1}\frac{1}{\left(\gamma-1\right)^{2}}\delta\left(\xi-\xi_{1}\right),$
(9)

where we have introduced the Lorentz factor $\gamma\equiv\sqrt{1+p^{2}}$ and
$\xi_{1}$ is defined by:

$\xi_{1}=-\sqrt{\frac{\gamma-1}{\gamma+1}}.$ (10)

Using Eq. (9), Eq. (6) reduces to

$\left.\frac{dn_{RE}}{dt}\right|_{av}=2\pi n_{e}n_{RE}cr^{2}_{e}\int
dpp^{2}\frac{v}{\left(\gamma^{2}-1\right)\left(\gamma-1\right)^{2}}P\left(p,\xi_{1}\right).$
(11)

Noting that the right-hand side of Eq. (11) is directly proportional to the RE
density $n_{RE}$, this implies an exponentially growing solution with a growth
rate given by:

$\tau_{c}\gamma_{av}=\frac{1}{2\ln\Lambda}\int
dp\frac{v}{\left(\gamma-1\right)^{2}}P\left(p,\xi_{1}\right),$ (12)

where $\xi_{1}$ is defined by Eq. (10). Thus, once $P\left(p,\xi\right)$ has
been evaluated, the avalanche growth rate can be directly inferred from Eq.
(12).

A caveat when evaluating Eq. (12) is that while this integral formally extends
to $p\to\infty$, the RPF is evaluated assuming a finite $p_{max}$. The error
induced by this approximation can be estimated by considering the magnitude of
the integrand in Eq. (12) across the range of momenta used in this paper (see
Fig. 4, where $p_{min}\approx 0.2$ and $p_{max}\approx 10.74$). Here we have
taken $P=1$ inside the integrand of Eq. (12), such that Fig. 4 provides an
upper bound on the value of the integrand. Noting that the integrand has
decayed to a value of $\approx 3.5\times 10^{-4}$ at the upper boundary, this
implies a small contribution to the avalanche growth rate for secondary
electrons born with $p>p_{max}$, particularly for parameters where the system
is well above marginality.

Figure 4: The integrand $v/\left(\gamma-1\right)^{2}/\left(2\ln\Lambda\right)$
of Eq. (12) when $P=1$. The low and high energy bounds are $p_{min}\approx
0.2$ and $p_{max}\approx 10.74$, respectively. The Coulomb logarithm was taken
to be $\ln\Lambda=15$.

###  IV.2 Parametric Dependence of RE avalanche

Figure 5: Avalanche growth rate versus electric field for different values of
$Z_{eff}$ and $\alpha$. Panel (a) is for $\alpha=0$, panel (b) is for
$\alpha=0.1$, and panel (c) is for $\alpha=0.2$. The blue points in panel (b)
represent Monte Carlo simulations with $n_{e}=5\times 10^{21}$ m-3, and Eq.
(13) is shown as the dashed green curve. The Coulomb logarithm was taken to be
$\ln\Lambda=15$.

Using the RPF evaluated in Sec. III, Eq. (12) can be used to infer the
avalanche growth rate across the parameter space
$\left(E_{\|},Z_{eff},\alpha\right)$ (see Fig. 5). Here, the avalanche growth
rate increases approximately linearly with the electric field when $E_{\|}\gg
1$, with the slope and threshold sensitive to $Z_{eff}$ and $\alpha$.
Specifically, the value of $\alpha$ significantly impacts the avalanche
threshold (i.e. where $\gamma_{av}\approx 0$), but has a negligible impact at
large electric fields. In contrast, $Z_{eff}$ strongly impacts the avalanche
growth rate for all values of the electric field. A feature of Eq. (12) is
that since it indicates the number of REs generated via large-angle
collisions, it is thus positive definite. As a result, avalanche growth rates
predicted by Eq. (12) will not account for the decay of the RE population when
$E_{\|}<E_{av}$ and will instead asymptote to zero. While such behavior is
strictly correct when focusing solely on the avalanche growth mechanism, when
developing a RE model appropriate for coupling with a self-consistent MHD
solver it will be necessary to account for the decay of the RE distribution
when $E_{\|}<E_{av}$. The inclusion of this additional physics will be the
subject of future work.

It will be of interest to compare the predictions of the PINN with available
analytic theories. Considering the avalanche growth rate given by Ref.
Rosenbluth and Putvinski (1997):

$\tau_{c}\gamma_{RP}=\frac{1}{\ln\Lambda}\sqrt{\frac{\pi}{2\left(Z_{eff}+5\right)}}\left(E_{\|}-1\right),$
(13)

a comparison between the PINN predictions and Eq. (13) is shown in Fig. 5(b).
For reference, a small number of avalanche growth rates computed by a Monte
Carlo code McDevitt et al. (2019) using the complete Møller source are also
shown [the blue points on Fig. 5(b)]. Details on how this Monte Carlo data set
was generated are given in Sec. IV.3 below. It is apparent that both the PINN
and Eq. (13) yield results in reasonable agreement with the Monte Carlo
points, where the PINN yields more accurate results near threshold. In
particular, while Eq. (13) implies an avalanche threshold field of $E_{\|}=1$
regardless of the values of $\left(Z_{eff},\alpha\right)$, the PINN is able to
account for the variation of the avalanche threshold for non-zero values of
$\alpha$ and different $Z_{eff}$. In addition, at larger electric fields a
modest difference in the slope of the avalanche growth rate predicted by Eq.
(13) compared to the Monte Carlo data is evident. The predictions of the PINN,
in contrast, are able to recover the correct slope of the avalanche growth
rate, albeit with a modest offset in the magnitude as will be discussed
further in Sec. IV.3 below.

When well above marginality the avalanche growth rate is most conveniently
described by evaluating the amount of poloidal flux needed to increase the
amplitude of a RE seed by one order of magnitude Rosenbluth and Putvinski
(1997); Boozer (2018). This quantity can be evaluated by noting that for
$E_{\|}\gg 1$, the avalanche growth rate scales linearly with $E_{\|}$, i.e.

$\gamma_{av}=\gamma_{0}\left(\frac{E_{\|}}{E_{c}}-1\right)\approx\gamma_{0}\frac{E_{\|}}{E_{c}}\approx\frac{\gamma_{0}}{E_{c}}\frac{1}{R_{0}}\frac{\partial\psi}{\partial
t},$ (14)

where $\gamma_{0}$ is a constant that depends on
$\left(Z_{eff},\alpha,\ln\Lambda\right)$, $\ln\Lambda$ is the Coulomb
logarithm, $R_{0}$ is the major radius, and $\psi$ is the poloidal flux
function. After integrating Eq. (14) over the time interval $t_{f}-t_{i}$, the
number of exponentials of the RE population can be written in terms of the
change of poloidal flux, and a constant $\psi_{\text{exp}}$ that defines the
efficiency of the avalanche. In particular, integrating Eq. (14) over the time
interval $t_{f}-t_{i}$, yields

$N_{exp}=\int^{t_{f}}_{t_{i}}dt\gamma_{av}\approx\frac{\gamma_{0}}{E_{c}}\frac{1}{R_{0}}\int^{t_{f}}_{t_{i}}dt\frac{\partial\psi}{\partial
t}=\frac{\Delta\psi}{\psi_{\text{exp}}},$ (15)

where $\Delta\psi\equiv\psi\left(t_{f}\right)-\psi\left(t_{i}\right)$ and we
have defined $\psi_{exp}\equiv R_{0}E_{c}/\gamma_{0}$, which is related to the
amount of poloidal flux required to effect one exponential amplification of
the RE population. Thus, once $\gamma_{0}$ is inferred, the efficiency through
which the decay of the poloidal flux leads to an amplification of the seed RE
population can be evaluated. Further defining the quantity $\psi_{10}\equiv\ln
10\psi_{exp}$ (which delineates base ten amplifications), the efficiency of
the avalanche growth rate for a broad range of parameters is shown in Fig.
6(a). Here, the avalanche growth rate is most efficient for low values of
$Z_{eff}$ and increases by nearly a factor of two for $Z_{eff}=10$, agreeing
with previous results McDevitt et al. (2019). We have not indicated the
dependence of $\psi_{10}$ on $\alpha$ since for large electric fields, the
avalanche growth rate will be independent of $\alpha$.

Figure 6: (a) The value of $2\pi\psi_{10}/\mu_{0}R_{0}$ as a function of
$Z_{eff}$, with a Coulomb logarithm of $\ln\Lambda=15$. (b) Avalanche
threshold $E_{av}$ as a function of the synchrotron radiation strength
$\alpha$. The solid lines represent the PINN predictions, the dashed lines and
‘o’ markers represent Eq. (B15) and Monte Carlo results from Ref. McDevitt et
al. (2019), respectively, and WebPlotDigitizer was used to extract the Monte
Carlo values. The blue values represent $Z_{eff}=1$ and the red values
represent $Z_{eff}=5$. The Coulomb logarithm was taken to be $\ln\Lambda=20$.

In addition, the avalanche threshold is strongly impacted by the parameters
$\left(Z_{eff},\alpha\right)$. In particular, as $Z_{eff}$ or $\alpha$ are
increased, the threshold electric field $E_{av}$, where the avalanche growth
rate is zero, is increased. The dependence of the avalanche threshold on
$\left(Z_{eff},\alpha\right)$ is shown in Fig. 6(b), where the PINN
predictions are the solid lines. The Monte Carlo results are the ‘o’ markers,
and the dashed lines represent empirical fit formula given by Eq. (B15) of
Ref. McDevitt et al. (2019). Here, since the avalanche growth rates predicted
by the PINN are positive definite, we define the avalanche threshold as the
value of the electric field where $\tau_{c}\gamma_{av}=2\times 10^{-3}$ . From
Fig. 6(b) it is apparent that the PINN predicts $E_{av}$ particularly well
across $\alpha$ for $Z_{eff}=1$, with somewhat larger deviations evident for
$Z_{eff}=5$. The systematic under-prediction of the avalanche threshold is due
to the assumption in Eq. (12) that primary electrons have infinite energy.
Such an assumption overestimates the number of secondary electrons generated,
and perhaps more importantly, neglects that the primary electron population
itself will slowly decay to the thermal bulk when $E\approx E_{av}$.

###  IV.3 Verification of the Surrogate Model

In this section we will verify the PINN’s predictions of the avalanche growth
rate against first principle Monte Carlo simulations across the training
region ($E_{\|}$, $Z_{eff}$, $\alpha$). The RunAway Monte Carlo (RAMc) code is
employed (see Ref. McDevitt et al. (2019) for a detailed description), which
evolves the guiding center motion of relativistic electrons and includes
effects from small-angle collisions, large-angle collisions, and synchrotron
radiation. Large-angle collisions are evaluated using the full Møller source
Møller (1932). In order to avoid toroidal corrections to the avalanche growth
rate McDevitt and Tang (2019); Arnaud and McDevitt (2024), all REs are
initialized near $r=0$, and a large tokamak device was chosen with a minor
radius of $a=2\;\text{m}$ and major radius $R_{0}=6\;\text{m}$, in order to
render spatial transport negligible. A geometry with circular flux surfaces
was selected, with a safety factor profile taken to be
$q\left(r\right)=2.1+2\left(r/a\right)^{2}$.

Figure 7: (a) Avalanche growth rate comparison between the PINN and Monte
Carlo solver. The grey dashed line represents the coefficient of determination
of $R^{2}=1$. (b) Same as panel (a), but with the predictions of the avalanche
PINN multiplied by the factor $0.94886$. (c) Avalanche growth rate comparison
between Eq. (13) and the Monte Carlo solver. The avalanche growth rates were
evaluated across $E_{\|}\in$ (1,10), $Z_{eff}\in$ (1,10), and
$\alpha\in\left(2.8\times 10^{-3},0.2\right)$. The other parameters were
chosen to be $T_{e}$ = 10 eV and $n_{e}$ = 1021 m-3.

The Monte Carlo avalanche simulations are set up by initializing a small
population of electrons (8 for this analysis) at high momentum [$p\in$
(10,20)] and strongly aligned with the magnetic field [$\xi\in$ (-0.9,-1.0)].
They are then allowed to exponentially grow in time until a saturated growth
rate can be identified, generally resulting in the initial seed RE population
growing by several orders of magnitude. The plasma parameters chosen were a
density of $n_{e}=10^{21}$ m-3, a temperature of $T_{e}=10$ eV, and a toroidal
magnetic field strength that was varied to give the appropriate value of
$\alpha$. One caveat is that at a fixed density and temperature $\alpha\propto
B^{2}$, which in turn can lead to orbit drifts at the low boundary of
synchrotron radiation due to the weak magnetic field. These orbit drifts are
thus ensured to be negligible by enforcing a minimum toroidal magnetic field
of $B=$2 T, which corresponds to a lower bound ($\alpha\approx 2.8\times
10^{-3}$) in the training region for synchrotron radiation. Fifty randomly
selected samples of ($E_{\|}$, $Z_{eff}$, $\alpha$) are chosen as input
parameters for the Monte Carlo solver and the avalanche PINN, where cases
below the avalanche threshold are discarded, leaving forty-one samples. A
comparison between the predictions of the avalanche PINN and Monte Carlo
solver is shown in Fig. 7(a). Here, the grey dashed line represents equality
(y-axis equal to x-axis) between the PINN and Monte Carlo predictions. The
coefficient of determination Guilford and Guilford (1936) between the PINN
predictions and the Monte Carlo avalanche growth rate was evaluated to be
$R^{2}\approx 0.9849$, indicating that the PINN was able to accurately predict
the avalanche growth rate across a broad range of parameters. One noticeable
feature present in Fig. 7(a) is the systematic over prediction of the
avalanche growth rate by the PINN. This feature is expected, and is due to the
use of the simplified Rosenbluth-Putvinski secondary source term described by
Eq. (9), which assumes all primary electrons to have $p\to\infty$ and
$\xi=-1$. We note that this feature can be mitigated by introducing an order
unity factor to correct for the discrepancy between the PINN and Monte Carlo
predictions. After optimization, it was found that multiplying the predictions
of the PINN by a factor of 0.94886 significantly reduced the discrepancy
between the PINN and Monte Carlo predictions. After applying this factor [see
Fig. 7(b)], the resulting coefficient of determination increases to $R^{2}$ =
0.9988.

The PINN is also shown to robustly provide better avalanche growth rate
predictions than the analytical expression provided by Eq. (13) [compare Figs.
7(a) and (c)]. Here, the avalanche growth rate provided by Eq. (13) compared
with the Monte Carlo solver has a weaker correlation than that between the
PINN and Monte Carlo solver, where the coefficient of determination for Fig.
7(c) was $R^{2}\approx 0.9644$. The avalanche PINN thus improves on the
accuracy of the avalanche growth rate predictions, even in the absence of the
0.94886 factor.

##  V Summary and Conclusions

This work utilized a physics-informed neural network to evaluate the steady
state solution of the adjoint to the relativistic Fokker-Planck equation.
Noting that the PINN takes the physical parameters of the problem as inputs,
once trained the PINN provides an efficient surrogate model of the RPF. In
addition, while a comprehensive description of the avalanche growth rate
requires evaluating the primary electron distribution, a quantity not
evaluated in the present approach, by invoking the often employed
simplification that primary electrons have infinite energy and a pitch
$\xi=-1$ Rosenbluth and Putvinski (1997) the RE avalanche growth rate was
shown to be directly linked to an integral of the RPF [Eq. (12)]. While this
approximation to the secondary source of REs is known to lead to a modest
overestimate of the avalanche growth rate at modest values of the electric
field McDevitt et al. (2019), predictions from the PINN were shown to agree
with direct Monte Carlo simulations that utilized a complete Møller source
across a broad range of parameters with an offset of roughly 5% [Fig. 7(a)].
This small offset can be largely removed by multiplying predictions of the
PINN by a factor of order unity [compare Figs. 7(a) and 7(b)]. We note,
however, that such a correction will have little impact on the avalanche
threshold, and thus the present model predicts a smaller avalanche threshold
electric field compared to predictions using a full Møller source [see Fig.
6(b)], though it does recover the dominant dependencies on the physics
parameters $\left(Z_{eff},\alpha\right)$.

An additional aim of this paper was to provide a proof-of-principle
demonstration that physics constrained deep learning methods offer an
attractive avenue through which RE surrogate models can be developed. In
contrast to data driven deep learning approaches, the present method encodes
the underlying kinetic solution into the neural network, rather than just the
QoI (avalanche growth rate in this case), and thus allows for greater
interpretability and hence greater confidence in the accuracy of the
prediction. While the present paper has focused on an idealized description of
REs, generalization to more complete models of RE formation, incorporating
partial screening effects Hesslow et al. (2017) for example, can be
accomplished by modifying the collision coefficients in the adjoint equation
described by Eq. (3) and retraining the PINN. This extension will be the
subject of future work. We thus anticipate that the present approach provides
a flexible means through which RE surrogate models can be developed for a
broad range of plasma conditions.

###### Acknowledgements.

This work was supported by DOE OFES under award Nos. DE-SC0024634 and DE-
SC0024649. The authors acknowledge the University of Florida Research
Computing for providing computational resources that have contributed to the
research results reported in this publication.

## References

  * Hender et al. (2007) T. Hender, J. Wesley, J. Bialek, A. Bondeson, A. Boozer, R. Buttery, A. Garofalo, T. Goodman, R. Granetz, Y. Gribov, et al., Nuclear fusion 47, S128 (2007).
  * Matthews et al. (2016) G. Matthews, B. Bazylev, A. Baron-Wiechec, J. Coenen, K. Heinola, V. Kiptily, H. Maier, C. Reux, V. Riccardo, F. Rimini, et al., Physica Scripta 2016, 014070 (2016).
  * Hoppe et al. (2021) M. Hoppe, O. Embreus, and T. Fülöp, Computer Physics Communications 268, 108098 (2021).
  * Bandaru et al. (2021) V. Bandaru, M. Hoelzl, C. Reux, O. Ficker, S. Silburn, M. Lehnen, N. Eidietis, J. Contributors, J. Team, et al., Plasma Physics and Controlled Fusion 63, 035024 (2021).
  * Liu et al. (2021) C. Liu, C. Zhao, S. C. Jardin, N. M. Ferraro, C. Paz-Soldan, Y. Liu, and B. C. Lyons, Plasma Physics and Controlled Fusion 63, 125031 (2021).
  * Sainterme and Sovinec (2024) A. Sainterme and C. Sovinec, Physics of Plasmas 31 (2024).
  * Dreicer (1959) H. Dreicer, Physical Review 115, 238 (1959).
  * Kruskal and Bernstein (1962) M. Kruskal and I. Bernstein, Princeton Plasma Physics Lab. Report No. Matt-Q-20 (unpublished) (1962).
  * Connor and Hastie (1975) J. Connor and R. Hastie, Nuclear fusion 15, 415 (1975).
  * Hesslow et al. (2019a) L. Hesslow, L. Unnerfelt, O. Vallhagen, O. Embréus, M. Hoppe, G. Papp, and T. Fülöp, Journal of Plasma Physics 85, 475850601 (2019a).
  * Helander et al. (2004) P. Helander, H. Smith, T. Fülöp, and L.-G. Eriksson, Physics of Plasmas 11, 5704 (2004).
  * Smith and Verwichte (2008) H. M. Smith and E. Verwichte, Physics of plasmas 15, 072502 (2008).
  * Yang et al. (2023) M. Yang, P. Wang, D. del Castillo-Negrete, Y. Cao, and G. Zhang, arXiv preprint arXiv:2306.05580 (2023).
  * Rosenbluth and Putvinski (1997) M. Rosenbluth and S. Putvinski, Nuclear Fusion 37, 1355 (1997).
  * Martín-Solís et al. (2015) J. Martín-Solís, A. Loarte, and M. Lehnen, Physics of Plasmas 22, 092512 (2015).
  * Aleynikov and Breizman (2015) P. Aleynikov and B. N. Breizman, Physical review letters 114, 155001 (2015).
  * Hesslow et al. (2019b) L. Hesslow, O. Embréus, O. Vallhagen, and T. Fülöp, Nuclear Fusion 59, 084004 (2019b).
  * McDevitt (2023) C. J. McDevitt, Physics of Plasmas 30, 092501 (2023), ISSN 1070-664X, URL https://doi.org/10.1063/5.0164712.
  * SOKOLOV (1979) I. SOKOLOV, JETP Letters 29, 218 (1979).
  * Jayakumar et al. (1993) R. Jayakumar, H. Fleischmann, and S. Zweben, Physics Letters A 172, 447 (1993).
  * Martín-Solís et al. (2017) J. Martín-Solís, A. Loarte, and M. Lehnen, Nuclear Fusion 57, 066025 (2017).
  * Vallhagen et al. (2020) O. Vallhagen, O. Embreus, I. Pusztai, L. Hesslow, and T. Fülöp, Journal of Plasma Physics 86 (2020).
  * Lagaris et al. (1998) I. E. Lagaris, A. Likas, and D. I. Fotiadis, IEEE transactions on neural networks 9, 987 (1998).
  * Karpatne et al. (2017) A. Karpatne, G. Atluri, J. H. Faghmous, M. Steinbach, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, and V. Kumar, IEEE Transactions on knowledge and data engineering 29, 2318 (2017).
  * Karniadakis et al. (2021) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Nature Reviews Physics 3, 422 (2021).
  * Lusch et al. (2018) B. Lusch, J. N. Kutz, and S. L. Brunton, Nature communications 9, 1 (2018).
  * Wang et al. (2020) R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu, in _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ (2020), pp. 1457–1466.
  * Abadi et al. (2016) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., in _12th USENIX symposium on operating systems design and implementation (OSDI 16)_ (2016), pp. 265–283.
  * Paszke et al. (2017) A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer (2017).
  * Sun et al. (2020) L. Sun, H. Gao, S. Pan, and J.-X. Wang, Computer Methods in Applied Mechanics and Engineering 361, 112732 (2020).
  * McDevitt and Tang (2024) C. J. McDevitt and X.-Z. Tang, arXiv preprint arXiv:2402.08495 (2024).
  * Wang et al. (2022) S. Wang, X. Yu, and P. Perdikaris, Journal of Computational Physics 449, 110768 (2022).
  * Lu et al. (2021) L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, SIAM Review 63, 208 (2021).
  * Antonsen Jr and Chu (1982) T. Antonsen Jr and K. Chu, The Physics of Fluids 25, 1295 (1982).
  * Taguchi (1983) M. Taguchi, Journal of the Physical Society of Japan 52, 2035 (1983).
  * Fisch (1986) N. Fisch, The Physics of fluids 29, 172 (1986).
  * Karney and Fisch (1986) C. F. Karney and N. J. Fisch, The Physics of fluids 29, 180 (1986).
  * Liu et al. (2016) C. Liu, D. P. Brennan, A. Bhattacharjee, and A. H. Boozer, Physics of Plasmas 23, 010702 (2016).
  * Liu et al. (2017) C. Liu, D. P. Brennan, A. H. Boozer, and A. Bhattacharjee, Plasma Physics and Controlled Fusion 59, 024003 (2017), URL http://stacks.iop.org/0741-3335/59/i=2/a=024003.
  * Zhang and del Castillo-Negrete (2017) G. Zhang and D. del Castillo-Negrete, Physics of Plasmas 24, 092511 (2017).
  * Andersson et al. (2001) F. Andersson, P. Helander, and L.-G. Eriksson, Physics of Plasmas 8, 5221 (2001).
  * Decker et al. (2016) J. Decker, E. Hirvijoki, O. Embreus, Y. Peysson, A. Stahl, I. Pusztai, and T. Fülöp, Plasma Physics and Controlled Fusion 58, 025016 (2016).
  * Guo et al. (2017) Z. Guo, C. J. McDevitt, and X.-Z. Tang, Plasma Physics and Controlled Fusion 59, 044003 (2017).
  * McDevitt et al. (2018) C. J. McDevitt, Z. Guo, and X.-Z. Tang, Plasma Physics and Controlled Fusion 60, 024004 (2018).
  * Møller (1932) C. Møller, Annalen der Physik 406, 531 (1932).
  * Ashkin et al. (1954) A. Ashkin, L. A. Page, and W. Woodward, Physical Review 94, 357 (1954).
  * Boozer (2015) A. H. Boozer, Physics of Plasmas 22, 032504 (2015).
  * McDevitt et al. (2019) C. J. McDevitt, Z. Guo, and X. Tang, Plasma Physics and Controlled Fusion 61, 054008 (2019).
  * Boozer (2018) A. H. Boozer, Nuclear Fusion 58, 036006 (2018).
  * McDevitt and Tang (2019) C. McDevitt and X.-Z. Tang, EPL (Europhysics Letters) 127, 45001 (2019).
  * Arnaud and McDevitt (2024) J. Arnaud and C. McDevitt, arXiv preprint arXiv:2401.11291 (2024).
  * Guilford and Guilford (1936) J. P. Guilford and R. B. Guilford, The journal of Psychology 2, 109 (1936).
  * Hesslow et al. (2017) L. Hesslow, O. Embréus, A. Stahl, T. C. DuBois, G. Papp, S. L. Newton, and T. Fülöp, Phys. Rev. Lett. 118, 255001 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.255001.


# A reduced model for particle calcination for use in DEM/CFD simulations

Lucas Reineking  Automatic Control and System Theory, Ruhr University Bochum,
Germany  Torben Bergold  Energy Plant Technology, Ruhr University Bochum,
Germany  Enric Illana  Energy Plant Technology, Ruhr University Bochum,
Germany  Viktor Scherer  Energy Plant Technology, Ruhr University Bochum,
Germany  Martin Mönnigmann  Automatic Control and System Theory, Ruhr
University Bochum, Germany

###### Abstract

We treat the accurate simulation of the calcination reaction in particles,
where the particles are large and, thus, the inner-particle processes must be
resolved. Because these processes need to be described with coupled partial
differential equations that must be solved numerically, the computation times
for a single particle are too high for use in simulations that involve many
particles. Simulations of this type arise when the Discrete Element Method
(DEM) is combined with Computational Fluid Dynamics (CFD) to investigate
industrial systems such as quick lime production in lime shaft kilns.

We show that, based on proper orthogonal decomposition and Galerkin
projection, reduced models can be derived for single particles that provide
the same spatial and temporal resolution as the original PDE models at a
considerably reduced computational cost. Replacing the finite-volume particle
models with the reduced models results in an overall reduction of the reactor
simulation time by about 60% for the simple example treated here.

##  1 Introduction

Quick lime (CaO) is an important industrial commodity used for many
applications, such as wastewater treatment, flue gas desulfurization or steel
production. Quick lime is produced from lime stone ($\text{CaCO}_{3}$) by
calcination. The reactors commonly employed for quick lime production are
shaft kilns. Such kilns can reach heights of $20\,\mathrm{m}$, diameters of
$4\,\mathrm{m}$ and production capacities of $400\,\mathrm{t/d}$. Due to their
size, measurements are extremely challenging in these reactors because the
moving particle bed is densely packed and high temperatures above
$850\,^{\circ}\mathrm{C}$ are present. Therefore, numerical simulation of
these kilns is an attractive alternative.

Models and numerical tools used for lime shaft kiln simulation vary in
complexity. Many studies are based on one-dimensional heat and mass balances.
For example, Bes [5] investigated the operational conditions inside a kiln for
different fuels such as lean gas, natural gas and lignite neglecting gradients
in the radial direction of the kiln. Shagapov and Burkin [22] presented a
mathematical model to compute a one-dimensional profile of gas-phase
properties and solid temperatures along a kiln considering calcination of
limestone and burning of coke.

Other approaches use CFD simulations to resolve the fluid flow. The particles
are, however, typically treated in an approximate fashion. For example,
Mohammadpour et al. [19] carried out CFD simulations of a shaft kiln
representing the presence of the particles by a porous medium in the CFD
simulations. Similarly, Duan et al. [7] employed CFD simulations combined with
a shrinking core calcination model on a sub-grid level to investigate the
inner state of a shaft kilns [7, 8].

For a more detailed description of the three-dimensional processes in lime
shaft kilns, a method must be used that reflects the correct time and
temperature evolution of the limestone bed in the kiln. Such a method is the
Discrete Element Method (DEM). DEM tracks the movement of the individual
particles and their reaction. When combined with Computational Fluid Dynamics
(CFD), this allows for a detailed description of the thermochemical processes
in the kiln. Studies that combine DEM and CFD in this manner are scarce,
however. A single shaft kiln was investigated with DEM/CFD in [6] and [15].
The same method was applied to a two-shaft regenerative kiln in a subsequent
study [16]. The simulation of a shaft kiln in oxyfuel operation was treated
[13], where combustion takes place in a $\text{O}_{2}-\text{CO}_{2}$
atmosphere to allow for $\text{CO}_{2}$ separation at the kiln gas outlet.

DEM/CFD simulations for lime shaft kilns are very time-consuming as the number
of particles can easily reach more than one million, and particle and gas-
phase reaction must be described in addition to particle and fluid flow. The
present paper demonstrates that the models and simulations for the intra-
particle calcination processes can be replaced by reduced order models. In
contrast to the original PDE models, the reduced models consist of ordinary
differential equations (ODEs) that require only a fraction of the
computational time, while the error of the reduction can be controlled.
Several different model reduction methods have been explored for their
potential to reduce the overall computation time of DEM/CFD with complex
particle models, among them methods based on proper orthogonal decomposition
and projection both for particles and flow fields (see, e.g., [24]), methods
based purely on machine learning (see, e.g., [14, 18, 17]), and combinations
of these two classes of methods (see, e.g., [18]).

We use model reduction based on proper orthogonal decomposition and Galerkin
projection (see, e.g., [2]). This choice is motivated by their successful
application to DEM/CFD simulations of drying processes. These investigations
started with single particles [21] and proving they can be used for optimal
control [3] and state reconstruction purposes [4]. More recently, it has been
shown that DEM/CFD simulations for wood chip drying can be accelerated by a
factor of about three if the original single particle PDE models are replaced
with reduced models [20].

In contrast to drying, which can be modeled with coupled diffusion processes,
calcination involves reactions and, hence, exponentially depends on
temperature as long as transport does not limit the conversion. More
specifically, the conversion of limestone to quick lime inside the particles
is the challenge. The local conversion rate is determined by the local partial
pressure of $\text{CO}_{2}$, which is in turn governed by the local production
of $\text{CO}_{2}$ of the calcination and the transport of $\text{CO}_{2}$ in
the gas phase. The gas-phase transport in limestone is considerably different
from that in quick lime because quick lime has got a distinctly higher
porosity, which leads to a lower transport resistance. Effectively, a reaction
front is formed that moves from the particle surface to the particle core.

The intra-particle model for calcination is introduced in section 2.
Subsequently, the derivation of the reduced model is explained in section 3.
The sample reactor system, and the methods used to simulate it, are described
in section 4. Reference simulation results obtained with the finite-volume
method and the reduced model are presented and compared in section 5.

##  2 Particle model

The particle is modeled as a homogeneous isotropic porous medium. The solid
phase is composed of limestone ($\text{CaCO}_{3}$) which reacts to quicklime
(CaO) according to the endothermic calcination reaction

$\text{CaCO}_{3}(s)\,\rightarrow\,\text{CaO}(s)\,+\,\text{CO}_{2}(g)\,+\Delta^{R}H.$
(1)

The reaction progress $R$ is defined as the mass ratio between produced mass
of quicklime $m_{\text{CaO}}$ and the maximum mass of quicklime
$m_{\text{CaO},\max}$ that can be produced, i.e.,

$R=\frac{m_{\text{CaO}}}{m_{\text{CaO},\max}}=\frac{m_{\text{CaO}}}{m_{\text{CaCO}_{3},0}}\frac{M_{\text{CaCO}_{3}}}{M_{\text{CaO}}}$
(2)

where $m_{\text{CaCO}_{3},0}$ denotes the initially present mass of limestone,
and $M_{\text{CaCO}_{3}}$ and $M_{\text{CaO}}$ are the molar masses of
$\text{CaCO}_{3}$ and CaO, respectively.

The gas phase, a mixture of air and $\text{CO}_{2}$ from (2) with density
$\varrho$ and $\text{CO}_{2}$ mass fraction $y$, is transported through the
pores of the particle. This is modeled with the porosity approach, i.e., the
intra-particle pores are not resolved and all quantities are volume-averaged.
The conservation equations for mass and species are stated in terms of density
$\rho$ and mass fraction of carbon dioxide $y$. They read

$\displaystyle\frac{\partial}{\partial
t}(\varepsilon\varrho)+\nabla\cdot(\varepsilon\varrho\mathbf{v})$
$\displaystyle=M_{\text{CO}_{2}}\dot{n}$ (3a)
$\displaystyle\frac{\partial}{\partial t}(\varepsilon\varrho
y)+\nabla\cdot(\varepsilon\varrho y\mathbf{v})$
$\displaystyle=\nabla\cdot(\varepsilon\varrho D_{\mathrm{eff}}\cdot\nabla
y)+M_{\text{CO}_{2}}\dot{n}$ (3b)

with three-dimensional gas-phase velocity $\mathbf{v}$, molar reaction rate
$\dot{n}$ of (1), and molar mass $M_{\text{CO}_{2}}$ of $\text{CO}_{2}$. The
effective diffusivity $D_{\mathrm{eff}}$ takes the binary diffusion of
$\text{CO}_{2}$ in air and Knudsen diffusion into account. Details are given
in (22) - (24) in appendix A. The intra-particle porosity $\varepsilon$,
defined as the fraction of the volume occupied by the gas phase to the volume
occupied by both phases, depends on the reaction progress according to

$\varepsilon=R\varepsilon_{\text{CaO}}+(1-R)\varepsilon_{\text{CaCO}_{3}},$
(4)

where $\varepsilon_{\text{CaO}}$ and $\varepsilon_{\text{CaCO}_{3}}$ denote
the porosity of quicklime and limestone, respectively. Note that the porosity
of quicklime is much larger than the porosity of limestone (see Table 1).

The velocity vector $\mathbf{v}$ determines the advection of the gas phase
through the porous particle. We calculate it with Darcy’s law

$\mathbf{v}=-\frac{K}{\mu}\nabla p,$ (5)

where the permeability $K$ is a material constant and $\mu$ denotes the fluid
dynamic viscosity. The gas-phase pressure $p$ is calculated according to the
ideal gas law

$\varrho R_{m}T=Mp,$

where $M$ denotes the molar mass of the gas phase.

The reaction rate $\dot{n}$ depends on the local temperature $T$ and
concentration of $\text{CO}_{2}$ $c_{\text{CO}_{2}}$ according to the
Arrhenius equation

$\dot{n}=k_{0}A_{\mathrm{eff}}T\,\exp\bigg{(}-\frac{E_{A}}{R_{m}T}\bigg{)}Y_{\mathrm{TC}}\,(c_{\text{CO}_{2}}-c_{\mathrm{eq}}(T))$

where $k_{0}$, $A_{\mathrm{eff}}$, $E_{A}$, $R_{m}$ and $Y_{\mathrm{TC}}$
denote the constant preexponential factor (in m/s), the effective reaction
surface (according to the size of the control volume), the activation energy,
the universal gas constant and the temperature correction factor,
respectively. The equilibrium concentration $c_{\mathrm{eq}}(T)$ is calculated
according to the ideal gas law

$c_{\mathrm{eq}}=\frac{p_{\mathrm{eq}}(T)}{R_{m}T},$

where the equilibrium pressure $p_{\mathrm{eq}}(T)$ can be determined with the
empirical correlation (20) given in appendix A.

The balance equations (3) need to be combined with an energy balance, which is
stated in terms of the temperature $T$. We recall the temperature must be
resolved inside the particle, because the particles are large. The temperature
of the gas and solid phase are assumed to be equal. This approximation is
valid because due to the large surface of the interface between the two
phases. The PDE for the conservation of energy for the combined solid-gas
phases reads

$\frac{\partial}{\partial t}(\rho c_{p}T)=\nabla\cdot(\lambda\cdot\nabla
T)+M_{\text{CO}_{2}}\,\Delta^{R}H\,\dot{n},$ (6)

where the heat capacity $c_{p}$ and thermal conductivity $\lambda$ are mass-
averaged from the individual phases. They depend on the temperature and
composition of the gas and the solid phases. As the particle diameter stays
constant during calcination, the density $\rho$ can be computed from the sum
of masses of the phases. The loss of solid mass in the reaction (1) must be
compensated by an increase of porosity.

Fluid and particles exchange heat and mass flux across the particle surface.
Hence, the system of governing equations of (3) and (6) is completed by the
boundary conditions

$\displaystyle\rho c_{p}\frac{\partial
T}{\partial\mathbf{n}}|_{\partial\Omega}$
$\displaystyle=\alpha(T_{\infty}-T)|_{\partial\Omega}$ (7a)
$\displaystyle\varepsilon\varrho D_{\mathrm{eff}}\frac{\partial
y}{\partial\mathbf{n}}|_{\partial\Omega}$
$\displaystyle=\beta\varrho_{\infty}(y_{\infty}-y)|_{\partial\Omega}$ (7b)
$\displaystyle\varepsilon\varrho\frac{K}{\mu}\frac{\partial
P}{\partial\mathbf{n}}|_{\partial\Omega}$
$\displaystyle=C^{*}(p_{\infty}-p)|_{\partial\Omega}$ (7c)

for the convective heat flux (7a), the convective transport of $\text{CO}_{2}$
(7b), and the advective mass flux (7c), where $T_{\infty}$, $\rho_{\infty}$
and $p_{\infty}$ are the properties of the fluid from the CFD simulation. The
heat transfer coefficient, the convective mass transfer coefficient and the
advective transfer coefficient are denoted by $\alpha$, $\beta$ and $C^{*}$,
respectively. All transfer coefficients are computed from the empirical
correlations (25) - (28) given in appendix A.

The partial differential equations (3) and (6) are discretized with a 3D
finite volume scheme of tetrahedral cells. The validation of the model and its
sensitivity to the spatial discretization are discussed in appendix B.

##  3 Reduced particle model

The reduced model essentially replaces the PDEs (3), (6) and the finite volume
solver required to solve them with a set of spatial modes $\varphi^{(\cdot)}$
and time-dependent weighting factors $a^{(\cdot)}(t)$ for these modes. The
time-dependent linear combination of the modes $\varphi^{(\cdot)}$ then yields
the temperature $T(x,t)$, gas phase density $\rho(x,t)$ and $\text{CO}_{2}$
mass fraction $y(x,t)$ in the particle as a function of location $x$ and time
$t$ according to

$\displaystyle T(x,t)$
$\displaystyle\approx\bar{T}(x)+\kappa_{T}\sum_{i=1}^{r}\varphi^{(T)}_{i}(x)a^{(T)}_{i}(t)$
(8a) $\displaystyle\varrho(x,t)$
$\displaystyle\approx\bar{\varrho}(x)+\kappa_{\varrho}\sum_{i=1}^{r}\varphi^{(\varrho)}_{i}(x)a^{(\varrho)}_{i}(t)$
(8b) $\displaystyle y(x,t)$
$\displaystyle\approx\bar{y}(x)+\kappa_{y}\sum_{i=1}^{r}\varphi^{(y)}_{i}(x)a^{(y)}_{i}(t)$
(8c)

where $\bar{T}(x)$, $\bar{\rho}(x)$ and $\bar{y}(x)$ are the temporal means
and $\kappa_{T}$, $\kappa_{\rho}$ and $\kappa_{y}$ are scale factors
introduced for numerical reasons. The description of the calcination process
with (8) is much more computationally efficient than solving the PDE model
(3), (6) because the modes $\varphi^{(\cdot)}$ only have to be determined once
as a preparation, and the coefficients $a^{(\cdot)}$ are governed by ODEs that
replace the PDEs (3), (6). Moreover, a small number of ODEs, coefficients
$a_{i}^{(\cdot)}$ and modes $\varphi^{(\cdot)}$ suffices to achieve a high
precision. This number, which is denoted by $r$ in (8), amounts to $r=5$ for
the calcination particle model treated here (see section 5).

We briefly describe how to determine the modes $\varphi^{(\cdot)}_{i}$ and the
ODEs that govern the coefficients $a_{i}^{(\cdot)}$ in sections 3.1 and 3.2.
These sections are given in order for the paper to be self-contained.

###  3.1 Determining the modes

The modes $\varphi^{(\cdot)}_{i}(x)$, as well as $T(x,t)$, $\rho(x,t)$ and
$y(x,t)$, are denoted as continuous functions of $x$ and $t$ for convenience
in (8). These quantities are approximated on the same grid used in the finite
volume solver. Let $\Omega$, $N$ and $x_{j}$ refer to the spatial domain the
particle occupies, the number of cells used to discretize $\Omega$, and the
cell center of cell $j$, respectively. Assume the finite volume method
provides results at the $M$ time points $t_{k}=k\Delta t$, $k=0,\dots,M-1$.
The matrix that holds the results provided by the finite volume solver is
called snapshot matrix. For $T(x,t)$ this matrix reads

$\tilde{T}=\begin{bmatrix}|&\,&|\\\ T(t_{0},x_{j})&\cdots&T(t_{M-1},x_{j})\\\
|&\,&|\end{bmatrix}\in\mathbb{R}^{N\times M}.$

The temporal mean $\bar{T}$ required in (8) is the row average of $\tilde{T}$.
The scale factor $\kappa_{T}$ used in (8) is chosen to be

$\kappa_{T}=\max_{x_{j},\,t_{k}}|\tilde{T}-\bar{T}|.$

The snapshot matrices, temporal means and scale factors for density $\varrho$
and mass fraction $y$ are calculated analogously. We calculate the modes
$\varphi_{i}^{(T)}(x)$

$\begin{bmatrix}|\\\ \varphi_{i}^{(T)}(x_{j})\\\
|\end{bmatrix}\in\mathbb{R}^{N},\,i=1,\dots,r$

from the thin singular value decomposition of the scaled meanfree snapshot
matrix

$\displaystyle\frac{1}{\kappa_{T}}(\tilde{T}-\bar{T})$
$\displaystyle=\Phi_{T}S_{T}W_{T}^{\intercal}$
$\displaystyle\varphi_{i}^{(T)}$
$\displaystyle=\left(\Phi_{T}\right)_{\cdot,i}.$

where $\Phi_{T}\in\mathbb{R}^{N\times M}$, $S_{T}\in\mathbb{R}^{M\times M}$,
and $W_{T}\in\mathbb{R}^{M\times M}$, which implies
$\mathrm{rank}(\tilde{T}-\bar{T})=M$.

The first $r$ columns of $\Phi_{T}$ are known to be the optimal truncated
basis with $r$ elements for the column space of
$(\tilde{T}-\bar{T})/\kappa_{T}$ according to the Eckart-Young theorem [see
Algorithm 8.6.2, 10]. Consequently, these $r$ columns are the optimal choice
for the approximation (8).

In general, the spatial domain $\Omega$ is discretized with finite volume
cells of varying volume $V_{j}$. Therefore, a volume-weighted scalar product
is introduced for functions $f,g$ on the domain $\Omega$

$\langle
f(x),g(x)\rangle=\int_{\Omega}f(x)\,g(x)\,\mathrm{d}v\approx\sum_{j=1}^{N}f(x_{j})g(x_{j})V_{j}$
(9)

The modes $\varphi^{(\cdot)}(x_{j})$, which are defined for each finite volume
cell center $x_{j}$, can be understood as the spatial discretization of
$\varphi^{(\cdot)}(x)$ with a continuous $x\in\Omega$. The approximation (9)
therefore reads
$\langle\varphi^{(\cdot)}_{i},\varphi^{(\cdot)}_{j}\rangle=\langle\varphi^{(\cdot)}_{i}(x),\varphi^{(\cdot)}_{j}(x)\rangle=$

$\int_{\Omega}\varphi^{(\cdot)}_{i}(x),\varphi^{(\cdot)}_{j}(x)dv\approx\sum_{k=1}^{N}\varphi^{(\cdot)}_{i}(x_{k})\varphi^{(\cdot)}_{j}(x_{k})V_{k}$
(10)

for the modes. We use the left hand side of (10) whenever convenient. The
scalar product (9) is used to scale the modes $\varphi^{(\cdot)}_{i}$
according to

$\displaystyle\langle\varphi^{(\cdot)}_{i},\varphi^{(\cdot)}_{j}\rangle$
$\displaystyle=\frac{1}{3}\delta_{ij}$

where $\delta_{ij}$ is Kronecker’s delta.

###  3.2 Galerkin Projection

It remains to derive ODEs that govern the coefficients $a_{i}^{(\cdot)}$
required for the approximation (8). These ODEs are derived by substituting (8)
into the PDEs (3) and (6) and calculating the scalar products (9) of the
resulting equations with the modes $\varphi_{i}^{(\cdot)}$. The variant of
this procedure used here is called Galerkin projection. We demonstrate the
Galerkin projection with the mass fraction $y$ because the dependence of the
projected ODEs appears for this case.

Forming the scalar product (9) of the first $r$ modes $\varphi_{i}^{(y)}$ with
(3b) yields

$\displaystyle\langle\varphi^{(y)}_{i}(x),\frac{\partial}{\partial{t}}$
$\displaystyle(\varepsilon\varrho
y)\rangle=\>\langle\varphi^{(y)}_{i}(x),\nabla\cdot(D_{\mathrm{eff}}\varrho\nabla{y}-\varepsilon\varrho
y\mathbf{v})\rangle$ (11) $\displaystyle+$
$\displaystyle\langle\varphi^{(y)}_{i}(x),M_{\text{CO}_{2}}\dot{n}\rangle,\qquad
i=1,\ldots,r.$

Note that $r$ linearly independent equations are created. It is convenient to
treat the scalar products in (11) separately. Applying the chain rule to the
left hand side of (11) yields

$\langle\varphi^{(y)}_{i}(x),\frac{\partial}{\partial{t}}(\varepsilon\varrho
y)\rangle=\sum_{j=1}^{r}M^{(yy)}_{ij}\>\dot{a}^{(y)}_{j}(t)+M^{(y\varrho)}_{ij}\>\dot{a}^{(\varrho)}_{j}(t)$
(12a) with $\displaystyle M^{(yy)}_{ij}$
$\displaystyle=\kappa_{y}\int_{\Omega}\varphi^{(y)}_{i}(x)\,\varepsilon\varrho\varphi^{(y)}_{j}(x)\,\mathrm{d}v$
(12b) $\displaystyle M^{(y\varrho)}_{ij}$
$\displaystyle=\kappa_{\varrho}\int_{\Omega}\varphi^{(y)}_{i}(x)\,\varepsilon
y\varphi^{(\varrho)}_{j}(x)\,\mathrm{d}v,$ (12c)

where $M^{(yy)}_{ij},M^{(y\varrho)}_{ij}\in\mathbb{R}$. Using Gauss theorem
and substituting the approximations (8), the first scalar products on the
right handside of (11) results in

$\langle\varphi^{(y)}_{i}(x),\nabla\cdot(D_{\mathrm{eff}}\varrho\nabla{y}-\varepsilon\varrho
y\mathbf{v})\rangle=S^{(y)}_{i}-I^{(y)}_{i}$

where

$\displaystyle S^{(y)}_{i}$
$\displaystyle=\int_{\partial\Omega}\varphi^{(y)}_{i}(x)\>\big{(}D_{\mathrm{eff}}\varrho\nabla{y}-\varepsilon\varrho
y\mathbf{v}\big{)}\cdot\mathbf{n}\,\mathrm{d}a$ (13a) $\displaystyle
I^{(y)}_{i}$
$\displaystyle=\int_{\Omega}\nabla\varphi^{(y)}_{i}(x)\cdot\big{(}D_{\mathrm{eff}}\varrho\nabla{y}-\varepsilon\varrho
y\mathbf{v}\big{)}\,\mathrm{d}v$ (13b)

can be interpreted as the surface fluxes $S^{(y)}\in\mathbb{R}^{R}$ and the
innner fluxes $I^{(y)}\in\mathbb{R}$, respectively. This interpretation is
advantageous for two reasons. First, due to the surface integral in (13a) it
is possible to explicitly incorporate the boundary conditions (7b) and (7c)
into the reduced model. This yields the surface fluxes

$S^{(y)}_{i}=\int_{\partial\Omega}\varphi^{(y)}_{i}(x)\>\big{(}\beta\varrho_{\infty}(y_{\infty}-y)+yC^{*}(p_{\infty}-p)\big{)}\,\mathrm{d}a$
(14)

Second, only first derivatives appear in (13b). The remaining scalar product
in (11) reads

$\langle\varphi^{(y)}_{i}(x),M_{\text{CO}_{2}}\dot{n}\rangle=\int_{\Omega}\varphi^{(y)}_{i}(x)M_{\text{CO}_{2}}\dot{n}\,\mathrm{d}v=J^{(y)}_{i},$
(15)

which defines the $i$th element of the inner source terms
$J^{(y)}\in\mathbb{R}^{r}$.

Substituting (12a)-(15) into (11) yields the ODE for the mass fraction $y$

$M^{(yy)}\dot{a}^{(y)}(t)+M^{(y\varrho)}\dot{a}^{(\varrho)}(t)=S^{(y)}-I^{(y)}+J^{(y)}.$
(16)

Note that both $\dot{a}^{(y)}$ and $\dot{a}^{(\rho)}$ appear in (16). In this
sense, the ODEs for $a^{(y)}$ and $a^{(\rho)}$ are coupled.

The procedure that yields (16) is analogously applied to the PDEs (3b) and
(6). Combining the results for the Galerkin projections of all three PDEs
yields

$\mathbf{M}\mathbf{\dot{a}}=\mathbf{S}-\mathbf{I}+\mathbf{J}$ (17)

where

$\mathbf{M}=\begin{bmatrix}M^{(TT)}&0&0\\\ 0&M^{(\varrho\varrho)}&0\\\
0&M^{(y\varrho)}&M^{(yy)}\end{bmatrix}\in\mathbb{R}^{3r\times 3r}$ (18)

and

$\mathbf{\dot{a}}=\begin{bmatrix}\dot{a}^{(T)}\\\ \dot{a}^{(\varrho)}\\\
\dot{a}^{(y)}\end{bmatrix},\>\mathbf{S}=\begin{bmatrix}S^{(T)}\\\
S^{(\varrho)}\\\ S^{(y)}\end{bmatrix},\>\mathbf{I}=\begin{bmatrix}I^{(T)}\\\
I^{(\varrho)}\\\ I^{(y)}\end{bmatrix},\>\mathbf{J}=\begin{bmatrix}J^{(T)}\\\
J^{(\varrho)}\\\ J^{(y)}\end{bmatrix}\in\mathbb{R}^{3r}.$ (19)

The results of the reduced model for a single particle with constant boundary
conditions are given in Figure 10 in appendix C.

##  4 Sample system

We consider a rectangular reactor with a simple arrangement of particles as a
test case. The domain of interest measures $50\,\mathrm{mm}$ in width,
$50\,\mathrm{mm}$ in depth, and $125\,\mathrm{mm}$ in height. Five cylindrical
particles with a diameter of $12\,\mathrm{mm}$ and a length of
$36\,\mathrm{mm}$ are arranged vertically with $6\,\mathrm{mm}$ separation
distance between them ($18\,\mathrm{mm}$ axis to axis, see Figure 1). The
cylinders are discretized by 8000 tetrahedral control volumes each. The
particles are assumed to initially consist of pure limestone with an uniform
initial temperature of 900 ${}^{\circ}\text{C}$. The material properties for
limestone and quick lime are given in Table 1.

Table 1: Properties for lime and quicklime according to [12] property |  | unit | $\text{CaCO}_{3}$ | CaO
---|---|---|---|---
density | $\varrho_{\text{j}}$ | $\left[\frac{\text{kg}}{\text{m}^{3}}\right]$ | $2812.05$ | $3310$
heat capacity | $c_{\text{p,j}}$ | $\left[\frac{\text{J}}{\text{kgK}}\right]$ | $836.8$ | $753.1$
thermal conductivity | $\lambda_{\text{j}}$ | $\left[\frac{\text{W}}{\text{mK}}\right]$ | $2.26$ | $0.07$
porosity | $\epsilon_{\text{j}}$ | $[-]$ | $0.04$ | $0.543$
specific surface | $S_{\text{j}}$ | $\left[\frac{\text{m}^{2}}{\text{kg}}\right]$ | $16\,000$ | $7000$

Air enters the reactor from the bottom ($z=0$) with a superficial velocity
$v_{\text{in}}=0.5\,\mathrm{m/s}$, mass flow
$\dot{m}_{\text{in}}=0.375\,\mathrm{g/s}$ and temperature
$T_{\text{in}}=900\,\mathrm{{}^{\circ}C}$, and leaves the reactor at the top
end ($z=0.125\,\mathrm{m}$) where the ambient pressure is set to
$p_{\text{out}}=101325\,\mathrm{Pa}$. The particle Reynolds number amounts to
$\mathrm{Re}_{0}=228$ based on an equivalent sphere.

Figure 1: System consisting of 5 limestone particles distributed uniformly
over the height

DEM and CFD simulations are coupled with the porous medium approach, also
called averaged volume method (AVM). In this approach, particles are
represented in the flow field by a porosity field, which implies the detailed
particle shape is not considered in the CFD simulations. The porous media and
the CFD are coupled with source terms for momentum, energy and species from
DEM to CFD. Conversely, the CFD solution provides far-field boundary
conditions, which are uniform along the particle surface, to the porous medium
simulation. Since the purpose of the paper is to investigate the acceleration
that can be achieved with reduced models for the particles, we used a one-
dimensional model for the gas phase for simplicity. The domain is spatially
discretized with 7 control volumes in the $z$-direction ($\delta
z=18\,\mathrm{mm}$). Each particle is coupled only with the control volume in
which its center of gravity is located. A zero gradient condition is imposed
at the side boundaries for all quantities. All simulations are carried out
with an inhouse DEM code. OpenFOAM version 20.12 [1] is used for simulating
mass, species and energy transport in the fluid phase.111 All simulations were
conducted with a single core AMD(R) EPYC(R) 7443 CPU @ 2.85GHz and 256 GB RAM
under Ubuntu 20.04.5 LTS

Because of the small number of particles, an immersed boundary method (IBM)
[see, e.g., 11] would also be applicable to accurately describe particle shape
in the flow field. However, we deliberately have refrained from IBM as it will
be too costly for large-scale simulations of technical reactors with several
hundred thousands of particles, which are the motivation of our work.

##  5 Results

Figure 2 depicts the temporal evolution of the calcination degree,
temperature, $\text{CO}_{2}$ mass fraction and gas-phase density. All diagrams
in Figure 2 show the spatial average over the five particles. As the particles
have an initial uniform temperature of 900 K, the $\text{CO}_{2}$ production
driven by kinetics is initially very high in the whole particle, leading to a
peak in the $\text{CO}_{2}$ mass fraction. In parallel, the particle
temperature drops fast due to the endothermic calcination reaction initially
present in the whole particle. The fast $\text{CO}_{2}$ release at the
beginning also leads to an initial peak in density, as the $\text{CO}_{2}$
produced is initially stored in the particle and not enough $\text{CO}_{2}$
can leave the particle, i.e., an overpressure builds up in the particle. The
$\text{CO}_{2}$ mass fraction decreases very fast after the initial peak since
the advective mass flow, which is proportional to the intra-particle pressure,
is large leading to a high driving force. As soon as this initial phase has
settled to a minimum temperature and an equilibrated pressure at around 250 s,
the calcination progresses steadily inwards. This is accompanied by an
increase in particle temperature as heat is transferred by convection from the
gas to the particle and, in addition, the outer control volumes become fully
calcinated, i.e., the heat sink is then absent. From 250 seconds until
approximately 1000 seconds, the $\text{CO}_{2}$ mass fraction and gas density
increase again with moderate gradients as the particle temperature is still at
a low level, i.e., the reaction rate is slow. Nevertheless, $\text{CO}_{2}$
production is larger than the amount of $\text{CO}_{2}$ being diffused and
advected outside the particles, leading again to a storage effect of
$\text{CO}_{2}$ in the particle. After 1000 seconds, due to the reduction in
the reaction area of the reaction front, diffusion and advection become
predominant, and the $\text{CO}_{2}$ mass fraction as well as density decrease
to the values of the surrounding gas. In general, the reduced-order model with
only 5 modes and 5 ODEs is able to reproduce the results obtained with the
finite DEM/CFD simulation in a satisfactory manner. This holds for the
calcination degree in particular, which is the main parameter defining quick
lime quality.

Figure 2: Spatially averaged (over all particles) calcination degree,
temperature, $\text{CO}_{2}$ mass fraction and gas-phase density inside the
particle over time.

Figures 3–6 show vertical cross-sectional views through the center plane of
the reactor. Black vertical lines divide the diagrams into left and right
parts that show the results of the DEM/CFD simulation and the results obtained
with the reduced-order model, respectively.

The gas phase temperature is depicted in Figure 3. It is obvious that the gas-
phase temperature decreases from the gas inlet to the outlet. This is due to
the fact that the gas phase transfers heat to the particles along its path to
the outlet to cover the endothermic calcination reaction. For each of the
particles, the temperature in the center is lower than at the surface. In the
outer layers, calcination progress is already 1 (see Figure 4) and the heat
transferred from the gas phase is no longer required to cover the heat for
calcination but available for particle heating. The mean temperature of the
particle further downstream is at a lower level, as the progress of
calcination is lower there, i.e., more heat is still needed for the
endothermic calcination reaction. The reduced-order model results in slightly
higher temperatures in the interiors of the particles. This deviation, which
amounts to a few Kelvin, is consistent with the temperature values in Figure 2
at time $t=2500\text{s}$.

Figure 3: Temperature in the center cross section of the reaction volume at
time $t=2500\text{s}$. DEM/CFD results shown to the left and reduced-order
model results shown to the right of the black vertical line.

Figure 4 shows the reaction progress as defined in (2) for the particles. It
resembles the intra-particle $\text{CO}_{2}$ mass fraction (see Figure 6),
however, is not identical to it. Whereas the reaction progress describes local
conversion of $\text{CaCO}_{3}$ to CaO, the $\text{CO}_{2}$ mass fraction only
gives an information of the local $\text{CO}_{2}$ concentration in the
particle. This implies, for example, a cell can have a non-zero
$\text{CO}_{2}$ mass fraction even if no reaction takes place in it, because
of diffusion of $\text{CO}_{2}$ from an neighboring cell. The results obtained
with the DEM/CFD simulation and the reduced-order model are in good agreement
in Figure 4.

Figure 4: Reaction progress in the center cross section of the reaction volume
at time $t=2500\text{s}$. DEM/CFD results shown to the left and reduced-order
model results shown to the right of the black vertical line.

The gas phase density is depicted in Figure 5.

Figure 5: Gas phase density at the center cross section of the reaction volume
at time $t=2500\text{s}$. DEM/CFD results shown to the left and reduced-order
model results shown to the right of the black vertical line.

The density is determined by the balance of $\text{CO}_{2}$ released from the
limestone and the transport of the $\text{CO}_{2}$ away from the location it
is produced at. When the transport is slower than $\text{CO}_{2}$ production,
the local density (and the pressure) increases. This is reflected in Figure 5.
In the center of the particle, where the calcination reaction is active, the
pressure is larger than in outer layers, where reaction progress is already 1,
i.e., no $\text{CO}_{2}$ is produced any longer and transport to the particle
surface takes place. It is evident from Figure 5 that the results obtained
with the reduced-order model are in very good agreement for the
$\text{CO}_{2}$ concentration.

Figure 6 shows that the $\text{CO}_{2}$ released from the limestone is
accumulated in the gas phase, i.e., the gas phase $\text{CO}_{2}$
concentration increases in the direction of the fluid flow. The
$\text{CO}_{2}$ mass fraction in the particle is lower at the reactor entrance
than further downstream for two reasons. First, the accumulated
$\text{CO}_{2}$ in the fluid further downstream hinders the mass transfer of
$\text{CO}_{2}$ from the particle core to the particle surface, i.e., more
$\text{CO}_{2}$ remains in the particles that are located closer to the gas
outlet of the reactor. Second, as the calcination reaction has progressed
further for the particles at the entrance, the porosity of these particles is
already large due to the presence of the porous CaO in the outer layers of the
particle (see (4)), i.e., resistance for mass transport to the surface is
lower for the particles at the reactor entrance. In general, the
$\text{CO}_{2}$ concentration in the center of the particle is larger than in
outer particle layers.

The reduced-order model results in lower $\text{CO}_{2}$ mass fractions than
the DEM/CFD simulation. Specifically, these deviations appear at the surface.
When integrated over the surface, these discrepancies are compensated. As a
consequence, the sum of $\text{CO}_{2}$ source terms are similar for the
finite-volume and the reduced-order model results, which results in very
similar $\text{CO}_{2}$ mass fractions across the entire fluid phase.

Figure 6: $\text{CO}_{2}$ concentration at the center cross section of the
reaction volume at time $t=2500\text{s}$. DEM/CFD results shown to the left
and reduced-order model results shown to the right of the black vertical line.

For the simulation results shown here, the reduced model from section 3
required approximately 40% of the computational effort needed by the finite
volume method from section 2.

## Acknowledgments

This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) - Project-ID 422037413 - TRR 287. Gefördert durch die
Deutsche Forschungsgemeinschaft (DFG) - Projektnummer 422037413 - TRR 287.

##  Appendix A Closure terms for the calcination model

The continuum model stated in section 2 requires several closures which are
given for completeness. The equilibrium pressure $p_{\mathrm{eq}}(T)$ of the
calcination reaction only depends on the temperature $T$. It is given by [23]

$\begin{split}p_{\mathrm{eq}}&(T)=101325\,\,\exp\big{(}17.74\\\
&-0.00108T+0.332\ln(T)-\frac{22020}{T}\big{)}.\end{split}$ (20)

The enthalpy of reaction $\Delta^{R}H(T)$ of the reaction (1) can be derived
from Kirchhoff’s law for the reactants. This yields

$\begin{split}\Delta^{R}H(T)&=\Delta
H_{\mathrm{form}}(\theta)|_{\theta=298\text{K}}+\int_{298K}^{T}c_{p,{CO_{2}}}(\theta)\,M_{CO_{2}}\\\
&+c_{p,CaO}M_{CaO}-c_{p,{CaCO_{3}}}M_{CaCO_{3}}\text{d}\theta,\end{split}$
(21)

where $c_{p,i}$ and $M_{i}$ denote the heat capacity and the molar mass of
specie $i$, respectively. The formation enthalpy $\Delta H_{\mathrm{form}}$ is
a constant of the reaction.

The effective diffusivity $D_{\mathrm{eff}}$ is computed from the binary
diffusivity of $\text{CO}_{2}$ in air $D_{b}$ and the Knudsen diffusivities
$D_{\mathrm{Kn},CaO}$ and $D_{\mathrm{Kn},CaCO_{3}}$

$D_{\mathrm{eff}}=\frac{\varepsilon}{\tau^{2}}\cdot\bigg{(}\frac{R}{\frac{1}{D_{b}}+\frac{1}{D_{\mathrm{Kn},CaO}}}+\frac{1-R}{\frac{1}{D_{b}}+\frac{1}{D_{\mathrm{Kn},CaCO_{3}}}}\bigg{)},$
(22)

where $\tau$ denotes tortosity. The Knudsen diffusivities are calculated from

$D_{\mathrm{Kn},i}=\frac{2}{3}r_{\mathrm{Pore},i}\sqrt{\frac{8R_{m}T}{\pi
M_{CO_{2}}}}\qquad i=\text{CaO},\text{CaCO}_{3}.$ (23)

The binary diffusivity of $\text{CO}_{2}$ in air is given by [9]

$D_{b}=\frac{1.343}{p}\,\left(\frac{T}{273}\right)^{1.75},\quad
p\>\mathrm{in\>Pa},\>T\>\mathrm{in\>K}$ (24)

Further, the heat and mass transfer coefficients $\alpha$ and $\beta$, are
calculated from

$\displaystyle\alpha$
$\displaystyle=\mathrm{Nu}\,\,\frac{\lambda_{\infty}}{L}$ (25a)
$\displaystyle\beta$ $\displaystyle=\mathrm{Sh}\,\,\frac{D_{b}}{L},$ (25b)

where the characteristic length $L$ is set to the square root of the surface
area of the facet of the finite volume tetrahedron. Nusselt number
$\mathrm{Nu}$ and Sherwood number $\mathrm{Sh}$ are approximated by the
correlation

$\displaystyle\mathrm{Nu}$
$\displaystyle=0.664\,\,\mathrm{Re}^{1/2}\mathrm{Pr}^{1/3}$ (26a)
$\displaystyle\mathrm{Sh}$
$\displaystyle=0.664\,\,\mathrm{Re}^{1/2}\mathrm{Sc}^{1/3}.$ (26b)

The Schmidt number $\mathrm{Sc}$ is calculated from properties of the
surrounding fluid

$\displaystyle\mathrm{Pr}$ $\displaystyle=0.71$ (27a)
$\displaystyle\mathrm{Sc}$
$\displaystyle=\frac{\eta_{\infty}}{\varrho_{\infty}D_{b}}$ (27b)

The advective mass transfer coefficient is obtained from Darcy’s law. This
results in

$C^{*}=\frac{K\varrho_{\infty}}{\mu_{\infty}L^{*}}$ (28)

with the characteristic length $L^{*}$ set to the radius of the volume
equivalent sphere of the local tetrahedron.

##  Appendix B Validation of the calcination model

We varied the number of tetrahedral cells used in the finite-volume model from
$800$ to $10,700$ cells to determine an appropriate resolution. We evaluated
the calcination degree as a function of the cell number for this purpose. The
calcination degree is defined as the ratio of $\text{CO}_{2}$ mass released
from a lime particle to the total initial mass of $\text{CO}_{2}$ bound in
$\text{CaCO}_{3}$ in the particle. It can be calculated from

$D_{\mathrm{calcination}}=\frac{\int_{0}^{t}\dot{m}_{\text{CO}_{2},\mathrm{\text{out}}}\,\mathrm{d}\tau}{\int_{\Omega}\rho_{\text{CaCO}_{3}}-\rho_{\text{CaO}}\,\mathrm{d}v}$
(29)

Figure 7 depicts the calcination degree as a function of time for a gas-phase
temperature of 900 °C and $\text{CO}_{2}$ mass fraction in the gas-phase of 0
$\frac{kg}{kg}$. The results are very similar for a number of cells above
2100.

Figure 7: Impact of discretization on the calcination degree of the
cylindrical particles

As a further validation, we compared finite-volume simulation results obtained
with 1700 cells to experimental data for a single spherical particle (diameter
$10\text{mm}$) from the literature [12]. This comparison is shown in Figures 8
and 9, where the calcination degree is plotted as a function of time for
different ambient $\text{CO}_{2}$ concentrations at constant gas-phase
temperature of $1150\,\mathrm{K}$ (top), and for a set of gas-phase
temperatures at zero ambient $\text{CO}_{2}$ concentration (bottom). As
expected, a larger ambient $\text{CO}_{2}$ concentration (Figure 8) slows down
calcination reaction, and an increase in ambient gas-phase temperature
accelerates calcination (Figure 9). The model is able to reproduce these
effects with reasonable accuracy.

Figure 8: Comparison of simulation results obtained with the finite-volume
model to experimental data [12], part I. Calcination degree for a sphere for
various $CO_{2}$ mass fractions. Figure 9: Comparison of simulation results
obtained with the finite-volume model to experimental data [12], part II.
Calcination degree for a sphere for various temperatures.

##  Appendix C Validation of the reduced-order model

Before using it in the simulations summarized in section 4, we validated the
single-particle reduced-order model with $r=5$ modes and ODEs by comparing it
to finite-volume simulations with 8000 cells. The fluid properties were set to
$T_{\infty}=1175\,\text{K}$, $y_{\infty}=0$ and $p_{\infty}=101325\,\text{Pa}$
and kept constant over time for this validation. The initial particle
temperature was $1175\,\text{K}$ and the gas phase was assumed to be in the
equilibrium state, such that there is no calcination reaction in the initial
state.

The results for the particle temperature, $\text{CO}_{2}$ mass density and
calcination degree averaged over all 800 cells are shown in Figure 10. For the
first approximately $60\,\text{s}$, $\text{CO}_{2}$ is driven out of the
particle because of the difference of the $\text{CO}_{2}$ concentration in the
particle and the fluid. The calcination reaction then starts and the mean
temperature of the particle sharply declines. A minimum in temperature is
reached at approximately $300\,\text{s}$. Subsequently, the heat due
convective heat transfer from the fluid to the particle exceeds the heat
consumed by the calcination reaction. In this first period, the calcination
reaction is mostly located at the surface of the particle. We stress that the
reduced model reproduces this first period particularly well. This is evident
from Figure 10.

At around $400\,\text{s}$ a second period begins which is marked by an
approximately constant average mass fraction of $\text{CO}_{2}$ inside the
particle. The calcination reaction takes place at a moving reaction front that
travels from the particle surface inwards. The reaction front can be observed
also in Figure 3, where a sharp gradient of the reaction progress can be seen.
While this moving reaction front persists, an equilibrium state is reached in
the sense that approximately the same amount of $\text{CO}_{2}$ is released at
the reaction front as is transferred to the fluid at the particle surface.
This indicates that the reaction is limited by the transport of
$\text{CO}_{2}$ in the particle. We note that the reduced model shows
differences in the prediction of the temperature and the gas density. As a
result, the reaction front travels effectively slower in the reduced model,
resulting in longer time of the total calcination reaction of the particle.
The slower movement of the reaction front in the reduced model is reflected in
a slower increase of the calcination degree in the reduced model between
approximately $1000\,\text{s}$ and $3000\,\text{s}$ compared to the finite-
volume model.

The last period is reached in the finite-volume model at around
$2700\,\text{s}$ and with the reduced model at around $3600\,\text{s}$. At
this time the moving reaction front reached the core and the state variables
converge to the respective fluid properties.

Figure 10: Spatially averaged (over all particles) calcination degree,
temperature, $\text{CO}_{2}$ mass fraction and gas-phase density in the
particle over time used for validation.

## References

  * [1] OpenFOAM, free cfd software, the openfoam foundation.  https://openfoam.org/.  Accessed: 2013-11-27.
  * Benner et al. [2015] P. Benner, S. Gugercin, and K. Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  _SIAM Review_ , 57(4):483–531, 2015.
  * Berner et al. [2020a] M. O. Berner, V. Scherer, and M. Mönnigmann.  Controllability analysis and optimal control of biomass drying with reduced order models.  _Journal of Process Control_ , 89:1–10, 2020a.
  * Berner et al. [2020b] M. O. Berner, V. Scherer, and M. Mönnigmann.  An observer for partially obstructed wood particles in industrial drying processes.  _Computers & Chemical Engineering_, 141:107013, 2020b.
  * Bes [2006] A. Bes.  _Dynamic process simulation of limestone calcination in normal shaft kilns_.  PhD thesis, Otto-von-Guericke-Universität Magdeburg, Universitätsbibliothek, 2006.
  * Bluhm-Drenhaus et al. [2010] T. Bluhm-Drenhaus, E. Simsek, S. Wirtz, and V. Scherer.  A coupled fluid dynamic-discrete element simulation of heat and mass transfer in a lime shaft kiln.  _Chemical Engineering Science_ , 65(9):2821–2834, 2010.
  * Duan et al. [2022a] S. Duan, B. Li, and W. Rong.  Numerical simulation study of mixed particle size calcination processes in the calcination zone of a parallel flow regenerative lime kiln.  _Materials_ , 15(13):4609, 2022a.
  * Duan et al. [2022b] S. Duan, B. Li, and W. Rong.  Study on gas-solid heat transfer and decomposition reaction of calcination process in an annular shaft kiln based on the finite volume method.  _Processes_ , 10(4):648, 2022b.
  * Fuller et al. [1966] E. N. Fuller, P. D. Schettler, and J. C. Giddings.  New method for prediction of binary gas-phase diffusion coefficients.  _Industrial & Engineering Chemistry_, 58(5):18–27, 1966.
  * Golub and Van Loan [2013] G. H. Golub and C. F. Van Loan.  _Matrix computations_.  JHU press, 2013.
  * Gorges et al. [2024] C. Gorges, M. Brömmer, C. Velten, S. Wirtz, E. I. Mahiques, V. Scherer, K. Zähringer, and B. van Wachem.  Comparing two ibm implementations for the simulation of uniform packed beds.  _Particuology_ , 86:1–12, 2024.
  * Hills [1968] A. Hills.  The mechanism of the thermal decomposition of calcium carbonate.  _Chemical Engineering Science_ , 23(4):297–320, 1968.
  * Illana et al. [2023] E. Illana, H. Merten, S. Wirtz, and V. Scherer.  DEM/CFD simulations of a generic oxy-fuel kiln for lime production,.  _Thermal Sciences and Eng. Progress_ , submitted, 2023.
  * Jendersie et al. [2023] R. Jendersie, A. Mjalled, X. Lu, L. Reineking, A. Kharaghani, M. Mönnigmann, and C. Lessig.  Neuropnm: Model reduction of pore network models using neural networks.  _Particuology_ , 86:239–251, 2023.  ISSN 1674-2001.
  * Krause et al. [2015] B. Krause, B. Liedmann, J. Wiese, S. Wirtz, and V. Scherer.  Coupled three dimensional DEM–CFD simulation of a lime shaft kiln—calcination, particle movement and gas phase flow field.  _Chemical Engineering Science_ , 134:834–849, 2015.
  * Krause et al. [2017] B. Krause, B. Liedmann, J. Wiese, P. Bucher, S. Wirtz, H. Piringer, and V. Scherer.  3d-dem-cfd simulation of heat and mass transfer, gas combustion and calcination in an intermittent operating lime shaft kiln.  _International Journal of Thermal Sciences_ , 117:121–135, 2017.
  * Mjalled et al. [2023a] A. Mjalled, R. Namdar, L. Reineking, M. Norouzi, F. Varnik, and M. Mönnigmann.  Parametric 3d convolutional autoencoder for the prediction of flow fields in a bed configuration of hot particles, 2023a.
  * Mjalled et al. [2023b] A. Mjalled, E. Torres, and M. Mönnigmann.  Reduced-order modeling framework using two-level neural networks.  In _Proc. in Applied Mathematics and Mechanics_ , 2023b.
  * Mohammadpour et al. [2022] K. Mohammadpour, A. Chitsazan, and E. Specht.  The CFD simulation of reactive flow in parallel flow regenerative shaft kilns using porous media model.  _Thermal Science_ , 26(2 Part A):1175–1183, 2022\.
  * Reineking et al. [2024] L. Reineking, J. Fischer, A. Mjalled, E. Illana, S. Wirtz, V. Scherer, and M. Mönnigmann.  Convective drying of wood chips: Accelerating coupled dem-cfd simulations with parametrized reduced single particle models.  _Particuology_ , 84:158–167, 2024.
  * Scherer et al. [2016] V. Scherer, M. Mönnigmann, M. O. Berner, and F. Sudbrock.  Coupled DEM-CFD simulation of drying wood chips in a rotary drum - baffle design and model reduction.  _Fuel_ , 184:896–904, 2016.
  * Shagapov and Burkin [2008] V. S. Shagapov and M. Burkin.  Theoretical modeling of simultaneous processes of coke burning and limestone decomposition in a furnace.  _Combustion, Explosion, and Shock Waves_ , 44:55–63, 2008\.
  * Silcox et al. [1989] G. D. Silcox, J. C. Kramlich, and D. W. Pershing.  A mathematical model for the flash calcination of dispersed calcium carbonate and calcium hydroxide particles.  _Industrial & engineering chemistry research_, 28(2):155–160, 1989.
  * Sommer et al. [2023] K. D. Sommer, L. Reineking, Y. P. Ravichandran, R. Skoda, and M. Mönnigmann.  Estimating flow fields with reduced order models.  _Heliyon_ , 9(11):e20930, 2023.  ISSN 2405-8440.