Modelling of Complex Dynamics in Fixed Bed Reactors

(1)

M A T S A P P E L B L O M

Master's Degree Project

Stockholm, Sweden 2005

(2)

Abstract

Chemical reactors are a part of modern industry and the catalytic tubular ﬁxed bed reactor examined in this work is an important reactor for chemicals production.

In this work two different types of models for the reactor are studied; a pseudohomo-geneous model and a heteropseudohomo-geneous model. The goal is to find differences in behaviour between these two types of reactor models and explain these.

In a real reactor there exists two phases, a solid catalyst and a ﬂuid reactant. Both these phases are in the pseudohomogeneous model treated as a single phase, a pseudoﬂuid. In the heterogeneous model the two phases are treated separately.

When comparing these types of models a few structural differences exist, void fraction, heat exchange between two phases, and heat dispersion in the phases, and all of these will affect the behaviour of the models differently.

The models are studied using bifurcation analysis and linear analysis. Bifurcation theory is used to ﬁnd and track diﬀerent solutions depending on a certain parameter and to get a good overall picture of a system’s solutions and their type, steady state or sustained oscillation.

Linear analysis is used to study linearization around a speciﬁc solution and to determine stability and frequency dependency.

It is found that the concept of void fraction in the reactor model affects the behaviour only as a time scaling, while the concept of interfacial heat exchange affects the stability. The distribution of heat dispersion between phases has a significant impact on the reac-tion behaviour. Feedback is determined as the main cause for instabilities and oscillative solutions.

(3)

1 Introduction

A tubular fixed bed reactor is a reactor in which a chemical reaction transforms one substance to another. A reactant (gas or fluid) is fed into the reactor. The reactant then passes through the fixed bed of solid catalyst where the desired reaction takes place and the desired product comes out at the other end.

One popular way to model this kind of reactor is to make a pseudohomogeneous model where the fluid reactant and the solid catalyst together are treated as one single phase, a pseudofluid, which inherits properties from both the fluid and the solid, see [6], [8], [9], and [10].

Another way to model the reactor is to consider the two diﬀerent phases as separate phases. This will yield a model called a heterogeneous model, see [8].

When a model is made, different assumptions will yield different models and different types of models for the same system. In the example with the catalytic reactor studied in this work, a lot of interesting behaviour can be found in certain models of the reactor. Some of these behaviours can increase the effectiveness of the reactor by a considerable amount.

Examining such possibilities in the models will yield important information to design-ing future reactors. The models used are derived from physical considerations and their equations and, since the reactor itself is a truly complex system, a number of assump-tions and simpliﬁcaassump-tions have to be made in order to get a model simple enough to be of practical use.

When these assumptions are made, one has to consider where it is realistic to include or omit a certain physical property and where it is reasonable to do a simplification. This work aims at examining how different assumptions in the modelling will affect the qualitative behaviour of the model.

This work will focus on two families of models introduced above: the

pseudohomo-geneous model and the heteropseudohomo-geneous model. Moreover this work will further focus

on regions or modes of work where the reactors are assumed to work more eﬃciently by using physical feedback to increase the eﬀect. Even further, to modes where oscillating solutions are present.

In many applications a stable working mode is desirable over an oscillating mode, but in the case with this type of reactor recent work has shown that it is possible to extract a higher conversion rate from the reactor in oscillating modes, thus making the reactor more eﬃcient.

1.1 Problem Deﬁnition

The goal with this work is to examine how diﬀerent modelling assumptions aﬀect the predictions of the reactor model.

It is important to know how the modelling assumptions affect the models for a number of reasons. It may give a hint of what examinations on one type of model may produce and the effects that does not appear in that model. It may tell which assumptions that have a great influence on the models and which assumptions that does not. It may point out different things that has to be considered when designing a real reactor.

The main focus of this examination will be on the differences that comes from the choice of modelling the reactor as a single phase (pseudofluid) pseudohomogeneous model or that of a multi phase (fluid/solid) heterogeneous model which takes into consideration a more detailed description of the different constituting parts of the reactor.

(6)

The analysis will primarily focus on how stability is aﬀected by modelling assumptions, the oscillative behaviours and how they appear.

The models will be examined by means of bifurcation analysis and the results from that will in turn be examined in a linear feedback analysis.

Bifurcation analysis will be used to determine the differences between the models and give information about the model differences and how they are governed by the different assumptions.

The linear feedback analysis uses more traditional control theory to examine the dif-ferences found in the bifurcation analysis and explain why the diﬀerences occur.

To be able to compare the diﬀerent kinds of reactor models a heterogeneous model has be derived based on already existing pseudohomogeneous models.

(7)

1.2 Report outline

Section 2 will give a physical description of the reactors used in this work. It will be described what they are used for and how they work. This section will also describe methods used to increase eﬀectiveness of the reactor and the problems this is related with. Section 3 will consider the models used; the pseudohomogeneous and the heterogeneous models. A derivation of these models will be presented and it will be shown that these models may give the same results under certain operating conditions. In this section the method used for spatially discretizing the model will also be introduced.

Section 4 will deal with the issues concerning the bifurcation analysis part of the work. An introduction to bifurcation theory will ﬁrst be presented. The tests and test results from the bifurcation analysis will then be presented and discussed.

Section 5 will present the linear analysis. Here the methods used for examining the system properties with respect to a linear open loop model of the system will be treated. The results from this analysis will also be discussed in this section.

Section 6 is the conclusion of the report and it summarizes the results found and discuss them further. This section will also present possibilities for future work with obtained results as a starting point.

(8)

2 Fixed Bed Reactor

In the physical modelling of a reactor, a number of properties are accounted for. First the structure of the reactor. The reactor studied here is a so called tubular ﬁxed-bed

catalytic reactor. The reactor consists of a cylindrical tube where one end is the inlet

and the other end is the outlet. Into this tube a gas or fluid1 called the reactant is fed. The purpose is to make this fluid react and make some specific chemical compound in the reactor which can be collected at the outlet.

This chemical reaction is triggered by a so called catalyst which is already present in the reactor. The catalyst in most cases consists of a tightly packed bed of ceramic pellets, or a net of some metal, hence the name fixed-bed. The fluid reactant flows through this porous bed of the solid catalyst and the reactions take place.

Except from the presence of the reactant and the catalyst, a certain amount of heat needs to be added to the system in order to trigger the reaction. This is mainly solved by heating the gas or ﬂuid before it is fed into the reactor.

This type of reactor is common in many applications, for example it is used to purify the exhaust gases from combustion engines. One of its largest areas of use is to produce chemicals used to make crop fertilizer which is a vital part of all modern agricultural food production.

The reaction itself can be either endothermal or exothermal, meaning it either consumes heat or produces heat as it is converts the chemical compounds. The processes studied here are all exothermal, so they produce heat at the same time as they produce the desirable chemicals. This exothermal property of the reaction can be proﬁted upon. The heat the process releases can be used to heat the ﬂuid reactant at the inlet. This can be done in a number of ways but quantitatively the results are similar whatever method is used. Here the simple solution of letting a part of the warm outlet gas be fed back and mixed with the input gas and this combination is then fed into the inlet.

The normal use of these reactors is to use them in a steady state mode, where the temperature and concentration profiles are constant in time. The profiles for temperature and concentration are not flat throughout the length of the reactor due to the fact that the reaction takes place along the whole length of the reactor. Positive feedback may however make the steady state of the reactor unstable. It may also give raise to oscillations in the solution or in some cases even chaotic behaviour.

A common model used for this kind of reactor is the pseudohomogeneous model. This model is however quite simplistic and does not consider the eﬀects a two phase heteroge-neous model may introduce.

When feedback causes oscillative behaviour, the use of a multi phase model may predict diﬀerent behaviour than the pseudohomogeneous single phase model of the reactor.

1_{Either a gas or a ﬂuid can be used and both will show similar properties. For the remainder of this}

(9)

3 Modelling

The derivation of the two kinds of reactor models will be presented here. The goal is to get two models that are as similar as possible with respect to the mathematical structure, in order to facilitate comparison. As a starting point for this an already existing model for a pseudohomogeneous model is used. This model comes from [3].

To be able to get the heterogeneous model as structurally similar as possible to the pseudohomogeneous model the derivation of the heterogeneous model will be based on the derivation of the pseudohomogeneous model.

The reason for including a derivation of the pseudohomogeneous model as well as the derivation of the heterogeneous model is that the derivation procedures then can be compared but it is also to oﬀer some insight into the concept of pseudohomogeniety. The derivation of the pseudohomogeneous model presented here is an adaptation of the model in [3], reinforced with my own comments and explanations.

3.1 Derivation of a pseudohomogeneous model of a tubular

reac-tor.

The derivation of the pseudohomogeneous model begins with a heterogenous model, where the heat balances are treated separately for the two diﬀerent phases, the solid catalyst and the ﬂuid.

For this, three balance equations are needed: the mass balance(for the fluid only, the solid is fixed and do not change), the heat balance for the fluid and the heat balance for the solid. For simplicity this is done in one dimension only.

The models presented below will start with a very simple one. When this model is derived more and more properties will be added to the models.

3.1.1 Mass balance

The model assumes that the reactor is a cylinder with length L_R, radius R_R, and a cross sectional area A_R = R2π. Looking at a small volume of the reactor from x = x_i−1 to

x = x_i= x_i−1+ ∆x the change in mass can be described by equation (1), where C is the concentration of the fluid and F is the flow of the fluid.

∆V_idCi

dt = F (Ci−1− Ci) (1)

This is divided by ∆x and diﬀerentiated with respect to x and the following expression is obtained. dV dx ∂C ∂t =−F ∂C ∂x (2)

The volume derivative of equation (2) above is constant an can thus be replaced with the area A_R as shown in equation (3).

dV dx = AR

dx

dx = AR (3)

It is preferable to make the variables of the model dimensionless. First the time and space variables are made dimensionless by the introduction of the dimensionless time and space coordinates τ and z.

τ = F

V_Rt; z = x

(10)

It takes one dimensionless time unit τ for the ﬂuid to pass through the reactor and the length of the reactor is now one dimensionless length unit z. The new coordinates inserted in equation (2) give equation (5).

A_R F V_R ∂C ∂τ =− F L_R ∂C ∂z (5)

Which reduces to:

∂C

∂τ +

∂C

∂z = 0 (6)

3.1.2 Heat balance

For the solid in a small element ∆x of the reactor, the heat balance can be expressed as in equation (7).

∆V_siρ_sc_psdTsi

dt = 0 (7)

In the reactor both the solid and the ﬂuid has to share the volume provided by the dimensions of the reactor. this is done as in equation (8), where ε is a constant describing the distribution of the two phases.

V_s+ V_f= V_R; V_f = εV_R; V_s= (1− ε)V_R (8)

After dividing with ∆x and diﬀerentiating we get:

A_sρ_sc_psdTs

dt = 0 (9)

The heat balance equation for the ﬂuid is presented in equation (10). ∆V_{f i}(ρ_fc_pf)dTf i dt =−F (ρfcpf)(Tf (i−1)− Tf i) (10) A_f(ρ_fc_{f s})∂Tf ∂t =−F (ρfcpf) ∂T_f ∂x (11)

Adding the heat balance equations for the solid and substituting for equation (8) gives equation (12). εA_R(ρ_fc_pf)∂Tf ∂t + (1− ε)AR(ρscps) ∂T_s dt =−F (ρfcpf) ∂T_f ∂x (12)

Assuming that the heat transferring between the solid and the ﬂuid is inﬁnitely good, such that T_f = T_s= T we get equation (13).

A_F(ερ_fc_pf+ (1− ε)ρ_sc_ps)∂T

∂t =−F (ρfcpf) ∂T

∂x (13)

A coefficient called the Lewis number is introduced next. The Lewis number describes the heat capacity of the pseudofluid relative the flow speed. If only a fluid is examined a corresponding Lewis number would be one; the heat will flow with the fluid at the same speed. However when a solid material is present this solid material will absorb and contain some of the heat. Normally the solid has a much higher heat capacity than the fluid. This results in a scenario where the heat is transported much slower than the mass in the pseudofluid.

L_e= ερfcpf+ (1− ε)ρscps

(11)

Using the Lewis number presented above in equation (14) and the dimensionless time and space coordinates in equation (4), the resulting equation is equation (15).

L_e∂T

∂τ +

∂T

∂z = 0 (15)

3.1.3 Extending the model

The model can be extended with consideration to more physical details, yielding more complex models. The next step is to add the reaction term. The reaction is an exothermic reaction which will reduce the concentration of the ﬂuid and at the same time produce heat in the catalyst, or in the pseudohomogeneous ﬂuid.

Responsible for producing the desired substance and the heat is the reaction. Chemical reactions are complex processes that can be hard to describe and can be modelled as complex but often quite simple models will give accurate approximations of the reality.

r_A= Cnke−ERT (16)

The reaction in this work is modelled as an Arrhenius process and is given in equation (16), where r_A describes how much reaction there is in the reactor, and is dependant on both temperature and concentration. Excitation issues are not considered here, mostly because they are very small at the working temperatures, so reaction is always present no matter how small the concentrations or temperatures are in this model.

In the model the reaction will be present in both the heat balance and in the mass balance because the reaction aﬀects them both by reducing the concentration of the reac-tant while producing the desired substance and at the same time produce heat due to the exothermal nature of the reaction. The terms where the reactions appear are called the reaction term.

With the reaction term the mass balance equation (1) will become the mass balance equation in equation (17).

∆V_idCi

dt = F (Ci−1− Ci)− rA∆Vi (17)

This equation becomes equation (18) after diﬀerentiating and substituting as in the pre-vious case. ∂C ∂t + ∂C ∂x =−rA V_R F (18)

The heat from the reaction is produced in the solid catalyst so the heat balance equation for the solid, equation (7) will with the reaction term added turn into equation (19), where ∆H is the heat generated in the reaction. The heat balance equation for the ﬂuid, equation (10), is unchanged because no heat is produced in the ﬂuid.

∆V_siρ_sc_psdTsi

dt = ∆HrA∆Vi (19)

When the heat balance equations are added together and the variables are substituted this yields the pseudohomogeneous model for the heat balance presented in equation (20).

L_e∂T ∂τ + ∂T ∂z = V_R∆H F ρ_fc_pfrA (20)

(12)

3.1.4 Dimensionless concentration and temperature

In the pseudohomogeneous model described above in equation (17) and (20), the concen-tration C and the temperature T are not dimensionless. It is desirable to make them dimensionless so that the model easily can be adapted and used to examine a wide range of reactors of diﬀerent dimensions without the need for a speciﬁc model for each reactor.

The reactor model is supposed to work in the neighborhood of some nominal values for concentration and temperature. These are the values of the input ﬂuid and are called

C_{F 0}∗ and T_{F 0}∗ . Using the nominal values we make a new dimensionless variables for the concentration which is the deviation from the nominal value normalized with the nominal value as can be found in equation (21).

α = C

∗ F 0− C

C_{F 0}∗ (21)

This new dimensionless concentration will give the expression for the mass balance as:

∂α ∂τ + ∂α ∂z =−rA V_R C_{F 0}∗ F (22) Now introducing: Da = VR F₀C_{F 0}∗ r 0 A; r0A= rA(CF 0∗ , TF 0∗ ); RA= r_A r0_A (23)

With the use of the quantities described in equation (23) and (4) the mass balance equation (22) can be rewritten as in equation (24). The Damkholer number Da is a measure widely used in chemical modelling and hydrodynamics that describes the relation between ﬂuid speed and reaction rate.

∂α ∂τ +

∂α

∂z = (1− f)DaRA (24)

Making the temperature dimensionless is similar to making the concentration dimension-less as in equation (21), but some scaling is done to make the equations more tractable. This is found in equation (25).

θ = 1 β T_{F 0}∗ − T T_{F 0}∗ ; β = B T_{F 0}∗ ; B = ∆HC_{F 0}∗ ρ_fc_pf (25)

Considering this the heat balance equation for the pseudohomogeneous model becomes equation (26). L_e∂θ ∂τ + ∂θ ∂z = (1− f)DaRA (26) 3.1.5 Cooling

Heat is transported from the reaction to the walls of the reactor which are held at a constant temperature T_H. This phenomena can be described in the model by adding a term to the heat balance equation for the ﬂuid.

∆V_{f i}(ρ_fc_pf)dTf i

dt =−F (ρfcpf)(Tf (i−1)− Tf i) + U OH(TH− Tf i)∆x (27)

It is done in equation (27), where O_H is the circumference of the cross-section of the reactor and U is a heat transfer coefficient between the wall and the fluid. This coefficient

(13)

U is in this model assumed to be constant even though it in reality is strongly dependant

of the ﬂow speed.

Diﬀerentiating this equation and adding the heat balance for the solid will lead to a pseudohomogeneous heat balance equation, which is found in equation (28).

L_e∂T ∂τ + ∂T ∂z = V_R∆H F ρ_fc_pfrA+ U A_H F ρ_fc_pf(TH− T ) (28)

The coeﬃcients from the cooling term are collected in one constant δ. This is shown in equation (29).

δ = U AH

F₀ρ_fc_pf (29)

Substituting for all the coeﬃcients found in equations (4), (25), and (23) gives the heat balance equation found in equation (30).

L_e∂θ ∂τ +

∂θ

∂z = (1− f)DaRA+ (1− f)δ(θH− θ) (30)

3.1.6 Dispersion

Dispersion in concentration and heat is the next phenomena added to the model.

In a fluid mass and heat can be subjected to diffusion, in a solid only the diffusion of heat exists of course. Since a pseudofluid is used, and it behaves like a fluid, both mass diffusion and heat diffusion will be present. Apart from these two obvious ways of dispersion similar effects can be contributed to some turbulence in unaccounted directions in the actual system, these effects will not be modelled in the one dimensional model. These effects are accounted for in the corresponding diffusion terms and they are therefore called dispersion terms instead of diffusion since the terms describes all dispersion not only the diffusion.

Taking the mass balance equation in equation (17) and adding the dispersion term yields equation (31), where D_M is a dispersion coeﬃcient describing the dispersion of mass.

∆V_idCi

dt = F (Ci−1− Ci)− rA∆Vi−

D_MA_R

∆x (Ci−1− Ci) (31) If this equation is diﬀerentiated with respect to x and the time and space variables are substituted with the dimensionless ones in equation (4). The resulting equation is equation (32). ∂C ∂t + ∂C ∂x = D_MA_R L_RF ∂2C ∂z2 − rA V_R F (32)

The coefficients in the diffusion term are collected in a new coefficient that is called the Peclet number. The Peclet number is described in equation (87) and it is a measurement of how resistive the process is to dispersion. Large Peclet number means little dispersion.

P e_M = LRF

D_MA_R (33)

Now substituting the rest of the variables yields the mass balance equation (34).

∂α ∂τ + ∂α ∂z = 1 P e_M ∂2α ∂z2 + (1− f)DaRA (34)

(14)

Dispersion of the heat works in a similar way but with a diﬀerent Peclet number P e_H which has to do with the resistivity to heat dispersion. The heat balance equation for the pseudohomogeneous reactor is equation (35).

L_e∂θ ∂τ + ∂θ ∂z = 1 P e_H ∂2θ ∂z2 + (1− f)DaRA+ (1− f)δ(θH− θ) (35) 3.1.7 Boundary conditions

To get a set of solvable diﬀerential equations the boundary conditions have to be stated. These are stated separately for the mass and the heat balance equations and for the outlet and inlet as well.

For the mass balance at the inlet the conditions is as in equation (36), and based on the nominal input, the feed back, and the dispersion. The boundary condition for the outlet is based on the dispersion alone, and can be found in equation (37).

α(0, τ ) = (1− f)α_{F 0}+ f α(1, τ ) + 1 P e_M ∂α ∂z|z=0 (36) ∂α ∂z|z=1 = 0 (37)

The heat balance at the inlet and at the outlet are calculated based on the same principles as for the mass balance, but regarding heat instead of mass, so concentration is instead temperature et.c. The boundary conditions for the heat balance are found in equation (38), for the inlet, and equation (39), for the outlet.

θ(0, τ ) = (1− f)θ_{F 0}+ f θ(1, τ ) + 1 P e_M ∂θ ∂z|z=0 (38) ∂θ ∂z|z=1 = 0 (39) 3.1.8 Collected model

The ﬁnal heat and mass balance equations are here collected.

∂α ∂τ + ∂α ∂z = 1 P e_M ∂2α ∂z2 + (1− f)DaRA (40) L_e∂θ ∂τ + ∂θ ∂z = 1 P e_H ∂2θ ∂z2 + (1− f)DaRA+ (1− f)δ(θH− θ) (41)

(15)

3.2 Derivation of a heterogeneous model of a tubular ﬁxed bed

reactor

In the heterogeneous model the two diﬀerent phases are considered as separate phases; both the solid catalyst and the ﬂuid reactant. In this model only a heat balance will be considered for the solid.

In this model the volume of the reactor is divided between the solid and the ﬂuid according to the constant ε as shown in equation (42) where V is the volume and the indices_f,_s, and_R in order refer to ﬂuid, solid, and the total of the reactor.

V_R= V_f+ V_s; V_f = εV_R; V_s= (1− ε)V_R (42)

The catalyst in a real rector consists of a bed closely packed pellets through which the fluid flows. Because of that the area distribution between solid and fluid may not be the same in any two cross sections of the reactor. To be able to use the cross sectional areas of the solid and the fluid in computation the average areas over all possible cross sections is defined as A_f and A_s. The relations are shown in equation (43), A_R is constant throughout the whole reactor and L_Ris the total length of the reactor.

A_R= A_f+ A_s; A_f = εA_R; A_s= (1− ε)A_R; V_j = A_jL_R; j = f, s, R (43) The ε in the equations is the void fraction. This is a dimensionless variable that account for area and material properties and its easy to calculate this for the cooling where the areas and materials are known, in the heterogeneous case with the interfacial heat ex-change, the heat exchange coefficient is a lot more complex. Depending on the structure and composition of the solid material the interfacial area can be almost anything. This interfacial area can be very difficult to calculate in the different nodes so an assumption is made that the heat exchange coefficient is the same everywhere in the reactor, this factor may also be highly sensitive to turbulence and such.

Now it is time to state the heat balance equations for the solid and the ﬂuid. Starting with the solid this model contains three things that contribute to changes in the heat. These three are: the reaction, which produces heat, the dispersion of heat along the reactor due to conduction, and the heat exchange with the ﬂuid phase.

This equation is stated in equation (44), where ρ is the density, c_pthe heat capacity, T the temperature, ∆H the heat produced by the reaction, r_Athe amount of reaction, U_I the interfacial heat exchange coefficient (assumed constant), O_I the interfacial circumference (the total length of the average borders between the areas of the solid and the fluid in all possible cross sections of the reactor), D_Hthe heat dispersivity coefficient, and x the axial coordinate along the reactor.

∆V_siρ_sc_psdTsi dt = ∆HrA∆Vf i+ UIOI(Tf i− Tsi)∆x− DHsAsρscps dT_s dx|i−1− dT_s dx |i (44) In equation (44), above, the reason for using ∆V_{f i} in the reaction term is because the reaction r_A is dependent of heat and mass in the ﬂuid phase only.

The heat balance for the fluid looks a bit different. In this model the causes of change in heat for the fluid are: convection (due to mass transport), heat exchange with the solid, dispersion of heat along the reactor, and heat exchange with the surrounding shell (held at a constant temperature). this heat balance is shown in equation (45), where F i the flow trough the reactor, U_H is the heat exchange coefficient with the shell (assumed constant),

(16)

and O_H is the circumference of the cross section. ∆V_{f i}ρ_fc_pfdTf i dt =−F ρfcf s(Tf (i−1)− Tf i) + UIOI(Tsi− Tf i)∆x −DHfAfρfcpf dTf dx|i−1− dTf dx|i + U_HO_H(T_H− T_{f i})∆x (45) The two heat balance equations above, equations (44) and (45), can easily be diﬀerentiated by division by ∆x, and letting ∆x→ 0. This and the fact that dV_j/dx = A_jfor j = s, f, R yields the heat balance equations (46) and (47) below.

A_sρ_sc_ps∂Ts ∂t = ∆HrAAf+ UIOI(Tf− Ts) + DHsAsρscps ∂2T_s ∂x2 (46) A_fρ_fc_pf∂Tf ∂t = F ρfcf s ∂T_f ∂x +UIOI(Ts−Tf)+DHfAfρfcpf ∂2T_f ∂x2 +UHOH(TH−Tf) (47)

The next step in the derivation is to rescale the length and time coordinates and make the dimensionless as a method for making the model applicable to any reactor of this kind. The new length coordinate z is scaled so the total length of the rector is one. The new time coordinate τ is scaled so that the time it takes for mass to pass through the reactor is one. How this is done is shown in equation (48).

z = x

L_R; τ = F t

V_f (48)

The heat balance equations (46) and (47) are, using the dimensionless coordinates, rewrit-ten into equations (49) and (50).

F A_sρ_sc_ps V_f ∂T_s ∂τ = ∆HrAAf+ UIOI(Tf− Ts) + D_HsA_sρ_sc_ps L2_R ∂2T_s ∂z2 (49) F A_fρ_fc_pf V_f ∂T_f ∂τ =− F ρ_fc_pf L_R ∂T_f ∂z + UIOI(Ts− Tf) + D_HfA_fρ_fc_pf L2_R ∂2T_f ∂z2 + UHOH(TH− Tf) (50) Now the temperatures have to be made dimensionless. In order to do that we first define the temperature T_{F 0}∗ which is the nominal temperature of the input fluid. The rescaling to dimensionless temperatures uses this nominal temperature value both for the fluid and the solid. The parameter β is a scaling constant depending of the exothermicity B and will be properly defined later on. In equation (51) the dimensionless temperature is found.

θ_j =Tj− TF 0∗ βT_{F 0}∗ ; j = s, f, H; β = B T_{F 0}∗ =⇒ θj= T_j− T_{F 0}∗ B (51)

The heat balance equations can now be rewritten as in equations (52) and (53).

F A_sρ_sc_psB V_f ∂θ_s ∂τ = ∆HrAAf+ UIOIB(θf− θs) + D_HsA_sρ_sc_psB L2_R ∂2θ_s ∂z2 (52) F A_fρ_fc_pfB V_f ∂θ_f ∂τ =− F ρ_fc_pfB L_R ∂θ_f ∂z+UIOIB(θs−θf)+ D_HfA_fρ_fc_pfB L2_R ∂2θ_f ∂z2 +UHOHB(θH−θf) (53) Considering a nominal reaction rate r_A0, which is the reaction rate at the nominal tem-perature T_{F 0}∗ and nominal concentration C_{F 0}∗ a dimensionless reaction rate R_A can be obtained as in equation (54) where the exothermicity is also stated.

r_A0 = r_A(T_{F 0}∗ , C_{F 0}∗ ); R_A=rA

r_A0; B =

∆HC_{F 0}∗

(17)

Inserting this in the heat balance equations (52) and (53) and rearranging gives equations (55) and (56). ∂θ_s ∂τ = V_fR_Ar0_A F C_{F 0}∗ ρ_fc_pf ρ_sc_ps ε (1− ε)+ U_IO_IL_R F ρ_fc_pf ρ_fc_pf ρ_sc_ps ε (1− ε)(θf− θs) + D_HsA_f F L_R ∂2θ_s ∂z2 (55) ∂θ_f ∂τ =− ∂θ_f ∂z + U_IO_IL_R F ρ_fc_pf (θs− θf) + D_HfA_f F L_R ∂2θ_f ∂z2 + U_HO_HL_R F ρ_fc_pf (θH− θf) (56)

By introducing a few dimensionless coefficients like the Peclet number and the Damkohler number and a feedback number, as well as dimensionless heat exchange coefficients for the flow to the equations they become a lot simpler to grasp. The dimensionless parameters are presented in equation (57).

P e_j = F LR D_HjA_f; j = s, f ; f = F F₀; Da = V_fr0_A F₀C_{F 0}∗ ; δk= U_kO_kL_R F₀ρ_fc_pf ; k = I, H (57) Le = ερfcpf+ (1− ε)ρscps ρ_fc_pf (58)

This gives the concluding equations (59) and (60).

∂θ_s ∂τ = (1− f)DaRA ε Le− ε + (1− f)δI ε Le− ε(θf− θs) + 1 P e_s ∂2θ_s ∂z2 (59) ∂θ_f ∂τ =− ∂θ_f ∂z + (1− f)δI(θs− θf) + 1 P e_f ∂2θ_f ∂z2 + (1− f)δH(θH− θf) (60) Mass balance

The mass balance is only considered for the ﬂuid phase. For convenience C_f is therefore written C. ∆V_{f i}dCi dt = F (Ci−1− Ci)− rA∆Vi− DMAf dC_s dx |i−1− dC_s dx |i (61) Diﬀerentiating yields: A_f∂C ∂t =−F ∂C ∂x − rAAf− DMAf ∂2C ∂x2 (62)

Inserting dimensionless coordinates yields:

F A_f V_f ∂C ∂τ =− F L_R ∂C ∂z − rAAf+ D_MA_f L2_R ∂2C ∂z2 (63)

The dimensionless concentration (or rather lack of concentration):

α = C

∗ F 0− C

C_{F 0}∗ (64)

and rearranging yields:

∂α ∂τ =− ∂α ∂z − V_f F C_{F 0}∗ rA+ A_fD_M F L_R ∂2α ∂z2 (65)

substituting the rest of the dimensionless variables gives:

∂α ∂τ =− ∂α ∂z − ε(1 − f)DaRA+ ε P e_M ∂2α ∂z2 (66)

(18)

3.2.1 Boundary Conditions

The boundary conditions for the ﬂuid phase are stated exactly as in the pseudohomoge-neous case, and are found in equations (67),(68),(69), and(70).

α(0, τ ) = (1− f)α_{F 0}+ f α(1, τ ) + 1 P e_M ∂α ∂z|z=0 (67) ∂α ∂z|z=1 = 0 (68) θ_f(0, τ ) = (1− f)θ_{f F 0}+ f θ_f(1, τ ) + 1 P e_M ∂θ_f ∂z |z=0 (69) ∂θ_f ∂z |z=1= 0 (70)

For the solid phase the boundary conditions are based on prohibiting heat dispersion at the ends. The conditions are found in equation (71), and (72).

∂θ_s

∂z|z=0= 0 (71)

∂θ_s

∂z|z=1= 0 (72)

3.2.2 Collected model

The ﬁnal balance equations for the heterogeneous model are here collected.

∂θ_s ∂τ = (1− f)DaRA ε Le− ε+ (1− f)δI ε Le− ε(θf− θs) + 1 P e_s ∂2θ_s ∂z2 (73) ∂θ_f ∂τ = − ∂θ_f ∂z + (1− f)δI(θs− θf) + 1 P e_f ∂2θ_f ∂z2 + (1− f)δH(θH− θf) (74) ∂α ∂τ = − ∂α ∂z − (1 − f)DaRA+ 1 P e_M ∂2α ∂z2 (75)

(19)

3.3 Validation of heterogeneous model through condensation into

the pseudohomogeneous model

Under certain conditions the heterogeneous model and the pseudohomogeneous model will have the same asymptotical behaviour. The assumed conditions for this to happen are the following:

• Since both solid and ﬂuid phases are treated as one, the heat exchange between the

two phases has to be instant and absolute, i.e. the heat exchange coeﬃcient has to be inﬁnitely large (or very large when used in practical experiments).

• All of the volume has to be treated as a ﬂuid phase which implies that the void

fraction of the packed bed, ε, is set to one.

When those conditions are applied to the heterogeneous model, it can be shown by ex-periment that the behaviour is identical to that of the pseudohomogeneous model. It is however necessary to show that this is true, and that it can be validated through pure mathematical methods.

This validation starts from the equations for the heterogeneous model, as derived earlier with the time variable based on the mass residing time, presented below in equations (76), (77), and (78). ∂α ∂τ = − ∂α ∂z − (1 − f)DaRA+ 1 P e_M ∂2α ∂z2 (76) ∂θ_s ∂τ = (1− f)DaRA ε Le− ε+ (1− f)δI ε Le− ε(θf− θs) + 1 P e_s ∂2θ_s ∂z2 (77) ∂θ_f ∂τ = − ∂θ_f ∂z + (1− f)δI(θs− θf) + 1 P e_f ∂2θ_f ∂z2 + (1− f)δH(θH− θf) (78)

Since the pseudohomogeneous model only consists of one heat balance equation the main focus should be on reducing the number of heat balance equations for the heterogeneous model. The key to this is to consider the one term the two heat balance equations of the heterogeneous model, equations (77) and (78), have in common; namely the heat exchange term in equation (79) below (Notice the order of θ_sand θ_f compared to the original heat balance equations).

(1− f)δ_I(θ_s− θ_f) (79)

The heat exchange term (79) can be isolated in the two heat balance equations resulting in the equations (80) and (81), presented below.

(1− f)δ_I(θ_s− θ_f) = −Le− ε ε ∂θ_s ∂τ + (1− f)DaRA+ Le− ε ε 1 P e_s ∂2θ_s ∂z2 (80) (1− f)δ_I(θ_s− θ_f) = ∂θf ∂τ + ∂θ_f ∂z − 1 P e_f ∂2θ_f ∂z2 − (1 − f)δH(θH− θf) (81)

Equations (80) and (81) reveals a quite obvious way of combining the heat balance equa-tions from the heterogeneous model. Carried out, this combination result in equation (82). −Le−ε ε ∂θ∂τs + (1− f)DaRA+Le−εε P e1s ∂2θs ∂z2 = = ∂θf ∂τ + ∂θf ∂z −P e1f ∂2θf ∂z2 − (1 − f)δH(θH− θf) (82)

(20)

Returning to the assumptions, the assumption that the heat exchange between the two phases is inﬁnitely good implies that the temperatures of the two phases are equal, i.e.

θ_s = θ_f. Even though this knowledge is of no use as it is, it in turn implies a couple of important equalities stated in equation (83) below.

θ_s= θ_f = θ =⇒ ∂θs ∂u = ∂θ_f ∂u = ∂θ ∂u =⇒ ∂2θ_s ∂u2 = ∂2θ_f ∂u2 = ∂2θ ∂u2 (83)

Introducing the equalities obtained in (83) to equation (82) and rearranging it yields equation (84) below. Le ε ∂θ ∂τ =− ∂θ ∂z + (1− f)DaRA+ Le− ε εP e_s + 1 P e_f ∂2θ ∂z2 + (1− f)δH(θH− θf) (84)

The second assumption to consider is that the new pseudohomogeneous phase ﬁlls the total volume of the reactor. The ε parameter is set to one to satisfy that requirement. Doing this, equation (76) and (84) can be collected into the pseudohomogeneous model described by equation (85) and (86).

∂α ∂τ = − ∂α ∂z − (1 − f)DaRA+ 1 P e_M ∂2α ∂z2 (85) Le∂θ ∂τ = − ∂θ ∂z + (1− f)DaRA+ 1 P e_H ∂2θ ∂z2+ (1− f)δH(θH− θf) (86)

The Peclet number for this model can be one of the following in equation (87) depending on if the heat diﬀusion in the solid phase is neglected or not.

P e_H = _{P e}

f; No heat diﬀusion in solid

Le−1 P es + 1 P ef −1

; Heat diﬀusion in solid

(87)

This shows that the heterogeneous model gets not only the same behaviour but can be condensed into the same form if the required assumptions are made.

(21)

3.4 Discretization

The models obtained are all partial differential equations. To be able to use these partial differential equations in simulations they have to be discretized and approximated into ordinary differential equations. The most usual and logical way to do this is to discretize the spatial variable and to keep the time in the resulting ordinary differential equations. This is done with Taylor series expansion, as shown in equations (88) and (89).

u_i− u_i−1 ∆z = ∂u ∂z|i− ∆z 2 ∂2u ∂z2|i+ O(∆z 3₎ ₍₈₈₎ u_i+1− 2u_i+ u_i−1 ∆z2 = ∂2u ∂z2|i+ O(∆z 3₎ ₍₈₉₎

The third and higher order terms are eliminated and the equations are rearranged to get equations (90) and (91). ∂u ∂z|i= u_i− u_i−1 ∆z + ∆z 2 u_i+1− 2u_i+ u_i−1 ∆z2 (90) ∂2u ∂z2|i= u_i+1− 2u_i+ u_i−1 ∆z2 (91)

Using these the equations are transformed to:

dαi dτ = − αi − αi−1 ∆z − ε(1 − f)DaRA + ε P eM − ∆z 2

αi+1 − 2αi + αi−1

∆z2 (92) dθsi dτ = (1 − f)DaRA ε2 Le − ε+ (1 − f)δI ε Le − ε(θfi − θsi) + ε P es

θs(i+1) − 2θsi + θs(i−1)

∆z2 (93)

dθfi

dτ = −

θfi − θf(i−1)

∆z + (1 − f)δI (θsi − θfi) +

ε P ef − ∆z 2

θf(i+1) − 2θfi + θf(i−1)

∆z2 + (1 − f)δH (θH − θfi)(94)

If the dispersion terms in the discretized equations are closely examined it can be observed that there exists two dispersion terms in the balance equations where it prior to the discretization only existed one. One of these dispersion terms corresponds to the original dispersion whereas the other is introduced by the discretization and corresponds to the numerical dispersion. It can sometimes be convenient to chose a discretization level so that the numerical dispersion is equal to the actual dispersion and therefore cancel it out leaving the equations free of dispersion terms.

(22)

4 Bifurcation Analysis

The purpose of this section is to perform a bifurcation analysis of the models. First a brief introduction to what bifurcation analysis is and how it is applied in this work is oﬀered.

A description of the different tests that are carried out in the bifurcation analysis to find and verify the differences in the models as well as the actual testing and the results they give are also treated in this section.

4.1 Overview

Bifurcation analysis is a method to examine how the solutions of a system are dependent of one or more of the system parameters.

Many systems consisting of differential equations may have more than one solution, and these solutions may also differ in characterization; steady state solutions, periodic solutions, and in some cases even chaotic solutions. The steady state and periodic solutions can also be either stable or unstable, which will contribute to the diversity of different possible solutions.

Bifurcation analysis is a method for surveying and analyzing different solutions for a system of differential equations, with special respect to a specific, or a few specific parameters of the system.

Bifurcation analysis uses continuation methods to find solutions to a system, see [1]. The strength of these methods lies in the property of parameter dependency of the so-lutions. The number of solutions and their classification may in many cases be highly dependant of the parameters defining the system.

4.2 Parameter dependency

To yield some insight into how parameter dependency works a simple illustrative example is oﬀered.

A common feed-back loop is considered. It consists of the open loop system where the output signal is ampliﬁed and added to the input signal. Moreover, assume that we can change the gain of the feed-back as we like. The solution and the type of solution will change when the gain is altered. The feed-back gain is chosen to be the bifurcation

parameter.

Now at diﬀerent values for the gain the behaviour of the solutions will be diﬀerent. For a value of zero for the gain the solution will be that of the open loop system. When the gain is increased the solution will change. In many cases the solution will go from a steady state solution to a point where the solution splits into a self sustaining oscillation and an unstable steady state solution.

Now these solutions can be plotted in a bifurcation diagram for each value of the bifurcation parameters. Here at ﬁrst the lower values of the bifurcation parameter the steady state solution will be plotted as a curve. But at some point called a bifurcation

point the solution will bifurcate into two diﬀerent solutions, one unstable steady state

solution and one stable oscillation, before the stable oscillation in turn becomes an unstable oscillation.

A bifurcation point may give rise to many diﬀerent behaviours and are classiﬁed ac-cordingly. To illustrate this an example from this work is presented, a bifurcation diagram for the reactor is shown in Figure 1. The α(1) is the conversion at the outlet and the f is a parameter describing the amount of feedback in the reactor.

(23)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 f α (1)

Figure 1: Bifurcation diagram for the reactor.

Here we start to look at a low value for the bifurcation parameter f , for which it exists only one solution for the conversion rate α(1) and this is a stable steady state solution. If the line consisting of these stable steady state solution is followed we come to a point where the line turns back into an unstable steady state solution. The point where this happens is called a limit point. The steady state curve can be tracked further to another limit point completing a S-shaped curve. There are two other bifurcation points along this curve. These two are called Hopf bifurcations and at these points the solution curve is branching, or bifurcating, into a periodic solution. The periodic solution can also be tracked and here it is represented by the min and max values of the amplitude, thus the double curve represents one solution only.

4.3 How is bifurcation theory applied to this work

Since this work emphasizes diﬀerences between behaviour in a number of models and modelling assumptions, bifurcation analysis will work as a good mean to observe and examine these diﬀerences.

Bifurcation diagrams may give important information on the qualitative diﬀerences in the behaviour of the diﬀerent models, and the impact of some of the constituting parameters.

The periodic solutions are of special interest and the bifurcation analysis will strongly contribute to the information about this kind of behaviour since the frequency of said solutions is obtained through bifurcation analysis.

When oscillative behaviour is studied the focus will be on the Hopf bifurcation points. This is because once a periodic solution has appeared it then follows the same pattern in

(24)

amplitude and shape, and the period of the oscillation does not change signiﬁcantly over the course of the particular periodic solution.

The study of the Hopf bifurcation points will suﬃce for describing in what regions of the solution space that periodic solution is present, and will also give useful information of the frequency, or equivalently, period.

4.4 Numerics

The theory of bifurcation analysis is based on a number of numerical procedures to track the solutions of a system without the need to actually solve the system for all parameter values. The actual numerical procedures will not be considered here, they can be found in [1].

4.4.1 XPPAUT and AUTO

To cope with the substantial amount of numerical calculations used for bifurcation analysis a program named AUTO is used, see citem. To solve differential equations, examine different solutions and provide starting solutions for the bifurcation analysis a program named XPPAUT is used, see [12]. This program has a built-in version of AUTO which will, except for offering an easy transition for data between a solver and a bifurcation tool also provides me with a graphical interface and instant rendering of the bifurcation diagrams.

(25)

4.5 Results

The differences in structure between the pseudohomogeneous model and the heterogeneous model are, apart from the obvious fact that there is a different number of constituting equa-tions, that some new constants, and concepts, and corresponding terms in the equations are present in the heterogeneous model but not in the pseudohomogeneous one. Looking only at the constituting equations the differences can be summarized:

• The concept of a void fraction has been added to the heterogeneous model to handle

the distribution of the reactor volume between the solid and the ﬂuid phases.

• Since the heterogeneous model consists of two diﬀerent phases with separate heat

balance, a heat exchange term governs the heat distribution and exchange between the solid and the ﬂuid.

• Two diﬀerent phases also need two diﬀerent terms for handling the spatial heat

distribution in each phase separately, for this the heat dispersion terms are used. Except for the differences stated above all the other terms in the pseudohomogeneous model are also present in the heterogeneous model. These other terms are the same for the two different models and should produce no differences in behaviour.

With the structural diﬀerences in mind a number of possible tests can be devised to quantitatively examine the impact of modelling assumptions concerning the structural diﬀerences.

One preferable property of these tests are that they give results that can be compared to similar results in the diﬀerent models. It is also desirable to be able to examine one, or as few as possible, traits at a time, so the underlying reasons for diﬀering results are clear. Since it has been shown that the pseudohomogeneous model and the heterogeneous model will behave in the same way under certain conditions the pseudohomogeneous model will be the reference model in these tests.

Starting from a condition where the different models give rise to exactly the same behaviour, one single difference in the properties of the heterogeneous model can be studied and compared to the reference behaviour of the pseudohomogeneous model. At the point where the impact of the single traits are fully explored, multiple traits can be examined based on the reference and the results from the studies of the differences in single traits.

4.5.1 Parameter values used in simulations

parameter value parameter value

Da 0.06 n 1.5 Le 1000 β 0.9 P e_M 100 γ 12 P e_f 100 nodes 33 δ_H 0 4.5.2 Void fraction

One test is done for the impact of the void fraction. The void fraction determines how the ratio between the volume of the ﬂuid and solid phases, and may have an impact on the behaviour of the reactor models.

In the pseudohomogeneous case all of the reactor volume is filled with a pseudofluid consisting of both the solid and the fluid phase, behaving like a fluid with respect to mass

(26)

transport but has the combined heat capacity of the two phases. In the heterogeneous case the fluid and the solid are treated as two different phases where the solid behaves like a solid and the fluid like a fluid.

The solid phase has a much larger heat capacity than the fluid. This fact is accounted for in the pseudohomogeneous model when the Lewis number was introduced. Since the heterogeneous model also uses the Lewis number to determine the heat capacity in both phases, the heat capacity of both the solid and fluid phase together is the same as the heat capacity for the pseudofluid used in the pseudohomogeneous model, the heat capacity of the heterogeneous model is thus unaffected by the void fraction.

The ε in the reaction term for the solid phase comes from the fact that the amount of reaction is dependent of the amount of reactant present in the reactor, and the reactant is only present in the ﬂuid phase so the amount of reactant is proportional to the void fraction.

The time in the models are scaled so that it takes one time unit for the mass to pass through the reactor.

The void fraction for the reactor is examined separately. This is achieved by choosing a heat exchange coefficient δ_I that is very large. This has two effects: first, the heat exchange will be almost instantaneous as it is assumed to be in the pseudohomogeneous model, so the behaviour of the heterogeneous model and the effects of the void fraction can more easily be compared to the behaviour of the pseudohomogeneous model; second, since the heat exchange term consists of a product of ε and δ_I, this will allow this term to be relatively unchanged for a wide range of void fraction values (still very good and almost instantaneous heat exchange for ε > 0.01).

No heat conduction in the solid assumed in this step, and the whole term responsible for this is omitted. This too makes the behaviour close to the pseudohomogeneous model for comparison issues.

The bifurcation diagram for the model is obtained for different values of ε. The results are that the steady state solutions are the same for all values of ε, how these bifurcation diagram looks can be seen in figure 2. For values where the heat exchange can be considered instantaneous all the bifurcation points are located at the same place with respect to the bifurcation parameters for different values of the void fraction. The one thing that differs with ε is the frequency of the periodic oscillation. The oscillations are slower for low values of ε as expected, since the time constant is scaled with the void fraction. The period of the oscillation at the Hopf bifurcation point where the steady state solution loses stability and the periodic oscillation commences can be seen in figure (3). Data points are acquired at steps of 0.05 for ε from 0.05 to 1. No data points acquired below 0.05 in order to avoid the effects of a heat exchange no longer instantaneous, which will move the Hopf bifurcation point and redound to the period time of the oscillations. For illustrative purposes the best fit exponential function is added to the picture.

Here the presumed eﬀects of the heat contained in the solid can be seen. The process gets a slower behaviour with a small void fraction.

Since the ratio between reactor dimensions and flow is kept constant in the trans-formations to dimensionless variables and coefficients this phenomena occurs. Physically speaking this means that the same amount of fluid has to pass through a smaller hole in the same time as before resulting in a higher speed.

The void fraction will only work as a time scaling as was expected. The reason for this is that it reduces the amount of ﬂuid that will pass through the reactor, while the time scaling is based on the whole reactor volume and the ﬂow is constant.

(27)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 f α(1) δΙ=1000000

Figure 2: Bifurcation diagram for the void fraction examination. This diagram is the same for diﬀerent values of the void fraction ε.

4.5.3 Heat exchange

The effect of the heat exchange coefficient is another trait that is examined, and in this case tests can be designed so that all of the other parameters and conditions are kept constant and equal to the conditions that holds for the pseudohomogeneous model. When the heat exchange is very large the heterogeneous process should behave almost exactly as in the pseudohomogeneous case as shown in the validation. The differences in behaviour are expected to occur at values of the heat exchange coefficient that are relatively low.

The heat exchange coefficient determines the speed of the heat exchange in the hetero-geneous model. If this parameter goes to infinity the heat exchange is instantaneous and the temperature in the solid is equal to the temperature of the fluid at any given time.

It has been shown in section 3.3 that the behaviour of the heterogeneous model is the same as for the pseudohomogeneous model under certain conditions. To obtain these conditions, and similarity in behaviour the void fraction ε has to be one and the heat exchange coeﬃcient has to be inﬁnity (or a very large value).

Starting from these conditions a first analysis of the effects of the heat exchange coef-ficient can be made. During this analysis the void fraction is kept constant at ε = 0.5 to isolate the effect of the assumption that there are two different phases and there exists a transfer of heat between these two. No heat conduction in the solid phase is assumed.

As can be observed in figure 4 and 5 the Hopf bifurcation points change their positions when the heat exchange coefficient changes. For a lower value of the coefficient, and thus a lower heat exchange speed between the two phases, the Hopf bifurcation points move to a lower f , closing in on the upper limit point, and finally meet just above the upper limit

(28)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.5 1 1.5 2 2.5x 10 4 eps Period

Figure 3: Period of the oscillation at the upper Hopf bifurcation point for diﬀerent values of the void fraction.

point cancelling each other out making all periodic solutions disappear. The steady state solution curve in the bifurcation diagram is unchanged with the heat exchange coefficient. How the Hopf bifurcations moves and the region where there exists an oscillation is best analyzed in a two parameter diagram where, naturally the heat exchange coefficient δ_I is chosen as the one of the parameters and the feedback f as the other, since the feedback is used as main parameter in most bifurcation analysises in this work. These diagrams are found in figure 6 and figure 7. Here it can be observed that the distance between the Hopf points gets smaller as δ_I gets smaller and that there exists values for δ_I where no periodic solutions are present at all. By looking at the possible locations for the Hopf bifurcation points it can be observed that the curve has its lowest value at a point separated from the lowest allowed δ_I. This is because the lower Hopf bifurcation point passes along the steady state solution curve in the bifurcation diagram and the two Hopf bifurcations conjoin in the upper branch. As can be seen in figure 8. Apart from the earlier analysis of the isolated void fraction the heat exchange coefficient contribute to the behaviour of the model and has a stabilizing effect on the solutions compared to a pseudohomogeneous model, when a low heat exchange is present between the two phases.

The period of the oscillation of the solution at the upper Hopf bifurcation point is also monitored, this is shown in ﬁgure 9. Here we can observe that the period is dependent of at what feed back value f the Hopf point is located when the Hopf point is moved by changing the heat exchange coeﬃcient. It can be noted that in most of the range the

(29)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 f α(1) δΙ=1000000

Figure 4: Bifurcation diagram for a very large heat exchange coeﬃcient (δ_I = 1000000).

dependance is almost linear. The explanation for this eﬀect lies in the physical properties of the pseudohomogeneous model. A low heat exchange will stabilize the reactor by allowing oscillations for a smaller region of f the oscillations will also be slower.

4.5.4 Void fraction combined with heat exchange coeﬃcient

The results from the two tests above can with little extra eﬀort be collectively examined by looking whether a certain result in one of the tests will change independently of the results in the other test.

At diﬀerent values of δ_I and ε the behaviour of the solution is identical to when ε = 0.5 except for the frequency of the oscillations which is scaled according to the results in the analysis of the void fraction.

The conclusion is that the void fraction does not affect the stabilizing effect of the heat exchange coefficient examined above other than changing the time constant.

4.5.5 Heat conduction distribution between solid and ﬂuid phase

The process of examining the third big diﬀerence in the architecture of the models is a bit more tedious and needs changing in the equations since the heat conduction in the solid phase only can be examined at full extent if such a term is added in the equations. When no heat conduction is allowed in the solid the heterogeneous model will behave according to the results obtained from the earlier tests on the two other structural diﬀerences in the design. With a heat exchange however a few other things are possible as objects for examination.

(30)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 f α(1) δΙ=10

Figure 5: Bifurcation for a small heat exchange coeﬃcient(δ_I = 10). Note that no Hopf bifurcations exists and the whole of the upper branch is stable.

The distribution of the heat conduction between the two phases is the priority in this test. The same over-all heat dispersion that is used in the pseudohomogeneous model is preferably used. To be able to do this the Peclet number used in the pseudoﬂuid in the pseudohomogeneous model has to be expressed in the two new Peclet numbers used for the ﬂuid and solid respectively.

By introducing two different phases to the model the heat dynamics may look different from the pseudohomogeneous model. In the fluid phase as in the pseudofluid of the pseudo-homogeneous model we got heat transportation due to the mass flow through the reactor and heat dispersion due to conduction, mass dispersion and turbulence in unaccounted directions. In the solid phase only conduction is present.

In the pseudohomogeneous model the heat dispersion is made up from both the heat dispersion in the fluid and the heat conduction in the solid. In the heterogeneous model one simple condition is to let all the heat dispersion be accounted for in the fluid phase, thus dropping the diffusion term in the equations for the solid. This has the advantages that it makes the model less complex and makes its behaviour to be closer to the behaviour of the pseudohomogeneous model, which is good for comparative reasons. Another assumption is to allow conduction in both the fluid and the solid to get use of the extra amount of dynamical behaviour introduced by the heterogeneous model.

The analysis of the effect of heat conduction in the solid starts with some basic condi-tion. In the verification of the heterogeneous model an expression where the Peclet number from the pseudohomogeneous model could be expressed as a function of the Peclet numbers for the fluid and solid phase in the heterogeneous model was obtained. This expression is

Modelling of Complex Dynamics in Fixed Bed Reactors

M A T S A P P E L B L O M

Master's Degree Project

Stockholm, Sweden 2005

Abstract

Contents

1

Introduction

1.1

Problem Deﬁnition

1.2

Report outline

2

Fixed Bed Reactor

3

Modelling

3.1

Derivation of a pseudohomogeneous model of a tubular

reac-tor.

3.2

Derivation of a heterogeneous model of a tubular ﬁxed bed

reactor

3.3

Validation of heterogeneous model through condensation into

the pseudohomogeneous model

3.4

Discretization

4

Bifurcation Analysis

4.1

Overview

4.2

Parameter dependency

4.3

How is bifurcation theory applied to this work

4.4

Numerics

4.5

Results