CET vol 99


 
 
 
 
 
 
 
 
 
 
                                                                                                                                                                 DOI: 10.3303/CET2399088 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Paper Received: 21 November 2022; Revised: 25 March 2023; Accepted: 2 May 2023 
Please cite this article as: Cutillo E.A., Bizon K., Continillo G., 2023, Assessing the Effect of Variable Ambient Temperature on the Self-ignition 
of a Reaction-diffusion System Employing a Reduced Order Modelling Methodology, Chemical Engineering Transactions, 99, 523-528  
DOI:10.3303/CET2399088 
  

 CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL.  99, 2023 

A publication of 

 

The Italian Association 

of Chemical Engineering 

Online at www.cetjournal.it 

Guest Editors: Sauro Pierucci, Flavio Manenti 

Copyright © 2023, AIDIC Servizi S.r.l. 

ISBN 978-88-95608-98-3; ISSN 2283-9216 

Assessing the Effect of Variable Ambient Temperature on The 

Self-Ignition of a Reaction-Diffusion System Employing a 

Reduced Order Modelling Methodology 

Enrico A. Cutilloa*, Katarzyna Bizonb, Gaetano Continilloa  

aDipartimento di Ingegneria, Univeristà degli studi del Sannio, Piazza Roma 21, Benevento, Italia 
bFaculty of Chemical Engineering and TechnologyCracow University of Technology, Kraków, Poland 

cutillo@unisannio.com 

The system under study in this work is a self-igniting pile of solid material. To predict and understand the effect 

of steep changes of the state variables on such systems, a reaction-diffusion model is employed. These systems 

can exhibit complex oscillatory behaviour, and changes in ambient conditions over time may strongly impact the 

inherent oscillations. To simulate the unsteady evolution of the pile, both a classical numerical technique 

(method of lines) and a reduced order approach are employed in combination with a stiff ODE solver. To account 

for circadian fluctuations in temperature, time-variable boundary conditions are assumed upon formulating the 

problem. The reduced order model is introduced in view of understanding if an approximated formulation 

characterized by a much lower number of state variables can accurately predict the complex behaviour of the 

system even in the case of sudden, steep variations of the values of the state variables due to the phenomenon 

of self-ignition, intensified here by variable boundary conditions. The selected case studies have the goal of 

exploring the effect of stockpile properties on the self-ignition phenomenon. Numerical solutions show the 

anticipated coupling between the system intrinsic dynamics and the oscillating temperature imposed at the 

boundary. All of the analysed cases are accurately replicated by the reduced order model. 

1. Introduction 

Spontaneous ignition is a phenomenon whose danger must be carefully foreseen with a view to its prevention, 

to reduce the risk of material and economic losses in the storage of porous flammable materials. Today, storage 

of solid materials can be relevant to many industries. Traditional production processes such as energy 

production from coal, biobased processes where fossil sources are replaced by biomass (Thorenz et al., 2018), 

or waste treatment plants (Zambra et al., 2011), share the need of handling and storing porous flammable 

materials. The problem of spontaneous combustion and even explosion can also arise in the agricultural or food 

industry (Ramírezn et al., 2009). Non-linear dynamical systems are often employed to study and anticipate the 

occurrence of such unwanted phenomena via numerical simulations (Fu et al. 2021). In this study, the reaction-

diffusion model adopted in Continillo et al. (2000), is used to investigate the effects of ambient temperature 

variations. To describe the spatial and temporal evolution of these initial-boundary value problem, finite 

difference approximations of spatial derivatives may be utilized, which feature a high number of discretization 

nodes, that are required to maintain a high level of accuracy even in the presence of steep changes in the state 

variables, either spatially or temporally, as occurs in the self-combustion of a solid stack (Cutillo et al., 2023). In 

this context, Reduced Order Models (ROMs) are paramount to enable sound design and optimization activities 

of such high dimensional and computationally expensive systems (Jain et al., 2017). The methodology of model 

order reduction has been shown to be both fast and accurate in a variety of applications, ranging from 

aerodynamics (Krath et al., 2021) to chemical engineering problems (Bizon and Continillo, 2021). In this work it 

is shown how knowing non-trivial information about the system, such as the influence of the dimensionless 

activation energy γ upon the dynamics of the system, can help in deriving reduced models that are valid even 

outside of the parameter value domains for which they are built. In fact, if properly constructed, ROMs that are 

based on original physical models can predict system behaviour in areas of parameter space that were not 

523



explored during the model's development. In this study, a ROM based on Proper Orthogonal Decomposition 

(POD) and Galerkin projection is developed and utilized to simulate the impact of temperature changes on self-

ignition phenomena. 

2. Mathematical model and computational methods 

2.1 Mathematical model of a self-igniting pile 

The problem of self-ignition of a stockpile is formulated under the following assumptions (Continillo et al., 2000): 

The lower side of the pile is in contact with the ground, while the upper side is exposed to ambient conditions. 

Given that the physical and geometrical properties of the pile in the horizontal directions are assumed constant, 

a one-dimensional model is appropriate to describe the problem. To account for reactant transport inside the 

porous medium, Fickian diffusion of the gas is considered. The chemical reaction between the reactant and the 

solid fuel is modelled using a first-order kinetic. The temperature dependence of the reaction rate is described 

using the classical Arrhenius exponential expression. Gas- and solid-phase are assumed to be in local thermal 

equilibrium, resulting in a single energy balance equation. In addition, neglecting solid material consumption 

leads to the use of only one mass balance equation. Based on these assumptions, the problem is expressed 

using the following dimensionless model equations: 

   
 

   

2
2

2

c c γ
=Le -Φ cexp -

t x T
 (1) 

2
2

2

T T γ
βΦ exp

t x T
c

   
= +  − 

   
 (2) 

The associated boundary and initial conditions are: 

( ) ( ) ( ) ( )
x 1

c
T 0,t A sin ω , T 1,t 1, c 0,t 1, 0 for t 0

x
t

=


=    =     =    =      


 

 
(3) 

( ) ( )  ss ssT x,0 =T ,   c x,0 =c  for x 0,1         (4) 

Following the work of Zambra et al. (2011), who analysed a pile of composting material under variable ambient 
conditions, the naturally occurring temperature oscillations throughout the day and night are formulated in this 
work using a sinusoidal boundary condition, as reported in Equation 3. The above system of partial differential 
equations is solved applying the finite difference method. Particularly, the discretization nodes employed in this 
study are N=200. The resulting system of ordinary differential equations (ODEs) is called Full Order Model 
(FOM) and is solved with a stiff ODE solver available in MATLAB (ODE15s). 

2.2 Reduced order model formulation 

The methodology employed in this work to build the reduced order model is known as Proper Orthogonal 

Decomposition (Stanislav et al. 2012) and can be seen as a method to extract the most important information 

about the dynamics of the system and to use that information to build a set of empirical basis functions that 

allow an efficient approximation of the state variables of the system. This method is based on the collection of 

experimental observations that in the case of this work come from numerical simulations performed with the 

FOM and will be referred to as snapshots. In general, it is recommendable that the experimental or simulation 

data coming from the FOM be chosen as representative of all different observed dynamic behaviours.  

In this work, the set of snapshots chosen to build the ROM are sampled using a uniform sampling strategy. The 

snapshots are then collected in the following matrix: y(x,ti), x=[x1, x2, …,xn], i = 1, …, M (note that, for simplicity, 

the dependencies of y will not be indicated in the following). To derive the set of empirical basis functions, the 

following eigenvalue problem is solved (Holmes et al., 2012): 

   Φ=C λC  (5) 

where:  

TC =  y y  (6) 

is the correlation matrix and   denotes an arithmetic mean of the standard inner product of the snapshot 

matrix. Then, the truncated state variables are expressed as a linear combination of a sub-set of first K 

determined basis functions: 

524



( ) ( )
1

Φ
K

j i j

j

a t
=

= y x  (7) 

The POD method allows to extract N basis function, but only the most relevant are employed to build the reduced 

model. In fact, the basis function associated with the highest eigenvalues carries most of the information about 

the system dynamics. For this reason, it is possible to approximate the state variable with only K << N POD 

functions, leading to the desired reduction in the dimension of the original system. Following a previous work 

(Cutillo et al., 2023), the set of snapshots was sampled for values of the parameter γ equal to 11 and 13. These 

two values are located in the vicinity of a region of the parameter space where the system has chaotic behaviour, 

as earlier reported in Continillo et al. (2000). The choice of sampling the FOM solutions for these two parameter 

values stems from the fact that two fundamental pieces of information can be derived from the evolution of the 

state variables. The first is the characteristic behaviour of an ignition event with a subsequent stationary state, 

the second is the characteristic behaviour of an ignition event leading to a stable limit cycle. After sampling, the 

reduced model was built and, particularly, 11 modes were chosen for the approximation of the state variables. 

In this way, the original system with 400 differential equations was reduced to 22 differential equations. Like with 

the FOM, the reduced model was also solved using the method of lines, through the same stiff integrator in 

MATLAB (ode15s). 

3. Results and discussion 

3.1 Key parameter values 

Although the model employed in this work was originally derived for simulating coal stockpiles, its crucial kinetic 

parameters can be adjusted to forecast the self-ignition of a variety of materials used across multiple industries, 

including biomass. It is important to consider the potential implications of changing these parameters before 

using the model for new materials, for this reason the selection of parameter values used in the simulation was 

made with the goal of representing realistic materials and ambient conditions. The parameter values reported 

in Table 1 are typical of lignin-rich biomass such as the one studied in Bach et al. (2017). Lignocellulosic 

biomass, in particular, is gaining attention as it is one of the key components of second-generation biofuels (De 

Blasio et al., 2022). Meanwhile, the values used in this study for the combustion enthalpy were taken from the 

research conducted by Krigstin and Wetzel (2016). The frequency of the oscillating temperature at the top 

boundary has been chosen to represent circadian oscillation; being a dimensionless variable, its value depends 

on L (the pile height) and its value can be computed as ω=2𝜋L2 (86400 ∙ 𝛼)⁄ . For the other dimensionless 

variables, reference can be made to Continillo et al. (2000). 

Table 1: Dimensionless and dimensional parameters employed 

Parameter Value Parameter Value 

Φ2 [0.45, 6.5]·105 D 5·10-5 

β  [4.5-4.67] α 1.8·10-4 

γ  [13-17] ω 0.396 

Le 0.2722 ε 0.3 

ΔH [1.8, 2.0]·105 L [0.8, 3] 

Ea [3.25, 4.25]·105 A 20 

3.2 Simulation results 

The studied system can exhibit complex behaviour even without external variable forcing. In fact, Figure 1 shows 

the evolution of the system when a constant temperature value is applied to the outer boundary. The 

dimensionless parameter values used in the analyzed simulation are: φ2 = 5.9529·105, β = 4.5, γ= 17.0386, Le 

= 0.27, which are compatible with the following dimensional parameter values: ΔH = 1.88·105, Ea = 4.25·104, L 

= 2.9 , while the other values are reported in Table 1. As it can be observed, after the transient regime has 

extinguished, the system starts oscillating around what is called a chaotic attractor. This regime is characterized 

by aperiodic behaviour that is strongly dependent on the initial conditions. Although this behaviour has been 

already predicted in Continillo et al. (2000), in the cited work the simulations were carried out for different values 

of the parameters, that were representative of a coal stockpile, whereas in this work are characteristic of a 

stockpile of woody biomass. Because of the hypothesis made in deriving the model, such as: the consumption 

of the solid fuel is neglected; the properties of the pile are constant during the time of the simulation, while the 

concentration of oxidant is constant at the top of the domain; it is seen that the consumption of reactant caused 

by combustion is followed by a new build-up, and so on. This phenomenon can be found in the literature under 

525



the name of relaxation-oscillation regime, in which extremely slow build-up is followed by a sudden discharge 

(Strogatz, 2018). In systems characterized by relaxation-oscillation phenomena, the boundary conditions define 

the range of possible solutions that can be challenging in terms of numerical stability, and this of course can 

affect the applied method of the Proper Orthogonal Decomposition for the construction of the Reduced Order 

Model. As seen in the following results, a change in boundary conditions can also affect the amplitude and 

frequency of the oscillations. When the boundary conditions are changed, the system may switch from an 

oscillatory to a non-oscillatory state or can switch between oscillatory states characterized by different 

amplitudes and frequency. For this reason, time-dependent boundary conditions have a significant impact on 

the behaviour of dynamic systems that exhibit relaxation-oscillation phenomena.  

 

Figure 1: Chaotic orbit (a) and time series (b) representative of the evolution of the concentration and the 

temperature inside the porous matrix. The data is sampled at a point in the middle of the domain (x = 0.5).  

Figure 2 shows that, after forcing the system with sinusoidal boundary conditions, its dynamic behaviour 

switches from chaotic to more regular periodic oscillations, featuring much larger bursts of peak temperature, 

resulting in more violent ignition phenomena, as opposed to the low amplitude oscillations that can be found 

without the periodic forcing, shown earlier in Figure 1b. Figure 2a shows how the “stresses” accumulated during 

the slow build-up are released during the sudden “discharge”. Here, the “stress” is the concentration of reactant 

that accumulates within the porous matrix, resulting in a sudden “discharge” when the explosion happens. Figure 

2a shows that the build-up takes almost 3 days. The occurring phenomenon can be visualized better by looking 

at Figure 2b, where snaps of the spatiotemporal evolution of temperature and concentration are reported.  

 

Figure 2: Temperature evolution at the middle of the domain before and after oscillating boundary condition 

are considered, mimicking circadian temperature variations (a), and state variables profiles along the domain 

in four different salient moments. 

The first snapshot is taken at 567.5 h (or 23.64 days) from the beginning of time. Here the system has already 

gone through an ignition phenomenon (as it can be seen from Figure 2a): both temperature and concentration 

are decreasing inside the domain. This decrease is followed by a new, slow build up, as can be seen from the 

snap at 605.7 h (25.24 days). After this state, it takes almost 30 h for the system to achieve the right combination 

of temperature and concentration that allow for the combustion event to newly occur. At 635.3 h the combustion 

starts at the middle of the domain, and in less than 0.1 h (6 minutes), the front of the flame has spread throughout 

the whole domain, consuming the accumulated oxygen. Figure 2a reports the evolution of the system before 

and after the sinusoidal forcing is applied. After the forcing, the system is simulated both by the FOM and the 

526



ROM. In this case study, as emphasized by the overlapping curves, the ROM appropriately reproduces the 

dynamic behaviour analysed, which is remarkable if one thinks that this data-driven model has been built from 

solutions of the FOM derived with the hypothesis of constant boundary condition. Clearly, this accuracy arises 

from the fact that the Reduced Order Model (built upon the original first principles model solutions), maintains 

key information about the system dynamics.  

 
Figure 3:  Time series representative of the evolution of the state variables inside the porous matrix. Comparison 

of FOM results with the time-varying boundary condition (a) and FOM with ROM (b). The data is sampled in a 

point at the middle of the domain (x = 0.5). 

While keeping the oscillating boundary temperature, the behaviour of the system can drastically change as the 

values of the parameters of interest are varied. In fact, the results of the simulation reported in Figure 3 are 

obtained for the following values of the parameters: φ2 = 4.5302·104, γ= 12.9894, which are compatible with the 

following dimensional parameter values: ΔH = 1.88·105, Ea = 3.24·104, L = 0.8, while maintaining constant the 

other parameters. Figure 3 shows that, in this case, self-ignition phenomena happen more frequently: taking as 

reference the oscillating boundary condition (reported as the dashed curve in Figure 3a), during one period of 

the sinusoidal boundary condition there are between six and seven combustion events which are followed by 

low amplitude and high frequency temperature and concentration oscillations. Moreover, the most energetic 

combustion events take place when the boundary temperature reaches its minimum. This is caused by the fact 

that, in this case, the concentration of oxidiser in the first layers of the pile cannot drop below the critical level of 

self-ignition because the lower temperature of the “night” inhibits the consumption of reactant in the first layer of 

the domain; hence the reactant  is capable of filling the pile and self-ignite. During ignition, the combustion fronts 

travel along the domain. At the top of the domain the pile is already at low temperature, while in the bottom of 

the domain the pile is still violently burning. The concentration of the reactant keeps building up in the first layers 

and, as the combustion event is ended, a new one can start. Thus, multiple ignitions are observed during one 

night. As the temperature rises (during the day) the concentration of oxygen no longer has the possibility to 

accumulate inside the pile because it is consumed in the first layers and, for this reason, impulsive combustion 

cannot be observed. With these parameter values the ROM appears to be less accurate than the previous case, 

and in fact the curves are not overlapped anymore but, as clearly shown in Figure 3b, the behaviour of the 

system is still well depicted in terms of frequency and amplitude of the oscillations of the state variables. The 

error between FOM and ROM is due to a phase shift, that is characteristic of ROMs built with POD method 

when they are used for the simulation of periodic regimes (Bizon et al., 2008). 

4. Conclusions 

The size chosen for the system under study does not require high computational costs for its numerical 

resolution, and this allowed to evaluate the performance of the ROMs under various conditions. Particularly, the 

reduced model was formulated in view of understanding what would happen if time-dependent boundary 

conditions were simulated, mimicking circadian phenomena, on top of relaxation-oscillation phenomena 

occurring with constant boundary conditions. This in some cases enhanced extreme ignition phenomena that 

might affect numerical stability: in this respect, the MATLAB ODE solver, in association with the model reduction 

technique developed, proved robust and reliable. The system dynamics was shown to be clearly influenced by 

varying boundary conditions, and several parameter values, chosen as characteristic of second-generation 

biofuels, were investigated, and the system behaviour was described quantitatively (through the simulations) 

and physically (giving a physical explanation of the observed phenomenon). Particularly, the relaxation-

oscillation behaviour has been found and described in the neighbourhood of chaotic solutions. Simulations 

results show that the non-stationary character of burning piles (i.e. braziers, stored biomasses, charcoal) already 

527



emerges from a simple reaction-diffusion model that neglects most of the real system complexity. Thanks to the 

ROM, it was possible to accurately reproduce the results of the classical model applied in the literature, with a 

reduction of the systemic order of magnitude (going from 400 equations to 22) with obvious savings in CPU 

time. The work will serve as a basis for extending the analysis to more complex systems, such as 2D or 3D 

domains, which may represent the fuel stack more accurately. It is expected that, in addition to the reduction in 

the system, significant savings in CPU time will allow the model to be applied to solve computationally onerous 

problems, such as parametric analysis, optimization or control studies, with a view to predicting and reducing 

the risks these systems pose in terms of economic or material losses.  

Nomenclature

c – dimensionless reactant concentration,  - 

T – dimensionless temperature,  - 

φ – thermal Thiele modulus,  -  

Le – Lewis number,  - 

β – dimensionless heat reaction,  - 

γ – dimensionless activation energy,  - 

t – dimensionless time,  - 

L – layer thickness, m 

D – mass diffusivity, m2/s 

α – thermal diffusivity, m2/s 

ΔH – enthalpy of reaction, J/mol 

Ea – activation energy, J/mol 

ε – porosity,  - 

cp – specific heat, kJ/(kg·K) 

A – amplitude of circadian oscillations,  - 

ω – angular velocity of circadian oscillations,  - 

y – snapshots matrix, - 

x – discretized coordinate vector, - 

M – number of time samples, - 

Φ – set of basis functions, - 

aj – modal coefficients, - 

λ – eigenvectors associated with basis functions, - 

N – number of discretization nodes, - 

K – truncation order of the state variables, - 

C – correlation matrix,  - 

References 

Bach Q.V., Tran K.Q. and Skreiberg Ø., 2017, Combustion kinetics of wet-torrefied forest residues using the 

distributed activation energy model (DAEM), Applied Energy, Vol. 185, pp.1059-1066. 

Bizon K., and Continillo G., Russo, L., Smuła, J., 2008, On POD reduced models of tubular reactor with periodic 

regimes. Computers & Chemical Engineering, Vol. 32, pp. 1305-1315. 

Bizon K., and Continillo G., 2021, Efficient optimization of a multifunctional catalytic fixed-bed reactor via 

reduced-order modeling approach. Chemical Engineering Research and Design, Vol. 165, pp.214-229. 

Continillo G., Faraoni V., Maffettone P.L., Crescitelli S., 2000, Non-linear dynamics of a self-igniting reaction-

diffusion system, Chemical Engineering Science, Vol 55, 303-309. 

Cutillo E. A., Petito G., Bizon K., Continillo G., 2023, Analysis of an innovative sampling strategy based on k-

means clustering algorithm for POD and POD-DEIM reduced order models of a 2-D reaction-diffusion 

system, Combustion Theory and Modelling, 1-28. 

De Blasio C., Salierno, G. and Özdenkci, K., 2022, Lignocellulosic Biomass for Biofuels Production, an 

Integrated Approach, Enzymes in the Valorization of Waste, pp. 23-49. 

Fu S.C., Haixiang W., Sidhu S. D., Luangwilai H. S., Shu Y., 2013, Numerical Study on Effect of Ambient 

Humidity Variation on Self-heating and Spontaneous Ignition of the Eucalyptus Bark Pile, Fire Technology, 

Vol 57, 1803-1825. 

Holmes P., Lumley J.L., Berkooz G. and Rowley C.W., 2012, Turbulence, coherent structures, dynamical 

systems and symmetry. Cambridge university press. 

Jain S., Tiso P., Rutzmoser J.B., Rixen D.J., 2017, A quadratic manifold for model order reduction of nonlinear 

structural dynamics. Computers & Structures, Vol. 188, pp.80-94. 

Krath E.H., Houchens B., Marian, D. Pol, S. and Westergaard C., 2021, Multivariate Design and Optimization 

of AeroMINE Internal Turbine Blade. In AIAA Propulsion and Energy 2021 Forum. 

Krigstin S., and Wetzel S., 2016, A review of mechanisms responsible for changes to stored woody biomass 

fuels, Fuel, Vol. 175, pp.75-86. 

Ramírez A., García-Torrent J., Aguado P.J., 2009, Determination of parameters used to prevent ignition of 

stored materials and to protect against explosions in food industries, Journal of Hazardous Materials, Vol 

168, 115-120. 

Strogatz S.H., 2018, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and 

engineering, CRC press. 

Thorenz A., Wietschel L., Stindt D., Tuma A., 2018 , Assessment of agroforestry residue potentials for the 

bioeconomy in the European Union, Journal of Cleaner Production, Vol 176, 348-359. 

Zambra C.E., Moraga N.O., Escudey M., 2011, Heat and mass transfer in unsaturated porous media: Moisture 

effects in compost piles self-heating, International Journal of Heat and Mass Transfer, Vol 54, 2801-2810. 

528


	383cutillo.pdf
	Assessing the Effect of Variable Ambient Temperature on The Self-Ignition of a Reaction-Diffusion System Employing a Reduced Order Modelling Methodology