Dynamic analysis on methanation reactor using a double-input-multi-output linearized model☆

Xingxing Li ,Jiageng Li ,Bolun Yang ,*,Yong Zhang

1 Department of Chemical Engineering,State Key Laboratory of Multiphase Flow in Power Engineering,Xi'an Jiaotong University,Xi'an 710049,China

2 Northwest Research Institute of Chemical Industry,Xi'an 710049,China

Keywords:Coke oven gas methanation Synthetic natural gas Double-input-multi-output linearized system Dynamic analysis Limit cycle behavior

A B S T R A C T A double-input-multi-output linearized system is developed using the state-space method for dynamic analysis of methanation process of coke oven gas.The stability of reactor alone and reactor with feed-effluent heat exchanger is compared through the dominant poles of the system transfer functions.With single or double disturbance of temperature and CO concentration at the reactor inlet,typical dynamic behavior in the reactor,including fast concentration response,slow temperature response and inverse response,is revealed for further understanding of the counteraction and synergy effects caused by simultaneous variation of concentration and temperature.Analysis results show that the stability of the reactor loop is more sensitive than that of reactor alone due to the positive heat feedback.Remarkably,with the decrease of heat exchange efficiency,the reactor system may display limit cycle behavior for a pair of complex conjugate poles across the imaginary axis.

1.Introduction

The research on alternative fuel is highly attractive with the concern of environmental protection and growing energy demand.Natural gas,mainly containing methane,is considered as one of the most promising alternative fuels because of its excellent combustion characteristics such as high octane number,good antiknock property and less emission[1,2].How ever,it is reported that natural gas will be exhausted in 40-60 years,so a sustained effort has been made to convert various coal-based and biomass feedstock into synthetic natural gas[3,4].As a by-product in the coke production,the coke oven gas(COG),most of which is either discharged directly into the atmosphere or committed to flames,is considered as an important hydrocarbon resource to the synthesis of natural gas through methanation process[5].

The methanation of COG most likely occurs in adiabatic fixed bed reactors with nickel-based catalysts.The inherent distribution nature in the reactor(variables change with axial position),interaction of heat,mass transport and reaction,and thermal inertia of catalysts may give great challenge to the study on dynamic behavior,where significant dynamic lags,multistability,and sudden loss of stability could be induced[6,7].The COG methanation reactor is also sensitive to disturbances of CO concentration and temperature,because COG from different coking furnaces varies in components and enormous reaction heat will bring about a large temperature rise in the adiabatic bed.Furthermore,the feed gas is always preheated with the product gas from the view point of energy cascade utilization,which could introduce a feedback mechanism and make the reactor stability worse and the anti-interference ability poor[8].Therefore,it is necessary to analyze these dynamic characteristics in the methanation process of COG.

Up to now,many methods have been adopted to address the dynamic model of distributed parameter system.With the finite difference method,Kordabadi and Jahanmiri[9],and Rahimpour[10]carried out a dynamic simulation on a dual-stage reactor for methanol synthesis in the face of catalyst deactivation.Baldea and Daoutidis[11]computed the singular perturbation of an autothermal fixed bed reactor for methane steam reforming,where reactor extinction may arise.With the orthogonal collocation method,Margarida and Rosa[12]simulated the start up and the wrong way behavior of a fixed bed reactor with two distinct zones for methanol oxidation.Continillo et al.[13]characterized the dynamic behavior of a self-ignition reactor by considering the diffusion effect,where the period-doubling bifurcation could befound.Other methods such as polynomial approximation,wavelet-like collocation method and piecewise linearization method are also developed to deal with the distributed parameter system[14-18].These methods can describe the transient behavior of the fixed bed reactors well,while lacking of a complete system identification that characterizes the input-output relationship of a reactor loop.Furthermore,double perturbation may also occur,but the counteraction and synergy effects caused by simultaneous variation of concentration and temperature were seldom taken into account in above researches.

In this paper,a double-input-multi-output(DIMO)linearized system based on a simple dynamic model is developed for dynamic analysis of a methanation reactor with feed-effluent heat exchanger(FEHE).The stability of the reactor loop and the reactor alone is compared through the system transfer functions,and further demonstrated by the responses of reactor variables when the inlet temperature and CO concentration are perturbed separately or simultaneously.In addition,the essential behavior,including fast concentration response,slow temperature response,wrong way behavior and limit cycle behavior,is also revealed for the study of counteraction and synergy effects through the concentration and temperature variations.The main purpose of our work is not to predict the industrial case with great numerical accuracy,but rather to yield qualitative and further insight into possible phenomena,which could contribute to the system identification and optimal control for the methanation process of COG.

2.Modeling

2.1.Basic assumptions

Except for hydrogen,carbon monoxide and carbon dioxide,there is also some methane and a little nitrogen in COG[19].The following independent reactions are involved in the methanation system of COG

Fig.1 show s the methanation system,which consists of three adiabatic fixed bed reactors with interstage cooling.The diameter of the reactor is 1 m,and the lengths of the three beds are 2.8,1.9 and 2.5 m.Besides,a feed-effluent heat exchanger is added for surplus heat utilization,and its heat exchange area is 11.51 m2.

To develop the dynamic model,w e make a few assumptions.

(1)In the COG methanation system,methanation of CO2may be restrained in the presence of CO[20,21],so methanation of CO mainly takes place in the first stage and continues in the middle stage,and methanation of CO2mainly occurs in the last stage.The CO2concentration is relatively low(about 2%)and the reaction heat of CO2methanation is less than that of CO methanation,so only the CO methanation in the first stage is considered.

Fig.1.Reactor system of COG methanation.

(2)The axial dispersion and radial gradients in the reactor are negligible,the gas mixture moves with a uniform velocity in the reactor and Ni/Al2O3catalyst is uniformly distributed.Catalyst specifications are listed in Table 1.

(3)The methanation reaction of CO is controlled by the rates of intrinsic kinetics and external diffusion,which can be described as

where ηeis the effectiveness factor,Da is the Damkohler number expressing the ratio of reaction rate to mass transfer rate,K is the adsorption coefficient,and k0and kgrepresent the reaction rate constant and mass transfer coefficient,respectively.The intrinsic kinetic parameters are taken from Xu and Froment[22].kgis evaluated as follows[23]

Here Sh is the Sherwood number and is written as

where Dg=8.34×10-6T1.75/p is the diffusion coefficient of CO.

2.2.Nonlinear dynamic model

Based on the mass and heat balances,the dynamic performance of the reactor can be described by tw o partial differential equations

where v is the gas velocity,ε is the voidage of catalyst bed,Cpgand Cpsrepresent the specific heat of gas mixture and catalyst particles,respectively,and ΔH is the reaction heat of CO methanation.

Table 1Physical parameters of Ni/Al2O3 catalyst in the methanation reactor

With the feed-effluent heat exchanger,the thermal efficiency εhis defined as[24]

where Tinand Toutare the inlet and outlet temperature of the reactor,and Tfstands for the temperature of the feed gas.

2.3.Local linearization of the nonlinear dynamic model

The partial derivatives in nonlinear Eqs.(12)and(13)are discretized along reactor length by a first-order finite difference approximation,which is well suited for the smooth reactor profile and successfully used in many studies[25,26],so the dynamic model can be converted to a set of ordinary differential equations as follows,with the catalyst bed divided into 28 cells and 56 difference equations generated.

The nonlinearity arises in Eqs.(15)and(16)due to the Arrhenius dependency of the reaction rate on temperature.Expanding these nonlinear terms in a Taylor series and truncating high order terms,following linear ordinary differential equations are obtained

where Cisand Tisare the CO concentration and temperature along the reactor,which must be the steady state values.Letting the time variation terms in Eqs.(15)and(16)to be zero,the New ton iteration method can be adopted to calculate these steady state values.In this case,changes of concentration and temperature along the reactor caused by tiny disturbances around the steady state can be monitored using Eqs.(17)and(18).

2.4.Double-input-multi-output linearized model

Rewrite Eqs.(17)and(18)with vectors,a standard state-space model[27]is developed as follow s

where state vector x is the changes of CO concentration and temperature along the reactor,input vector u consists of inlet CO concentration and temperature perturbations,output vector y represents the changes of CO concentration and temperature at one quarter(1/4 N),half(1/2 N),three quarters(3/4 N)and the outlet(N)of the reactor,matrices A,B,C and D are the state,input,output and feed forward matrix,respectively.Through this state-space form,a double-input-multi-output system is attained to reveal the nature of the methanation reactor by indicating the “input-state-output”relationship,where the input-output transfer function G(s)is calculated by Eq.(24)[28],with G(s)to be a 8×2 order matrix.With the heat recycle,the state-space model can be reformed using the MATLAB program:[Ah,Bh,Ch,Dh]=cloop(A,B,C,D,8,2),and the schematic diagram of the DIMO system is illustrated in Fig.2.The method for calculating the system transfer function is the same as Eq.(24),and the transfer function for each variable could be described by Eq.(25).

Fig.2.Schematic diagram of the double-input-multi-output linearized system.

3.Results and Discussion

3.1.Operating point under steady state condition

Fig.3 show s the profiles of reactor temperature and CO concentration under steady state condition.The results are in fair agreement with the data from the industrial plant,where the reactor temperature ranges from 250 °C to 550 °C,and the mole faction of CO at the reactor outlet is within 1%.

Fig.3.Pro fi les of CO concentration and reactor temperature under steady state condition.●temperature;■CO concentration.

3.2.Analysis on poles and zeros

Fig.4 show s the poles in transfer function Gl,k′(s)of the DIMO system(roots of the denominator in the transfer function),where the thermal efficiency of the FEHE is 0.425Notice that the map divides the poles into two distinct groups.Specifically,some are further left in the complex plane,with the remainder close to the imaginary axis,which are amplified in the top right corner of Fig.4.Every pole has negative real part to confirm the stable system,according to Lyapunov stability theory[29].Moreover,the further the dominant poles(closest to the imaginary axis)are away from the imaginary axis,the more stable the reactor runs.Table 2 is the comparison for the dominant poles of the transfer function without and with heat feedback.It show s that the former has better stability than the latter,which will be identified by the subsequent dynamic responses to different disturbances.opposite direction firstly to the inlet temperature disturbance,and then return tow ard the direction of the inlet temperature disturbance.How ever,because the even number of positive real zeros exist in the transfer function G5,2′,the response trend of temperature at one quarter(1/4 N)of the reactor is consistent with the direction of the inlet temperature perturbation,and the same dynamic behavior occurs when the transfer function has no positive real zeros(minimum phase system).Based on these facts,simultaneous increase of inlet temperature and CO concentration may eliminate the undershoot in the temperature response(such as G8,1′and G8,2′),while simultaneous increase of inlet temperature and decrease of CO concentration may exacerbate the undershoot phenomenon.These analysis results can be verified by subsequent dynamic responses to different disturbances.

Fig.4.Poles map of the transfer function with heat feedback.

Table 2Dominant poles of the transfer functions without and with heat feedback

3.3.Analysis on dynamic behavior

In this study,the temperature and CO concentration at the reactor inlet are varied for creating different disturbances,with the purpose focused on three aspects.Firstly,we compare the response speed in the reactor loop and in there actor alone to certify the impact of heat recycle on system stability.Secondly,we are interested in the transient behavior with single or double perturbation of temperature and CO concentration,which may relate to the positive real zeros of the linearized system and provide further understanding of the counteraction and synergy effects through the concentration and temperature variations.Finally,the instability in the reactor must be eliminated to prevent sudden accident,so the critical state that continuous oscillations appear requires particular attention.Different behaviors based on above evaluation is discussed below.

3.3.1.Change in inlet CO concentration

The first disturbance is a step change of the inlet CO concentration,owing to the changeability of COG component from different coke ovens.Fig.5 presents the variation of CO concentration versus time with the inlet CO concentration slightly perturbed from 13.80 to 14.80 mol·m-3.State variables at one quarter(1/4 N),half(1/2 N),three-quarters(3/4 N)and the outlet(N)of the reactor are analyzed.In the reactor without thermal feedback,we find that immediately after the additional increase of carbon monoxide,within the gas residence time in the catalyst bed,a carbon monoxide concentration profile is established rapidly in about 5 s.How ever,there action heat causes the temperature profile to develop gradually in about 1000 s,as show n in Fig.6,much more slowly than the change of CO concentration.These effects are referred as fast concentration response follow ed by slow temperature response[30].There are a couple of reasons for the phenomenon.First,the change of CO concentration travels through the reactor with substantially the same velocity as the gas mixture.Second,an appreciable amount of heat is absorbed by the catalyst,the heat capacity of which is far greater than that of gas mixture.Therefore,the temperature equilibrium between the catalyst and gas phase is established much more slowly.

As presented in Figs.7 and 8,the changes of CO concentration and temperature inside the reactor with heat exchanger show the same trend as that without heat exchanger.The stable time of the concentration profile changes little(Figs.5 and 7),but the temperature in the reactor loop takes much more time(5000 s)to reach stable than that in the reactor alone(Figs.6 and 8).Theoretically,it is mainly because the

Table 3 shows the positive real zeros of the transfer functions with heat feedback,which could be related to the dynamic responses(similar to the transfer functions without heat feedback).The positive real zeros in the transfer function indicate that the system could be regarded as the non-minimum phase system.If the number of positive real zeros is odd,undershoot will appear initially in the step response(G6,2′,G7,2′and G8,2′).In this case,the responses of temperature at half(1/2 N),three-quarters(3/4 N)and the outlet(N)of the reactor will act in theinsertion of the feed-effluent heat exchanger makes the dominant poles of the system transfer function shift toward the imaginary axis,which shall result in lower level stabilization,as shown by the above stability analysis.The obvious fluctuations in the temperature response curves(Fig.8)also support the conclusion.

Table 3Positive real zeros of the transfer functions with heat feedback

Fig.5.Fast concentration response in reactor alone with the increase of inlet CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.____1/4 N;----1/2 N;...........3/4 N;N.

Fig.6.Slow temperature response in reactor alone with the increase of inlet CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.____1/4 N;----1/2 N;...........3/4 N;N.

Fig.7.Fast concentration response in reactor with heat fe_e_d_b_ack to the increase of inlet CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.1/4 N;----1/2 N;...........3/4 N;N.

Fig.8.Slow temperature response in reactor with heat feedback to the increase of inlet CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.____1/4 N;----1/2 N;...........3/4 N;N.

3.3.2.Change in inlet temperature

From the viewpoint of system control,it is important to learn how the reactor behaves after the disturbance of inlet temperature.Therefore,an increase of the inlet temperature(523 K→543 K)is imposed on the methanation system.As for the reactor without feed-effluent heat exchanger,the variation of reactor temperature is illustrated in Fig.9.The reactor could be divided into two regions with the split point locating probably at one quarter of the reactor.The temperature in the upstream part(less than one quarter of the reactor)increases directly,because higher inlet temperature increases the reaction rate and heat generation rate correspondingly.Conversely,the temperature in the downstream part of the reactor decreases initially and increases later.The interpretation is that more carbon monoxide is converted to methane in the upstream part of the reactor for the increased reaction rate,while the reaction heat is absorbed by the catalyst in this part and cannot arrive at the downstream part promptly.As a consequence less carbon monoxide reaches the downstream part where reduced reaction heat causes a drop in temperature.This state lasts until the downstream catalyst bed is heated enough by the upstream catalyst bed so that the effect of the increase in temperature on reaction rate overshadow s the competitive effect of the decreased CO concentration on the reaction rate.This counterintuitive phenomenon is referred as“wrong way behavior”or reverse response behavior[31].

Fig.9.Wrong way behavior in reactor alone with the increase of inlet temperature from 523 K to 543 K.____1/4 N;----1/2 N;...........3/4 N;

For the reactor with feed-effluent heat exchanger,the rising of temperature inside the reactor would encounter difficulties with the inlet temperature disturbance,as show n in Fig.10.The increasing trend becomes slow and even a slide may occur.It is something like“secondary response”,and also vividly described by the name of“snowball effect”[32].This can be explained with reference to the presence of an inverse response due to the inlet temperature disturbance combined with the positive heat feedback caused by the feed-effluent heat exchanger.The increase of the inlet temperature,which makes an initial decrease of temperature at reactor outlet,will suffer a mitigation effect through the heat exchanger,so several wavelike bumps in the temperature/time curves can be found.

Fig.10.Wrong way behavior in reactor with heat feedback to the increase of inlet temperature from 523 K to 543 K.____1/4 N;----1/2 N;...........3/4 N;N.

3.3.3.Simultaneous change in inlet temperature and CO concentration

Fig.11.Temperature responses in reactor with heat feedback to simultaneous increase of inlet temperature from 523 K to 543 K and CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.____1/4 N;----1/2 N;...........3/4 N;N.

It is probable that misfortunes may not come singly.The coupling effect of CO concentration and temperature disturbances on dynamic behavior of the methanation system is discussed.As show n in Fig.11,the wrong way behavior disappears with simultaneous increase of the inlet temperature(from 523 K to 543 K)and CO concentration(from 13.80 mol·m-3to 14.80 mol·m-3).This is mainly because the gas mixture moves fast in the reactor,consequently the increased CO arriving at the downstream part of there actor makes up for the consumption of CO due to the increase of inlet temperature.On the contrary,simultaneous increase of inlet temperature and decrease of CO concentration will force the split point of the reverse response in the reactor to move upstream(less than 1/4 N)and exacerbates there verse response behavior,as illustrated in Fig.12.These temperature responses are consistent with that in the zeros analysis of the DIMO system.

Fig.12.Temperature responses in reactor with heat feedback to simultaneous increase of inlet temperature from 523 K to 543 K and decrease of inlet CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3.____1/4 N;----1/2 N;...........3/4 N;N.

3.3.4.Limit cycle behavior

Due to long period running,the smudge and corrosion may influence the efficiency of the feed-effluent heat exchanger.With the thermal efficiency εhfurther decreasing to 0.345,the inlet and outlet temperature of the reactor will decline.Fig.13 is the CO concentration and temperature pro files along the catalyst bed,where the inlet temperature is 444.8 K and the outlet temperature is 724.4 K.Simultaneously,the dominant poles of the system transfer function with heat feedback is found locating at the imaginary axis(±0.0106),so under this critical state,temperature and concentration oscillations could be observed with the temperature and concentration disturbances.

Fig.13.Pro fi les of CO concentration and reactor temperature under critical state.●temperature;■CO concentration.

Fig.14.Phase diagram of temperature and concentration oscillations with the increase of inlet temperature from 444.8 K to 464.8 K under critical state.

Fig.14 is the phase diagram of temperature and concentration oscillations for a sudden increase of the inlet temperature(444.8 K→464.8 K).Firstly the temperature in the upstream of the reactor increases[Fig.14(a)]and that in the downstream declines[Fig.14(b)-(d)],then through the positive heat feedback,the decrease of the outlet temperature due to the inverse response can even neutralize the increase of inlet temperature,from the evidence that the minimum temperature variation locating at one quarter of the catalyst bed is below 0 K[Fig.14(a)],so the oscillations will turn up.Interestingly,in the migration process of the temperature throughout the reactor,the temperature synergy action also exists,which makes the temperature variation locating at about half of the reactor reach the maximum[Fig.14(b)].

It is worth mentioning in Fig.15 that the temperature and concentration oscillations still occur for a simultaneous increase of CO concentration and inlet temperature,despite the addition of CO may make the reverse response go away.This explanation perhaps lies that the methanation process runs at a relatively low inlet temperature under such a critical state condition.When the inlet temperature jumps suddenly,the change of reaction rate is more sensitive at lower temperature than at higher temperature.Therefore,the same effect that more carbon monoxide is depleted in the upstream part of the reactor occurs,which cannot be compensated by the increase of CO concentration.In other words,the wrong way behavior still exists.Nevertheless,this refers to a slight disturbance of CO concentration implemented,since dramatic changes are not suitable for the developed linear model.

4.Conclusions

A double-input-multi-output system on the basis of the state-space strategy is developed for the dynamic simulation of methanation reactor with feed-effluent heat exchanger.The stability of the reactor loop is more sensitive than the reactor alone due to the positive heat feedback,indicating that the dominant poles of the transfer function are useful for the investigation of the system stability.The transient behavior of temperature and CO concentration in the reactor is also exhibited,which are closely related to positive real zeros of the system transfer function.Step change in inlet CO concentration will lead to fast concentration response follow ed by slow temperature response.Increase in inlet temperature will lead to reverse response,which will be intensified by a simultaneous decrease of inlet CO concentration but will disappear by an increase of inlet CO concentration.Furthermore,with the decrease of heat exchange efficiency,temperature and concentration oscillations are induced for a pair of complex conjugate poles across the imaginary axis,which may result in process instability or even reactor extinction.Therefore,corresponding control strategy is highly desirable for safe operation of reactor.

Nomenclature

A state matrix

a specific area of catalyst particle,m2·m-3

B input matrix

C output matrix

CiCO concentration at different positions of the reactor,mol·m-3

CisCO concentration at different positions in steady state,mol·m-3

Cpgspecific heat of gas mixture,J·kg-1·K-1

Cpsspecific heat of catalyst particle,J·kg-1·K-1

D feed forward matrix

Dgdiffusion coefficient,m2·s-1

Da Damkohler number

dpdiameter of catalyst particle,m

E activation energy,kJ·mol-1

G(s) transfer function without feedback

ΔH reactor heat of CO methanation,kJ·mol-1

I identity matrix

Fig.15.Phase diagram of temperature and concentration oscillations with simultaneous increase of inlet temperature from 444.8 K to 464.8 K and CO concentration from 13.80 mol·m-3 to 14.80 mol·m-3 under critical state.

K adsorption coefficient,Pa

Kl,ksystemic gain

k reaction rate constant of CO methanation,m3·(kg cat)-1·s-1

k0pre-exponential factor,m3·(kg cat)-1·s-1

kggas phase mass transfer coefficient,m·s-1

L reactor length,m

N number of sections that the reactor is divided into

p operating pressure,MPa

R reactor diameter,m

Rgmolar gas constant,J·mol-1·K-1

Re Reynolds number

Sc Schmidt number

Sh Sherwood number

Tftemperature of feed gas,K

Titemperature at different positions of the reactor,K

Tininlet temperature of reactor,K

Tistemperature at different positions in steady state,K

Toutoutlet temperature of reactor,K

v superficial velocity of gas mixture,m·s-1

Δz reactor infinitesimal,m

ε voidage of catalyst bed

εhthermal efficiency

μ viscosity,Pa·s

ρgdensity of gas mixture,kg·m-3

ρsdensity of catalyst particle,kg·m-3

Subscripts

g gas mixture

i different position

s catalyst particles