An Efficient Discretization Approach for Partial Differential Equations describing Chemical Reaction Systems

19th European Symposium on Computer Aided Process Engineering – ESCAPE19 J. JeĪowski and J. Thullie (Editors) © 2009 Elsevier B.V. All rights reserved.

901

An Efficient Discretization Approach for Partial Differential Equations describing Chemical Reaction Systems Jan C. Schöneberger, Harvey Arellano-Garcia, Günter Wozny Chair of Process Dynamics and Operation, Berlin Institute of Technology, Sekr. KWT9, Str. des 17. Juni, D-10623 Berlin, Germany, E-mail: [email protected]

Abstract In this work, an efficient solution approach based on the orthogonal collocation on finite elements is proposed for PDE-Systems describing chemical reaction systems. For this purpose, three discretization methods are compared to each other in terms of convergence, accuracy, and computation time. First, the simplest and most applied discretization method is the method of lines (ML). Here, the differential terms in the equation system are replaced by linear difference approximations. The second one is the multiple shooting (MS), which makes use of the boundary conditions, and thus, the PDE is reformulated to a root finding problem on ODE. Finally, the orthogonal collocation (OC) is a Runge-Kutta type discretization method that leads to the highest possible error order by collocating the discrete points at their optimal positions. These three approaches are tested for different numbers of discrete points on a PDE of the type which results from modelling fixed bed reactors (FBR). A test system with an existing analytical solution was chosen for an impartial comparison. In addition, two further case studies are considered to point out the advantages of the proposed discretization approach.

Keywords: Partial Differential Equations, Discretization, Collocation, Multiple Shooting, Fixed Bed Reactor

Orthogonal

1. Introduction Catalytic fixed bed reactors have a wide range of industrial applications. The needs for the different heterogeneous catalyzed gas phase reactions led to different geometrical designs for these reactors e.g. tubular reactor, radial reactor, mono line reactor, etc. Due to the unknown precise geometry of the catalyst filling and the complex nature of the chemical reactions, it is generally not possible to apply CFD tools. However, the effects inside these reactors can be modeled very detailed by the consequent application of conservation equations, i.e. mass, energy and momentum balances, combined with transport equations. The resulting equation system contains partial differential equations (PDE), ordinary differential equations (ODE) and algebraic equations (AE) which are coupled in a complex manner. Even with today’s high CPU power it is still a challenging task to solve these kinds of problems. Moreover, the chosen discretization approach is decisive for success or failure. It affects considerably the robustness and the convergence velocity of the solving algorithm. The choice of the discretization method also affects the possibilities of physical modelling. The orthogonal collocation (OC) e.g. allows but does not require the application of a boundary condition for second derivatives. This second boundary condition often can not be declared physically, and in the case of ML, it additionally lowers the numerical accuracy.

J.C. Schöneberger et al.

902

2. Problem Statement Reaction systems can be modeled using the conservation equations such as component balances (1), energy balance (2), and momentum balance [1]. The momentum balance is not considered here, because it is impossible to solve it without the knowledge of the exact geometry. Thus, based on the assumptions taken, the Ergun equation is used here for the pressure drop [2]. NR dci & = div (Di ∇ci ) − ∇(wci ) + ¦ Ȟi,j rj − n PI a PI dt j

NC

¦cc i

i

P ,i

dT = div (λ ∇ T ) − dt

NC

¦cc i

P ,i

(1)

& ∇ (w T ) +

i

NR

¦

− Δ h R , j r j − q PI a PI

(2)

j

From this type of PDE the generically test equation (3) can be derived for e.g. the Temperature T and one spatial dimension x.

dT ∂ 2T ∂T + k3 2 = k 0 + k1T + k 2 dt ∂x ∂x

(3)

Typical boundary conditions for these type of equations are a Dirichlet bound at the inlet (x = 0) and a Neumann bound at the outlet (x = L), given in Eq. (4).

T ( x = 0 ,t ) = T0

§ ∂T · =0 ¸ ¨ © ∂x ¹ x = L ,t

( Dirchlet )

(4)

( Neumann )

The parameters for the test equation are given in Tab. 1. The advantage of this equation is that it can be solved analytically for the steady state case. This makes comparable the numerical solutions with the different discretization approaches. Due to the large terms of the analytical solution it is not posed here. Table 1: Parameter Set for the test PDE

k0

k1

k2

k3

T0

L

1 J/s

0.001 J/(K s)

0.01 Jm/(K s)

0.01 Jm²/(K s)

1000 K

1m

3. Discretization Approaches 3.1. Method of Lines Due to its simplicity, the Method of Lines is the most common discretization approach [3]. The spatial differential terms are replaced by difference terms as shown in Eq. (5). This transforms the PDE into an ODE. But, the inclusion of the boundary conditions leads to additional AE, thus, resulting in a DAE system. This affects the numerical behavior and the computational effort as well. The error order of the ML is 2.

§ ∂ 2T · Tk +1 − Tk −1 ML ¨¨ 2 ¸¸ ⎯⎯→ x Δx 2 ∂ © ¹ x = xk

T −T § ∂T · ML ¨ ¸ ⎯⎯→ k +1 k −1 2Δx © ∂x ¹ x = xk

(5)

Moreover, the boundaries lead to an additional problem, namely, for the inclusion of the Neumann boundary, the central difference has to be replaced by the backward

An Efficient Discretization Approach for Partial Differential Equations

903

difference in the last element reducing the error order to 1. This effect can be seen in the resulting numerical solutions. 3.2. Multiple Shooting The approach of the multiple shooting is to apply common ODE solvers for the solution of the PDE. For this purpose the boundary conditions have to be transformed to initial conditions. In the single shooting method, it is done by introducing a new parameter s. This is determined by the solution of the root finding problem given in Eq. (6). For a set of equations this method is instable, but it can be stabilized by introducing intervals (i.e. discrete points) in which the single shooting is applied fulfilling new boundary conditions that guarantee a continuous solution. This approach is called multiple shooting [4]. The error order of the MS depends on the error order of the applied ODE solver.

§ ∂T · ¨ ¸ =s © ∂x ¹ x =0

· § § ∂T ( s ) · ¸ − 0 ¸¸ = 0 ¹ © © ∂x ¹ x= L

Φ ( s ) = ¨¨ ¨

(6)

3.3. Orthogonal Collocation The main idea of the proposed discretization approach based on the OC is that the numerical solution of the PDE is given as a polynomial, Eq. (7), [5]. The degree of the polynomial is determined by the number of discrete points NK, and dn represents the polynomial coefficients. NK

T ( t , x ) = ¦ d n ( t )x n −1

(7)

n =1

With this assumption, the differential terms in Eq. (3) can now be replaced as shown in Eq. (8) and (9). NK § ∂ 2T · OC ¨¨ 2 ¸¸ ⎯⎯→ d n ( n − 1 )( n − 2 )xk( n−3 ) ¦ ∂ x n =1 © ¹ x = xk

NK § ∂T · ⎯OC ⎯→ ¦ d n ( n − 1 )xk( n −2 ) ¨ ¸ © ∂x ¹ x = xk n =1

(8)

(9)

This is an approach that includes also the boundaries. If there is no need for the outlet boundary, the polynomial can be chosen in a different way excluding this point. Due to the extrapolation of the polynomial, a solution will be obtained anyway. The discrete point’s xk are neither equidistant nor free, but determined in the way that the error order is maximized. This leads to an error order of 2NK, which denotes the maximum order for Runge-Kutta type discretization methods. The inclusion of the Neumann boundary can be implemented properly without affecting the error order by making use of the features of a polynomial. Similar to the multiple shooting approach, the space can be divided into finite intervals, in which the polynomials can be determined. This is related to the Orthogonal Collocation on Finite Elements (OCFE), which is used here in order to solve the resulting ODE after the spatial discretization, because it allows step size

J.C. Schöneberger et al.

904

control algorithms. Hence, it can be used as a common ODE solver. For this purpose, two and three OC points are used. The OCFE approach has the advantage over commercial ODE solver because derivative calculations can be included simply with the method of internal numerical differentiation. This also enables its use for the solution of optimization problems [6].

4. Comparison of the Approaches In Fig. 1 the analytical and the numerical solutions for 6 discrete points are shown. MS and OC give similar results but ML is far away from the solution. The error order, which reduces the influence of the Neumann boundary, affects also the accuracy of the solution at the inner points.

Figure 1: Trajectories of the Numerical Solutions with 6 Discrete Points.

The results for a varying number of discrete points are summarized in Tab. 2. To reach a comparable accuracy with the ML, 100 discrete points are necessary. In order to reach the same error with MS, only 4 points are used, and the computation time is reduced by the factor of 25. The OC with 4 points gives even better results and is 12 times faster than the MS. However, it should be noted that for few points the accuracy of the first derivative is higher with the MS. The results are mean values of 100 calculations. For the MS approach the Matlab® solver ode15s was employed, because it was found to be the most accurate solver for the examined system [7]. Table 2: Results for Different Numbers of Discrete Points

Method Discrete Points

ML 6

ML 100

MS 4

MS 6

OC 4

OC 6

Error in T [K]

5.85

0.06

0.06

0.03

5¨10-3

2¨10-6

Error in dT/dx [K/m]

23.89

0.28

0.02

0.08

Computational Effort [s]

0.08

3.18

0.12

7¨10-3 0.29

5¨10-5 0.02

0.01

Comparable results were obtained for the dynamic simulation of the problem. The additional dimension allows the combination of different approaches. These results are too extensive so as to be presented in detail here. However, the most efficient approach is the combination of the OC in spatial direction with the OCFE in time direction.

An Efficient Discretization Approach for Partial Differential Equations

905

5. Case Studies 5.1. Consecutive Chemical Reaction An isothermal consecutive reaction system with high differences in the reaction rates is analyzed. This case study is used to demonstrate the ability of the OC to solve extremely stiff equation systems. The results are compared to the solutions calculated with the commercial ODE solvers provided by Matlab®. 5.2. Two-Phase Model of a Catalytic Fixed Bed Reactor A catalytic fixed bed reactor is modeled using the conservation equations for the fluid phase and the solid phase. The two phases are coupled by complex heat- and mass transfer correlations. The time and space variant system involves a set of 15 state variables. This gives for 40 discrete points (that are a minimum for the ML) a set of 570 ODE and 150 AE. Here, the advantage of the few discrete points, which are necessary for the OC, is demonstrated.

6. Conclusions Although the ML is the most common and easiest to apply approach it is not giving satisfactory results for chemical reaction systems. Especially in multi phase systems the OC is of advantage because it manages a high accuracy with few discrete points which leads to much smaller algebraic equation systems. The MS also gives good results and is even more accurate than the OC when describing the gradients. Anyhow, to reach the same accuracy in the numerical solution more calculation effort is necessary. These positive effects are multiplied when the methods are applied to more dimensions, e.g. the time. Finally it can be concluded that the OC combined with the OCFE is the best approach for the considered PDE describing chemical reaction systems, especially catalytic fixed bed reactors.

Acknowledgements The authors gratefully acknowledge support from the Max-Buchner Forschungsstiftung.

References [1] H.D. Baehr, K. Stephan, Heat- and Mass Transfer, Springer, 1998 [2] J.C. Schöneberger, H. Arellano-Garcia, S. Körkel, H. Thielert, G. Wozny, Computer Aided Chemical Engineering, Vol. 25 (2008) [3] F. Shakeri, M. Dehghan, Computers & Mathematics with Applications Vol. 56, (2008), P. 2175 [4] R. England, R. Lamour, J. López-Estrada , Applied Numerical Mathematics, Vol. 42 (2002), P. 117 [5] B.A. Finlayson, Nonlinear Analysis in Chemical Engineering, McGraw-Hill, 1980 [6] L.T. Biegler,,Computers & Chemical Engineering, Vol. 8 (1984), P. 243 [7] L.F. Shampine, M.W. Reichelt, SIAM J. on Scientific Computing, Vol. 18 (1997), P. 1