An inverse problem in estimating simultaneously the non-linear reaction rates for a fixed-bed reactor

Accepted Manuscript An Inverse Problem in Estimating Simultaneously the Non-linear Reaction Rates for a Fixed-Bed Reactor Cheng-Hung Huang, Bo-Yi Li PII: DOI: Reference:

S0307-904X(14)00516-2 http://dx.doi.org/10.1016/j.apm.2014.10.032 APM 10185

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

11 March 2013 19 June 2014 13 October 2014

Please cite this article as: C-H. Huang, B-Y. Li, An Inverse Problem in Estimating Simultaneously the Non-linear Reaction Rates for a Fixed-Bed Reactor, Appl. Math. Modelling (2014), doi: http://dx.doi.org/10.1016/j.apm. 2014.10.032

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An Inverse Problem in Estimating Simultaneously the Non-linear Reaction Rates for a Fixed-Bed Reactor

by

*Cheng-Hung Huang and Bo-Yi Li Department of Systems and Naval Mechatronic Engineering National Cheng Kung University Tainan, Taiwan 701, R.O.C.

*: corresponding author e-mail: [email protected] Keywords : Nonlinear Problem, Inverse Problem, Fixed-Bed Reactor, Reaction Rate Estimation, Exothermic Reaction

ABSTRACT The Conjugate Gradient Method (CGM) is adopted in this study to estimate simultaneously the unknown conversion and temperature-dependent reaction rates for a fixed-bed reactor by using interior measurements of conversion and temperature distributions. The validity and accuracy of this inverse fixed-bed reactor problem was examined using the simulated exact and inexact conversion and temperature measurements in the numerical experiments. Numerical simulations indicate that the estimation of the conversion and temperature-dependent reaction rates can be obtained quickly on an Intel Xeon Core 2 2.00 GHz personal computer and the estimations are still reliable when measurement errors are included since in test case 1 the maximum average errors for the estimated R1(r,z) and R2(r,z) are obtained as 0.85 %, and 1.4 %, respectively, in test case 2 they are 2.58 % and 3.98 %, respectively.

NOMENCLATURE A1

: reciprocal of Peclet number for effective radial mass transfer

A2

: reciprocal of Peclet number for effective heat transfer

Bi

: Biot number

ERR1,ERR2 : relative errors J

: functional defined by equation (3)

J1' , J2'

: gradient of functional defined by equation equations (15) and (17)

P1 , P2

: direction of descent defined by equations (5a) and (5b)

r

: dimensionless coordinate

R1 (x1,x2,θ) : dimensionless reaction rate R2 (x1,x2,θ) : dimensionless reaction rate T

: dimensional temperature

x1,x2

: estimated dimensionless conversions

x

: dimensionless coordinate

Y1,Y2

: measured dimensionless conversions

Y3

: measured dimensionless temperature 1

z

: dimensionless coordinate

Greek Letters β1,β2

: search step sizes

γ1,γ2

: conjugate coefficients

λ1,λ2,λ3 : Lagrange multipliers defined by equation (12) λ4,λ5,λ6 : Lagrange multipliers defined by equation (16) θ(r,z)

: estimated dimensionless temperature

∆x1 , ∆x 2 , ∆θ : sensitivity functions defined by equation (7) ∆x
, ∆x
, ∆θ
: sensitivity functions defined by equation (8) 1

2

η

: convergence criteria

σ

: standard deviation for measurement

ω

: random number

Superscript n

: iteration index

INTRODUCTION The design of tubular type fixed bed catalytic reactors is generally based on a one-dimensional model. In this model it is assumed that conversion and temperature gradients occur only in the axial direction, and that the only transport mechanism operating in this direction is the overall flow itself, i.e. only the plug flow type is considered. This one-dimensional model leads to average values for the temperatures and conversions and provides no information concerning excessive temperatures along the axis. However, in many cases the inclusion of radial temperature gradients is unavoidable since the simulated temperature and conversion profiles in fixed bed catalytic reactors are very sensitive to the radial heat transfer. For this need Froment [1] derived a two-dimensional model to calculate conversion and temperature distributions in packed beds based on the continuity and energy equations with terms accounting for the reaction. This leads to a system of nonlinear second-order partial differential equations. 2

However, two non-linear reaction rate functions appeared in this twodimensional model, and their values are difficult to determine since for different kinds of reactions, both reaction rates varied. If the above mentioned two non-linear reaction rate functions are not available, the direct fixed-bed reactor problems can not be solved. The direct fixed-bed reactor problem is utilized to calculate the distributions of conversion and temperature distributions at interior points of a region when the other computational conditions are all specified. In contrast, the inverse fixed-bed reactor problem examined here involves the determination of the unknown conversion and temperature-dependent reaction rates in a two-dimensional fixed-bed reactor based on the information of conversion and temperature measurements taken within the domain. For the traditional inverse heat transfer problems, the conjugate gradient equations (CGM) [2] has been shown to be a very powerful tool to compute the inverse problems and applied to many different research fields. For instance Orlande and Ozisik [3] applied CGM in a reaction-diffusion parabolic problem to determine the reaction function. Huang and Yan [4] applied the CGM to predict the temperaturedependent thermal conductivity and heat capacity simultaneously. Huang and Chen [5] used boundary element method together with CGM to estimate the boundary heat fluxes for an irregular domain. Huang and Chin [6] used CGM in a 2-D inverse problem to determine the thermal conductivity in a non-homogeneous medium. Huang [7] applied the CGM in a nonlinear inverse vibration problem to estimate the unknown external forces. However, the coupled inverse problem to calculate the unknown conversion and temperature-dependent functions is very limited in the literature. Huang et al. [8] utilized the CGM in a non-linear inverse problem to estimate simultaneously the heat and mass production rates for a chemically reacting fluid. The purpose of this work is to utilize the CGM in the present non-linear inverse fixed-bed reactor problem and to estimate the conversion and temperature-dependent reaction rates.

3

The numerical experiments for this study with two different set of reaction rates and measurement errors will be illustrated in this study to show the validity of using CGM in the present two-dimensional inverse fixed-bed reactor problem.

I. DIRECT PROBLEM The following physical problem [1] is utilized to derive the expressions for use to estimate simultaneously the reaction rates for a fixed-bed reactor. The case considered here is of a rather complex nature, namely:

(1a) This reaction model represents the gas phase air oxidation of o-xylene into phthalic anhydride on V2O5 catalysts. A represents o-xylene, B phthalic anhydride, and C the final oxidation products to CO, CO2 and H2O. Therefore, this reaction model can be re-written as:

(1b) This model was occurring on a vanadium catalyst and the resultant steady-state material and energy balance partial differential equations can be obtained as followings [1]:  ∂ 2 x (r, z) 1 ∂x (r, z)  ∂x1 (r, z) 1 1 = A1  +  + B1R1 (x1, x 2 , θ) , 2 ∂z r ∂r   ∂r

in 0 < r <1, 0 < z <1  ∂ 2 x (r, z) 1 ∂x (r, z)  ∂x 2 (r, z) 2 2 = A1  +  + B1R 2 (x1, x 2 , θ) , 2 ∂z r ∂r   ∂r

(2a)

in 0 < r < 1, 0
(2b) 4

 ∂ 2θ(r, z) 1 ∂θ(r, z)  ∂θ(r, z) = A2  +  + B2R1 (x1, x 2 , θ) + B3R 2 (x1, x 2 , θ) , 2 ∂z r ∂r   ∂r in 0 < r <1, 0 < z <1

(2c)

∂x1 (r, z) ∂x 2 (r, z) ∂θ(r, z) = 0 ; = 0 ; = 0, at r = 0 ∂r ∂r ∂r ∂x1 (r, z) ∂x 2 (r, z) ∂θ(r, z) = 0 ; = 0 ; = Bi ( θ − θ w ) , at r = 1 ∂r ∂r ∂r x1 (r, z) = 0 ; x 2 (r, z) = 0 ; θ(r, z) = θ0 ,

at z = 0

(2d) (2e) (2f)

Where r and z represent the dimensionless radial and axial coordinates, respectively. x1 and x2 indicate the fractional conversion to B and C, respectively, while θ is the dimensionless temperature. R1 and R2 are two unknown reaction rate functions and need be determined. When computing conversions x1(r,z), x2(r,z) and temperature θ(r,z), i.e. giving R1(x1,x2,θ) and R2(x1,x2,θ) to calculate x1(r,z), x2(r,z) and θ(r,z), the direct problem (2) is nonlinear since the rates of reaction are function of temperature and conversions, therefore an iterative technique is needed to solve the direct problem with a suitable finite difference method. At this stage, R1(x1,x2,θ) and R2(x1,x2,θ) can not be replaced by R1(r,z) and R2(r,z), since x1(r,z), x2(r,z) and θ(r,z) are unknown before the direct problem calculations. However, when x1(r,z), x2 (r,z) and θ(r,z) are calculated under given R1(x1,x2,θ) and R2(x1,x2,θ) and some specified boundary conditions; the values of R1(x1,x2,θ) and R2(x1,x2,θ) at any position, (r,z), should be known and fixed because conversions x1(r,z), x2(r,z) and temperature θ(r,z) are known and fixed at any point (r,z). This implies that R1(x1,x2 ,θ) ≡ R1(r,z) and R2(x1,x2,θ) ≡ R2(r,z) can be used. For the inverse calculations considered here, the measured conversions and temperature are assumed to be known either from numerical simulations or from experimental measurements. Once x1(r,z), x2(r,z) and θ(r,z) are obtained, there exists some unknown but fixed reaction rate functions R1 and R2 that their values R1(r,z) and R2(r,z), at any specific position, (r,z), must satisfy equation (2) to match with the existing conversion and temperature distributions, x1(r,z), x2(r,z) and θ(r,z). Therefore 5

R1(x1,x2,θ) ≡ R1(r,z) and R2(x1,x2,θ) ≡ R2(r,z) can be used in the inverse calculations even when the problem is nonlinear. Under this consideration the original coupled direct problems can be de-coupled and become three independent problems. The implicit finite difference method with iteration can be used to solve the present direct problem.

II. INVERSE PROBLEM

For the present inverse problem, the reaction rates R1(x1,x2,θ) and R2(x1,x2,θ) are regarded as being unknown, but everything else in equation (2) is known. Besides, the measured conversion and temperature distributions are assumed available. The measured conversions at position (r,z) can be denoted as Y1(r,z) and Y2(r,z) and measured temperature is given as Y3(r,z). The present inverse problem can be stated as follows: by utilizing the above mentioned measured conversions and temperature data Y1(r,z), Y2(r,z) and Y3(r,z), estimate the unknown reaction rates R1(x1,x2,θ) and R2(x1,x2,θ) over the specified domain. Once the reaction rates R1(r,z) and R2(r,z) are determined, they can be transformed back to R1(x1,x2,θ) and R2(x1,x2,θ). The following functional J is defined and is to be minimized in the present inverse problem, zf

1

2 2 2 ∫  (x1 − Y1) + (x 2 − Y2 ) + (θ − Y3 )  rdrdz z=0 r=0

J  R1 ( r, z ) , R 2 ( r, z )  = ∫

(3)

Here x1(r,z) and x2(r,z) are the estimated (or computed) conversions and θ(r,z) is the estimated (or computed) temperature. They can be computed from Equation (2) by using the estimated reaction rates, R1(r,z) and R2(r,z). The traditional measurement techniques cannot be utilized in this study since there are very huge numbers of measurements. Alternative non-invasive techniques need be considered to measure the conversions and temperature. Zhang and Ethier [9] used a specially designed cell and associated optical equipment to measure the local refractive index and light extension coefficient of a polymer mixture. This non-invasive technique can measure the concentration in polymer mixtures and can be utilized to 6

measure the conversions x1 and x2 in this work. Besides, the temperature θ can be measured by infrared camera.

III. CONJUGATE GRADIENT METHOD FOR MINIMIZATION

An iterative process is now applied to determine R1(r,z) and R2(r,z) by minimizing the above functional J  R1 ( r, z ) , R 2 ( r, z )  : R1n +1 (r, z) = R1n (r, z) − β1n P1n (r, z)

n=0,1,2,..

(4a)

R 2n +1 (r, z) = R 2 n (r, z) − β2 n P2 n (r, z)

n=0,1,2,..

(4b)

where β1n and β 2 n are the search step sizes in going from iteration n to iteration n+1, and P1n (r, z) and P2 n (r, z) are the directions of descent (i.e. search directions) given by P1n (r, z) = J1′ n (r, z) + γ1n P1n −1 (r, z)

(5a)

P2 n (r, z) = J′2 n (r, z) + γ 2 n P2 n −1 (r, z)

(5b)

which are the conjugation of the gradient directions J1′ n (r, z) and J′2 n (r, z) at iteration n and the directions of descent P1n −1 (r, z) and P2 n −1 (r, z) at iteration n-1. The conjugate coefficients can be determined from zf

1

n 2 ∫ [J1′ (r, z)] drdz

∫

γ1n = z = 0 r = 0 zf

1

∫

∫

0 with γ1 = 0

(6a)

0 with γ 2 = 0

(6b)

[J1′ n −1 (r, z)]2 drdz

z =0 r =0 zf 1

n 2 ∫ [J′2 (r, z)] drdz

∫

γ 2n = z = 0 r = 0 zf

1

n −1 2 ∫ ∫ [J′2 (r, z)] drdz z=0 r=0 The step sizes β1n and β 2 n and the gradients of the functional J1′ n (r, z) and

J′2n (r, z) need be calculated first before the iterative equation (4) can be executed. In

order to determine R1(r,z) and R2(r,z), two sensitivity problems and two adjoint problems are constructed as described below. IV. THE SENSITIVITY PROBLEMS AND SEARCH STEP SIZES

7

It is assumed that when R1(r,z) undergoes a variation ∆R1(r,z), x1(r,z), x2(r,z) and θ(r,z) are perturbed by ∆x1 ( r,z ) , ∆x 2 ( r,z ) and ∆ θ ( r,z ) , respectively. Replacing in the direct problem x1 by x1+ ∆x1 ( r,z ) , x2 by x2+ ∆x 2 ( r,z ) and θ by θ + ∆ θ ( r,z ) , subtracting the resulting expressions from the direct problem and neglecting the second-order terms, the sensitivity problems with variables ∆x1 ( r,z ) , ∆x 2 ( r,z ) and ∆ θ ( r,z ) can be obtained.

 ∂ 2 ∆x ( r,z ) 1 ∂∆x ( r,z )  ∂∆x1 ( r,z ) 1 1 =A1  +  + B1∆R1 (r,z) , 2 ∂z r ∂ r  ∂r 

in 0 < r < 1，0 < z < 1

(7a)

 ∂ 2 ∆x ( r,z ) 1 ∂∆x ( r,z )  ∂∆x 2 ( r,z ) 2 2 = A1  +  , in 0 < r < 1，0 < z < 1 2 ∂z r ∂r   ∂r

(7b)

 ∂ 2 ∆θ ( r,z ) 1 ∂∆ θ ( r,z )  ∂ θ ( r,z ) = A2  +  + B2∆R1 (r,z) , ∂z r ∂r  ∂r 2 

in 0 < r < 1，0 < z < 1

(7c)

∂∆x1 ( r,z ) ∂∆x 2 ( r,z ) ∂∆θ ( r,z ) = = =0, ∂r ∂r ∂r ∂∆x1 ( r,z ) ∂∆x 2 ( r,z ) ∂∆θ ( r,z ) = = 0 and = − Bi∆θ ( r,z ) , ∂r ∂r ∂r ∆x1 ( r,z ) = ∆x 2 ( r,z ) = ∆θ ( r,z ) = 0 ,

at

r = 0 (7d)

at

r = 1 (7e)

at

z =0

(7f)

Similarly, when R2(r,z) undergoes a variation ∆R2(r,z), x1(r,z), x2(r,z) and θ(r,z) are perturbed by ∆x
( r,z ) , ∆x
( r,z ) and ∆θ
( r,z ) , respectively. Finally the following 1

2

problem for the sensitivity functions ∆x
1 ( r,z ) , ∆x
2 ( r,z ) and ∆θ
( r,z ) can be obtained.  ∂ 2∆x
( r,z ) 1 ∂∆x
( r,z )  ∂∆x
1 ( r,z ) 1 1  , in 0 < r < 1，0 < z < 1 =A1  + 2 ∂z r ∂r   ∂r  ∂ 2∆x
( r,z ) 1 ∂∆x
( r,z )  ∂∆x
2 ( r,z ) 2 2 = A1  +  + B1∆R 2 (r,z) , 2 ∂z r ∂ r ∂r  

8

(8a)

in 0 < r < 1，0 < z < 1

(8b)

 ∂ 2∆θ
( r,z ) 1 ∂∆θ
( r,z )  ∂θ
( r,z ) = A2  +  + B3∆R 2 (r,z) , ∂z r ∂r   ∂r 2 

in 0 < r < 1，0 < z < 1

(8c)

∂∆x
1 ( r,z ) ∂∆x
2 ( r,z ) ∂∆θ
( r,z ) = = =0, ∂r ∂r ∂r ∂∆x
1 ( r,z ) ∂∆x
2 ( r,z ) ∂∆θ
( r,z ) = = 0 and = − Bi∆θ
( r,z ) , ∂r ∂r ∂r ∆x
( r,z ) = ∆x
( r,z ) = ∆θ
( r,z ) = 0 , 1

2

at

r=0

(8d)

at

r = 1 (8e)

at

z =0

(8f)

The above sensitivity problems can be solved by implicit finite difference method. The functional J  R1n +1 ( r, z ) , R n2 +1 ( r, z )  for iteration n+1 is obtained by   rewriting equation (3) as J  R1n +1 ( r, z ) , R 2n +1 ( r, z )  =   zf

∫

1

∫

{  x1 ( r, z; R1n − β1n P1n , R 2n − β2n P2n ) − Y1 ( r, z )

2

z =0 r =0

(

(9)

2

) 2 + θ ( r, z; R1n − β1n P1n , R n2 − βn2 P2n ) − Y3 ( r, z )  }rdrdz   +  x 2 r,z;R1n − β1n P1n , R 2n − β2n P2n − Y2 ( r, z )   

where R1n +1(r, z) and R 2n +1(r, z) were replaced by equations (4a) and (4b). If the estimated conversions and temperatures are linearized by a Taylor’s expansion, equation (9) takes the form: J  R1n +1(r, z), R n2 +1 (r, z)  =   zf

∫

1

∫

{  x1(r,z;R1n , R 2n ) − β1n ∆x1(P1n ) − β2n ∆x
1(P2n ) − Y1(r, z) 

2

z =0 r =0

(10)

2 +  x 2 (r,z;R1n , R 2n ) − β1n ∆x 2 (P1n ) − β2n ∆x
2 (P2n ) − Y2 (r, z)   

+ θ (r,z;R1n , R 2n ) − β1n ∆θ (P1n ) − β2n ∆θ
(P2n ) − Y3 (r, z)   

2

}rdrdz

where x1(r,z;R1n , R 2n ) , x 2 (r,z;R1n , R 2n ) and θ(r,z;R1n , R 2n ) are the solutions of the

direct problem by using the estimates for R1(r,z) and R2(r,z). 9

Equation (10) is differentiated with respect to search step sizes β1n and β 2 n , alternately, and then made equal to zero to obtain two equations. By solving these two equations for two unknowns β1n and β 2 n , finally the following expressions for search step sizes can be obtained as: (C C − C3C5 ) β1n = 2 4 (C1C2 − C3C3 ) (C C − C3C4 ) β2 n = 1 5 (C1C2 − C3C3 )

(11a) (11b)

Where zf

C1 = ∫

1

2 2 2 ∫ (∆x1 + ∆x 2 + ∆θ )rdrdz

(11c)

z=0 r=0 zf 1

C2 = ∫

2 2 2 ∫ (∆x
1 + ∆x
2 + ∆θ
)rdrdz

C3 = ∫

∫

C4 = ∫

∫  ∆x1(x1 − Y1) + ∆x 2 (x 2 − Y2 ) + ∆θ(θ − Y3 )  rdrdz

(11f)

C5 = ∫

∫  ∆x
1(x1 − Y1 ) + ∆x
2 (x 2 − Y2 ) + ∆θ
(θ − Y3 )  rdrdz

(11g)

z =0 r =0 zf 1

z = 0 r =0 zf 1

(11d)

( ∆x1∆x
1 + ∆x 2∆x
2 + ∆θ∆θ
)rdrdz

(11e)

z =0 r =0 zf 1

z =0 r =0

The sensitivity problem considered here is to calculate the sensitivity functions using ∆R1(r,z) and ∆R2(r,z), therefore, the implicit finite difference method can be used to solve the present sensitivity problems without iteration since ∆R1(r,z) and ∆R2(r,z) are not function of ∆x ( r,z ) , ∆x ( r,z ) , ∆ θ ( r,z ) , ∆x
( r,z ) , ∆x
( r,z ) and ∆θ
( r,z ) . 1

2

1

2

V. ADJOINT PROBLEMS AND GRADIENT EQUATIONS

To obtain the adjoint problem, equations (2a), (2b) and (2c) are multiplied by the Lagrange multipliers (or adjoint functions) λ1(r,z), λ2(r,z) and λ3(r,z), respectively, the resultant expressions are sum up together to become one term and then integrated over the correspondent space domains. The result is then added to the right hand side of equation (3) to yield a new expression for the functional. Firstly, the variation ∆J1 is obtained by perturbing R1 by R1+∆R1, x1 by x1+ ∆x1 ( r,z ) , x2 by x2+ ∆x 2 ( r,z ) and θ by θ + ∆ θ ( r,z ) in the new expression for the functional, subtracting the resulting expression from the original equation and 10

neglecting the second-order terms. Finally the boundary conditions of the sensitivity problem are utilized to perform the integration by parts. The vanishing of the integrands leads to the following adjoint problem for the determination of λ1(r,z), λ2(r,z) and λ3(r,z):  ∂ 2 λ ( r,z ) 1 ∂λ ( r,z )  ∂λ1 ( r,z ) 1 1 + A1  +  +2(x1 − Y1) = 0 , in 0
∂λ1 ( r,z ) ∂r ∂λ1 ( r,z ) ∂r

= =

∂λ 2 ( r,z )

=

∂r ∂λ 2 ( r,z )

∂λ3 ( r,z ) ∂r

= 0 and

∂r

=0, ∂λ3 ( r,z ) ∂r

=Bi λ3 ( r,z ) ,

λ1 ( r,zf ) = λ 2 ( r,z f ) = λ3 ( r,zf ) = 0 ,

at

r=0

at

r=1

at

z =zf

(12d) (12e) (12f)

The adjoint problem differs from the standard initial value problems in that the final condition at z = z f is specified instead of the traditional initial condition. However, this problem can be transformed back into an initial value problem by using the transformation variables z∗ = zf - z. The techniques of implicit finite differences method can then be used to solve the above adjoint problem. The only term which is not equal to zero in the Lagrange functional variation is zf

∆J1[R1(r,z)] = ∫

1

∫ [ λ1(r, z)B1 + λ3 (r, z)B2 ]∆R1 (r, z)rdrdz

(13)

z =0 r = 0

From definition [2], the residual functional variation in the Hilbert space L2 can be presented as zf

∆J1[R1(r, z)] = ∫

1

' ∫ J1 (r, z)∆R1 (r, z)rdrdz

(14)

z=0 r=0

A comparison of equations (14) and (15) leads to the gradient of functional J1 ' :

11

J1'  R1 ( r, z )  = λ1 (r, z)B1 + λ3 (r, z)B2

(15)

Similarly, equations (2a), (2b) and (2c) are multiplied by the Lagrange multipliers (or adjoint functions) λ4(r,z), λ5(r,z) and λ6(r,z), respectively. When R2(r,z) undergoes a variation ∆R2(r,z), x1(r,z), x2(r,z) and θ(r,z) are perturbed by ∆x
1 ( r,z ) , ∆x
( r,z ) and ∆θ
( r,z ) , respectively. Finally the adjoint problem for determining the 2

adjoint functions λ4(r,z), λ5(r,z) and λ6(r,z) can be obtained as  ∂ 2 λ ( r,z ) 1 ∂λ ( r,z )  ∂λ 4 ( r,z ) 4 4 + A1  +  +2(x1 − Y1 ) = 0 , in 0
∂λ 4 ( r,z ) ∂r ∂λ 4 ( r,z ) ∂r

= =

∂λ5 ( r,z ) ∂r ∂λ5 ( r,z ) ∂r

=

∂λ6 ( r,z ) ∂r

= 0 and

=0, ∂λ 6 ( r,z ) ∂r

=Bi λ 6 ( r,z ) ,

λ 4 ( r,zf ) = λ5 ( r,z f ) = λ6 ( r,z f ) = 0 ,

at

r=0

(16d)

at

r=1

(16e)

at

z =zf

(16f)

and the gradient of functional J 2 ' can be obtained as

J 2'  R 2 ( r, z )  = λ5 (r, z)B1 + λ6 (r, z)B2

(17)

VI. STOPPING CRITERION

When the problem contains no measurement errors, the following check condition can be used to stop the iteration:

(18)

J[R1 (r, z), R 2 (r, z)] < η

12

where η is the small numbers. However, the measured conversions and temperature data may contain measurement errors. Therefore, it is not expected that the functional equation (3) will be equal to zero at the last iteration step. Following the rigorous justification by Alifanov [2], the following discrepancy principle was chosen as the stopping criterion, i.e. it is assumed that the residuals for conversions and temperature can all be equal to σ  x1 ( r, z ) − Y1 ( r, z )  =  x 2 ( r, z ) − Y2 ( r, z )  = θ ( r, z ) − Y3 ( r, z )  ≈ σ (19)

here σ is the standard deviations of the measurements, which are assumed to be constant. By substituting equation (19) into equation (3), the following expressions are obtained for η : η = 3σ 2 z f

(20)

The stopping criterion is now given by equation (18) with η determined from equation (20).

VII. COMPUTATIONAL PROCEDURE

The computational procedure for the solution of the present inverse problem can be summarized as follows: Suppose

R1n (r, z) and R 2n (r, z) are available at iteration n.

Step 1. Solve the direct problem given by equation (2) for x1(r,z), x2(r,z) and θ(r,z). Step 2. Examine the stopping criterion. Continue if not satisfied. Step 3. Solve the adjoint problems given by equations (12) and (16) for λ1(r,z), λ2(r,z), λ3(r,z), λ4(r,z), λ5(r,z) and λ6(r,z). Step 4. Compute the gradient of the functional J1 '[R1 (r, z)] and J 2 '[R 2 (r, z)] from equations (16) and (18), respectively.

Step 5. Compute the conjugate coefficients γ1n and γ 2 n and the direction of descent P1n (r, z) and P2n (r, z) from equations (6) and (5), respectively. 13

Step 6. Set ∆R1(r,z) = P1n (r,z) and ∆R 2 (r,z) = P2n (r,z) , and solve the sensitivity problems given by equations (7) and (8) for ∆x1 ( r,z ) , ∆x 2 ( r,z ) , ∆ θ ( r,z ) , ∆x
1 ( r,z ) , ∆x
2 ( r,z ) and ∆θ
( r,z ) .

Step 7. Compute the search step sizes β1n and β 2 n from equations (11). Step 8. Compute the new estimation for R1n+1(r, z) and R 2n+1(r, z) from equations (4a) and (4b) and return to step 1.

VIII. RESULTS AND DISCUSSIONS

The objective of the present work is to show the validity of the CGM to estimate the unknown conversion and temperature-dependent reaction rates, R1(x1,x2,θ) and R2(x1,x2,θ) simultaneously, for a fixed-bed reactor undergoes a reaction regarding the gas phase air oxidation of o-xylene into phthalic anhydride on V2O5 catalysts, with no prior information on the functional form of the unknown reaction rates. In order to illustrate the accuracy of the CGM to determine R1(x1,x2,θ) and R2(x1,x2,θ) simultaneously with inverse analysis, based on the measured conversion and temperature distributions, the following two specific examples were considered for the estimation of the reaction rates R1 (x1,x2,θ) and R2(x1,x2,θ). In all the test cases considered here, the initial guesses of R1(r,z) and R2(r,z) are both chosen as R10(r,z) = R20(r,z) = 0.0. In order to compare the results for situations involving random measurement errors, a normally distributed uncorrelated error with zero mean and constant standard deviation will be assumed. The simulated inexact measurement data Y1, Y2 and Y3 can be expressed as Y1 = Y1,dir + ω σ

(21a)

Y2 = Y2,dir + ω σ

(21b)

Y3 = Y3,dir + ω σ

(21c)

where Y1,dir, Y2,dir and Y3,dir are the solutions of the Equation (2) with exact R1(x1,x2,θ) and R2(x1,x2,θ); σ is the standard deviation of the measurements; and ω is a

14

random variable that generated by subroutine DRNNOR of the IMSL [10] and is between -2.576 to 2.576 for a 99% confidence bound. Two numerical experiments in determining R1(x1,x2,θ) and R2(x1 ,x2,θ) by the inverse analysis are now presented below. (A). Numerical test case 1

The parameters for the direct problem are given as follows: A1 = 5.76 ; A2 = 10.94 ; Bi = 2.5 ; B1 = 5.106 ; B2 = 3.144 ; B3 = 11.16 ; θw = 1.0; θ0 = 1.0; zf = 1.3 Besides, the space increments used in numerical calculations for r and z directions are taken as ∆r = 0.01 (i.e. 100 grid points in r direction) and ∆z = 0.01 for zf = 1.0 (i.e. i.e. 100 grid points in z direction), respectively. Therefore totally of 13231 unknown discreted values for both R1(r,z) and R2(r,z) are to be estimated in this study. The exact conversion and temperature-dependent reaction rates, R1(x1,x2,θ) and R2(x1,x2,θ), for a fixed-bed reactor problem are assumed as R1 = k1 (1 − x1 − x 2 ) − k 2 x1

(22a)

R 2 = k 2 x1 + k 3 (1 − x1 − x 2 )

(22b)

where   1  ki = exp αi + γi 1 −   ,  θ  

i = 1, 2, 3

(23)

here α1 = −1.74; α2 = −4.24; α3 = −3.89; γ1 = 2.16; γ2 = 2.51; γ3 = 2.29 It should be noted that when generating simulated measurements for conversions Y1(r,z), Y2(r,z) and temperature Y3(r,z), exact R1(x1,x2,θ) and R2(x1,x2,θ) are used in Equation (2) and thus the problem is coupled and nonlinear, therefore iterative technique is needed for its solutions. However, in the inverse problem, the reaction rates exist in the forms of R1(r,z) and R2(r,z), therefore the problem becomes decoupled and linear. For this reason the estimated conversions and temperature can be calculated directly. In all the test cases considered here, the initial guess R10(r,z) = R20(r,z) = 0.0 are used. However, the exact values of R1(r,zf) and R2(r,zf) are not equal to zero in the 15

present study, it is therefore the singularity at final position z = zf will occur in the present test case and the estimated values for R1(r,zf) and R2(r,zf) must be the same as the initial guess values, i.e. R1(r,zf) = R2(r,zf) = 0.0. However, if the inverse calculations are extended to z = 1.3 and reserve the inverse solutions only to z = 1.0, the estimations for R1(r,z) and R2(r,z) become reliable under this strategy. The inverse calculation is first considered by assuming no measurement errors, i.e. σ = 0.0. By setting η = 10-7, after 65 iterations, the inverse estimations can be obtained. The exact conversions Y1 and Y2 and temperature,Y3, are shown in Figure 1, the exact reaction rates, R1(r,z) and R2(r,z), are shown in Figure 2, while the estimated x1, x2 and θ are illustrated in Figure 3 and the estimated R1(r,z) and R2(r,z) are shown in Figure 4, respectively. Due to the singularity near zf = 1.3, the estimated values, for just up to z = 1.0, have been illustrated in Figures 3 and 4. The average errors for the estimated R1(r,z) and R2(r,z) are calculated as ERR1 = 0.37 % and ERR2 = 0.69 %, respectively, where the average errors for the estimated R1(r,z) and R2(r,z) can be obtained by  101 101 R (m, n) − Rˆ 1(m, n)  ERR1 %=  ∑ ∑ 1  ÷ 10201×100% R1(m, n)   m =1 n =1  101 101 R (m, n) − Rˆ 2 (m, n)  ERR2 %=  ∑ ∑ 2  ÷ 10201×100% R 2 (m, n)  m =1 n =1 

(24a)

(24b)

here m and n represent the index of discreted grid in r and z directions, respectively, while R1(m,n) and R2(m,n) denote the exact values while Rˆ1 ( m, n ) and Rˆ 2 ( m, n ) represent the estimated values of reaction rates. It is clear from Figures 1 and 3 that there is an excellent agreement between the exact and estimated conversions and temperatures. From the above calculated average errors as well as Figures 2 and 4 it can be learned that the estimated reaction rates match very well to the exact values. It is also interesting if the issue regarding the influence of the measurement errors on the inverse solutions is examined. The measurement error for the conversions and temperature is taken as σ = 0.0003. After 55 iterations, the estimated reaction 16

rates R1(r,z) and R2 (r,z) can be obtained and are shown in Figure 7. Here the stopping criterion can be obtained by discrepancy principle which is given in Equation (21). The average errors for R1(r,z) and R2(r,z) are calculated as ERR1 = 0.44 % and ERR2 = 1.07 %, respectively. Next, the measurement error is further increased to σ = 0.001, after 11 iterations, the estimated reaction rates R1(r,z) and R2 (r,z) can be obtained and are shown in Figure 8. The average errors for R1(r,z) and R2(r,z) are calculated as ERR1 = 0.85 % and ERR2 = 1.4 %, respectively. This implies that reliable inverse solutions can still be obtained when measurement errors are considered. (B). Numerical test case 2

In the second test case, the parameters for the direct problem are chosen the same as those used in test case 1, but the exact reaction rates, R1(x1,x2,θ) and R2(x1,x2,θ), are now assumed as R 1 = k1 (1 − x1 +x 2 ) + exp(k 2 x1 )

(25a)

R 2 = exp(k 2 x1 ) + 0.3  k 3 ( − x1 − x 2 ) 

(25b)

where   1  k1 = cos α1 + γ1  1 −    θ     1  k 2 = sin  α 2 + γ 2 1 −    θ  

(26)

  1  k 3 = cos  α3 + γ 3  1 −    θ   where α1 = −1.74; α2 = −4.24; α3 = −3.89; γ1 = 2.16; γ2 = 2.51; γ3 = 2.29 First, assuming σ = 0.0 and setting η = 3 × 10-3, after 247 iterations, the inverse estimations can be obtained. The exact reaction rates, R1(r,z) and R2(r,z), are shown in Figure 7, while the estimated R1(r,z) and R2(r,z) are shown in Figure 8, respectively. The average errors for the estimated R1(r,z) and R2(r,z) are calculated as ERR1 = 1.98 % and ERR2 = 1.62 %, respectively. It is clear from Figures 7 and 8 that good match between the exact and estimated reaction rates can be seen.

17

It is then continue to discuss the influence of the measurement errors on the inverse solutions. The measurement error for the conversions and temperature is taken as σ = 0.06 and then increased to σ = 0.1. After 138 and 38 iterations, respectively, the estimated reaction rates R1(r,z) and R2 (r,z) can be obtained and are shown in Figures 9 and 10, respectively. The average errors for R1(r,z) and R2(r,z) are calculated as ERR1 = 2.52 % and ERR2 = 2.59 % for σ = 0.06 and ERR1 = 2.58 % and ERR2 = 3.98 % for σ = 0.1, respectively. From these two test cases it can be concluded that an inverse fixed-bed reactor problem to determine the conversion and temperature-dependent reaction rates, R1(x1,x2,θ) and R2 (x1,x2,θ) simultaneously, is now completed. Reliable estimations for R1(x1,x2,θ) and R2(x1,x2,θ) can be obtained when considering both exact and error measurements. IX. CONCLUSIONS

The Conjugate Gradient Method (CGM) was applied successfully for the solution of the inverse fixed-bed reactor problem to estimate the unknown conversion and temperature-dependent reaction rate functions, R1(x1,x2,θ) and R2(x1,x2,θ), by utilizing simulated conversion and temperature readings obtained from different measured positions. Two test cases involving different forms of reaction rates and measurement errors were considered. The results show that in test case 1 the maximum average errors for the estimated R1(r,z) and R2(r,z) are obtained as 0.85 %, and 1.4 %, respectively, in test case 2 they are 2.58 % and 3.98 %, respectively. This implies that the inverse solutions obtained by CGM are still reliable even when the measurement errors are considered. ACKNOWLEDGMENT

This work was supported in part through the National Science Council, R. O. C., Grant number, NSC-100-2221-E-006-011-MY3. REFERENCES

1. G. F. Froment, Fixed Bed Catalytic Reactors, Industrial and Engineering Chemistry, 59 (1976) 18-27. 18

2. O. M. Alifanov, Inverse Heat Transfer Problem, Springer-Verlag, Berlin, 1994. 3. H. R. B. Orlande and M. N. Ozisik, Determination of the reaction function in a reaction-diffusion parabolic problem, J. Heat Transfer, 116 (1994) 1041-1044. 4. C. H. Huang and J. Y. Yan, An Inverse Problem in Simultaneously Measuring Temperature Dependent Thermal Conductivity and Heat Capacity, Int. J. Heat and Mass Transfer, 38 (1995) 3433-3441. 5. C. H. Huang and C. W. Chen, A Steady Inverse Problem of Estimating Boundary Condition in An Irregular Domain with Statistical Analysis, Numerical Heat Transfer. 33 (1998) 251-268. 6. C. H. Huang and S. C. Chin, A Two-Dimensional Inverse Problem in Imaging the Thermal Conductivity of a Non-homogeneous Medium, Int. J. Heat and Mass Transfer. 43 (2000) 4061-4071. 7. C. H. Huang, A Nonlinear Inverse Vibration Problem of Estimating the External Forces for A System with Displacement-Dependent Parameters, J. Sound and Vibration, 248 (2001) 789-807. 8. C. H. Huang, C. Y. Yeh and Helcio R. B. Orlande, A Non-Linear Inverse Problem in Simultaneously Estimating the Heat and Mass Production Rates for A Chemically Reacting Fluid, Chemical Engineering Science, 58 (2003) 3741-3752. 9. W. Zhang, C. R. Ethier, Method for non-invasive concentration measurement in polymer mixtures: Tests with dextran and cytochrome-c, Chemical Engineering Communications, 189 (2002) 101-113. 10. IMSL Library Edition 10.0. User's Manual: Math Library Version 1.0, IMSL, Houston, TX, 1987.

19

Figure 1a. The exact conversion distribution, x1(r,z), in case 1.

Figure 1b. The exact conversion distribution, x2(r,z), in case 1.

Figure 1c. The exact temperature distribution, θ(r,z), in case 1. 20

Figure 2a. The exact reaction rate, R1(r,z), in case 1.

Figure 2b. The exact reaction rate, R2(r,z), in case 1.

21

Figure 3a. The estimated conversion distribution, x1(r,z), in case 1 with σ = 0.0.

Figure 3b. The estimated conversion distribution, x2(r,z), in case 1 with σ = 0.0.

Figure 3c. The estimated temperature distribution, θ(r,z), in case 1 with σ = 0.0. 22

Figure 4a. The estimated reaction rate, R1(r,z), in case 1 with σ = 0.0.

Figure 4b. The estimated reaction rate, R2(r,z), in case 1 with σ = 0.0.

23

Figure 5a. The estimated reaction rate, R1(r,z), in case 1 with σ = 0.0003.

Figure 5b. The estimated reaction rate, R2(r,z), in case 1 with σ = 0.0003.

24

Figure 6a. The estimated reaction rate, R1(r,z) , in case 1 with σ = 0.001.

Figure 6b. The estimated reaction rate, R2(r,z) , in case 1 with σ = 0.001.

25

Figure 7a. The exact reaction rate, R1(r,z), in case 2.

Figure 7b. The exact reaction rate, R2(r,z), in case 2.

26

Figure 8a. The estimated reaction rate, R1(r,z), in case 2 with σ = 0.0.

Figure 8b. The estimated reaction rate, R2(r,z), in case 2 with σ = 0.0.

27

Figure 9a. The estimated reaction rate, R1(r,z), in case 2 with σ = 0.06.

Figure 9b. The estimated reaction rate, R2(r,z), in case 2 with σ = 0.06.

28

Figure 10a. The estimated reaction rate, R1(r,z), in case 2 with σ = 0.1.

Figure 10b. The estimated reaction rate, R2(r,z), in case 2 with σ = 0.1.

29