Parametric Sensitivity Analysis for the Industrial Case of O-Xylene Oxidation to Phthalic Anhydride in a Packed Bed Catalytic Reactor

Abstract

The operation of packed bed tubular reactors, for exothermic catalytic reactions, presents special challenges provoked by hotspot development. Their potential safety risk can be assessed using different mathematical methodologies, among them, the so-called parametric sensitivity analysis (PSA). This study deals with the identification of safe operational conditions (e.g., feed temperature conditions) for the catalytic oxidation of o-xylene. Three different reaction networks, with different degrees of complexity, were analyzed. Thus, the critical values of the operating parameters, allowing us to define run-away and stable operation conditions, were provided for different reactive configurations. The obtained results were compared with the data reported by various authors who used similar reaction rate laws. The purpose of the present study is to illustrate the peculiarities of the PSA and its application for the design, analysis, and operation of o-xylene multitubular catalytic reactor.

1. Introduction

It is with great respect and profound admiration that we dedicate this work to Professor Hugo de Lasa on the occasion of honoring and recognizing his outstanding career contributions to the fields of heterogeneous catalysis, photocatalysis and catalytic reaction engineering. In the subject of modelling and analysis of catalytic chemical reactors, it is not unjust to claim a starring place for Professor Hugo de Lasa. Indeed, from the intelligent use of mathematical models during his early contributions (e.g., Soria Lopez and de Lasa, 1981 [1]; de Lasa, 1983 [2]; Arandes and de Lasa, 1995 [3]), Professor de Lasa exemplified his remarkable modeling and analytical skills. Professor de Lasa has set inspiration and charted the course between the oversimplification and clouding details for parametric sensitivity analysis of fixed bed catalytic reactors. So far, his intellect, wisdom and research contributions continue to influence us deeply.

Packed Bed Tubular Reactors (PBRs) are considered the workhorses of the chemical and petrochemical industries. Indeed, they are usually the first choice for large-scale production. PBRs, in their multi-tubular scheme, are preferred for highly exothermic reactions, such as the partial oxidation of different hydrocarbons [4]. Any effort to quantify the phenomena taking place in a PBR must consider the mass and heat transport between the reactive fluid and the catalyst, the pressure drop for flow through the packed bed, and the energy balance on the heat-exchange fluid. The inherent two-way coupling between concentration and temperature causes chemical reactors to show exceptional behaviors, such as a slight change in one or more inlet parameter, which can enormously alter the reactor effluent conditions. Bilous and Amundson [5] called this phenomenon parametric sensitivity. When a PBR operates in the parametrically sensitive region, its performance becomes unreliable. For instance, a slight increase in the inlet reactive-fluid temperature can lead to a large rise in the reaction rate (due to the exponential temperature dependence of the rate constant). The increased reaction rate results in more heat generation which furthermore increases the reaction temperature. This feedback mechanism results in the so-called reactor runaway (when the rate of heat generation far exceeds the rate of heat removal). Such increase in the reaction temperature can damage the reactor, can generate safety hazards, can cause catalyst deactivation, and can propitiate undesirable side reactions [6]. Thus, it would be of great value to reactor designers and operators to be able to predict parametric sensitivity conditions.

The solution of the rigorous model describing non-isothermal PBRs has proven to be extremely challenging. Additionally, the level of complexity increases when including mathematical tools for parametric sensitivity analysis. Fortunately, recent advances in computational resources have enabled their simultaneous simulation. During the last three decades, the use of the following four methodologies for parametric sensitivity analysis has become widespread: (i) Temperature-partial pressure phase plane [7]; (ii) The sensitivity indices method [8]; (iii) The divergence methods, based on chaos theory [9] and (iv) The methods of trajectory extension [10]. This rough classification implies that, although the runaway phenomenon is well known, a unique analysis approach does not exist. One of the purposes of the present work is to compare the results obtained using the temperature-partial pressure phase plane (e.g., by de Lasa, [2]) vs. those obtained using the sensitivity indices method for the case of the catalytic (with V₂O₅) oxidation of ortho-xylene (o-xylene) for producing phthalic anhydride. This is a classic example of a fast oxidation reaction which, if not controlled, produces CO and CO₂ (CO_x) instead of valuable products. It has been pointed out in the open literature that a multitubular PBR is indispensable and that procedures for obtaining nearly complete conversion without causing a runaway must be guaranteed [11]. Industrially, the reactive process proceeds at atmospheric pressure, in the temperature range of 360–400 °C, reaching almost complete o-xylene conversion and with a phthalic anhydride selectivity between 70–75% [12]. Ortho-tolualdehyde (o-tolualdehyde) and phthalide are the main intermediates. However, depending on operation conditions, maleic anhydride can also be generated (Figure 1). It is evident that the o-xylene oxidation process is highly complex.

Different approaches and assumptions, proposed by several researchers in the literature, yielded some models for o-xylene disappearance. The only probable agreement between the different authors is that the reactions occur through a “redox” mechanism.

Table 1. Reaction steps considered by different authors for kinetic studies or reactor modelling.

Mechanism	Reference	Reaction Step
		1	2	3	4	5	6	7	8	9	10	11	12
1	[7]						X
2	[2,16]					X	X					X
3	[17]	X	X	X			X					X
4	[12]	X	X	X						X		X
5	[14]	X	X	X		X	X					X
6	[18]	X	X	X	X				X			X
7	[19]	X	X	X			X	X		X		X
8	[20]	X	X	X		X	X	X	X	X	X	X
9	[21]	X	X	X		X	X	X	X	X	X	X	X

X: Reaction step from Figure 1.

includes some of the reaction steps considered by different investigators for kinetic studies or reactor modelling. For practical purposes, more than just the reaction rate of o-xylene to phthalic anhydride (step 6 in Figure 1) is needed. To model the highly temperature-sensitive system, the competition between phthalic anhydride and complete oxidation must be followed at least, which causes large heat effects [14,15].

To illustrate the effect of o-xylene oxidation mechanism on reactor performance prediction, three different models were considered in Figure 2 (data used for simulation will be presented later on). Regardless of the difference of the kinetic model used for simulation, a common behavior of the temperature profile occurs: a hotspot is predicted at some position in the PBR. However, the magnitude of the hotspot depends on the reaction kinetics. It is also interesting to note that, for the conditions used in these calculations, the predicted o-xylene conversions are significantly different.

The influence of the PBR parameters (parametric sensitivity), such as operation conditions and reaction kinetics on o-xylene oxidation, with constant physicochemical properties, was analyzed using the temperature-partial pressure phase plane by several authors as follows:

• van Welsenaere and Froment, (1970) [7] introduced, for an ideal, one-dimensional PBR model, with constant wall temperature, two criteria for runaway based on an intrinsic property of the reactive system. They transposed the peak temperature and the conditions of the inflexion points for the prediction of the critical values for the operating variables. Their work also treats the related problem of a hotspot which has to be limited for reasons other than runaway.
• Soria-López et al. (1981) [1], as an extension of van Welsenaere and Froment’s [7] work, analyzed an ideal, one-dimensional PBR model, with co-current external cooling, under two types of operation modes: (i) the temperature along the reactor axis increases monotonically; that is, without a hotspot can occur under certain conditions (referred to as pseudo-adiabatic operations, PAO), and (ii) the temperature shows a maximum at a finite axial reactor position (MFARP). Limiting operation conditions were defined between MFARP temperature curves and PAO curves.
• de Lasa (1983) [2], continuing Soria-López et al.’s work [1], studied the peculiarities of PAO and its significance for the design and operational conditions of o-xylene oxidation in multi-tubular reactors using reaction networks models of different complexity.
• Akella and Lee (1983) [22] studied, for an ideal, one-dimensional PBR model, a counter-current operation. The runaway conditions were derived from a phase plane of feed and coolant inlet temperatures.
• Hosten and Froment (1986) [23] extended the van Welsenaere and Froment (1970) [7] approach to co-currently cooled PBRs with wall temperature variation with reactor length.
• Later, Arandes and de Lasa (1995) [3] extended the PAO analysis to the following operation regimens in PBRs: isothermal, decreasing temperature profile, cold spot, hotspot–cold spot, and cold spot–hotspot.

From another perspective of analysis, Henning and Perez (1986) [16] defined, for o-xylene oxidation reaction, a criterion for runaway in PBRs based on the behavior of the sensitivity indices along the reactor, when the temperature variation in a co-current cooling medium is considered. Their analysis was not restricted to simple kinetics.

It is interesting to observe that all of the above-mentioned works neglected the pressure drop for flow through a packed bed. However, in gas phase reactions, the concentration of the reacting species is proportional to the total pressure. Therefore, a proper account of the effect of pressure variation on the PBR is a key factor in the success of reactor operation.

This work deals with the definition of suitable inlet reactor conditions (e.g., temperature and pressure leading to safe o-xylene oxidation reactor operation) applying a sensitivity indices method to the analysis of a non-isothermal PBR model, which includes: mass and heat balances inside (packed bed) and outside (shell side: in co-current or countercurrent flow) the reactor tubes and the pressure drop balance (Ergun equation). For comparison, three different reaction networks, with different degrees of complexity, are analyzed: (i) a pseudo-first-order kinetic model, used by van Welsenaere and Froment [7], which includes a single reaction (step 6 in Figure 1); (ii) a three-step reaction model (steps 5, 6 and 11 in Figure 1), as used by de Lasa [2], and (iii) the Chandrasekharan and Calderbank [17] kinetic model involving five steps of reaction (steps 1, 2, 3, 6 and 11 in Figure 1). Thus, the previously reported a priori runaway conditions for hotspot operation are compared with the new obtained criteria. Additionally, the boundary between safe and unsafe operation regions will be established.

2. Reactor Model

Previous work [24,25] indicated that a one-dimensional, pseudo-homogeneous, PBR model can be considered adequate for the analysis for o-xylene oxidation. Thus, the model used in this work is based on the following assumptions:

• The PBR is operating at steady state conditions.
• Molar flows, pressure and temperature gradients only occur in the axial direction. The only transport mechanism operating in this direction is the overall flow itself, and this is considered to be of the plug flow type.
• The parameters M, u, ρ_g, ${\mathrm{C}}_{{\mathrm{p}}_{\mathrm{g}}}$ ρ_c, ${\mathrm{C}}_{{\mathrm{p}}_{\mathrm{c}}}$ w_c, U, G, Re_p are assumed to be independent of temperature.
• There are neither temperature nor concentration differences between the fluid and the catalyst particle.

Thus, the mass, heat, and pressure balances can be written as follows [26]:

$$\begin{array}{}\mathrm{(1)}& \mathrm{Mass} \mathrm{balance}& \frac{\mathrm{d}{\mathrm{F}}_{\mathrm{j}}}{\mathrm{d}\mathrm{z}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{r}}_{\mathrm{i}\mathrm{j}}}\hfill \mathrm{(2)}& \begin{array}{c}\mathrm{Heat} \mathrm{balance}\\ \mathrm{inside} \mathrm{the} \mathrm{reactor}\end{array}& \frac{\mathrm{d}\mathrm{T}}{\mathrm{d}\mathrm{z}}=\frac{{\mathsf{\rho }}_{\mathrm{b}}}{\mathrm{u}{\mathsf{\rho }}_{\mathrm{g}}{\mathrm{C}}_{{\mathrm{p}}_{\mathrm{g}}}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{q}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}\left(-{\mathrm{r}}_{\mathrm{i}\mathrm{j}}\right)\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{\mathrm{i}}}}}-\frac{4\mathrm{U}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{C}}\right)}{\mathrm{u}{\mathsf{\rho }}_{\mathrm{g}}{\mathrm{C}}_{{\mathrm{p}}_{\mathrm{g}}}{\mathrm{D}}_{\mathrm{t}}}\hfill \mathrm{(3)}& \begin{array}{c}\mathrm{Pressure} \mathrm{drop}\\ \mathrm{equation}\end{array}& \frac{\mathrm{d}\mathrm{P}}{\mathrm{d}\mathrm{z}}=-\frac{{\mathrm{G}}^2\left(1-\mathsf{\epsilon }\right)}{{\mathsf{\rho }}_{\mathrm{g}}{\mathrm{D}}_{\mathrm{p}}{\mathsf{\epsilon }}^3}\left(\frac{150\left(1-\mathsf{\epsilon }\right)}{{\mathrm{Re}}_{\mathrm{p}}}+1.75\right)\hfill \mathrm{(4)}& \begin{array}{c}\mathrm{Heat} \mathrm{balance}\\ \mathrm{outside} \mathrm{the} \mathrm{reactor}\end{array}& \begin{array}{l}\frac{\mathrm{d}{\mathrm{T}}_{\mathrm{C}}}{\mathrm{d}\mathrm{z}}={\frac{{\mathrm{t}}_{\mathrm{n}}\mathsf{\pi }{\mathrm{D}}_{\mathrm{t}}\mathrm{U}\left(\mathrm{T}-{\mathrm{T}}_{\mathrm{C}}\right)}{{\mathrm{C}}_{{\mathrm{p}}_{\mathrm{C}}}{\mathrm{w}}_{\mathrm{C}}}}_{}{\to }_{}\mathrm{Co}-\mathrm{current} \mathrm{flow}\\ \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{C}}}{\mathrm{d}\mathrm{z}}={\frac{{\mathrm{t}}_{\mathrm{n}}\mathsf{\pi }{\mathrm{D}}_{\mathrm{t}}\mathrm{U}\left({\mathrm{T}}_{\mathrm{C}}-\mathrm{T}\right)}{{\mathrm{C}}_{{\mathrm{p}}_{\mathrm{C}}}{\mathrm{w}}_{\mathrm{C}}}}_{}{\to }_{}\mathrm{Countercurrent} \mathrm{flow}\end{array}\hfill \end{array}$$

The parameter values used for simulation runs are listed in

Table 2. Data used for the o-xylene oxidation simulation [7,12,27].

Molecular weight of gaseous mixture	$\mathrm{M}$	$29480 \frac{\mathrm{k}\mathrm{g}}{\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}}$
Superficial velocity	$\mathrm{u}$	$3600 \frac{\mathrm{m}}{\mathrm{h}}$
Total pressure	$\mathrm{P}$	$101.325 \mathrm{k}\mathrm{P}\mathrm{a}$
Heat capacity of the gaseous mixture	${\mathrm{C}}_{{\mathrm{p}}_{\mathrm{g}}}$	$0.2498 \frac{\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}}{\mathrm{k}{\mathrm{g}}_{} \mathrm{K}}$
Tube diameter	${\mathrm{D}}_{\mathrm{t}}$	$0.025 \mathrm{m}$
Number of tubes	${\mathrm{t}}_{\mathrm{n}}$	$3000$
Gas mixture density	${\mathsf{\rho }}_{\mathrm{g}}$	$1.293 \frac{\mathrm{k}\mathrm{g}}{{\mathrm{m}}^3}$
Oxygen’s partial pressure	${\mathrm{P}}_{{\mathrm{O}}_2}^0$	$21.07 \mathrm{k}\mathrm{P}\mathrm{a}$
Bulk density of the fixed bed	${\mathsf{\rho }}_{\mathrm{b}}$	$1300 \frac{\mathrm{k}\mathrm{g}}{{\mathrm{m}}^3}$
Heat capacity of the cooling medium	${\mathrm{C}}_{{\mathrm{p}}_{\mathrm{c}}}$	$0.3105 \frac{\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}}{\mathrm{k}{\mathrm{g}}_{} \mathrm{K}}$
Cooling medium flow rate	${\mathrm{w}}_{\mathrm{c}}$	$72000 \frac{\mathrm{k}\mathrm{g}}{\mathrm{h}}$
Overall heat transfer coefficient	$\mathrm{U}$	$82.7 \frac{\mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}}{{\mathrm{m}}^2 {\mathrm{h}}_{} \mathrm{K}}$
Total reactor length	$\mathrm{L}$	$2 \mathrm{m}$
Particle diameter	${\mathrm{D}}_{\mathrm{p}}$	$0.022 \mathrm{m}$
Bed porosity	$\mathsf{\epsilon }$	$0.5$
Particle Reynolds number	${\mathrm{Re}}_{\mathrm{p}}$	$121$

.

3. Kinetic Models

In the context of the present study, three different o-xylene oxidation models were selected: (i) the single reaction model used by van Welsenaere and Froment [7]; (ii) the three-step oxidation process reported by de Lasa [2]; and (iii) the Chandrasekharan and Calderbank [17] kinetic model, which includes five steps of reaction. The rate law and kinetic parameters for each of these models as well as the heat of reaction of each step are presented in

Table 3. Kinetic models for o-xylene oxidation considered in this work.

van Welsenaere and Froment [7]
$$-{\mathrm{r}}_{\mathrm{A}}={\mathsf{\rho }}_{\mathrm{L}}{\mathrm{k}}_1{\mathrm{e}}^{\left(-\frac{\mathrm{a}}{\mathrm{T}}\right)}{\mathrm{P}}_{\mathrm{A}}$$ (5)
${\mathsf{\rho }}_{\mathrm{L}}={0.5203}_{} \frac{\mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}}}{\mathrm{m}}$ $\mathrm{a}={13636}_{} {\mathrm{K}}^{-1}$ ${\mathrm{k}}_1=21.07{\mathrm{e}}^{19.837} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$	$\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_1}=-306707 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$
de Lasa [2]
$$-{\mathrm{r}}_{\mathrm{A}}={\mathsf{\rho }}_{\mathrm{L}}\left({\mathrm{k}}_{21}+{\mathrm{k}}_{23}\right){\mathrm{P}}_{\mathrm{A}}$$ (6)
$${\mathrm{r}}_{\mathrm{B}}={\mathsf{\rho }}_{\mathrm{L}}\left({\mathrm{k}}_{21}{\mathrm{P}}_{\mathrm{A}}-{\mathrm{k}}_{22}{\mathrm{P}}_{\mathrm{B}}\right)$$ (7)
$${\mathrm{r}}_{\mathrm{C}}=8{\mathsf{\rho }}_{\mathrm{L}}\left({\mathrm{k}}_{23}{\mathrm{P}}_{\mathrm{A}}+{\mathrm{k}}_{22}{\mathrm{P}}_{\mathrm{B}}\right)$$ (8)
${\mathrm{k}}_{21}=21.07{\mathrm{e}}^{\left(10.59-\frac{13587.68}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{22}=21.07{\mathrm{e}}^{\left(11.62-\frac{15801.97}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{23}=21.07{\mathrm{e}}^{\left(9.73-\frac{14392.88}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$	$\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{21}}=-307122 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{22}}=-7837000 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{23}}=-1090822 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$
Chandrasekharan and Calderbank [17]
$$-{\mathrm{r}}_{\mathrm{A}}={\mathsf{\theta }}_{\mathrm{v}}\mathsf{\alpha }{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{t}}\left({\mathrm{k}}_{31}+{\mathrm{k}}_{22}+{\mathrm{k}}_{33}\right){\mathrm{P}}_{\mathrm{A}}$$ (9)
$${\mathrm{r}}_{\mathrm{B}}={\mathsf{\theta }}_{\mathrm{v}}\mathsf{\alpha }{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{t}}\left({\mathrm{k}}_{32}{\mathrm{P}}_{\mathrm{A}}+{\mathrm{k}}_{35}{\mathrm{P}}_{\mathrm{E}}\right)$$ (10)
$${\mathrm{r}}_{\mathrm{C}}=8{\mathsf{\theta }}_{\mathrm{v}}\mathsf{\alpha }{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{t}}{\mathrm{k}}_{33}{\mathrm{P}}_{\mathrm{A}}$$ (11)
$${\mathrm{r}}_{\mathrm{D}}={\mathsf{\theta }}_{\mathrm{v}}\mathsf{\alpha }{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{t}}\left({\mathrm{k}}_{31}{\mathrm{P}}_{\mathrm{A}}-{\mathrm{k}}_{34}{\mathrm{P}}_{\mathrm{D}}\right)$$ (12)
$${\mathrm{r}}_{\mathrm{E}}={\mathsf{\theta }}_{\mathrm{v}}\mathsf{\alpha }{\mathsf{\rho }}_{\mathrm{b}}{\mathrm{F}}_{\mathrm{t}}\left({\mathrm{k}}_{34}{\mathrm{P}}_{\mathrm{D}}-{\mathrm{k}}_{35}{\mathrm{P}}_{\mathrm{E}}\right)$$ (13)
${\mathsf{\theta }}_{\mathrm{v}}=\frac{0.026}{0.026+\left({\mathrm{k}}_{31}+3{\mathrm{k}}_{32}+6.5{\mathrm{k}}_{33}\right){\mathrm{P}}_{\mathrm{A}}+{\mathrm{k}}_{34}{\mathrm{P}}_{\mathrm{D}}+{\mathrm{k}}_{35}{\mathrm{P}}_{\mathrm{E}}}$ $\mathsf{\alpha }=3.2832\times 10^{-3} {\mathrm{h}}_{} {\mathrm{m}}^2/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ ${\mathrm{F}}_{\mathrm{t}}=\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}/\mathrm{h}$
${\mathrm{k}}_{31}=137.808{\mathrm{e}}^{\left(-\frac{7381.11}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{32}=46.62{\mathrm{e}}^{\left(-\frac{6552.43}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{33}=12.888{\mathrm{e}}^{\left(-\frac{6153.13}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{34}=20.016{\mathrm{e}}^{\left(-\frac{5585.44}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$ ${\mathrm{k}}_{35}=114.84{\mathrm{e}}^{\left(-\frac{6962.56}{\mathrm{T}}\right)} \mathrm{k}\mathrm{m}\mathrm{o}{\mathrm{l}}_{} {\left({\mathrm{h}}_{} \mathrm{k}{\mathrm{g}}_{\mathrm{c}\mathrm{a}\mathrm{t}} \mathrm{k}\mathrm{P}\mathrm{a}\right)}^{-1}$	$\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{31}}=-108999 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{32}}=-295999 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{33}}=-550000 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{34}}=-99001 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$ $\mathsf{\Delta }{\mathrm{H}}_{\mathrm{r}\mathrm{x}{\mathrm{n}}_{35}}=-87800 \mathrm{k}\mathrm{c}\mathrm{a}\mathrm{l}/\mathrm{k}\mathrm{m}\mathrm{o}\mathrm{l}$

.

4. Computation of Sensitivity Indices

Sensitivity analysis can be divided as “global” and “local”. Global sensitivity analysis describes the effect of the simultaneous large variation of all parameters on the dependent variables. On the other hand, the local sensitivity method provides information on the effect of a model output when only one input factor is changed at a time while all other input factors are held fixed (at their nominal value). Local sensitivity analysis is essentially linear, and it is the one used in this work. Following the ideas proposed by Varma et al. [8], Equations (1) to (4) can be described by a vector y of the n dependent variables (e.g., F_j, T, P, and T_c), which changes in reactor length (z) rendering to the following general differential equation:

$$\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}\mathrm{z}}=\mathbf{f}\left(\mathbf{y};\mathsf{\phi };\mathrm{z}\right)\to \mathbf{y}=\mathbf{y}\left(\mathrm{z};\mathsf{\phi }\right)$$

(14)

where the initial conditions y(0) = y₀ and φ representing the vector comprising the m system input parameters. It is assumed that f is continuously differentiable in all its arguments (as denoted at the right-hand side of the arrow in Equation (14)). If a small alteration of any input parameter occurs (e.g., from φ_j to φ_j + Δφ_j), the corresponding solution for the dependent variables becomes $\mathbf{y}=\mathbf{y}\left(\mathrm{z};{\mathsf{\phi }}_{\mathrm{j}}+\mathsf{\Delta }{\mathsf{\phi }}_{\mathrm{j}}\right)$. In this way, n first-order local sensitivity indices (S) of the vector of dependent variables, y, concerning the same input parameter, ${\mathsf{\phi }}_{\mathrm{j}}$, follows the form:

$$\mathbf{S}\left(\mathbf{y};{\mathsf{\phi }}_{\mathrm{j}}\right)=\frac{\partial \mathbf{y}\left(\mathrm{z};{\mathsf{\phi }}_{\mathrm{j}}\right)}{\partial {\mathsf{\phi }}_{\mathrm{j}}}=\underset{\mathsf{\Delta }{\mathsf{\phi }}_{\mathrm{j}}\to 0}{\mathrm{L}\mathrm{i}\mathrm{m}}\left(\frac{\mathbf{y}\left(\mathrm{z};{\mathsf{\phi }}_{\mathrm{j}}+\mathsf{\Delta }{\mathsf{\phi }}_{\mathrm{j}}\right)-\mathbf{y}\left(\mathrm{z};{\mathsf{\phi }}_{\mathrm{j}}\right)}{\mathsf{\Delta }{\mathsf{\phi }}_{\mathrm{j}}}\right)$$

(15)

The local sensitivity indices can be evaluated using the direct differential method (Varma et al. [8]), which implies differentiate both sides of Equation (14) leading to:

$$\frac{\mathrm{d}\left(\raisebox{1ex}{$\partial \mathbf{y}$}\!\left/ \!\raisebox{-1ex}{$\partial {\mathsf{\phi }}_{\mathrm{j}}$}\right.\right)}{\mathrm{d}\mathrm{z}}=\frac{\mathrm{d}\left(\mathrm{S}\left(\mathbf{y};{\mathsf{\phi }}_{\mathrm{j}}\right)\right)}{\mathrm{d}\mathrm{z}}=\frac{\partial \mathbf{f}}{\partial \mathbf{y}}\frac{\partial \mathbf{y}}{\partial {\mathsf{\phi }}_{\mathrm{j}}}=\frac{\partial \mathbf{f}}{\partial \mathbf{y}}\mathrm{S}\left(\mathbf{y};{\mathsf{\phi }}_{\mathrm{j}}\right)+\frac{\partial \mathbf{f}}{\partial {\mathsf{\phi }}_{\mathrm{j}}}$$

(16)

Its initial conditions can also be obtained by differentiating the initial condition of Equation (14). Thus, depending on which input parameter in the vector φ is chosen, it is possible to obtain:

$${\mathrm{S}\left(\mathbf{y};{\mathsf{\phi }}_{\mathrm{j}}\right)|}_{\mathrm{z}=0}=\{\begin{array}{c}0{,}_{}{\mathsf{\phi }}_{\mathrm{j}}\ne \mathbf{y}\left(0\right)\\ 1{,}_{}{\mathsf{\phi }}_{\mathrm{j}}=\mathbf{y}\left(0\right)\end{array}$$

(17)

Thus, the simultaneous solution of Equation (14) (e.g., Equations (1) to (4)) and sensitivity Equation (16), along with their corresponding initial conditions, obtain the value of the dependent variable and the corresponding local sensibility index, as a function of the independent variable (z). For comparison with the above-mentioned works [7,12,13,16,28], sensitivity indices with respect to the inlet temperature of the reactants (T_o) are evaluated for Fj, T and Tc as follows:

$$\frac{\mathrm{d}\mathrm{S}\left({\mathrm{F}}_{\mathrm{i}},{\mathrm{T}}_{\mathrm{o}}\right)}{\mathrm{d}\mathrm{z}}={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{M}}\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial {\mathrm{F}}_{\mathrm{j}}}\mathrm{S}\left({\mathrm{F}}_{\mathrm{j}},{\mathrm{T}}_{\mathrm{o}}\right)+}\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial \mathrm{T}}\mathrm{S}\left(\mathrm{T},{\mathrm{T}}_{\mathrm{o}}\right)$$

(18)

$$\frac{\mathrm{d}\mathrm{S}\left(\mathrm{T},{\mathrm{T}}_{\mathrm{o}}\right)}{\mathrm{d}\mathrm{z}}={\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{M}}\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial {\mathrm{F}}_{\mathrm{j}}}\mathrm{S}\left({\mathrm{F}}_{\mathrm{j}},{\mathrm{T}}_{\mathrm{o}}\right)+}\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial \mathrm{T}}\mathrm{S}\left(\mathrm{T},{\mathrm{T}}_{\mathrm{o}}\right)+\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial {\mathrm{T}}_{\mathrm{C}}}\mathrm{S}\left({\mathrm{T}}_{\mathrm{C}},{\mathrm{T}}_{\mathrm{o}}\right)$$

(19)

$$\frac{\mathrm{d}\mathrm{S}\left({\mathrm{T}}_{\mathrm{C}},{\mathrm{T}}_{\mathrm{o}}\right)}{\mathrm{d}\mathrm{z}}=\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial \mathrm{T}}\mathrm{S}\left(\mathrm{T},{\mathrm{T}}_{\mathrm{o}}\right)+\frac{\partial {\mathrm{f}}_{\mathrm{i}}}{\partial {\mathrm{T}}_{\mathrm{C}}}\mathrm{S}\left({\mathrm{T}}_{\mathrm{C}},{\mathrm{T}}_{\mathrm{o}}\right)$$

(20)

with the initial conditions:

$$\mathrm{z}=0, {\forall }_{\mathrm{i}} \mathrm{S}({\mathrm{F}}_{\mathrm{i}}, {\mathrm{T}}_{\mathrm{o}})=0, \mathrm{S}(\mathrm{T}, {\mathrm{T}}_{\mathrm{o}})=1, \mathrm{S}({\mathrm{T}}_{\mathrm{C}}, {\mathrm{T}}_{\mathrm{o}})=0$$

(21)

5. Numerical Simulation

It is worthy of note that, depending on the reaction rate expressions included in the mathematical model, the final form and the number of parametric sensitivity equations will be different for each case analyzed in this work. Consequently, each reaction problem will correspond to a different numerical challenge, because each mechanism must fulfill the operating sensitivity limit (S (T; To) = 1) at a single point for the interval z = (0.2]. Their solution strategies are summarized in Section 5.1. and Section 5.2.

5.1. For Constant Temperature or Co-Current Flow in the Shell Side

These cases imply the integration of a set of, initial value, ordinary differential equations (ODE). Their solution is based on the execution of the ode15s MatLab® function, linked to the following condition: ∃! z [0,2]: S(T; T_o) = 1. This means that it is necessary to search, within the solution range of the sensitivity index of the reactive temperature, relative to the inlet temperature, a single condition different from the initial one, such that S(T; T_o) = 1. Thus, if the inlet temperature is set, the inlet o-xylene pressure will be varied until the criterion will be met (using additionally the fsolve MatLab® function with a tolerance of 1 × 10⁻⁶).

5.2. For Countercurrent Flow in the Shell Side

The solution of this case infers, in addition to the sensitivity index restriction, the fulfillment of the boundary conditions imposed by the fluid in the shell side. Thus, it is necessary to include two more restrictions: i) z = 2 → Tc_o = T_o and ii) z = 2 → S(T_C; T_o) = 0; generating a quite challenging numerical problem which requires modifying, simultaneously, three joint variables and a trial-and-error procedure. To achieve this, the shooting method was used where the inlet pressure or temperature (depending on the case), the outlet temperature of the reactive fluid and the sensitivity index for the inlet temperature, relative to the reactive temperature at the reactor outlet, were set. Then, to minimize the objective functions, a restricted nonlinear least-squares minimization method (using the lsqnonlin MatLab® function) allows for the limiting of the solutions to those corresponding to logical intervals of system operational conditions. For very demanding cases, such as the one including five reactions, numerical difficulties were overcome by an orthogonal collocation method (using the bvp4c Matlab® function).

6. Results

To analyze the behavior of o-xylene conversion and the reaction temperature, and to give a clear interpretation of the temperature sensitivity index, the reactor performance was examined on the basis of each of the three kinetic models. Moreover, the influence of the two types of heat exchange operation (e.g., variable fluid temperature with co-current or countercurrent operation) was also compared.

6.1. Constant Cooling Temperature (T_C = T₀ = 625 K)

At first, critical feed pressure values, under constant cooling temperature of 625 K, were calculated according to the approach proposed in this manuscript as well as to the models previously presented by van Welsenaere and Froment [7]; Hlavacek [28]; and Henning and Perez [16] and are resumed in

Table 4. Critical feed pressures at different service fluid configurations. Comparison between different criteria obtained from different kinetic models.

Kinetic Model	Constant Tc = To = 625 K		Co-Current Tco = To = 625 K		Countercurrent Tco = To = 625 K
	Author Criterion ${\mathbf{P}}_{\mathit{\text{o-x}}\mathit{y}\mathit{l}\mathit{e}\mathit{n}\mathit{e},\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}^0$ kPa	This Work ${\mathbf{P}}_{\mathit{\text{o-x}}\mathit{y}\mathit{l}\mathit{e}\mathit{n}\mathit{e},\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}^0$ kPa	Author Criterion ${\mathbf{P}}_{\mathit{\text{o-x}}\mathit{y}\mathit{l}\mathit{e}\mathit{n}\mathit{e},\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}^0$ kPa	This Work ${\mathbf{P}}_{\mathit{\text{o-x}}\mathit{y}\mathit{l}\mathit{e}\mathit{n}\mathit{e},\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}^0$ kPa	This Work ${\mathbf{P}}_{\mathit{\text{o-x}}\mathit{y}\mathit{l}\mathit{e}\mathit{n}\mathit{e},\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{a}\mathit{l}}^0$ kPa
One reaction	1.677 [7] 1.778 [16] 1.815 [28]	1.863	1.422 [1]	1.733	1.417
Three reactions	0.867 [16] a	1.187	0.867 [16] a 1.014 [2]	1.106 a 1.097	0.847
Five reactions	-	0.692	0.800 [2]	0.692	0.689

^a for a cooling flow of 79,200 kg/h.

. Notice that these simplified models [7,16,28] did not contemplate the pressure drop equation. For the single reaction kinetic model, the application of our approach allows us to use the higher value of ${\mathrm{P}}_{\mathit{\text{o-x}}ylene}^0$, implying less conservative operation conditions. Notice that, in addition, the value of critical inlet o-xylene pressure depends on the complexity of the kinetic model: it decreases with an increase in the number of considered reactions. This is probably due to the high exothermicity of additional reactions.

6.2. Variable Cooling Temperature According to Co-Current or Countercurrent Pattern (T_Co = T_o = 625 K)

Next, critical feed pressure values, obtained under variable cooling temperature, for both co-current or countercurrent pattern, were calculated using the proposed model, and are also presented in Table 4 for the three kinetic models. In the case of co-current flow pattern, the comparison of the values obtained using the proposed approach with these previously reported by Soria-López et al. [1], de Lasa [2] and Henning and Perez [16] shows significant differences. They can be related to the effect of pressure drop balance used in our model (approach). Both for the single step and three-step reaction kinetic model, higher critical feed pressure values are predicted using our model. However, the contrary occurs for the kinetic model involving five steps of reaction. These results prove the importance of both the reactor model and the kinetic model for the selection of safe operation conditions. Figure 3 shows the resulting conversion, reaction temperature and temperature sensitivity index profiles, illustrating the sensitivity of the system with respect to the kinetic model and cooling flow configuration. The co-current case was also analyzed by Soria-López et al. [1] and de Lasa [2], assuming a negligible pressure drop in the catalytic bed. Total o-xylene conversion was predicted only for the five-reaction kinetic model (Figure 3a,b). An inspection of temperature profiles shows how the hotspot position and intensity (ca. 675 K) are similarly predicted from one- or three-reaction kinetic models for the co-current case. However, the more rigorous kinetic model predicts a faster and higher temperature increase in the PBR (Figure 3c,d). A change in the cooling pattern, to the countercurrent one, amplifies the maximum temperature for one- and three-reaction kinetic models (maximum of temperature of 695 K and 691 K, respectively), while it remains constant for the five-reaction mechanism. The corresponding S(T, To)-z trajectories are presented in Figure 3e,f). As defined previously, the zone where the derivative of the sensitivity index with respect to z is positive (for z ≠ 0) is of special interest. Thus, for the operation conditions presented in Table 2, the critical inlet o-xylene pressure is summarized in Table 4 for different kinetic models and criteria. As in the constant cooling temperature case, our criterion is less conservative than these previously reported in the literature.

6.3. Variable Cooling Temperature According to Co-Current or Countercurrent Pattern (P_Ao = 0.9322 kPa)

This case is concerned with the determination of the maximum permissible inlet temperature of the reactor fluid, as a function of cooling flow pattern and the kinetic model, for a given inlet o-xylene partial pressure (the selected value of P_Ao = 0.9322 kPa was justified elsewhere [27,29]). Here, as defined in Section 5.2, for each kinetic model, the solution must fulfill both the sensitivity index restriction and the boundary conditions imposed by the fluid in the shell side. Figure 4 presents the obtained results for both cooling flow configurations. It illustrates the sensitivity of inlet temperature with respect to both the kinetic model and cooling flow configuration. In all cases, the lower feed temperatures, for safe operation conditions, is suggested by the five-reaction kinetic model (e.g., (i) van Welsenaere and Froment [7] = 655 K and 650 K, (ii) de Lasa [2] = 633 K and 619 K and (iii) Chandrasekharan and Calderbank [17] = 598 K and 598 K for co-current and countercurrent patterns, respectively).

Finally, the critical operation values are found readily by imposing the conditions of Equation (21), as in Figure 3e,f, Figure 4e,f. Thus, safe operation behavior is fully determined by inspecting a very large number of numerical solutions of the model (reactor and sensitivity index) for a wide variation of a specific inlet parameter. Figure 5 illustrates the runaway diagram for the o-xylene oxidation, specifically the effect of o-xylene partial pressure feed, in the range of 0.0046 to 0.0184, according to [29], for the three analyzed kinetic models. The curves in Figure 5 define a sepatrix that limits two regions: if the operating conditions are such that they lead to a coordinate in Figure 5 above the corresponding curve, extreme parametric sensitivity and runaway is expected; but if it is located under the curves, the reactor is insensitive to small fluctuations. Notice that only the five-step reaction kinetic model defines a single runaway region, regardless of the cooling flow pattern. As expected, the simplest kinetic model predicts the smaller runaway region.

7. Conclusions

The parametric sensitivity of a Packed Bed Catalytic Reactor depends on the cooling flow pattern and kinetic model. The required mathematical model, involving a one-dimensional reactor model and a sensitivity index formulation, was analyzed to predict an a priori estimation of safe reactor operation conditions. An effective numerical routine for the solution of the countercurrent flow pattern was also presented. The o-xylene oxidation was analyzed as an example of a complex kinetic model and it was successfully solved. It provided less conservative operational conditions compared to the models previously proposed, even in the case of low inlet pressures of a key reactant. The benefits of detailed kinetic models include the possibility of analyzing product selectivity as additional criteria for the selection of reactor operation conditions (to be explored in a future contribution). The presented methodology is generic and can accommodate rigorous kinetic, reactor and sensitivity models, which could provide significant results for reactor design and operation.