Thermal Stability of Gas Oil Hydrotreating Processes: Numerical Issues of the Matrix-Eigenvalue Approach

Résumé — Stabilité thermique de procédés d’hydrotraitement des gazoles : aspects numériques de l’approche par valeurs propres matricielles — Les procédés qui mettent en œuvre des réactions très exothermiques nécessitent une attention particulière afin d’éviter l’augmentation non contrôlée de la température connue comme emballement thermique. Une analyse de stabilité thermique est nécessaire afin d’établir les conditions d’opération sûre et productive des procédés exothermiques. L’hydrotraitement de gazoles met en œuvre principalement des réactions d’hydrogénation ; l’hydrotraitement de charges très insaturées comme les gazoles light cycle oil peut être fortement exothermique. Pour cette raison, ce procédé fait l’objet d’une étude de stabilité thermique. La théorie des perturbations a déjà été appliquée pour effectuer une analyse de stabilité de ce procédé en conditions dynamiques. Cette méthode consiste à perturber le modèle et à résoudre le modèle perturbé sous forme d’un problème


NOMENCLATURE a
Reaction order of hydrocarbon pseudo-component in the liquid phase A Heat exchange surface (m 2 /m 3 ) b Reaction order of hydrogen in the liquid phase C Concentration (mol/m 3 ) Cp Specific heat Capacity (J/kg/K) D ax Axial dispersion coefficient (m 2

INTRODUCTION
One of the security priorities involving highly reactive systems is the thermal runaway risk.The temperature increase for the reactions that follow an Arrhenius law induces the rise of the heat generation that further increases the reaction temperature; this situation may result in a thermal runaway.Some refining processes carry out highly exothermic reactions; some examples are Fischer-Tropsch synthesis, hydrotreating, hydrocracking and selective hydrogenations.Hydrotreating (HDT) is one of the most important processes in the oil refining industry.There are more than one thousand HDT units all around the world covering a capacity of about 45 millions of barrels per calendar day [1].Hydrotreating aims to purify oil the perturbed model solution as a standard eigenvalue problem or as a generalized eigenvalue problem are presented.The computation of the Jacobian by a numerical approach or with the analytical expressions is also carried out.In both cases, results are compared and their influence on the stability/instability results is presented.
fractions to produce fuels at commercial specifications.It is carried out at high hydrogen partial pressures and high temperatures with a solid catalyst.Hydrotreating of gas oil cuts is exothermic, the main reactions being aromatics and olefins hydrogenation, hydrodesulfurization (HDS) and hydrodenitrogenation (HDN).The exothermic level attained mainly depends on the content of unsaturated compounds in the gas oil processed.This content depends on the gas oil type.Light Cycle Oils (LCO) usually contain 70-80 weight percent of aromatics and 5-10 weight percent of olefins.Hydrotreating this kind of feeds is consequently highly exothermic.For safety reasons and their poor quality, LCO hydrotreating is carried out by dilution with less exothermic feeds.An alternative to improve diesel production can be the incorporation of higher contents of LCO in the HDT feed stream.Thus, it is necessary to determine the consequence on the process exothermicity.For this reason, this study deals with the thermal stability study of LCO hydrotreating.Thermal stability analysis is an engineer's support for reactor design and operation.Such a study must guarantee reliable operation and also ensure to attain the target selectivities and conversions.It is also essential to determine the regions of operating conditions where the reactor behavior becomes unreliable or dangerous (runaway regions).In practice, design of reactive systems is mainly based on a stationary criterion described by van Heerden [2].The heat generated by the reactions and the heat removal is represented as a heat removal efficiency diagram.The stationary stability is verified if the slope curve of heat generated by the reactions is lower than the slope of the heat transferred through the wall.If this condition is not verified, three possible steady states can be observed leading to unstable operation.Then, a hysteresis behavior is observed.
Another approach for thermal stability analysis under steady state conditions is parametric sensitivity.This method evaluates the sensitivity of the system's behavior following to changes in parameters of the system.In the parametrically sensitive regions the behavior of the reactor changes sharply with small variations of inlet parameters.This situation can lead to thermal runaway [3].This concept has been widely applied for chemical reactors thermal stability analysis [4][5][6][7][8][9][10][11][12][13][14].Both approaches, performed under stationary conditions, help to define stability/instability maps as a function of the most important parameters such as operating conditions and model parameters.However, as the thermal balance transient term is not taken into account, dynamic stability cannot be ensured in the boundary of stationary stability.The stationary stability is a necessary condition but not sufficient to ensure stability.Hence, a dynamic analysis is essential to confirm or enlarge the unstable regions determined by the stationary study; the work presented here is focused on the dynamic stability study.
The mathematical model consists of a system of Partial Differential Equations which is nonlinear, in particular due to the Arrhenius law of chemical reaction kinetics.By applying dynamical system theory (cf.[15,16]) to this system, steady equilibrium states are determined.The Jacobian of the nonlinear system at those stationary points defines the linear dynamical system for small perturbation: if the linearized system is stable, then small perturbations will ultimately decay and the stationary point is linearly stable.If however the linearized system is unstable, then even infinitesimal perturbations will amplify and linear stability theory hence provides sufficient conditions for instability.This dynamical system approach is commonly used in mechanics and it has also been applied in the field of chemical engineering.Aris and Amundson [17] used this linear approach to assess the stability of a reactor and the method was applied for control purposes to decide which parameter is more efficient to control an operating point (steady state).Perlmutter [18] also exposed the application of dynamic linear stability theory to chemical reactors.
While the application of the perturbation theory to chemical systems is not new, for most works reported in literature [19][20][21] the reaction systems adopted are in general rather simple, including simple reaction networks, and the reactors are frequently assumed to be homogeneous or pseudo-homogeneous.More recent works deal with the modeling of local hot spots in packed bed reactors [22][23][24].
Few studies of complex reactive systems, such as gas oils hydrotreating, are available.In a previous work, a dynamic thermal stability analysis of an HDT pilot plant has been undertaken [25,26].The local stability of several steady-state points was investigated using the linear stability approach.
The contribution of this work is to apply this theory to a fairly complex reactive system.The HDT model describes the packed bed upflow reactor.It involves 3 phases (gasliquid-solid), 6 lumped chemical families and 3 boiling point cuts.In this way, 18 Partial Differential Equations (PDE) and 18 algebraic equations constitute the model discretized in the axial axis.A regulation system was implemented to determine the steady-state solution [25].The Jacobian was approximated using finite differences and the algebraic equations were explicitly solved to be able to apply the standard eigenvalue problem formalism for the linear stability analyses.Several cases were tested.The spectral analysis of the eigenvalues indicating the stable/unstable behavior of the reactive system is compared with dynamic simulations.An excellent agreement was found between the simulations and the stability analysis.
The present work focuses on formalism and numerical assumptions for the spectral analysis of the Jacobian matrix.The first objective is to determine the impact of the stability conclusions if the problem is solved as a generalized eigenvalue problem instead of its corresponding standard eigenvalue problem.Another point addressed in this paper is the accuracy of the Jacobian and the consequence on the stability results.The numerical approach using finite differences to calculate the Jacobian is easy to implement for such a complex system.However, the order of magnitude of some variables can be very different (for instance pressure in Pa versus concentrations in mol/m 3 ).Therefore, the estimation of some derivatives may be less accurate.The analytical expressions resulting from the PDE and algebraic equations are developed and used to calculate the exact Jacobian (analytical).Results are compared to those obtained by the numerical approach.

HYDROTREATING REACTOR MODEL
Gas oils hydrotreating is the refining process studied in this work.Since main HDT reactions are exothermic (hydrogenations), it is important to study the stability of the reactive system.A reactor model has already been developed and is described elsewhere [25].The main features of this model are described as follows.The model represents the triphasic gasliquid-solid system.The reactor has a fixed bed of solid catalyst.The feed is injected in the reactor in upflow mode.Gas flow is considered to be in plug flow while the liquid phase flow is described by an axial dispersion model.The complex chemical composition and reactivity of gas oils is taken into account with 15 hydrocarbon lumps distributed in three boiling point cuts.The chemical lumps are: triaromatics (TRI), diaromatics (DI), monoaromatics (MONO), saturates (SAT), sulfur compounds (SULF) and olefins (OLEF).The boiling point cut ranges are: Initial Boiling Point (IBP) -200°C, 200-300°C and 300°C -Final Boiling Point (FBP).The reaction network is illustrated in Figure 1.Reactions taken into account are aromatics hydrogenation (R1, R2, R5, R6, R9 and R12), olefins saturation (R4, R8 and R10) and hydrodesulfurization (R3, R7 and R11).In addition hydrogen and hydrogen sulfide have also to be taken into account (H 2 S marked in gray since it is produced by hydrodesulfurization).This results in 17 lumps (Tab.1).Reaction kinetics is described via a power law type as shown by Equation ( 1), (lump i = 16 corresponds to hydrogen): (1) Parameters were fitted with pilot plant experiments carried out at industrial operating conditions and real LCO gas oil feeds.
The transient material balance for each lump is given by Equation ( 2).This means that there are 17 PDE concerning the material balance.The transient thermal balance of the reactor is indicated in Equation (3) (see Eq. 2, 3).

Lump index (i)
Reference One should note that gas velocity also changes all along the reactor.This is due to several facts: a) the temperature profile may change along the axis of the reactor, b) H 2 and H 2 S are respectively consumed and produced and c) hydrocarbon lumps are partially vaporized.As illustrated by Equation ( 4), the last depends on gas-liquid equilibrium (G-L-E) and also on temperature where the "Henry constants" H i are in fact variables since they are also temperature and pressure dependent.There are 17 algebraic equations describing the gas-liquid-equilibrium: (4) Moreover, as isobaric conditions are assumed, at each time step and along the reactor's axis the constant pressure balance must be satisfied.Hence, an additional algebraic relation is added to the reactor model: (5) In total, for the discretized system the number of PDE to be solved, resulting in 17 × nz material balances for each chemical lump and nz thermal balances, where nz is the number of discretization points along the axis z of the reactor.In addition, the algebraic relations associated with gas-liquid equilibrium have to be solved, which gives rise to 17 × nz G-L-E equations.Finally, the 1 × nz pressure balance equations have also to be taken into account.The integration along the z axis is carried out using an upwind scheme for the convective terms and a centered scheme for the diffusion-dispersion terms.An explicit 1st order Euler scheme was used for integration in time.Results presented in this work are all obtained for a discretization with nz = 50.

Stationary Analysis
The stability of operation of chemical reactors at stationary conditions is based on the van Heerden criterion [2].The principle is also described by Froment and Bischoff [27].This criterion imposes that the slope of heat generated by the reactions must be lower than the slope of the heat transferred through the wall.This can be expressed as: The comparison between both slopes is generally represented as a stationary stability diagram.Two representations are commonly used.One representation consists in plotting the reactor temperature (T) as a function of the operating parameter allowing to control T (for instance a wall cooling temperature or feed inlet temperature).To satisfy the steady state stability, for a given cooling or inlet temperature, only one reactor temperature must be obtained.For unstable reactors, a curve having an S shape will be obtained.Consequently, three reactor temperatures will be possible for a given cooling or inlet temperature.A hysteresis behavior is then observed.The stationary analysis is useful to identify a preliminary stable/ unstable region.However, as stated before, this condition is necessary but not sufficient to ensure stability.

Dynamic Analysis
The dynamic stability analysis carried out in this work consists in the following main steps.The system once discretized with respect to the spatial variable z, the resulting large algebraic-differential system is advanced in time and a regulator is applied and (possibly unstable) steady states are computed.Once an equilibrium state retrieved, the nonlinear system is linearized leading to the linear dynamical system for small perturbations involving the Jacobian matrix.The eigenvalues of the Jacobian are computed, the real parts of which determining whether the equilibrium state is stable or unstable.
The disturbance solution vector x, containing the perturbation at the discrete points z k along the reactor axis, is sought as an expansion in terms of the eigenvectors U j of the Jacobian at the stationary point, that is: From Equation (7), it is clear that if there are eigenvalues with positive real parts, any perturbations will grow when time increases.Conversely, the condition of asymptotic stability implies that all the eigenvalues must have a negative real part.Nonzero imaginary part of the eigenvalues will give rise to oscillatory behavior, the perturbation amplitude being attenuated or amplified, depending on the sign of the real part of the eigenvalues.
The question that remains is how to formalize the eigenvalue problem depending on the model.The HDT reactor model is a system of PDE-algebraic equations.The perturbed model can be solved as a standard or as a generalized eigenvalue problem and both approaches are exposed in Section 3.

DYNAMIC ANALYSIS: FORMULATION AS A MATRIX EIGENVALUE PROBLEM
For differential-algebraic systems, there are two ways to evaluate the stability of a steady-state.One is based on the standard eigenvalue problem and the second is referring to a generalized eigenvalue problem.These two methods were applied to the dynamic model for LCO hydrotreating.
A comparison of results obtained with both approaches is presented.The whole eigenvalues spectra presented in this work are computed using the EISPACK Fortran library.

Matrix-Eigenvalue Problem
The spatial differential operators in the variable z once discretized, one recovers an algebraic-differential system which may be written formally as: y being the solution vector containing the 17 liquid and gaseous phase concentration of lumps as well as the temperature and the gas velocity at the nz points of discretization.A steady state solution is therefore: f(y) = 0 (9) and superimposing a perturbation x, the general linearized perturbation model can be written as Ax ⋅ = Jx where A has the following structure (see Scheme 1).
The Jacobian matrix J, that is the partial derivatives of the nonlinear function f(y) evaluated at the steady state, is defined as: (10) According to the expansion indicated in Equation ( 7), the eigenvalues λ j and eigenvectors U j are solution of Equation ( 11): λ j AU j = JU j (11) This formalism corresponds to the generalized eigenvalue problem and this system can be solved as a standard eigenvalue problem only if the matrix A is invertible, which is not the case for the present HDT model.However, it is possible to formulate the stability of the system as a standard eigenvalue problem, by introducing the algebraic Equations ( 4) and ( 5) into the PDE Equations ( 2) and (3).This operation transforms the problem into a purely PDE system given by Equations ( 12) and ( 13) (see.Eq. 12, 13).
We obtain a new linearized model of perturbation: A' ⋅x ⋅ = J' ⋅x where A' is now an invertible matrix and J' is the Jacobian matrix of the new system.Then, this system is solved as a standard eigenvalue problem since x ⋅ = A' -1 ⋅J' ⋅x = M⋅x.Due to the reformulation of the system, the size of the Jacobian J' is reduced to the dimension [18 × nz] × [18 × nz].Of course, here a standard eigenvalue problem is recovered:

Stationary Stability Results
In order to evaluate the impact on the dynamic stability results, the perturbed model is solved as a standard or as a generalized eigenvalue problem and some reference cases were considered.For this purpose, a stationary stability curve was traced by model simulation.The following conditions were kept constant: feed composition Z = 100% LCO; liquid hourly spatial velocity LHSV = 2 h -1 ; inlet feed temperature T inlet = 200°C; heat transfer coefficient U wall = 3 W/m 2 /K; hydrogen/hydrocarbon feed ratio H 2 /HC = 1000 Nm 3 /m 3 and total pressure P t = 150 barg.To represent the stability curve, the average temperature (T mean ) along the reactor was obtained by imposing a constant cooling temperature (T c ) along the reactor (z axis).For each cooling temperature it is possible to calculate an average temperature (T mean ).In the case of multiplicity of steady states, the intermediate points cannot be determined since they are unstable.For these cases, it is necessary to apply a regulation that consists to artificially stabilize this point in order to determine the solution of this steady state.The regulation system is explained in Reference [25].
Figure 2 illustrates the curve of average temperature along the reactor (T mean ) as a function of the reactor wall cooling temperature T c .The gray zone inside the dotted lines corresponds to the unstable zone according to stationary van Heerden criterion.Some steady states are marked as case 1 to case 7.These seven steady states were used for the dynamic analysis, considering the formulations of the stability problem as both a standard and generalized matrix eigenvalue problem.Table 2 details the coordinates of the seven cases in the T mean -T c plane.The van-Heerden stationary criterion (Eq.6) was verified for all cases.Only case 1 and case 7 proved to be stable according to this criterion.

Comparative Results between Standard and Generalized Eigenvalue Problem
Both the standard and generalized eigenvalues were computed for the seven cases illustrated in Figure 2 and   2. For all cases the highest real parts of eigenvalues calculated by both approaches are compared.Figure 3 illustrates the results for all cases.The abscissa cuts at a zero value, hence positive real parts and negative real parts can be clearly distinguished.The standard eigenvalue approach indicates that cases 4 and 5 are unstable while the generalized eigenvalue approach results indicate instability for cases 2, 3, 4, 5 and 6.For case 2, instability is obtained with a very low real part (7 × 10 -7 ) of the eigenvalue which indicates almost marginal instability.Indeed, as shown in Figure 2, cases 2 and 3 are very close to the limit point corresponding to the lower branch of the stability curve.This limit point is the boundary between the stable and unstable part of the curve and this behavior is retrieved by the stability analysis.
Hence, the stability results obtained with both formulations slightly differ, the generalized eigenvalue problem seeming to be more reliable with regard to the van Heerden criterion.The differences between both approaches may however be within the numerical uncertainty, the eigenvalues having a real part very close to zero in case 3. Figure 4 shows the forty eigenvalues with the highest real parts obtained for case 3 by the standard and generalized problem.Although the unstable (or least stable) eigenvalue is slightly different, both spectra are very similar.This confirms that it is difficult to conclude between both approaches for cases having eigenvalues very close to zero.
Differences between both approaches may be explained by the reformulation of the algebraic-differential system to a pure differential one.This step requires the introduction of the algebraic equations into the differential ones.In particular, in Equation ( 12) the temporal derivative of the Henry constant divided by the temperature has to be linearized.This operation may affect the accuracy of the results.
Hence, for algebraic-differential systems it is preferable to keep the generalized eigenvalue approach instead of the Highest eigenvalue real parts for the seven cases studied.
Comparison between standard and generalized eigenvalue problem.
standard one.Moreover, one has to keep in mind that transforming the generalized eigenvalue problem into a standard one necessitates some manipulations and this operation may not be easily achieved for more complex models.

Numerical Jacobian
The impact of the Jacobian matrix calculation on the accuracy of the stability results is also evaluated in this work.The calculation of the elements of the Jacobian matrix requires the values of the derivatives of the model function as a function of each variable, at the stationary points.In our case both the numerical and analytical Jacobians are computed and stability results obtained with both Jacobians are compared.The numerical Jacobian is computed column wise, approximating the partial derivatives by finite differences with the value ε = 10 -6 . (15)

Analytical Jacobian
The analytical expressions of the reactor model to compute the Jacobian are derivated.This task, apparently simple, is in fact rather intricate.In the reactor model (Eq.1-5) the complex dependency between material balances, thermal balances and ε ε algebraic equations is not immediately apparent.In addition, derivatives have to be expressed as a function of the discretized scheme at the discrete points z k , k = 1, …, nz, along the reactor's axis.Due to the boundary conditions, one has to distinguish the reactor inlet coordinate z 1 and outlet coordinate z nz from the interior positions z k , k = 2, …, nz-1.
Only a brief description of the most relevant analytical derivatives of the discretized scheme is carried out in this section, by considering only the interior points.
For reading ease purposes, the material balance given in Equation ( 2) will be referred as F 1 .The thermal balance in Equation ( 3) will be noted as F 2 .Gas-liquid equilibrium (G-L-E) in Equation ( 4) will be referred as F 3 and pressure balance (PB) given in Equation ( 5) as F 4 .The reaction rates expressions Equation ( 1) are noted as F 5 .

Analytical Derivatives for Material Balance
The quantities at the positions z k are simply denoted by The derivative of gas concentration with respect to the material balance in Equation ( 2) at the k and k-1 axial positions along the reactor (for k ∈[2 to nz-1]) is calculated as follows: (16) (17) with z k+1z k = Δz the discretization width along the axis.Note that, as mentioned before, first and second derivatives with respect to z for any quantity y are approximated respectively by first order and second order finite differences with: For the liquid concentration in the material balance, the derivative of liquid concentration with respect to the material balance for the k-1, k and k+1 ∈[2 to nz-1] positions can be written as:  Spectrum highest eigenvalues for case 3. Comparison between standard and generalized eigenvalues, both calculated with a numerical Jacobian.(22) Reactions are considered to be carried out only in the liquid phase in contact with the catalyst.Hence, from Equation (1), for each j reaction, the general expression of the derivative for the reaction rates (F 5 ) with respect to liquid concentration of hydrocarbon lumps can be written as: (23) In the same way, for each j reaction, the derivative of the reaction rates (F 5 ) with respect to hydrogen liquid concentration becomes: (24) The dependency of material balance as a function of temperature is taken into account through the source derivative: (25) The source derivative as a function of temperature is expressed in Equation (26).As for liquid concentration dependency in the source term, the contribution from all reactions (R1 to R12) has to be accounted for in Equation (25) (not detailed here): Finally, since gas velocity also changes along the reactor, material balance also depends on this variable: (27) (28)

Analytical Derivatives for Thermal Balance
The analytical derivatives for thermal balance are more involved than the expressions obtained for material balance.For example, heat capacities, thermal conductivity and the source term are all function of the liquid concentration.The derivative of liquid concentration with respect to the thermal balance for the k and k-1 ∈[2 to nz-1] positions can be written as: (29) where: The derivative of temperature with respect to the thermal balance for the k-1, k and k+1 ∈[2 to nz-1] positions can be resumed as: (see Eq. 32-33).
(33) where: At positions k-1 and k+1 one gets: The derivative of gas concentration with respect to the thermal balance for the k ∈[2 to nz-1] positions can be resumed as: (see Eq. 34).
The derivative of gas velocity with respect to the thermal balance for the k ∈[2 to nz-1] positions is: (37)

Analytical Derivatives for Gas-liquid Equilibrium
Gas-liquid equilibrium (G-L-E) is an algebraic expression Equation ( 4).The analytical derivatives for gas-liquid equilibrium are illustrated in this section.It should be noted that Henry "constants" (H i ) are function of total liquid concentration, temperature and total pressure.The derivative of G-L-E with respect to liquid and gas concentrations for the k ∈[2 to nz-1] positions are respectively: The derivative of G-L-E with respect to temperature for the k ∈[2 to nz-1] positions is: (40)

Analytical Derivatives for Pressure Balance
Pressure balance is also an algebraic expression (eq.5).Pressure balance depends only on gas concentration of lumps and temperature.Hence, both analytical derivatives for pressure balance are given by Equations ( 40) and ( 42):

Results: Comparison between Numerical and Analytical Jacobian
The seven cases already discussed in Sections 3.2 and 3.3 are evaluated.The generalized eigenvalue problem is applied in this section.Jacobians of the seven stationary points corresponding to each case are computed with numerical and analytical approach.Figure 5 illustrates the highest eigenvalue real part for each case.Again the abscissa cuts at a zero and positive real parts can be clearly distinguished from negative real parts.As illustrated, cases 1 and 7 are stable whether the eigenvalues are calculated by numerical or analytical Jacobian.Stability conclusions are also in agreement for the unstable cases 3, 4, 5 and 6.Different conclusions are obtained for case 2: the computation with the numerical Jacobian predicts instability while the calculation with the analytical Jacobian indicates that this operating point is stable.However, since eigenvalues are very close to zero the results for case 2 are at the stability threshold.
A comparison between non zero elements of both Jacobians (numerical and analytical approach) corresponding to case 2 is also carried out.The relative error of each no zero element is calculated, taking as reference the analytical Jacobian.Two numerical Jacobians are computed with two different absolute derivative steps ε = 10 -6 and ε = 10 -7 .Results are shown in Figure 6.For a step of 10 -6 , relative Highest eigenvalue real parts for the seven cases studied.
Comparison between numerical and analytical Jacobians.
errors are comprised within ±0.06%.Calculations with a lower step of 10 -7 lead to a higher relative error range (±1.3%).In that case ε is too small and cancellation occurs in the numerator of the difference quotient (15) due to machine precision.The whole eigenvalues spectrum obtained for case 2 with analytical and numerical Jacobians is illustrated in Figure 7.It can be noticed that both spectra perfectly superimpose.
As concluding remark, results indicate that the two approaches are very close.It may be noticed that the numerical Jacobian approach has larger implementation advantages over the exact analytical approach, which indeed may not be feasible any more for more complex systems.

CONCLUSIONS
The reliable and safe operation of chemical reactors is an essential goal for industry.A safety concern for highly exothermic processes is thermal runaway.This risk must be assessed and reduced with safe designs and with identification of safe/unreliable/dangerous operating regions.Hence, it is essential to carry out a thermal stability analysis to identify a priori stable and unstable regions.Hydrotreating is a refining process that mainly consists on hydrogenation reactions which are highly exothermic.A thermal stability study of gas oils hydrotreating is therefore tackled in this work.A dynamic model that accurately represents a hydrotreating pilot plant was employed for this study.It takes into account three phases, the complex chemical gas oil composition and the most representative reactions contributing to heat generation.
The dynamic stability theory, capable of predicting the stability properties of the system close to stationary points, is applied in this work.The tangent system close to the stationary point is linearized and the stability condition states that all eigenvalues of the linearized system must have a negative real part.However, some questions arise concerning the calculation of eigenvalues.Indeed, the problem can be formulated as a standard or generalized eigenvalue problem.Also, the accuracy of the Jacobian, whether calculated by a numerical approach or obtained through analytical expressions, may have some on the reliability of the results.Both issues are addressed in this work by means of seven test cases, that is seven steady state points for which stationary stability were also determined.
Concerning the standard and generalized eigenvalue problem, spectra obtained are relatively similar.However, in the neighborhood of limit points, the generalized eigenvalue approach is in agreement with van Heerden criterion, whereas there are some differences when using the standard eigenvalues.The reformulation of the algebraic-differential system to a pure differential system diminishes the size of the system but at the same time this approach introduces new terms, which apparently slightly deteriorate accuracy.Hence, for algebraic-differential systems it should be preferred to use the corresponding generalized eigenvalue approach for stability.Also, for the present model it was indeed possible to eliminate the algebraic part of the system to recover a pure dynamical formulation.This rewriting of the system could however be much more involved or even impossible for more complex reactor models.Relative error between numerical and analytical Jacobians for non zero elements (Case 2).Numerical Jacobians calculated with two absolute derivative steps (epsilon).
-  The results may also depend on the way the Jacobian is computed.As discussed above, the partial derivatives with respect to the components of the state vector may be written explicitly, or alternatively they can be evaluated by divided differences.For the seven cases considered here, stability conclusions are equivalent for six cases.For the case where conclusion differs, it must be pointed out that the stationary point is located in the limit point neighborhood.Since results are very similar between both approaches, it is preferable (and much easier in general) to implement the numerical calculation of the Jacobian.In this case, attention must be paid to the absolute derivative step (ε in Eq. 15) which has to be chosen appropriately for a given machine precision.

Figure 2
Figure 2Average temperature along the reactor axis (T mean ) as a function of wall cooling temperature (T c ) at steady state Z = 100%; LHSV = 2 h -1 ; T inlet = 200°C; U wall = 3 W/m 2 /K; H 2 /HC ratio = 1000 Nm 3 /m 3 and Pt = 150 barg.Gray zone inside dotted lines: unstable zone according to stationary van Heerden criterion.a) Seven steady states evaluated (Tab.2), b) zoom close to a limit point.

Figure 7 Eigenvalues
Figure 7 Eigenvalues spectra obtained from numerical and analytical Jacobians (Case 2).

TABLE 2 Reference
Schweitzer et al. / Thermal Stability of Gas Oil Hydrotreating Processes: Numerical Issues of the Matrix-Eigenvalue Approach J-M Schweitzer et al. / Thermal Stability of Gas Oil Hydrotreating Processes: Numerical Issues of the Matrix-Eigenvalue Approach 781