Review of Modelling of Pyrolysis Processes with CFD-DEM

In a pyrolysis reactor, organic polymers from biomass or plastic waste are thermally decomposed into volatile gases, condensable vapours (tar or bio-oil) and solid residues (char). Since these products may serve as building blocks for downstream chemical refinement or form the basis of bio-derived fuels, pyrolysis is thought to be instrumental in our progress towards a circular economy. A pyrolysis reactor constitutes a multiphase reactive system whose operation is influenced by many chemical and physical phenomena that occur at different scales. Because the interactions and potential reinforcements of these processes are difficult to isolate and elucidate experimentally, the development of a predictive modelling tool, for example, based on the CFD-DEM (discrete element method) methodology, is attracting increasing attention, particularly for pyrolysis reactors operated with biomass as feedstock. By contrast, CFD-DEM descriptions of plastic pyrolysis remain a challenge at present, mainly due to an incomplete understanding of their melting behaviour. In this article, we provide a blueprint for describing a pyrolysis process within the scope of CFD-DEM, review modelling choices made in past investigations and detail the underlying assumptions. Furthermore, the influence of operating conditions and feedstock properties on the key metrics of the process, such as feedstock conversion, product composition and residence time, as determined by past CFD-DEM analyses is surveyed and systematised. Open challenges that we identify pertain to the incorporation of particle non-sphericity and polydispersity, the melting of plastics, particle shrinkage, exothermicity on part of the gas-particle chemistry and catalytic effects.

Heat of vaporisation of water (J kg −1 ) h i Specific enthalpy of species i (J kg −1 ) k evap Modelling constant in equilibrium model for drying (kgm 3 mol −1 s −1 ) L Characteristic length of a particle (m) l(r ) Distance between two particles at radial distance r (m) l min Minimum conduction distance (m) n Refraction index (

Introduction
Pyrolysis is the decomposition of organic matter, both naturally occurring and synthetic, into a range of light weight compounds through thermal treatment in the absence of oxygen. The liberated compounds can be used as energy sources or building blocks for further chemical processing (Sadaka and Boateng 2009). In addition to being an end-to-end process on its own, pyrolysis is also encountered as a preceding step of other processes, such as combustion and gasification (Grønli and Melaaen 2000). Within the scope of coal combustion, for instance, pyrolysis can be employed as a means of stabilising the flame, as loose volatiles are released before the oxidation commences. This initial devolatilisation has been found to have a considerable influence on the subsequent combustion behaviour of coal particles (Solomon and Fletcher 1994;Zhang et al. 2022b).
While the reactors carrying out pyrolysis may be operated with any organic polymer as feedstock, there is an apparent interest in biomass or plastics for this role. By virtue of its wide availability and carbon-neutrality, biomass has certainly asserted dominance over other renewable energy sources during the past decade. In fact, bio-energy is the largest renewable energy source, accounting for two-thirds of the total renewable energy throughput. Currently, it occupies about 14% of the global energy production (World Bioenergy Association 2020). Pyrolysis, alongside other conversion routes, such as liquefaction, fermentation, and enzymatic conversion, enables the extraction of bio-oil from biomass. Once upgraded, this bio-oil has the potential to serve as an alternative to fossil fuels (Matayeva et al. 2019). Plastics, on the other hand, face several operational challenges during their conventional conversion procedure, i.e., mechanical recycling. These challenges include the intrinsic irreversibility of thermosetting plastics and the need for sorting and separation based on plastic type prior to recycling-a process that is naturally labour-intensive (Eze et al. 2021). Consequently, from 6300 Mt of plastics that had been produced until 2015, only 9% has been recycled, whereas an alarming 79% has ended up as landfill (Geyer et al. 2017). A recycling strategy based on pyrolysis, by contrast, does not suffer from the aforementioned limitations concerning recycling; it provides a strategy to convert thermosetting plastics to reusable substances and can readily be implemented on unsorted feedstocks, e.g., municipal plastic waste. Nowadays, plastics are even allowed in the same pyrolyser as biomass, in pursuit of energy optimisation and better product quality (Wang et al. 2021c).
The extensive list of pyrolysis products can be classified into three main groups: pyrolysisoil, pyrolysis-gas and char residue. Pyrolysis-oil, which is the desired product of most applications, is recovered by condensing the vapour collected during the process. The noncondensable fraction of this vapour is termed the pyrolysis-gas and can be utilised for heat and power generation or as a feedstock in manufacturing industrial chemicals, such as methanol and ammonia. The solid residue, which remains after the reactions, is rich in carbon and has excellent absorbent qualities. It is important for the metallurgical industry as a reducing agent, or for households as a cooking fuel or as a soil amender with the ability to sequestrate carbon (Grønli and Melaaen 2000). Clearly, the yield and the composition of pyrolysis product groups strongly depend on the composition of the feedstock. For example, a particle of woody biomass mainly encompasses three distinct polymeric components, namely cellulose (C 6 H 1 O 5 ) n , hemicellulose, and lignin. As listed in Table 1, this lignocellulosic biomass displays a different composition depending on the wood type from which it is sourced (Gellerstedt et al. 2009). The compositions may also fluctuate depending on the seasons and during the prepossessing of the feedstock (Williams et al. 2016). Due to its comparatively shorter chain lengths, the decomposition of hemicellulose commences at the lowest temperature Table 1 Compositions (mass %) of wood species used in pyrolysis (Gellerstedt et al. 2009 (200 • C), followed by cellulose at 240 • C and lignin at 280 • C. During the initial decomposition, both hemicellulose and cellulose produce more condensable vapour, while lignin predominantly yields char residue. The extracted bio-oil from hemicellulose mainly contains acids and furfural derivatives. Levoglucosan (C 6 H 10 O 5 ) is the most abundant species in biooil, which is obtained mostly from the conversion of cellulose (Rabinovich et al. 2009). The average molecular weight of the bio-oil is seen to increase with the amount of lignin in the feedstock, since the resulting vapour contains heavier compounds (Fahmi et al. 2008).
Most of our municipal plastic waste consists of polypropylene (PP). PP along with high and low density polyethylenes (HDPE and LDPE), polyethylene terephthalate (PET), and polystyrene (PS) are the most common plastic types which are being pyrolysed (Vijayakumar and Sebastian 2018). A forecast of plastic pyrolysis products is possible based on the proximate analyses that break-down a plastic into four components: moisture, fixed carbon, volatiles, and ash. Generally, high volatile and low ash contents are known to result in high yields of pyrolysis-oil (Anuar Sharuddin et al. 2016). Pyrolysis of HDPE, LDPE, PP, and PS are characterised by a higher pyrolysis-oil yield since they are rich in volatile matter (Vijayakumar and Sebastian 2018). PET, on the other hand, converts mostly into pyrolysis-gas, while its oil is unfeasibly acidic with benzoic acid (Çepelioǧullar and Pütün (2013)). Although pyrolysis seems to be a convenient pathway to mitigate the issues of plastic recycling, it also possesses limitations, mainly due to the chemical compositions of certain plastics. The appearance of polyvinyl chloride (PVC) in pyrolysis feedstock raises health concerns as it decomposes into toxic hydrochloric fumes when undergoing pyrolysis. Furthermore, the high sulphur contents found in the liquid fuel obtained from pyrolysis of plastics render it unsuitable for deployment in existing combustion engines without any down-stream purification (Eze et al. 2021).
Based on the operating conditions, pyrolysis can be categorised as either slow, fast or flash pyrolysis. Slow pyrolysis is characterised by slow heating rates 0.1−1 • C s −1 and long particle and vapour residence times (450−550s). These conditions promote the devolitilisation to occur in favour of the char residue. Also, if the residence time is long, then the produced vapour will inevitably partake in secondary decomposition reactions, further reducing the yield of pyrolysis-oil. On the contrary, fast pyrolysis requires much faster heating rates 10−100 • C s −1 , and the residence time (0.5−10s) of particles and the liberated vapour is minimised by employing rapid gas flow rates. In this heating regime, the devolatilisation produces more condensable vapours. Due to their rapid release from the reactor, these vapours are largely prevented from further decomposition. The flash pyrolysis features even faster heating > 1000 • Cs −1 and particle and vapour entrainment rates (< 0.5s). Consequently, the product will largely contain pyrolysis-oil and only a little amount of char residue (Uddin et al. 2018).

Fig. 1
Chemical and physical processes and interactions on different scales in a pyrolyser. The pyrolysis system is affected by the chemical composition at the molecular scale, by particle properties at the particle scale, and by operating conditions at the reactor scale A generic pyrolysis reactor epitomises a complex multiphase reactive system, wherein a decomposing particulate phase continuously interacts with a moving gas that, too, contains reactive species. In some applications, the solid phase may be a combination of the feedstock and inert bed material that acts as a heat capacitance and improves mixing and heat transfer inside the reactor. As illustrated in Fig. 1, pyrolysis can be viewed as a multi-scale problem, where attributes at different scales co-dependently co-influence the efficacy and outcome of the overall process, measured in terms of the final conversion, product composition, particle residence time, etc. Apart from the feedstock composition, the existence of catalytic species in the system also tremendously influences the rates and paths of decomposition reactions at the molecular level. For instance, the naturally existing potassium in biomass has proven catalytic effects on thermochemical conversion processes, including pyrolysis (Safar et al. 2019;Fan et al. 2021). At the particle scale, the physiothermal properties of a particle, such as its size, shape, density, heat capacity, and conductivity, govern the intra-particle distributions of temperature and composition and control the shrinkage behaviour. Lastly, at the reactor scale, operating conditions, such as temperature and inert flow rate, largely influence the interphase heat, mass, and momentum interactions. Moreover, these scale-attributed phenomena are inherently linked to functions on other scales as well. For example, the temperature and the velocity of the carrier gas contribute to the intra-particle temperature distribution, which in turn affects its local decomposition rates.
The most common choices for pyrolysis reactor systems include fluidised bed reactors, fixed bed reactors, auger reactors, and entrained flow reactors (Uddin et al. 2018). Moreover, pyrolysis may be operated as an auto-thermal process using a circulating fluidised bed setup, where the energy demand for thermal decomposition in the primary bed is supplied by the heat extracted from the combustion of char residue inside a secondary bed through a circulating heat carrier (Liu et al. 2017b). Recently, novel technologies, such as microwave pyrolysis, grinding pyrolysis, and solar pyrolysis have also been researched in this field ). Due to their large heat and mass transfer rates and mixing degrees, fluidised bed reactors are often preferred to carry out pyrolysis (Ku et al. 2019). However, fluidised beds, too, are not free of drawbacks. The backmixing effect and premature exits of particles, which are characteristic of fluidised beds, are absolutely undesired in pyrolysis (Wang et al. 2021b). Furthermore, a fluidised bed may not be a viable option for the conversion of plastics, as they melt rather easily causing adhesion in particles, resulting in de-fluidisation (Zhang et al. 2022a). Many experimental studies are available on pyrolysis of various feedstocks (biomass (Guizani et al. 2017;Septien et al. 2012;Shen et al. 2009), plastics (Mastral et al. 2002;Chaudhary et al. 2023;Stelmachowski 2010;Lin and Yen 2005)), whose main focus has been analysing the effect of operating conditions. Despite being essential in the process of developing and optimising real-life reactors, they offer the barest insight into the physiochemical interactions taking place inside the system. Additionally, building an experimental setup is expensive and time-consuming. Its operation is also fairly challenging because of the accumulation of tar, that ends up causing blocks or even erosion of equipment (Chen et al. 2020). Consequently, researchers have shifted to employing numerical simulations, or what is known as "virtual experiments", employing computational fluid dynamics (CFD) to investigate the behaviour of pyrolysis reactors across scales. In particular, CFD-based tools have permitted the analysis of many-particle systems without sacrificing intricate features thanks to the exponential growth in computing power over the past few decades.
In the past, pyrolysis reactors have been modelled using two main approaches, the two-fluid model (TFM) and the CFD framework combined with the discrete element method (DEM), on which we focus here. The two approaches differ in terms of the physical description of the particulate feedstock and possess distinct advantages and challenges. In the TFM, the particulate feedstock and the carrier gas phase are considered as two interpenetrating and coexisting continua that may exchange mass, momentum and heat. While these exchanges are driven by physical and chemical processes that live on the single particle level, the particulate phase is conceptually described on the level of an entire particle population in terms of the particle number or mass density, the mean particle velocity (or momentum density) and an enthalpy density. Strictly, these densities represent means that are defined or computed as averages over all particles in a control volume about the present spatial location and can be related to statistics of the so-called particle property distribution. The property distribution quantifies how the particles within a control volume differ in terms of property values, for example, the particle size, velocity or temperature. The evolution in space and time of the property distribution can be described by the population balance equation (PBE) (Rigopoulos 2010), but its solution is practically challenging owing to its high dimensionality. For the particulate phase, the two-fluid equations follow from the PBE as the equations for the zeroth and first order moments of the property distribution and include two distinct types of closures that require modelling. The first type encompasses all assumptions that have led to the formulation of the PBE, for example, a continuity hypothesis on the distribution of the particles in physical and property space (Hulburt and Katz 1964) or assumptions on twoparticle statistics influencing collisions (Marchisio and Fox 2013, Chapter 6). The second type originates from the moment reduction of the PBE and requires assumptions on how the particle-averaged behaviour of the particulate phase is modified by property-polydispersity. In this regard, the closure of collisional interparticle momentum transfer, for instance, is achieved with the aid of a pseudo-stress model (Lun et al. 1984). Despite these closure challenges, the TFM has been very popular for the investigation of pyrolysis reactors, as it has many practical advantages. First, since the TFM is a fully Eulerian description, it supports solution through standard CFD tools. Second, the number of fields describing the particulate phase is independent of the number of physical particles, rendering the model particularly economical for applications with a large number of physical particles. Third, particle-particle interactions are accounted for statistically and, therefore, do not require the computationally expensive contact search and tracking algorithms of CFD-DEM.
In the context of pyrolysis, most of the past contributions to the modelling of gas-particle flows have been based on the TFM, analysing the flow in fluidised beds (Lathouwers and Bellan 2001), spouted beds (Hooshdaran et al. 2021), entrained-flow reactors (Xiu et al. 2008), and auger reactors (Aramideh et al. 2015). The investigation of the effects of operating conditions on biomass pyrolysis in fluidised beds by Lathouwers and Bellan (2001) is considered as one of the seminal works. They concluded that the temperature of the fluidising gas is the most influential parameter for the yield of bio-oil. Ranganathan and Gu (2016) compared the product yields predicted by global, semi-global, and detailed chemical kinetic models of biomass pyrolysis using TFM and demonstrated that the detailed scheme produces the most accurate results. Xiong et al. (2016) coupled TFM with a distributed activation energy model (DAEM) to account for the varying bond strengths across molecules of the same component. Considering a multi-fluid formulation, Liu et al. (2017a) and Xue and Fox (2015) were able to compute particle size distributions and the density of pyrolysis feedstock by adopting the direct quadrature method of moments (DQMOM) to solve the PBE.
Within the scope of the CFD-DEM approach, on the other hand, the dispersed particles are considered as separate individual entities, possessing distinct sets of degrees of freedom. As mass, momentum, and heat are exchanged between a particle and the carrier medium, the particle's degrees of freedom evolve in time. Both the number of degrees of freedom per particle and their physical meaning depend on the fidelity of the particle-level model. For example, a single particle temperature is sufficient if the heat conduction within the particle is much faster than the rate at which heat is supplied to the particle, while several temperature degrees of freedom are needed to parameterise the intra-particle temperature profile of thermally thick particles. For very small particles, temperature could even be entirely removed from the list of degrees of freedom if the particle and surrounding gas equilibrate (nearly) instantaneously to a common temperature. From the perspective of the continuous carrier phase, the DEM particles constitute point sources of mass, momentum, and enthalpy and are, therefore, often referred to as point-particles. In particular, the flow around the particles is not spatially resolved; this is in line with population balance models and the TFM and causes CFD-DEM to realise a meso-scale description. The dynamics of DEM pointparticles follow from single-particle mass, momentum, and enthalpy balances, including all particle-particle interactions. By design, the CFD-DEM framework can accommodate multiparticle interactions and particle-particle polydispersity without top-level approximations, and amendments of the point-particle model are readily incorporated. Furthermore, the CFD-DEM approach enables the prediction of residence times of individual particles, a feature that is absent in TFM ). This information is instrumental during the design and optimisation cycle of a continuous pyrolysis system to prevent both premature exits and accumulation of fully converted feedstock particles.
While less common than TFMs, the CFD-DEM methodology has been applied to many pyrolysis systems over the course of the past decade, including fluidised beds (Bruchmüller et al. 2012;Wang et al. 2021b), packed beds (Mahmoudi et al. 2016b), entrained flow reactors (Gao et al. 2021a), auger reactors (Qi and Wright 2020), batch reactors (Cordiner et al. 2018), and rotary drums . Papadikis et al. (2008Papadikis et al. ( , 2009b pursued a hybrid TFM-DEM modelling approach to model the pyrolysis of biomass in fluidised beds with inert bed material. Here, the inert particles were modelled as a continuous phase using a TFM description, while the feedstock was represented by DEM point-particles. Beyond pyrolysis, the applications of CFD-DEM contain many other multi-phase reactive systems such as gasifiers (Ku et al. 2015a), furnaces (Wei et al. 2019) and shaft kilns (Krause et al. 2015(Krause et al. , 2017. Pyrolysis is a maturing technology with the potential to have a positive impact on society and the environment, as it is a building block for a circular economy based on the recycling and sustainable upgrading of spent resources. In line with our objective to realise these benefits optimally at an industrial scale, CFD-DEM as a computer-based design and analysis tool offers valuable and explanatory insights into the pyrolysis processes that are difficult (or impossible) to gain through experimental investigations. Although many successful attempts at studying pyrolysis of biomass with CFD-DEM are reported in the literature, to the best of our knowledge, it is yet to be applied to the pyrolysis of plastics. We believe that a succinct review of the CFD-DEM framework with regard to pyrolysis systems, and a reevaluation of the insights provided by past studies will motivate researchers to not only analyse more systems in depth in the future, but also to develop effective tools to enhance the usability of CFD-DEM in this context. Our work contrasts with many reviews done by researchers on related topics, which include chemical kinetics of plastic pyrolysis (Armenise et al. 2021) and biomass pyrolysis (Papari and Hawboldt 2015), application of DAEM in biomass pyrolysis (Cai et al. 2014), co-pyrolysis of biomass (Abnisa and Wan Daud 2014), particle-scale models for biomass pyrolysis (Luo et al. 2022), biomass pyrolysis with CFD (Xiong et al. 2017), and MFiX-based CFD simulations of biomass pyrolysis (Lu et al. 2021a). There are several conversion methods, such as torrefaction and hydrothermal carbonisation, which are closely related to pyrolysis but which we do not discuss here because they have hardly been investigated using CFD-DEM tools.
This review article is structured as follows. In Sect. 2, we summarise the fundamental principles of the CFD-DEM framework while presenting the governing equations of both gas and solid phases encountered in a typical pyrolyser. Several potential simplifications for the descriptions of each phase as well as for interphase interactions are also detailed based on the characteristics of the targeted systems. Section 3 reviews the numerous physical mechanisms occurring in a pyrolyser, their influence on the overall process, and possible modelling techniques that lie within the scope CFD-DEM. Section 4 is reserved for a discussion of the impact of operational parameters on the results of pyrolysis, as obtained from past CFD-DEM studies on the subject. In Sect. 5, we identify open challenges with regard to modelling of pyrolysis using CFD-DEM, before conclusions are drawn in Sect. 6.

Carrier Gas
In the CFD-DEM framework, the continuous gas phase is described in a Eulerian fashion based on field variables that vary both along the spatial coordinates x and in time t. The mechanical field variables encompass the gas velocity u g and the mechanical pressure p and are, strictly, obtained by applying a spatial averaging or filtering operator to the true, spatially completely resolved fields. For conciseness, the averaging operation is not indicated notationally here, but we emphasise that its application is implied throughout. In practice, the averaging length scale is often associated with the grid on which the governing gas phase equations are solved using a direct discretisation scheme. Alternatively, the averaging length scale may be chosen independent of a discretisation in such a way that the gas volume fraction ε g obtained upon filtering is insensitive (or hardly sensitive) to the length scale. Conceptually, the averaging results in a removal of small-scale gradients from the describing fully resolved fields, resulting in averaged fields that vary smoothly on the filter scale. Within the scope of a CFD-DEM description, these averaged fields are the principal descriptors of the continuous phase. Particularly in the presence of turbulence, the influence of the residual scales that were eliminated by the filtering on the behaviour of the averaged fields may be accounted for through an amendment of the single phase constitutive laws that govern, for example, viscous stresses or heat conduction. In this review, however, we exclude both turbulence and the explicit modelling of residual scales other than through gas-particle mass, momentum or heat exchanges.
The filtering approach mentioned above was pioneered by Anderson and Jackson (1967) and leads, upon adaptation to a DEM description of the particulate phase, to the following mass and momentum balance equations governing the gas velocity u g and the mechanical pressure p, The right hand side of Eq. (1) includes the mass S mass j that is released to the gas by particle j per unit of time. The Dirac δ-distribution δ(x − x j (t)) indicates that the mass is released at the present location x j (t) of particle j and that the associated mass density is, therefore, concentrated entirely at x j . Similarly, the momentum source terms on the right hand side of Eq. (2) encompass the momenta S mom j imparted to the carrier flow by the dispersed particles j = 1, . . . , N p . Lastly, g (in [N/kg g ]) indicates a volume force that may be instrumented to accommodate gravity, for instance, while τ denotes the viscous stress tensor. Frequently, the gas phase is assumed to behave like a Newtonian fluid such that τ is given by where μ g is the dynamic viscosity and I denotes the identity tensor. The gas density ρ g that appears in Eqs.
(1) and (2) is a constitutive property that defines the material make-up of the gas and is linked to the gas temperature T g and the thermodynamic pressure P through an equation of state. In the modelling of pyrolysis reactors, the gas phase is often considered to behave like an ideal gas, whence the three thermal state variables obey where R is the universal gas constant and W denotes the mean molecular weight of the gas molecules. If the gas velocity is much smaller than the speed of sound, then acoustic pressure perturbations propagate almost instantaneously through the reactor, causing the thermodynamic pressure P in Eq. (4) to be nearly homogeneous (but possibly time-dependent), P = P(t). In this case, the mechanical pressure p in Eq.
(2) remains distinct from P and takes on the role of a kinematic variable akin to a multiplier that enforces the continuity equation (Eq. (1)). In reactors that are open to the environment or that are operated at constant pressure by means of a movable piston, P is prescribed as the ambient pressure. If a reactor is operated at a prescribed volume V (t), by contrast, P can be computed in terms of the total reactor volume from Eq. (4). Physically, the assumption that the thermodynamic pressure P is homogeneous and independent of p is associated with a low Mach number, Ma 0.3. In practice, the low Mach number assumption is very advantageous because it permits the thermal state (ρ g , T g , P) of the gas to be defined independent of the gas motion. For example, in a constant pressure reactor and for given W , the gas density (and, hence, the entire thermal state) is only controlled by T g .
At larger Mach numbers, the propagation of acoustic pressure perturbations is tied in a much stronger fashion to the gas flow field and the almost instantaneous, omnidirectional equidistribution of pressure that is so characteristic of small Mach numbers ceases. The thermodynamic pressure P then coincides with the mechanical pressure p, P = p(x, t), which acquires the meaning of a kinetic field variable. In the model formulations for pyrolysis reactors we review in this article, the low Mach number assumption is invoked throughout and we, therefore, maintain the simplifications it entails below.
In Eqs.
(1) and (2), the constriction of the flow domain due to the presence of pyrolysing particles is accounted for through the gas volume fraction ε g . Since ε g is the complement of the particle volume fraction ε p , it can be linked to the volumes V j (t) presently occupied by all particles j = 1, . . . , N p according to If the gas temperature T g and the mean molecular weight W are specified and if we leave aside the particle phase variables for a moment, then Eqs. (1) through (3) are completely closed in terms of the main gas phase fields u g , p and ρ g . In pyrolysis reactors, however, external heating, conductive heat losses to the particles, radiation, and gas phase chemical reactions induce temperature heterogeneities that evolve through conductive and convective heat transfer. Jointly, these effects enter an evolution equation for a calorific variable that quantifies both the sensible heat and the energy stored chemically within the gas. In open or constant pressure systems, this calorific variable is commonly chosen as the gas enthalpy h g (in [J/kg g ]), where h i denotes the specific enthalpy of gaseous species i, Δh 0 f ,i is its enthalpy of formation at the reference temperature T ref and C p,i represents its specific heat capacity at constant pressure. In Eq. (6), the chemical composition of the gas is represented in terms of the mass fractions Y g,i of all gaseous species i = 1, . . . , N g . Note that, by definition, the species mass fractions sum to unity.
The evolution equations for the species mass fractions follow from a material balance of the species-specific masses upon localisation and read where S mass−reac g,i is the rate at which species i is produced (S mass−reac g,i > 0) or consumed (S mass−reac g,i < 0) due to gas phase chemical reactions. Similarly, S mass j,i represents the mass of species i released to the gas phase by particle j; jointly, the source terms S mass j,i associated with all species yield the gas-particle mass exchange term S mass j in Eq. (1), We emphasise that the works reviewed in this article involve reactor setups that are operated exclusively in pyrolysis conditions, where the injected gas flows exclusively comprise inert species. Therefore, the term S mass−reac g,i in Eq. (7) is a consequence only of the secondary thermal decomposition of gas species produced and released by the decomposing feedstock and does not account for combustion or partial oxidation reactions that are possible in actual pyrolysis reactors, which may contain oxidising species in minute quantities.
Strictly, the species-specific diffusion velocities U i with follow from the binary diffusion coefficients through the Maxwell-Stefan equations. In computational practice, however, the solution of these equations for U i is not only expensive, but also difficult to take advantage of because U i drives an advective term whose discretisation may elicit unphysical artificial diffusion unless a costly Riemann solver is invoked. Poinsot and Veynante (2005) summarise two common approximations with different physical degrees of fidelity that result in a conversion of the advective U i -term into a smoothing diffusion term with far less stringent resolution requirements. The first approximation is based on the assumption that all binary diffusion coefficients are identical and equal to D g , yielding In a more general approach, Hirschfelder et al. (1964) assumed a single binary diffusion coefficient per species, yet, irrespective of the second species it diffuses into. This admits a species-unique coefficient of diffusion, given by D g,i , and results in the diffusion velocity approximation with (Wilke 1950) where D ki is the binary mass diffusion coefficient pertinent to diffusion of species k into species i, while X g,i denotes the mole fraction of species i. Unlike Fick's law in Eq. (10), the Hirschfelder and Curtiss approximation violates Eq. (9), causing Eq. (7) to be inconsistent with Eq. (1). This mass error can be eliminated by amending the right-hand side of Eq. (11) by a correction velocity u c g , where W i is the molecular weight of species i. In an alternative method, Eq. (7) is solved only for N g − 1 species, and the mass fraction of the last species is computed from the condition that all mass fractions sum to unity. In contrast to the approach based on u c g , this workaround results in an aggregation of the entire mass conservation error on the last species and should only be used if the last species is abundant (Poinsot and Veynante 2005).
If the gas enthalpy h g of Eq. (6) is employed to measure the chemical and sensible energy storage of the carrier gas, then Eq. (7) is complemented by the balance law where S heat j denotes the heat released by particle j. The flux vector q that appears in Eq. (14) accounts for the transport of heat across the bounding surface of a material gas volume; it is traditionally modelled as and encompasses two contributions. The first term in Eq. (15) corresponds to Fourier's law and represents the conductive heat transfer through the bulk gas that is, physically, mediated by molecular collisions and molecular diffusion. The intensity of this molecular heat migration is measured by the heat conductivity λ g . The second contribution serves to amend the heat flux by the net species-bound enthalpy transport. The third term on the right-hand side of Eq. (14) represents the heating of the gas by viscous tractions. In low Mach number flows, however, the viscous heating term can be neglected and the mechanical pressure in Eq. (14) may be replaced by a distinct thermodynamic pressure P = P(t), as in the equation of state above (Eq. (4)) (Poinsot and Veynante 2005, Section 1.2.1). Introducing these simplifications along with Eqs. (1) and (6) into Eq. (14), the enthalpy balance can be recast in terms of the gas temperature, yielding where C p,g denotes the gas' heat capacity at constant pressure,

Particulate Phase
In the applications we target in this review, the carrier gas is charged with dispersed particles that represent the pyrolysis feedstock and, possibly, inert bed material. Contrary to the gas, the particulate phase is observed from a Lagrangian perspective (DEM) and modelled as a collection of evolving point-particles. Specifically, the motion of particle j, translating and rotating with velocities v j and ω j , obeys Newton's second law (Fitzpatrick 2021), where x j , m j and I j are the position, mass and the inertia tensor of particle j. The symbols F c j , F d j , and F b j , moreover, represent the forces experienced by the particle due to particleparticle or particle-wall contacts, drag, and buoyancy, respectively. Commonly, the torque τ j is computed solely based on particle-particle contact interactions, disregarding the effects of forces exerted by the gas phase. This is tantamount to assuming negligible variations in the flow properties, particularly density and velocity, over the length scale of a particle. Note that in Eq. (18c) the shape of the particle is assumed to be spherical. In the case of non-sphericity, the term −ω j × I j ω j needs to be added on the right-hand side to account for the particle's directional variation of inertia and the orientation of the particle must be tracked, for example, through a set of quaternions as additional degrees of freedom (Diebel 2006;Zhao and van Wachem 2013).
In pyrolysis, the change of particle mass is controlled by two phenomena, drying and devolatilisation. Drying describes the evaporation of contained liquid water and requires the supply of heat through exothermal intra-particle reactions or by external heat transfer from the ambience. Devolatilisaion, by contrast, is the process of liberation of gaseous species from particles, produced as a result of the thermal decomposition of the feedstock. In the following, we denote the rates of drying and devolatilisation byṁ evap j andṁ devol j , respectively, These aspects are revisited along with their corresponding modelling approaches in subsequent sections of this article, specifically, drying in Sect. 3.3 and devolatilisation as part of chemical reactions in Sect. 3.4. Besides the mass m j , the thermochemical state of a particle is characterised in terms of the particle temperature T j and a chemical composition Y j = Y j,s | s ∈ {1, . . . , N s } . Here, Y j,s is the mass fraction of the sth particulate species, which represents either a solid phase educt or product or inherent moisture. From the definition Y j,s = m j,s /m j , the evolution of a species mass fraction in a particle is obtained upon time differentiation as where m j,s is the mass of species s contained in particle j and dm j,s /dt describes the rate-of-change of m j,s due to drying or intra-particle chemical reactions.
In a lumped description based on T j , a particle is assumed to possess a homogeneous internal temperature distribution. By an enthalpy balance, T j changes in time according to where denotes the specific heat capacity of the particle andQ j accounts for the radiative and conductive heat transfer to the particle. The modelling approaches that relateQ j to the particle's thermochemical state and the ambient conditions are reviewed in Sect. 3.1 below.
In the investigation of pyrolysers, the mass transfer from particles to the continuous phase is commonly assumed to be rate-limited by chemical kinetics, that is, the rate at which gaseous species are liberated from the free surface of the particles, rather than by mass transfer. In view of the generally small size of pyrolysis feedstock and the formation of cracks in the fuel particles undergoing pyrolysis, Rabinovich et al. (2009) emphasised that this assumption is well suited. By Spalding's mass transfer theory (Spalding 1953), the mass flux emitted from a particle in the continuum regime decreases inversely proportionally to the particle diameter, whence mass transfer limitations apply to large particles. Additionally, the formation of cracks and the small particle size promote a low intra-particle mass transfer resistance, facilitating the transport of gaseous species within the particle. Accordingly, all formed gaseous species are instantaneously released into the carrier gas and the interphase mass transfer term can be formulated without reference to a mass transfer coefficient, whereṁ devol j,i is the production rate of gas species i = {H 2 O} in particle j following its devolatilisation. Implicit in Eq. (23) is the assumption that water vapour is not produced by chemical reactions. If this assumption is invalid, thenṁ evap j needs to be complemented by a termṁ devol j,H 2 O that accounts for the chemical production of water vapour.

Discretisation and Gas-Particle Coupling
To simulate a flow laden with particles while accommodating the physical and chemical interaction mechanisms, a discretisation strategy for the partial differential equations and a strategy for coupling the Eulerian fluid fields and the states of the Lagrangian particles must be chosen. An important aspect of the simulation framework for predicting particle-laden flows is how the spacing of the grid points on which the discretisation of the fluid phase equations is based compares to the size of the particles. This aspect is commonly quantified in terms of the ratio of particle size to fluid mesh spacing. Currently, there exists a substantial body of research focusing on fully resolved simulations of particle-laden flows. In these frameworks, the particles are still described in a Lagrangian framework, but spatially resolved. Similarly, the fluid flow and the fluid-particle boundary layers are completely resolved using Direct Numerical Simulation (DNS). This approach is frequently referred to as particle-resolved (PR) DNS or PR-DNS (Tenneti and Subramaniam 2014). The advantage of employing such an approach is that it minimises the reliance on empirical parameters since all significant length and time scales are resolved through the methodology. On the downside, a considerable computational cost is incurred, even when dealing with a small number of particles in the flow, rendering the method only viable for low to moderate Reynolds number flows. As a result, the PR-DNS framework is not suitable for predicting the behavior of most realistic pyrolysis processes. Nonetheless, it can still be employed to investigate highly specific information within very small, typically academic, cases of the pyrolysis of a few particles.
In most particle-laden flow frameworks, the particle surfaces are not fully resolved by the Eulerian mesh, and the interaction between the fluid and the surfaces of the individual particles must be modelled. Here, the Eulerian mesh cells are larger than the individual particles. A commonly used framework (Eaton 2009;Kuerten 2016) is the particle-sourcein-cell Eulerian-Lagrangian (PSIC-EL) framework (Crowe et al. 1977), where the particles exchange information only with the cell in which they are present. The common assumption is made, that the Eulerian mesh cells are at least 10 times larger than the particle diameter, although even then errors of up to 10% can occur (Evrard et al. 2021). While various filtering techniques have been developed recently for a more accurate coupling between the Eulerian mesh and the Lagrangian particles (Capecelatro and Desjardins 2013;Evrard et al. 2019), these have not been used in the context of modelling pyrolysis processes with CFD-DEM.

Physical and Chemical Mechanisms
The CFD-DEM framework is kinetically driven by different chemical and physical processes that control not only the gas-particle mass, momentum, and heat exchanges, but also cause intra-particle alterations or changes on part of the carrier gas. At the end of this section, in Table 6, we summarise past CFD-DEM investigations of pyrolysis, with an emphasis on the modelling choices made by the investigators. The models referenced in Table 6 are discussed in detail in this section.

Heat Transfer
Accurate modelling of heat interactions, within and between the phases, is invaluable for an accurate prediction of the overall process for several reasons. The accuracy of many physical Table 2 Commonly used Nusselt number correlations in modelling interphase heat transfer during pyrolysis Ranz and Marshall (1952) Tavassoli et al. (2015) Nu j = 7 − 10ε g + 5ε 2 Particle Reynolds number Re j = ρ g u j L j μ g and Prandtl number Pr j = C p,g μ g k g phenomena, such as the drying and chemical kinetic rates, strongly depends on the accuracy of the heat transfer rates and the resulting temperature profiles. Also, the enthalpy budget of a reactor is strongly tied to the system's energy consumption, which is, in turn, mediated by the heat transfer in the system. Fundamentally, three modes of heat transfer can be distinguished in the system: conductive, convective, and radiative heat transfer. Following this distinction, the external heat transfer to a particle, which appears on the right hand side of Eq. (21), is decomposed intoQ with the termsQ cond j ,Q conv j , andQ rad j denoting the conductive, convective, and radiative contributions, respectively.
Conduction describes the heat transfer through direct contact from an element of high temperature to another at a lower temperature, without an apparent motion of matter. In fluids, this is facilitated by collisions and diffusion at the molecular level, whereas in solids, it proceeds through lattice vibrations and the movement of free electrons. The conduction in the gas phase is accounted for by the third term on the right hand side of Eq. (16), and needs no further detailing. On the other hand, the interparticle conductive termQ cond j encompasses two distinct mechanisms: the direct heat transfer in colliding particles, and the indirect heat transfer through the gas layer between non-colliding particles. The rates of these transfers, with reference to a particle j, are denoted byQ Considering the heat transfer from particle j 2 to another particle j 1 , the direct component is commonly computed with the model of Batchelor and O'Brien (1977), while the indirect or particle-gas-particle conduction is modelled according to the scheme introduced by Rong and Horio (1999) as follows, Here, r c is the contact radius and λ j , T j are the thermal conductivity and the temperature of particle j ∈ { j 1 , j 2 }, respectively. Additionally, r in and r out define the radial bounds of the corresponding heat transfer region and l(r ) is the thickness of the gas layer separating the Here, r c , l(r ), l min , r in , r out are the contact radius, normal gap, minimum conduction distance, and lower and upper radial bounds of heat transfer regions, respectively. Further, r j is the radius of particle j ∈ { j 1 , j 2 }, and θ is the thickness of the thermal boundary layer of the comparatively bigger particle particles at the radial distance r , while l min is a minimum (threshold) conduction distance. The threshold term in the denominator of Eq. (27) guards against a singularity of the heat transfer rate as l(r ) approaches zero and the particles touch. Schematics of both direct and indirect heat transfer mechanisms are shown in Fig. 2. The conduction term is sometimes neglected in the particle enthalpy balance equation, stating its inferiority compared to the convective term (Bruchmüller et al. 2012;Ku et al. 2019;. In this regard, the CFD-DEM study by Hu et al. (2019) reveals that the interphase convective heat transfer is approximately twice more effective than the interparticle conductive heat transfer. Lu et al. (2020) pointed out that artificially low spring constants are generally implemented in DEM to ensure numerical stability. Consequently, the overlaps between colliding particles are severely overestimated and, therefore, also the direct heat transfer rates in Eq. (26). In order to avoid this error, a correction may be imposed on the overlap prior to heat transfer calculations. Studying the impact of shrinkage patterns and operating parameters on biomass pyrolysis, Hu et al. (2019) implemented such a correction based on the scheme proposed by Zhou et al. (2010). The convective heat transfer, in contrast, is assisted by a bulk flow of matter that transports heat from a point of high temperature to a point of low temperature. Convective heat transfer can occur in two ways, either through free convection, where fluid flow is caused by density variations in a stagnant medium, or through forced convection, where fluid flow is controlled and directed. The termsQ conv j in Eq. (24) and the interphase heat transfer term S heat j in Eq. (16) originate from this same principle, however, with a change in direction by definition. Adopting the definition of an interphase heat transfer coefficient, denoted by α j for particle j, the heat transfer rates are given by where A j is the surface area of particle j. Here, α j is an indirect measurement of the thickness of the thermal boundary layer between the particle and the surrounding gas and may be related to the dimensionless Nusselt number Nu j of particle j according to Re j ≥ 1000 (43b) Kunii and Levenspiel (1977) with L j representing the characteristic length of the particle. Although Eq. (28) is very often applied in the context of pyrolysing DEM particles, it is, strictly, only valid if there is no mass exchange between the particle and the surrounding gas (S mass j,i = 0). Evaporation or devolatilisation, by contrast, induce a net gas-particle mass transfer that is accompanied by an outward-pointing Stefan flow, affecting both the species and temperature profiles across the gas-particle boundary layer. For spherical particles, the heat transfer rate in the presence of gas-particle mass transfer is obtained as (Turns 2012, Chapter 10) whereT indicates an estimated mean film temperature and may, for example, be computed using the 1/3-rule asT = T j + (T g (x j ) − T j )/3 (Abramzon and Sirignano 1989). Many empirical correlations exist and are referred in this context for the definition of the Nusselt number, as listed in Table 2 . Despite not having a volume fraction correction, the Ranz-Marshall model (Ranz and Marshall 1952) is seemingly the preferred choice of most researchers (Papadikis et al. 2009c;Mahmoudi et al. 2014;Cordiner et al. 2018;Ku et al. 2019;Lu et al. 2020) in studying the pyrolysis of spherical particles. On the other hand, the Tavassoli model (Tavassoli et al. 2015), which is a modified version of the Gunn model (Gunn 1978), is primarily employed for systems of non-spherical particles (Gao et al. 2021b;Lu et al. 2021d). In this case, the characteristic length is replaced by the diameter of a volume-equivalent sphere.  adopted different approaches for evaluating the heat transfer coefficient depending on the position of the particle inside the fluidised bed. The Ranz-Marshall correlation was applied for a particle elutriated by the fluidising gas, while the heat transfer coefficient of a particle in the dense bed was found by applying the penetration theory proposed by Kuipers et al. (1992). Simulating pyrolysis of cylindrical biomass particles using a glued-sphere method, Lu et al. (2021c) found an insignificant effect of the choice of Nusselt number correlation on the particle temperatures and the product compositions. Gao et al. (2021b) also reported a similar behaviour, having studied slow pyrolysis of cubic particles, and further argued that the system had been operated in a chemical kinetics dominated regime.
Radiative heat transfer enables a transfer of heat between elements that are not in direct contact, as it is mediated through electromagnetic waves. Both particles and gas molecules constantly emit and absorb radiation of heat energy to and from the system. The extent of radiative emission is regulated by their absolute temperatures and emissivities, whereas the absorption is controlled by the absorptivities and irradiation. A radiative heat transfer occurs if there is a disparity between the extents of these transactions; the temperature of an element increases following a surplus in absorption, and decreases if emission dominates. Often in fluidised beds, the gas phase is characterised as an optically thin medium, meaning that it transmits radiation of all wavelengths without any absorption. This idealisation makes it immune to heat transfer by radiation. The Stefan-Boltzmann law states that the radiative heat flux is proportional to the fourth power of a body's absolute temperature. According to this relation, radiative heat has an elevated effect in high-temperature applications, but is usually negligible if the system is operated at a moderate temperature. Stressing this fact, many researchers, also within the context of pyrolysis (Wang et al. 2021a;Hu et al. 2019;Lu et al. 2021d;Gao et al. 2021b), refrain from accounting for the effect of heat transfer by radiation, not only on the gas phase but also on part of the particles. In the following, we discuss how radiative heat transfer to a particle can be formulated and approximated with suitable models.
If G j represents the incident radiation at the location of particle j, then the rate of radiative heat transfer to the particle readṡ where σ s and a j are the Stefan-Boltzmann constant and the particle absorptivity, respectively. The exact magnitude of G j is derived from the radiative intensity field obtained as solution of the radiative transfer equation (RTE) (Modest and Mazumder 2021, Section 10.4). Conditional on effects of emission, absorption, and scattering in the medium, the RTE describes the propagation of radiative intensity through space. From a computational standpoint, solving the complete RTE is far too tedious and also expensive, as it involves the direction and wavelength as additional independent coordinates. To convert the RTE to a more practicable form, one can introduce a simplification to the description of the radiative intensity field, as in the methods of spherical harmonics (P N -approximation) (Jeans 1917) or discrete ordinates (S N -approximation) (Chandrasekhar 1960). A preferred practice, even in the context of pyrolysis modelling (Ku et al. 2019;Chen et al. 2020;Cordiner et al. 2018), is to simplify the radiative intensity term in the RTE with first order spherical harmonics (Modest and Mazumder 2021, Section 16.5). This is the so-called P 1 model, the simplest of its kind, that yields the incident radiation as solution of the following transport equation where a g , σ g , and n are the absorption and scattering coefficients, and the refraction index of the gas, respectively. The coefficient C accounts for the anisotropic scattering effects. The term S G on the right hand side of Eq. (35) can be instrumented to include the influence of radiation sources in the system. Following a low emissivity assumption for the gas phase, Bruchmüller et al. (2012Bruchmüller et al. ( , 2013) took a different approach to determine the incident radiation G j in Eq. (34). Their approximation reads withT being the arithmetic mean temperature of the particles found in near proximity of particle j. The eligible region when considering a particle for the calculation ofT was composed of the same and adjacent CFD cells as the cell that contained particle j at the time of measurement.

Momentum Transfer
The exchange of momentum, both between particles and between the particles and the gas phase, plays a significant role in the desired operation of pyrolysis reactors. Fluidised bed reactors frequently contain inert particulate material to improve the intensity and uniformity of momentum transfer from the carrier gas to the feedstock particles. This supports the suspension of the feedstock particles, allowing for an enhanced heat and mass transfer and promoting thermal decomposition while ensuring a uniform conversion across the pyrolysis reactor. Furthermore, a well-regulated momentum transfer is vital to maintain the balance between conversion and removal rates, in order to achieve a proper and efficient utilisation of the feedstock. Both the linear and angular momentum of a particle are affected when it collides with other particles in the system. Physically, a collision lasts for a small but measurable time span, during which repulsive and frictional forces act on the colliding particles, inhibiting their compression in the normal direction and the relative motion in the tangential direction, respectively. Correspondingly, the contact force F c j is typically decomposed into a normal component F n j and a tangential component F t j and evaluated separately based on the physics underlying each. Note that the term F c j in particle j's linear momentum balance equation (Eq. (18b)) represents the net collision force of all the pair-wise computed collisions in which the particle is involved at the current point in time. Similarly, τ j is the sum of the collisiondriven torques pertinent to particle j. For a contact between particles j 1 and j 2 , the torque exerted by particle j 2 on particle j 1 reads where F c j 1 ← j 2 is the contact force exerted by particle j 2 on particle j 1 and x c j 1 j 2 represents the point of contact. In DEM, the particle-particle contact force is typically modelled using a softsphere approach, whereby the points of first contact of both particles are interconnected by two spring-dashpot pairs, one along the normal direction and one in the tangential direction. The tangential force is additionally limited in magnitude by the Coloumb friction condition that specifies a maximum value in terms of the instantaneous normal force and the friction coefficient μ f . Accordingly, the normal and tangential forces F n j 1 ← j 2 and F t j 1 ← j 2 on particle j 1 due to its collision with particle j 2 are given by and respectively. Here, k j 1 j 2 , v j 1 j 2 , δ j 1 j 2 andδ j 1 j 2 represent the spring constants, damping coefficients, overlaps, and collision speeds, respectively, and the superscripts 'n' and 't' indicate the normal and tangential contact directions. The normal unit vectorn points from the centre of particle j 1 to that of particle j 2 . The tangential unit vectort, on the other hand, points along the direction of the initial tangential velocity of contact (Cundall and Strack 1979).
As mentioned in Sect. 2.2 above, the flow of gas is commonly assumed to leave the angular momentum of a particle unaffected. However, the gas flow influences the linear momentum of a particle via two distinct forces, the buoyancy force and drag. The buoyancy force F b j results from the spatial change in the mechanical pressure across the particle's extent and can be computed according to Conversely, the drag force results from tangential traction or viscous shear exerted by the gas on the particle surface, and is driven by the relative motion between a particle and the carrier gas. The main constitutive parameters affecting the shear-based momentum transfer are the local gas density and its dynamic viscosity. On part of the gas phase, drag gives rise to the momentum source S mom j in Eq.
(2) as the momentum transferred from a particle to the gas per unit of time. In practice, the drag force is commonly related to an interphase momentum transfer coefficient β (Anderson and Jackson 1967), that may either be calibrated experimentally in terms of the particle shape and Reynolds number or determined from spatially fully resolved calculations of the momentum boundary layer about the particle, Table 3 shows a summary of the drag correlations that have been used in the context of pyrolysis in the past. All of these drag models are particular to spherical particles and are based on the assumption that the DEM point-particles are locally homogeneously distributed, that is, ∇ε g ≈ 0 near particle j. In reality, the feedstock particles of biomass are non-spherical, in general. To address this non-sphericity, Lu et al. (2021d, c) coupled the drag correction of Ganser (1993), featuring a measure of particle sphericity, with the drag model of Gidaspow et al. (1991). In recent years, researchers have also implemented hybrid drag models that acknowledge the contrasting fluidisation behaviours of different particles as stipulated by their properties. According to the Geldart particle classification (Geldart 1973), biomass is regarded as type A and inert sand as type B by many research groups for the applications in fluidised beds. In the works of Lu et al. (2020) and Houston et al. (2022), the drag on Geldart B sand particles was computed using the Syamlal-O'Brien model (Syamlal et al. 1989), while the filtered drag model of Gao et al. (2018) was used to describe the drag on biomass particles. The drag model of Gao et al. (2018) pertains to the class of heterogeneous drag models and incorporates the effects of unresolved velocity scales on the drag experienced by Tar k 4 -4 .28 × 10 6 108.0 −42 a particle. In the context of CFD-DEM of pyrolysis reactors, heterogeneous drag laws were also employed by Wang et al. (2021a) and Gao et al. (2021b), who adopted, respectively, the energy-minimisation multiscale (EMMS) model of Yang et al. (2003) and the Di Felice-Hölzer/Sommerfeld hybrid drag model (Di Felice 1994; Hölzer and Sommerfeld 2008).

Drying
Moisture in the pyrolysis feedstock and its subsequent evaporation inside the reactor has a variety of effects on the pyrolysis process. First, from a product quality point of view, any moisture level is highly undesired as it ultimately contaminates the pyrolysis products (Mahmoudi et al. 2014). The emergence of water vapour increases the gas flow through the reactor, decreasing the ambient temperature and thereby curtailing both heat transfer and reaction rates. Further, a particle's fluidisation behaviour is also affected by its evaporation, following the change in mass. The moisture content of biomass feedstocks, widely used in pyrolysis and gasification processes, ranges between 10−20% by mass (Basu and Kaushal 2009). On the other extreme, many plastics, according to their proximate analyses, show moisture levels of below 1% (Abnisa and Wan Daud 2014). By consequence, the effect of drying on the pyrolysis of plastics is not as strong as in the pyrolysis of biomass, even to the extent that it can be excluded from the overall physical description. In a particle of biomass, water exists in one of two forms. Free water resides in and flows through the pores of a particle, while bound water is held in the cell walls by hydrogen bonds formed with the chemical constituents in wood. There is a thermodynamic phase equilibrium between the water vapour at the particle surface and the free liquid water inside. As a result, the surface vapour pressure of water coincides with its saturation vapour pressure at the current particle temperature (Machmudah et al. 2020). Evaporation occurs if the free-stream vapour pressure is lower than the saturation vapour pressure at the particle surface, and water vapour diffuses from the surface into the ambient gas. The migration of water vapour is accompanied by a diffusive flow of other species towards the particle that is balanced, at the semi-permeable gas-solid interface, by an outward-pointing Stefan flow (Akkus 2020). Here, the evaporation is driven both by diffusion and convection and, hence, mass transfer limited. A kinetic limitation, on the other hand, prevails if the vapour flow causes the vapour concentration at the surface to decrease faster than liquid water vaporises and phase equilibrium can be re-established. Physically, the main cause for a kinetic limitation of vaporisation may be an insufficient supply of heat to the particle, for example, due to slow conductive or radiative heating or insufficient exothermicity on the part of intra-particle chemical reactions. Conversely, if the ambient vapour pressure is larger than the saturation vapour pressure at the particle surface, water condenses onto the particle surface, releasing the heat of vaporisation. The released heat may either be transferred to the gas or drive endothermal chemical reactions within the particle. A larger gas temperature is accompanied by a larger heat supply to the particle, potentially mitigating any kinetic limitation and causing the particle temperature to rise. In turn, a larger particle temperature promotes evaporation because the saturation vapour pressure increases with temperature (Bellais 2007).
The influence of drying is introduced to the rest of the system through a drying rate, which quantifies the rate at which water is removed from a particle and added to the flow domain. Typically, it is computed by the heat sink model (Bruchmüller et al. 2013;Mahmoudi et al. 2016a;Ku et al. 2019), the chemical kinetic model (Gao et al. 2021b;Lu et al. 2021b), or the equilibrium model (Mahmoudi et al. 2016b). In the heat sink model, also referred to as the constant evaporation model, drying is assumed to be controlled by the heat transfer from the local surrounding to the particle, and any potential limitation by mass transfer is omitted. Thus, this approach is particularly suited for fast heating conditions where the diffusive evaporation below the boiling temperature is short-lived. The rate of evaporation of particle j, which takes effect at the boiling temperature denoted by T evap , is derived by the thermodynamic balance between the transferred heat and the enthalpy of the produced water vapour as (Bruchmüller et al. 2012) where Y j,H 2 O represents the current mass fraction of liquid water, while H evap and δt are the heat of vaporisation of water and the computational time step, respectively. One disadvantage of the formulation in Eq. (47) is that the evaporation rate is not differentiable in the particle properties, both at T evap and when the min function switches from returning the first argument to the second one. It also enforces a stagnant effect on the particle phase properties, which can be observed in the works of Ku et al. (2019) and Mahmoudi et al. (2014), in terms of momentary plateaus in both particle temperature and conversion profiles, in regions corresponding to drying in particles. Besides, unlike the models to follow, Eq. (47) does not account for drying below the boiling temperature of water (Bellais 2007).
In the chemical kinetic model, the conversion rate of liquid water to water vapour is assumed to obey a first order Arrhenius relation. In particular, the drying rate can be expressed asṁ where A evap and E evap are the pre-exponential factor and the activation energy of the phase change reaction H 2 O(l) −→ H 2 O(g), respectively. Herein, the Arrhenius coefficients are deliberately set to generate a steep increase in the evaporation rate as T j approaches T evap (Mehrabian et al. 2012). Compared to the heat sink model above, this model is both numerically well-behaved and easy to implement. In fact, it can be readily appended to the particle phase list of reactions in the chemical reaction model, along with the corresponding reaction heat, i.e., the heat of vaporisation of water (Bellais 2007).
The phase equilibrium between the liquid water (free and bound) and the water vapour residing in the pores of the particle underlies the formulation of the equilibrium model. The evaporation rate in this case is driven by the difference between the equilibrium and the present concentration of water vapour available at the particle surface (Wurzenberger et al. 2002).
Here, k evap is a modelling constant, whereas P * H 2 O and W H 2 O denote the temperature dependent saturation vapour pressure and the molecular weight of water, respectively. Mahmoudi et al. (2016b) adopted the equilibrium model in their numerical investigation of biomass pyrolysis and noticed that drying levels as high as 90% of the total moisture can be achieved at a temperature of 330 K. As opposed to the previous methods, this approach can capture the possible condensation of water vapour inside particles or parts of particles (if particles feature internal temperature distributions, see Sect. 3.6) with sufficiently low temperature and moisture content (Bellais 2007).

Chemical Reactions
Chemical reactions lie at the core of pyrolysis as they mediate the conversion of a used or little-valued feedstock into (re)usable and value-enhanced products. At incipient pyrolysis, the polymer chains within the feedstock are cut and reduced in length following the absorption of heat imparted by the ambience. This initial step is referred to as devolatilisation or primary decomposition. As a result, an array of organic vapours and a solid residue is produced, whose specific compositions depend on the chemical composition of the feedstock and the positions of the bond cleavages. The residual solid phase or char is largely non-reactive. The produced vapours may either be instantaneously liberated from a feedstock particle or spread through its pores and be gradually released over time. The fraction of the vapour which is condensable at room temperature is termed pyrolysis-oil or tar. Depending on the conditions to which the feedstock is exposed, particularly the heating rate and vapour residence time inside the reactor, a further disintegration of tar, or what is known as secondary decomposition, is also possible during pyrolysis. Lastly, the non-condensable fraction of vapour constitutes the pyrolysisgas. Generally, pyrolysis is considered to be an endothermal process overall. Yet, many experimental studies have hinted at the existence of exothermal reactions that may become significant in certain circumstances. In particular, this exothermicity has been observed in the form of central temperature overshoots in large feedstock particles on occasions where the reactors were operated at high pressure or temperature or with a high feedstock loading (Park et al. 2010;Basile et al. 2014;Di Blasi et al. 2014). One of the reasons for this effect is the exothermal nature of the decomposition of some components in the feedstock, for example, lignin in biomass (Di Blasi et al. 1999). On the other hand, this exothermicity has also been ascribed to a set of heterogeneous reactions of tar inside the particle's pores or at the particle surface that produce a secondary char (Anca-Couce and Scharler 2017; Di Blasi et al. 2017). A supporting evidence is the study of cellulose pyrolysis by Mok and Antal (1983) who observed a reduction of the energy required for pyrolysis with increasing production of char. Furthermore, the effect was found to be more pronounced in conditions that promote gas phase chemical reactions, such as high-pressure and long vapour residence times.
An evaluation of possible reaction mechanisms and their respective kinetic parameters is imperative for a quantitative description of pyrolysis. Additionally, since pyrolysis leads to many different product species (product groups), a plausible description of chemical reactions is required when optimising the process using simulations to maximise the yield of the most desired product. However, a precise determination of all possible reactions is not only difficult experimentally, but also to benefit from during modelling if the proposed mechanism is overly extensive. The inefficiency caused by the comprehensiveness is rather amplified in CFD-DEM, where particles are assessed individually and both temperature and composition may vary across particles. A general chemical reaction of pyrolysis is Educt ν 1 Product 1 + ν 2 Product 2 + ..., where a single educt is converted into either one or several products through an irreversible reaction. In Eq. (50), ν k with k = 1, 2, . . . represent the stoichiometric coefficients. Also, the chemical reactions are assumed to proceed at rates that follow first-order kinetics. In particular, the rate R n (in [mol m −3 s −1 ]) of reaction n is proportional to the current concentration of the educt species whose decomposition is described by the respective reaction. Accordingly, the rate of a particle phase reaction (primary decomposition) is given by where Y j,s is the mass fraction of species s in particle j, whose decomposition is described by reaction n, and k n, j represents the rate constant of reaction n in particle j. Similarly, the rate of a gas phase reaction (secondary decomposition) is described by with where X g,i is the mole fraction of gas species i participating in reaction s as the educt.
With the aid of the reaction rates in Eqs. (51) and (53), the rate of change of mass of species s in particle j is expressed aṡ and the rate of formation of gas species i aṡ where ν s,n and ν i,n are the stoichiometric coefficients of species s and i of the considered reaction n, while N pd denotes the number of primary reactions considered in the model. Furthermore, the species mass source/sink term S mass−reac g,i in Eq. (7) is given by where N sd gives the number of secondary decomposition reactions according to the employed chemical reaction model. In the subsequent parts of this section, we will focus on pyrolysis reaction models that have been utilised in CFD-DEM frameworks.
In the simplest approaches, the overall transformation of feedstock is described by a global chemical reaction. The feedstock is described in terms of a single species and the products are also lumped into groups on the basis of their final aggregate states as follows (Mahmoudi et al. 2014;Ku et al. 2019).

Feedstock
Tar + Char + Gas A slightly advanced scheme, such as the one proposed by Thurner and Mann (1981), includes either two or three competing reactions to describe the primary decomposition. An underlying assumption in the mechanisms of Eqs. (58) and (59) is that the characteristic time scales of gas phase reactions are much larger than the vapour residence time in the reactor, such that the released gas appears to be unreactive. However, a competing scheme still improves upon the global scheme in Eq. (58) as it permits capturing the shifts in yields of different product groups when the operating conditions change. In some competing reaction schemes, the primary decomposition of the feedstock is assumed to be preceded by an activation step that does not involve any change in composition or temperature (Beetham and Capecelatro 2019;Cordiner et al. 2018). Also incorporating the secondary decomposition of tar, such a pyrolysis reaction mechanism with lumped species reads Since the feedstock is represented irrespective of its precise chemical composition in the mechanisms in Eqs. (58) to (60b), they are only applicable to the specific feedstock for which the kinetic parameters have been calibrated. While the chemical make-up of plastics may be inferred from their types (Kayacan and Dogan 2008;Martín-Gullón et al. 2001), it is not that straightforward with biomass, where the composition varies significantly across its sources. Bradbury et al. (1979) first developed a kinetic model for the pyrolysis of cellulose, the most abundant component in woody biomass, including feedstock activation and its primary decomposition, given by two competitive reactions (char and gas are produced simultaneously). Since their model is based on a reduction of the earlier formulation of Broido (1976), it is now known as the Broido-Shafizadeh model. It has been incorporated into CFD-DEM models by Papadikis et al. (2009a) who investigated the pyrolysis of cellulose and by Hu et al. (2019) who, for simplicity, considered the biomass feedstock as cellulose. Later, Di Blasi (1994) extended the Broido-Shafizadeh model by the secondary decomposition of  tar. Miller and Bellan (1997) were the first to propose a constituent-based chemical reaction scheme for biomass pyrolysis by superposing the decomposition kinetics of its three major components cellulose, hemicellulose, and lignin, thereby extending the model by Di Blasi (1994). The kinetic parameters of the Miller-Bellan reaction mechanism are summarised in Table 4, with reference to the mechanistic model shown in Eq. (60a). Herein, the reaction heats of hemicellulose and lignin are assumed to be identical to those of cellulose. Furthermore, the char production is exothermal so as to reflect the experimentally observed increase in exothermicity with char yield. The Miller-Bellan reaction scheme has been employed in many CFD-DEM investigations of biomass pyrolysis (Bruchmüller et al. 2012;Rabinovich et al. 2009;Beetham and Capecelatro 2019;Wang et al. 2021a) due to its capability to describe the thermochemical conversion of any biomass feedstock with a lignocellulosic characterisation. Bruchmüller et al. (2012) simulated the pyrolysis of biomass sourced from red oak (cellulose-46.9%, hemicellulose-31.8%, and lignin-21.3%), using a Miller-Bellan scheme without secondary tar tracking. However, the yields of product groups were reasonably well predicted in three separate test cases that differed in terms of the operating temperatures (699 and 758 K) and particle moisture contents (7.0% and 16.1%). On all three occasions, tar and char showed a higher and a lower yield, respectively, compared to the experimental results. The deviation of the gas yield from the experimental measurements changed sign when the temperature of the fluidising inert gas was decreased from 758 to 699 K, while all other conditions remained fixed.
In recent years, more detailed chemical reaction models have been integrated into the CFD-DEM framework. Particularly the mechanistic models developed by the Chemical Reaction Engineering and Chemical Kinetics (CRECK) modelling group in Milan have been applied. An example is the reaction mechanism developed by Debiagi et al. (2018) which describes the primary decomposition in terms of 32 first-order chemical reactions involving 29 solid and 29 gaseous species. Here, the biomass feedstock is taken to consist of cellulose, three types of hemicellulose, three types of lignin, and two extractives. Product groups of tar and gas are composed of 22 and six species, respectively. In deriving the kinetic parameters, the authors examined 80 different biomass samples using thermogravimetric analyses (TGA) and fluidised bed experiments, subjecting the feedstock to different heating rates. Several CFD-DEM investigations of pyrolysis have incorporated this scheme to describe the thermochemical decomposition of biomass (Gao et al. 2021a;Lu et al. 2021c, b). Lu et al. (2021b) compared the yields of pyrolysis products predicted by the detailed kinetic scheme of Debiagi et al. (2018) for Loblolly pine and wood waste as feedstock with experimental values. Based on the percent composition of each main group (tar, char, and gas), the simulation results were comparable with those of the experiments. Particularly, the product yields of waste wood were nearly identical to the measured values at four different operating temperatures (400, 450, 500, and 550 • C). Gao et al. (2021a) used the same chemical reaction mechanism to simulate the pyrolysis of a single cylindrical particle (8 mm diameter, 19 mm length) as well as a pilotscale entrained flow reactor with polydispersed cylindrical particles (length-to-diameter ratio 5, diameters in the range 140 to 1791 µm), comparing the predictive fidelities of isothermal and non-isothermal particle descriptions. For the single particle, the non-isothermal model produced more accurate results, while in the case of the pilot-scale reactor the isothermal characterisation was better suited. Since, by design, the reaction model incorporated intraparticle transport, Gao et al. (2021a) surmised that coupling a particle-scale model in the second instance may have caused intra-particle effects to be accounted for twice. For us, however, this justification is hard to accept because it would then apply to the single particle test case as well. Recently, Houston et al. (2022) combined the primary decomposition mechanism of Debiagi et al. (2018) with the full secondary decomposition scheme developed by the CRECK group for gas phase kinetics (Ranzi et al. 2014;Debiagi et al. 2016;Pelucchi et al. 2019). At its furthest level of detail, the secondary decomposition model contains 509 species and 20,240 reactions. In order to reduce the associated expenses in the context of CFD-DEM, Houston et al. (2022) eliminated reactions that were either three-body or falloff reactions or involved radicals or reactants that are not present in their application. The reduced mechanism included 66 secondary reactions and 59 gaseous species.
In addition to the chemical kinetic models discussed above, there are many other models available in the literature that describe the devolatilisation process of biomass, individual plastics (Al-Salem and Lettieri 2010; Faravelli et al. 2001) and binary (Jing et al. 2013) and tertiary (Costa et al. 2010) plastic mixtures. However, these models have not yet been incorporated into CFD-DEM descriptions of pyrolysis reactors. The distributed activation model (DAEM) is an interesting modelling approach that has been employed in the context of biomass pyrolysis (Cai et al. 2014), aiming to address the observation that chemical bonds may vary in strength among molecules of the same species. This is done by assuming that the activation energy pertaining to a reaction is no longer constant, but normally distributed.
Particularly for biomass, there exists another branch of chemical kinetic models, called network models, that incorporates the chemical structure of the polymeric species in the determination of kinetic parameters and reaction stoichiometry. There are three main instances of such models, namely the Chemical Percolation Devolatilisation (CPD) model (Fletcher et al. 2012), the Flash-Distillation Chain-Statistics (FLASHCHAIN ® ) model (Niksa 2000), and the Functional Group-Depolymerisation, Vaporisation and Cross linking (FG-DVC) model (Solomon et al. 1988). Incidentally, all aforementioned models were invented in the context of coal pyrolysis and have then been extended to describe the pyrolysis behaviour of lignocellulosic biomass with suitable modifications. In the CPD model, the polymeric feedstock is assumed to comprises clusters of carbon rings adjoined through a labile bridge which acts as the active site during the decomposition. Within the scope of the CPD model of biomass (bio-CPD), three kinetic schemes have been derived for each representative species, cellulose, hemicellulose, and lignin, separately according to the reaction scheme

Bridge
Activated Bridge 2 Side chains 2 Light gas Char bridge + 2 Light gas (61) by assuming their decomposition behaviours to be independent. These are then superposed to arrive at an overall model which accounts for the variable compositions across feedstocks. The kinetic parameters relevant to the model in Eq. (61) are determined based on five parameters that describe the structure of the feedstock in addition to experimental measurements of pyrolysis products. Similar to the assumption on the structure of coal particles in the original work (Grant et al. 1989), the chemical structure of lignin is also represented by a Bethe lattice and the corresponding lattice statistics have been used in the model development. For cellulose and hemicellulose, the anomeric carbon, which connects to two oxygen atoms, is considered as the active site in the absence of aromatic rings. Furthermore, hemicellulose is modelled as an average between its major components, xylan and glucomannan. All the reactions in Eq. (61) follow first-order kinetics and involve Arrhenius-type rate constants which are evaluated with normally distributed activation energies (Fletcher et al. 2012). In the FLASHCHAIN model for biomass (bio-FC), biomass is represented by a chain copolymer of cellulose and lignin-like compounds, while hemicellulose is neglected. The solid reactant matrix is assumed to compris units of this copolymer with different chain lengths, ranging from monomers to infinite chains. The reactions proceed via a depolymerisation of the copolymers into smaller fragments, which are categorised based on size either as reactant, intermediate or metaplast. The reactant group is susceptible to further scission, producing either intermediates or metaplasts in the process, or may decompose to char by releasing light gases. The intermediates may similarly break down into metaplasts or condense into char links. The phase equilibrium between the tar species in the metaplast and the interacting gas phase drives the production of tar components, which are thought to be removed by a convective flow emanating from the solid matrix. Additionally, the metaplast groups may partake in recombination reactions and regress to intermediate fragments. All reactions in the proposed model are characterised by distributed activation energies and the reaction parameters are determined by fitting experimental measurements of cellulose and lignin pyrolysis. Furthermore, the stoichiometric coefficients are adjusted to the composition of lignin-based components in biomass (Niksa 2000).
The FG-DVC model comprises two sub-models; the FG model describes the release of gas as various functional groups in the polymeric feedstock decompose, while the DVC model governs the formation of tar and char as a consequence of depolymerisation, internal and external mass transfer, and cross-linking induced by heat treatment. The approximation of the chemical structure of the polymer involves 14 parameters (Solomon et al. 1988) and their calibration requires measurements from experimental techniques such as nuclear magnetic resonance (NMR) spectroscopy and field ionisation mass spectroscopy (FIMS). In its extension to describe the devolatilisation of biomass, an average macromolecule is characterised based on a van Krevlen diagram (plot of H/C versus O/C atomic ratios) to match the atomic make-up of the feedstock (Chen et al. 1998). This differs from the approach adopted in bio-CPD where the overall devolatilisation kinetics are a result of superposing the kinetics of individual representative components. However, similar to previous network models, the reaction rates are formulated based on normally distributed activation energies.

Shrinkage
As feedstock particles dry or decompose chemically and emit volatiles, their volumes may decrease by up to 70% of the original volumes (Davidsson and Pettersson 2002). Many experimental analyses of pyrolysis have affirmed that this volume reduction is influenced by several factors such as the type, size, and structure of the feedstock as well as the heating rate which may produce uneven rates of contraction along different directions of a particle (Davidsson and Pettersson 2002;Kumar et al. 2006). Correspondingly, an accurate description of how the geometry of a pyrolysing particle changes in time is extremely difficult. Further complications arise from other intriguing observations, for example, delayed shrinking or initial swelling of the feedstock particles (Zhou et al. 2013;Paulauskas et al. 2015;Caposciutti et al. 2019). Irrespective of the precise nature of volume changes, any change in the size or shape of a particle affects its fluidisation behaviour as well as heating and reaction rates (Wang et al. 2021b). Intra-particle processes will also be influenced by changes in the particle size or shape, through density and porosity . In the context of DEM, these dependencies may be accounted for through a sub-model that describes the shrinkage behaviour of a pyrolysing particle.
The existing CFD-DEM models for pyrolysis that we are aware of disregard the influence of particle drying on their shrinkage behaviour, and, further, the possibilities of swelling or delayed shrinking. Most model formulations employ either a constant size model (Papadikis et al. 2009b;Rabinovich et al. 2010) or a constant density model for this purpose. The constant size model appears to be the simplest and least realistic, as it features no apparent change in particle size. Considering a spherical particle with fixed diameter d j , the particle density ρ j adjusts to the current particle mass according to Conversely, the density is held constant in the constant density model and the diameter of a spherical particle undergoing pyrolysis is computed from the current mass m j (Cordiner et al. 2018;Chen et al. 2020), Ku et al. (2015b) pointed out that neither Eq. (62) nor Eq. (63) are physical as both the diameter and density of a particle simultaneously change during its decomposition. Bruchmüller et al. (2012Bruchmüller et al. ( , 2013 attempted to address this through a mass-proportional shrinkage model, where the particle diameter was calculated by Eq. (63), using, however, a particle density that varies linearly with the mass fraction Y j,B of virgin biomass in the particle j, Here, ρ B and ρ C are the bulk densities of the biomass feedstock and char, respectively. Di Blasi (1996) proposed a more flexible approach for modelling shrinkage during the pyrolysis of biomass and argued that the volume V j of any particle j consists of a solid volume V solid j and a pore volume V pore j , The solid phase is made up of virgin biomass as well as char that has been produced during the decomposition of the virgin biomass. The volume of the solid phase is assumed to vary linearly with the conversion of the virgin biomass according to where f 1 is a kinetic parameter, m B and m C are the current masses of virgin biomass and char, respectively, and m j,0 and V solid j,0 denote the initial mass and the solid volume of the feedstock particle. The pore volume of the particle, on the other hand, is controlled by two effects. First, it decreases as the feedstock is consumed, corresponding to the overall shrinkage of the particle. In particular, the pore volume is assumed to change linearly with the conversion of the particle, that is, the reduction of m B , from the initial pore volume V pore j,0 , to the final volume f 3 V pore j,0 at complete conversion (m B (t) = 0). Second, the pore volume increases as the solid phase disappears and the volume it occupied is filled with gases. Combining both effects, the pore volume of the particle is computed as where f 2 represents the fraction of volume recovered by particle internal pore space following the consumption of the solid phase. Like f 1 and f 3 , f 2 lies within the range from zero (total disappearance of particle as m B → 0) to unity (no shrinkage). In order to complete the specification of Di Blasi's model, the parameters f 1 , f 2 , and f 3 as well as the initial pore volume V pore j,0 must be calibrated based on the observed shrinkage behaviour of the biomass under consideration. Note that V pore j,0 may be deduced from the particle's initial porosity. Table 5shows that all shrinkage schemes discussed above can be formulated as special cases of the Di Blasi shrinking model. Additionally, the Di Blasi model parameters determined by several researchers for the feedstocks they considered are listed in Table 5.
Several CFD-DEM studies have been conducted, specifically in fluidised beds, with the objective of assessing the effect of the applied shrinkage model on pyrolysis behaviour. The single particle simulations of Papadikis and Gu (2009) indicated that a particle attains a higher temperature, both at its surface and the center, if less severe shrinking conditions are imposed via a Di Blasi scheme. This rather counter-intuitive observation may be explained by the larger surface areas of more slowly shrinking particles. Consequently, for a given time course of particle mass and composition changes, a slower size reduction results in larger heat transfer rates and a larger particle temperature. Based on simulations of the pyrolysis of thermally thin biomass particles in a fluidised bed, Ku et al. (2015b) observed that constant size shrinking leads to faster devolatilisation rates and greater residence times of particles compared to constant density shrinking. However, physically, their observation of the influence of the shrinkage model on the particle residence time is difficult to explain. Since the mass and heat transfer rates as well as the momentum transfer rate due to Stokes drag scale linearly with the particle size, we would expect a particle with a given mass to devolatilise quicker and be more susceptible to entrainment and removal if its size remained constant than if the density were unchanged. In a more recent work, Hu et al. (2019) investigated the pyrolysis of biomass with four different shrinking patterns and concluded that the residence time, in fact, increases with the severity of imposed shrinkage. Studies of both Hu et al. (2019) demonstrated only a minor effect by the employed shrinkage pattern on the product yields of biomass pyrolysis. By manually calibrating the Di Blasi parameters,  were able to achieve an acceptable agreement with the experimental particle size evolution, following the configuration f 1 = 0.3, f 2 = 0, f 3 = 0.3. Upon further investigation, the product yields were also found to be reasonably well predicted. However, we emphasise that this parameter set is not universal and pertinent only to a specific feedstock as well as, possibly, to a specific set of operating conditions.

Table 6
Summary of the CFD-DEM investigations of pyrolysing particles and pyrolysis reactors we review as part of this work   Sample intra-particle distributions of temperature (T ) and feedstock mass fraction (Y ) predicted by a an isothermal model, b a 1D continuous model, and c a 1D interfacial model for a spherical biomass particle undergoing pyrolysis. Here, T 0 and Y 0 are the initial temperature and feedstock mass fraction of the particle, respectively

Intra-particle Behaviour
When describing the evolution of a particle's degrees of freedom, specifically its composition and temperature in Eqs. (20) and (21), we considered the particle to be a homogeneous entity. This characterisation is corroborated by the assumption that a particle features a negligible resistance to intra-particle heat conduction and mass diffusion. Here, any amount of heat transferred to the particle from the outside is transported across its volume within an infinitesimal time span. As a result, the particle does not exhibit any internal temperature gradients and its decomposition will follow an identical route at a similar rate everywhere inside the particle, as shown in Fig. 3a). However, if the external rate of heat transfer exceeds the conductive heat transfer rate within the particle, the near-isothermality ceases and thermal gradients emerge in the outward direction of the particle. The spatially inhomogeneous temperature distribution may also induce intra-particle inhomogeneity on part of the chemical composition through drying and chemical reactions. On the other hand, compositional intraparticle heterogeneity may also result if the intra-particle species diffusion is slower than the species removal or absorption at the surface, a case that has not been considered within the scope of pyrolysis to our awareness. For heat transfer, the amenability to intra-particle temperature gradients is commonly assessed by analysis of the particle's Biot number a dimensionless quantity that corresponds to the ratio of the particle's internal thermal resistance to that at its surface. According to Bryden and Hagge (2003), a particle is characterised as thermally thick or non-isothermal if 0.2 < Bi j < 10. In this regime, the isothermal assumption leads to an overestimation of the particle's overall temperature and conversion, as revealed by many numerical investigations on pyrolysis (Bellais 2007;Luo et al. 2020;Ku et al. 2019) also with detailed chemical kinetic schemes (Gao et al. 2021b;Lu et al. 2021d;Gao et al. 2021a). In order to account for potential intra-particle inhomogeneities, the description of the point-particles can be amended by an intra-particle model. Based on the way of spatial discretisation, these models may be categorised as either 1D models (layer discretisation) or 3D models (control volume discretisation). The 1D models are limited to regular particle shapes and rely on the assumption that the particle's heat conductivity is isotropic. Two variants can be encountered in the context of pyrolysis modelling: the continuous model and the interfacial model. The continuous models are based on the one-dimensional heat conduction equation governing the temperature profile in a long plate (n = 0), a long cylinder (n = 1), or a ball (n = 2), respectively. Here, S heat−reac j,r describes the heat source/sink term due to chemical reactions and r corresponds to the thickness coordinate of the plate, the radial coordinate of the cylinder, or the radius of a sphere. Equation (67) is supplemented by boundary conditions that account for the symmetry at the particle centre and the external heat transfer rates at its surface, An example for the intra-particle distributions of temperature and feedstock mass fraction for a biomass particle is shown in Fig. 3b). Mahmoudi et al. (2014) employed the 1D continuous model to simulate the pyrolysis of biomass particles in a packed bed. The shapes of the particles were simplified as spheres having identical volumes, and, yet, a good agreement with the experiments was achieved, based on the system's mean temperature and mass loss, for a range of inlet temperatures (573−803K). The single particle simulations conducted by Gao et al. (2021b) with spherical and cylindrical particles, produced nearly identical results to that of the experiments, including the centre temperature overshoots. They further investigated the pyrolysis of cubic particles with Bi = 1.4 in a packed bed using both isothermal and 1D continuous intra-particle models. Intriguingly, the results suggested that the isothermal model required a considerably longer time for the complete conversion of the entire batch of particles, causing the total runtime to exceed the one required by the 1D models.
In the 1D interfacial model, by contrast, the volume of a particle is split into three main layers. Within the scope of biomass pyrolysis, these layers are referred to as the wet wood layer, the dry wood layer, and the char layer, in the order of appearance from the centre to the surface of a particle. The boundary between the wet and dry wood layers is called the drying front, while the devolatilisation front separates the dry wood and char layers. Each layer can be further discretised into sub-layers or grid points that may be used to resolve intra-layer temperature gradients. Formally, the main differences between the interfacial model and the continuous model lie with the distinct boundary conditions that apply to each of the main layers. For additional information on the assumptions and the governing equations of the 1D interfacial model, we refer to the pioneering work of Ström and Thunman (2013). Figure 3c) illustrates the partitioning of a particle into layers and indicates the intra-particle profiles of the temperature and feedstock mass fraction for an interfacial model. Wang et al. (2021b) adopted the interfacial method to simulate the pyrolysis of 0.75 mm poplar sawdust particles inside a three-chambered fluidised bed. Although the experimental feedstock contained ∼ 7% moisture, the particles in the model were considered to be completely dry and the interfacial model could, therefore, be reduced to having only two layers: an unreacted core and a reacted shell. Using a global kinetic model, Ku et al. (2019) analysed the effect of isothermal and 1D interfacial models on the pyrolysis behaviour of a single biomass particle (Bi = 0.67) suspended in a fluidised bed of inert particles. The results indicated that the devolatisation commences earlier and lasts longer in the interfacial model compared to the isothermal approximation. The authors also observed that the trajectories of the particles are strongly affected by the choice of the intra-particle model, and that the interfacial model results in a longer residence time.
Three-dimensional models, on the other hand, permit a more comprehensive description of intra-particle variability, including anisotropy, and allow for irregular particle shapes. In mesh-based approaches, the heat conduction or species transport equations on the volume occupied by a particle are spatially discretised using a direct discretisation technique. An alternative is the glued-sphere approach that can be initialised by replacing the cells in the mesh-based model with equal-volume-spheres. The evolution of the temperature of a composite sphere based on its position of appearance in a particle j is given by whereQ cond−intra j is the intra-particle conduction term evaluated by summing the heat transferred to the current composite sphere from all other spheres inside particle j with which it is in contact. Lu et al. (2021d) implemented the glued-sphere representation first on spherical and cylindrical particles to validate the model and, later, to investigate the pyrolysis behaviour of irregularly shaped biomass particles subjected to different operating conditions. Since the extent of heat transfer depends on the surface area perpendicular to the heat flow direction, they explicitly recorded the surface areas of the mesh cells composing the non-spherical particles during the initialisation of the glued spheres. Instead of the areas rendered by the generated composite spheres, the recorded areas were considered in heat transfer calculations. As expected, the difference between the surface and the centre temperatures increased with the size of a particle. Comparing with the results from an isothermal approach, it was found that the intra-particle variations are negligible below Bi = 0.41, where the isothermal assumption would become valid.
As pointed out by Papadikis et al. (2009b, c), the intra-particle models have a high computational demand and produce large output files in comparison to an isothermal model. A workaround to alleviate this loss of efficiency in the thermally-thick regime is by incorporating a corrected-isothermal model. Through a process of calibration, a corrected-isothermal model lumps the effect of intra-particle transport onto the external heat transfer and reaction rates. Then these modified rate expressions are used to compute a homogeneous temperature and composition for a physically thermally-thick particle. Luo et al. (2020) devised a corrected-uniform model with the aid of volume-averaged temperature and reaction rates predicted for a particle by a 1D continuous model. The derived model predicted the temperature and conversion of a 9.5 mm particle more accurately than a usual isothermal model, and its results were comparable with the 1D continuous model. Lu et al. (2020) calibrated the Nusselt number and reaction rates of a reactor-scale model by minimising the deviations of predictions from the results of a benchmark particle-resolved simulation.

The Effect of Operating Conditions on the Pyrolysis Efficacy and Outcome
In the past, the CFD-DEM framework has been instrumented by many investigators to assess the influence of operating parameters, such as the gas temperature, fluidisation velocity (Houston et al. 2022), reactor type (Lu et al. 2021d), or shape (Wang et al. 2021a), and the feedstock properties, for example, the particle size (Ku et al. 2019), shape (Lu et al. 2021d), and moisture content (Bruchmüller et al. 2012), on the overall pyrolysis behaviour or to optimise reactors with respect to product yields. In this section, we summarise key observations that may guide the design or operation of pyrolysis reactors. Operating Temperature: Rabinovich et al. (2009) analysed the effects of many operating parameters, including temperature, in the process of optimising the yield of bio-oil attained from pyrolysing biomass in an inert bed. After fixing the velocity of inert gas through the reactor at 1.1 times the hovering velocity of char particles produced during the process, they achieved a maximum yield of bio-oil of 72% for a temperature, particle size, and size of bed particles of 776 K, 440 µm, and 940 µm, respectively. They further conducted a sensitivity analysis of these optimum conditions to realise that the yield of bio-oil is more sensitive to an increase in the temperature than to a decrease, as such conditions accelerate the breakdown of heavier components in bio-oil. Having comprehensively validated their CFD-DEM model against results from two rather contrasting experiments, Ku et al. (2019) studied the impact of inert gas temperature on biomass pyrolysis in fluidised beds. Individual simulations were conducted with temperatures of 908, 1098, and 1198 K, while all other parameters influencing the system were fixed. They reported a rapid particle entrainment following faster heating and decomposition of particles, as the temperature was increased.
Similarly, Houston et al. (2022) also observed a reduction of the particle residence time when the inert gas temperature increased from 500 to 600 • C in a fluidised bed operated at three times the minimum fluidisation velocity. They emphasised that this behaviour was prompted by the increase in the devolatilisation rates of the particles whose peak rate doubled. An analysis of the product composition variation was also reported in their work. Accordingly, a monotonic decrease in char from 25.1 to 21.8% by mass in the final product was recorded as the temperature increased from 500 to 600 • C. This indicates that the temperature increase has intervened in the progression of decomposition reactions through the route of char production. The bio-oil, too, decreased by about 3%, while the gas yield increased by 13% following the temperature change. The observations by Houston et al. (2022) show that an increase in the operating temperature may cause the overall oil production to decrease, as the increase in the oil yield from the primary decomposition is outweighed by its secondary decomposition.
Similarly, the simulations of Chen et al. (2020), too, revealed an increase in gas and a decrease in char yields, each by 8% by mass, when the operating temperature in the fluidised bed was raised from 500 to 700 • C. Interestingly, the mass fraction of oil in the final product remained unchanged, implying that the contrasting effects of temperature on the severity of primary and secondary decompositions may have been balanced. The detailed chemical reaction mechanism employed by the authors further permitted an investigation into compositional changes of the bio-oil with operating temperature. Specifically, a temperature increase caused the low-molecular-weight species to be more abundant in the oil phase and the portion of high-molecular-weight species to be smaller. Although elevating the operating temperature accelerates the pyrolysis process through faster decomposition and quicker particle entrainment, Lu et al. (2021d) noticed that the heating of particles across the reactor would become non-uniform if only the temperature is increased. As a result, the routes of decomposition may differ from particle to particle, and, further, lead to undesired retention trends of particles in the reactor. This could be in the form of either prolonged retention of fully converted particles or premature exits of under-utilised feedstock.
Particle size: In order to ascertain the impact of the size of feedstock particles on the biooil yield, Rabinovich et al. (2009) carried out another sensitivity analysis of the established optimum operating conditions, varying the particle size from 400 to 1000 µm. While the change in the bio-oil yield was only minor, corresponding to about 1-2%, the pyrolysis time was found to increase notably. Employing a 1D interfacial model for intra-particle transport, Ku et al. (2019) compared the pyrolysis behaviour of biomass particles with four different sizes, 0.5, 1, 2, and 3mm. Prior to the study, their CFD-DEM model was fully validated against data from two experiments of rather large-sized particles of beech (20 mm) and poplar (9.5 mm) wood. Naturally, with the increase in size, the difference between the core and surface temperatures of the particle followed a rising trend. Furthermore, the time until a particle's core attains its surface temperature increased exponentially, with each duration reading approximately 0.4, 1.1, 3.2, and 9.0 s, respectively. Aware of how local temperature impacts the decomposition paths in pyrolysis, they highlighted the fact that an increase in particle sizes may promote intra-particle heterogeneity and render the pyrolysis rates more difficult to control.
Following an extensive validation, both against single particle and reactor-scale packed bed experiments, Mahmoudi et al. (2016b) conducted CFD-DEM simulations to investigate the influence of the size of beech wood particles on their pyrolysis behaviour. Three identical packed bed reactor setups that differed only in terms of particle sizes were compared. Setups 1 and 2 contained mono-dispersed particles of sizes 10 and 2.5 mm, respectively, while setup 3 encompassed a combination of particles with these sizes in layers. Due to the comparatively large surface area to volume ratio of small particles, the layer of particles that first contacts the flow of hot gas in setup 2 displayed a faster heating, contrary to that in setup 1. The small particles also have a higher surface area in total and larger heat transfer coefficients, compared to the large particles. Hence, the gas flow in setup 2 lost most of its thermal energy rather quickly, leading to a much narrower hot zone within the particle bed. Mahmoudi et al. (2016b) further examined the pyrolysis behaviour of a 2.5 mm particle in comparison with a 10 mm particle, drawn from a location downstream of the former from setup 3. Despite being downstream of the bigger particle, both drying and devolatilisation in the smaller particle started relatively earlier. Furthermore, the compositions of products from setups 1 and 2 affirmed that bigger particles produce more char compared to smaller particles. Chen et al. (2020) also observed a slight increase in the char yield as the size of biomass particles was increased in the range of 450-850 µm, and attributed it to the low heating rates of larger particles. Although the heating and mass loss rates were notably higher with small particles, the ultimate conversion of feedstock remained nearly constant, around 84%, unaffected by the change in particle sizes. In a coarse-grained CFD-DEM study of biomass pyrolysis, Lu et al. (2020) analysed the residence time and conversion distributions of three mono-dispersed particulate systems, of sizes 0.278, 0.344, and 0.426 mm. Unsurprisingly, the average residence time increased monotonically, from 4.07 to 5.62 s, with the size of the particles. In the condition they tested, the time-averaged conversion of particles also exhibited a growing trend with the particle size. The distributions of conversions of particles leaving the reactor uncovered an important concern pertaining to the choice of particle size in pyrolysis. In particular, about 15% of the smallest particles in the reactor leave the fluidised bed with a conversion below 80%. On the contrary, the fractions corresponding to this threshold in the other two systems are relatively insignificant. This behaviour evidences an improper utilisation of feedstock as its size is reduced. However, their study lacks a discussion on the variation of product compositions consequent to changes in feedstock size. In determining an optimum size for pyrolysis feedstock, such information must also be factored in.
Particle shape: Using CFD-DEM, Lu et al. (2021d) simulated the pyrolysis of biomass particles of several shapes, including spheres, cylinders, and one irregular shape. In order to resolve the shapes and accommodate intra-particle distributions of temperature and composition, a glued-sphere approach was implemented. The irregular shape contained a smooth tip at its top, resembling a locality of high surface area to volume ratio compared to the rest of the particle. In comparison to a cylindrical particle of identical volume and surface area, the irregularly shaped particle involved a wider intra-particle temperature profile along with a larger peak temperature owing to the appearance of a hot zone at its tip. With view towards a simplified geometrical representation of irregularly shaped particles, they advocated cylinders that possess the same volumes and surface areas as the irregular target particles instead of volume-equivalent spheres which are most commonly used in practice. Apart from the work of Lu et al. (2021d), Gao et al. (2021b) and Lu et al. (2021c) investigated the pyrolysis of non-spherical particles adopting either cubic shapes, formulated as superquadrics, or an overlapping glued-sphere approach, respectively. Unfortunately, the authors did not report any analysis on the influence of the particle shape on the pyrolysis behaviour.
Inert flow rate: Based on the three packed bed setups briefly described above, Mahmoudi et al. (2016b) also inspected the impact of the gas flow rate on the pyrolysis of beech wood particles. Upon increasing the flow rate from 2.5×10 −3 to 5.0×10 −3 m 3 , a 9% reduction of char yield was observed in the first and third setups, whereas in setup 2, this reduction was more compelling at 24%. From these observations, they concluded that the char yield becomes more sensitive to the changes in the flow rate if the feedstock comprises smaller particles. In a later investigation, Hu et al. (2019) monitored the changes in the composition of products formed during the pyrolysis of biomass as the gas velocity was adjusted to 2.2, 2.6, and 3.0 times the minimum fluidisation velocity. Here, the oil yield increased and the gas yield decreased, while the fraction of char in the pyrolysis product remained unaffected by the flow rate. By reviewing the progression of the multistep chemical reaction model that the authors had adopted, they concluded that the secondary decomposition of oils decreased drastically because of fast removal of produced vapour from the reactor at high flow rates.
For a constant inlet temperature of 550 • C, Houston et al. (2022) varied the inlet velocity of their fluidised bed simulations between two and six times the minimum fluidisation velocity. As the flow rate increased, the gas yield strictly decreased by 2.8%, while the char yield was found to increase by 4.1%. Oil, on the other hand, showed a weak dependency on the flow rate, featuring changes of only about 1%. Moreover, the distributions of the residence times also became more uniform as the inert flow velocity increased. In a similar vain, Lu et al. (2021d) experienced a reduction of the standard deviation of particle temperature in their pyrolysis system when the inlet flow rate of the fluidised bed reactor was increased. This effect was brought about by the enhanced heat transfer to feedstock particles at high ambient flow speeds.
Reactor type and configuration: Lu et al. (2021d) analysed the heating and mass loss rates of irregularly shaped particles that were subjected to pyrolysis in three different reactor systems, a fluidised bed, a spouted bed, and a packed bed. The three reactors were operated at the same average gas velocity of 5 ms −1 and the same temperature of 783 K. The fluidised bed delivered a rather homogeneous heating rate across all particles in the system as reflected by the low standard deviations in their temperature and density profiles. The fixed bed, on the other hand, yielded a higher mean temperature and a lower mean density for a considerable time span compared to the other two reactors. Concomitantly, the standard deviations of temperature and density were larger, suggesting a higher degree of non-uniformity in both the particle heating and reaction rates. The spouted bed offered the lowest particle heating rate, but had an appreciable uniformity, almost comparable to that of the fluidised bed. Wang et al. (2021a) employed a CFD-DEM framework to investigate how the fluidisation behaviour of pyrolysing biomass particles in fluidised beds changes with the shape of the reactor cross-section. To this end, they designed three pyrolysis reactors with similar inlet areas, but circular, square, and rectangular cross-sections. All the reactors were injected with nitrogen gas at 773 K and 0.36 ms −1 . The reactor with the circular cross-section produced the largest mean bed height, while the lowest one occurred with the rectangular cross-section. They also evaluated a bubble equivalent diameter at various heights for each reactor, accounting for the volumes of neighbouring CFD cells with a void fraction of at least 0.8. Nearly up to the top outlet, the circular bed possessed the largest equivalent diameter, while the rectangular reactor had the smallest at all measured heights. Both the average feedstock temperature and the resulting mass loss were slightly higher in the reactor with the circular cross-section compared to the other reactors. The yield of char appeared to be largely unaffected by the choice of the cross-sectional shape. Furthermore, to compare the erosion in differently shaped beds, they evaluated the cumulative erosion volume per unit area in each, adhering to the model proposed by Finnie (1960). After 10 s of operating time, the circular bed had endured the most wear and the square bed the least. Wang et al. (2021a) further analysed the fluidisation behaviour for different gas inlet areas of 25, 50, and 100% in separate circular-bottomed fluidised beds. The bed heights did not show a great dependency on the inlet area, while the frequency at which the heights fluctuated increased as the area of the inlet was reduced. The bed with the largest inlet area produced the lowest mass loss rate as well as the smallest bubble equivalent diameter. The severest degree of erosion was experienced in the reactor with 50% inlet area.

Open Challenges
During the last few years, CFD-DEM has become increasingly popular to predict the behaviour of pyrolysis processes. Different reactor setups have been studied using CFD-DEM, and attempts have been made to improve the level of detail with which pyrolysis systems are described in the scope of CFD-DEM. In terms of particle-side physics, researchers have begun to employ more pragmatic Di Blasi shrinkage schemes instead of the unfounded constant density or constant size shrinkage models. The incorporation of intra-particle models and internal thermal resistance has been demonstrated to impact the pyrolysis outcomes, particularly at high Biot numbers. Recently, a few authors have also examined the pyrolysis of non-spherical particles within the scope of CFD-DEM. Apart from a solitary attempt by Mahmoudi et al. (2016a) to spatially semi-resolve the gas phase in the pyrolysis reactor, the gas phase description remains largely unchanged over all works we reviewed. The interphase interactions have been approximated with many correlations available for the evaluation of the Nusselt number and the drag coefficient. A recent interest has been the use of distinct drag coefficients based on the characteristics of different particles that appear in the system (feedstock or inert bed particles), to closely represent differences in their fluidisation behaviours. The studies which feature analyses of the Nusselt number correlations on the outcome of pyrolysis have shown only a minor influence and have always highlighted the need for more accurate chemical kinetic models for a better prediction. Duly, the chemical kinetic models employed in CFD-DEM models of pyrolysis have become increasingly comprehensive by now. While early models were based on a single global reaction describing the overall transformation, the chemical kinetics now involve up to 32 and 60 reactions representing the primary and secondary decomposition of biomass, respectively.
Although a general framework exists and many chemical kinetic models have been developed, a complete CFD-DEM model describing the pyrolysis of plastics is yet to be formulated and validated. A particular challenge which this task entails is the modelling of the melting of plastic particles and its ensuing effects during the thermal treatment. For a molten feedstock, there are several efforts to model the pyrolysis of plastics (Flor-Barriga and Rodríguez-Zúñiga 2022; Yin et al. 2014) in which the volume-of-fluid method is coupled with CFD to spatially resolve the gas-liquid interface. Zhang et al. (2022a) implemented a CFD-DEM model to ascertain the influence of heat carrier loading and the rolling speed of a rotary furnace pyrolyser on the dynamics and heat transfer rates of polyethylene particles. However, this work hardly counts as an example of plastic pyrolysis modelling with CFD-DEM for two reasons. First, a chemical kinetic model was not included and, second, the melting of particles was omitted, although the system was operated at temperatures as high as 773 K.
In our opinion, one modelling aspect which needs to be improved is the representation of the feedstock particles. This includes not only their shape and size, but also their initial composition which may vary across the feedstock. The shapes of the particles of biomass feedstock often used in conversion processes, such as pyrolysis, differ drastically from spheres and are often closer to cylinders (Chaloupková et al. 2019). As highlighted by Lu et al. (2021d), a shape representation in terms of cylinders that possess the same volumes and surface areas as real particles could lead to improved CFD-DEM predictions. However, for particles of general shape, the geometrical representation becomes more difficult (Lu et al. 2015). Furthermore, the requirement of special collision models and models that approximate the interphase exchange, where the orientation of the particles is also of significant importance, are major challenges. These have been studied in recent contributions on the fundamentals of CFD-DEM for non-spherical particles (e.g., Wachs 2019; Zastawny et al. 2012;Cao and Tafti 2020) that mostly focus on 'regularly' shaped axisymmetric geometries, such as cylinders and ellipsoids. Although the CFD-DEM framework allows for the variation of individual particle properties, the interactions among particles as well as between the particles and the gas are complicated to model, and general, shape-independent interaction models do not exist at present. Moreover, the feedstock used for pyrolysis may contain particles of varying sizes and compositions, which are frequently represented by a mean value to simplify the simulations. Initialising the particles' properties based on experimentally measured size and composition distributions over the feedstock is a challenging but important step to allow for an accurate replication of the process.
The geometrical changes that a particle undergoes during its devolatilisation and the driving forces are yet to be fully understood, let alone to be incorporated into the CFD-DEM framework. In particular, there are experimental observations on the macroscopic behaviour of particle shrinkage, which are not yet reflected in CFD-DEM frameworks applied to pyrolysis reactors. For example, presently, a particle is thought to only shrink and the shrinkage commences at the same time as the mass consumption. Additionally, the shrinkage is also considered to be isotropic, regardless of the directional variations of heat transfer and reaction rates a feedstock particle may exhibit (Gentile et al. 2017). An easy but effective solution could be the calibration of the imposed shrinkage pattern based on a Di Blasi model (Di Blasi 1996), to match the experimental size evolution on the single particle level. This will not only improve the accuracy of the prediction of the particle size, but also all chemical and physical gas-particle interactions that are size-dependent. For a specific feedstock, one may even train the Di Blasi parameters by defining them as functions of variables of direct influence, such as feedstock conversion and heating and reaction rates, for a rather holistic approach. The anisotropy of shrinkage may be accounted for by integrating a 3D intra-particle model, and becomes more relevant if the particles are largely irregularly shaped or exposed to uneven heating rates from different directions. As particles collide, moreover, they may break up or form agglomerates. Although pyrolysis is a particle-centered process, particle breakage and agglomeration seem to have been overlooked and are not included in current CFD-DEM descriptions of pyrolysis reactors. In applications of fluidised beds, such phenomena may directly influence the fluidisation behaviour, in addition to modifying the heating, reaction, and exit rates of particles. In order to capture the effects of particle break-up (Bruchmüller et al. 2011;van Wachem et al. 2020) and agglomeration (Blum 2006), it is possible to enhance the existing CFD-DEM framework by corresponding sub-models.
Some topics that the prevailing chemical reaction models do not specifically address are the existence of heterogeneous tar reactions and the effect of catalytic species. These unaccounted reactions of tar, which are known to be exothermal, have an elevated importance in a predominately endothermal process such as pyrolysis. The heat release from exothermal reactions can promote endothermal reactions, particularly during primary decomposition, while lowering the overall energy demand of the process. Anca-Couce and Scharler (2017) extended the primary decomposition scheme of Debiagi et al. (2018) to account for this potential exothermicity; however, the extended chemical kinetics have not been deployed in a CFD-DEM model of pyrolysis yet. Although pyrolysis systems are operated in the presence of industrial catalysts for the purposes of enhancing the selectivity of desired species and minimising the energy demand, some species and products, such as ash, also possess natural catalytic effects which have been experimentally investigated (Yildiz et al. 2015;Di Blasi et al. 2018). In the absence of reliable chemical kinetic models, a study of a catalytic pyrolysis reactor remains an open challenge to date.

Conclusions
The computer-based modelling of pyrolysis systems permits an evaluation of effects and interactions in and across different scales that are difficult to isolate or observe experimentally. Although CFD-DEM models are computationally more expensive than fully Eulerian modelling approaches, they offer a direct resolution of particle-level phenomena. Within the scope of CFD-DEM, the carrier gas phase is described in a Eulerian fashion, while the dispersed feedstock and inert bed particles are tracked in a Lagrangian framework. A CFD-DEM description encompasses submodels that detail the rates at which mass, momentum, and heat are exchanged between the carrier gas and the particles, characterise the particles' shrinkage and drying behaviour, or treat particle-particle interactions. In the first part of this article, we summarise the existing closure models to describe these phenomena, and discuss their underlying physical assumptions or simplifications.
Subsequently, the article provides an overview of the influence of operating conditions on the pyrolysis behaviour, as reported in works that sought to optimise the process with the aid of CFD-DEM models. The optimisation of a pyrolyser is an attempt at maximising the degree of feedstock conversion and the yield of pyrolysis-oil, while minimising the residence time. In addition to accelerating the process, a high temperature is usually preferred for the operation of a pyrolyser, as it increases the selectivity of the pyrolysis-oil during the primary decomposition. On the downside, the secondary decomposition of the produced oils is promoted. Further, an incomplete utilisation of feedstock may also result through the inhomogeneous distribution of temperature across feedstock particles. On the other hand, increasing the gas flow rate through the pyrolysis reactor has been seen to homogenise both interparticle temperature distributions and particle residence times. Yet, a high flow rate can adversely affect the feedstock conversion as the particles may leave the reactor prematurely, particularly for excessively small sizes. Conversely, large-sized particles are not favoured in pyrolysis because they may feature inhomogeneous intra-particle temperature and composition distributions and lead to a larger production of char.
In the present article, we also discuss open challenges in the modelling of pyrolysis with CFD-DEM. Firstly, the extension of the CFD-DEM methodology to arbitrarily-shaped particles and size or shape distributions is a topic of fundamental CFD-DEM research and recent developments here are only beginning to be applied to pyrolysis processes. Secondly, almost all CFD-DEM studies in pyrolysis have been performed with biomass as feedstock and there are currently no CFD-DEM descriptions for the pyrolysis of plastic feedstock, although this is an increasingly important process. Thirdly, the feedstock particles change as they undergo pyrolysis; they typically shrink, but can also expand, break, melt or fuse. Although 'simple' shrinkage models do exist in the framework of CFD-DEM, such models are currently elementary, and do not treat the other changes in particle shape and size as the particles undergo pyrolysis. Finally, the complexity of the chemical reaction models currently applied in CFD-DEM descriptions of pyrolysis processes is limited, and further developments are required to include exothermicity and catalytic effects.