Matlab Ode Solver for Kinetic Modelling and Simulation of Claus Process Reactions

This paper is focused on a kinetic study and parameter analysis of a chemical reactions in a Claus process of a sulphur recovery unit at Mellita Gas plant. The chemical reactions that take place either in H 2 S conversion furnace or in a catalytic converters are studied by an analysis procedures that linked a kinetic modelling to a reactors performance. The kinetic data of a chemical reactions in a Clauss process were achieved a ProMax (software) simulator. A PROMAX process simulator is installed and connected to a Distribution Control System (DCS) through instrument analysers and transmitters of a Claus process in a sulphur recovery unit at Mellita Gas plant. A plant kinetic data were used to estimate a reactions parameters by a MATLAB code apply a nonlinear curve fitting regression. The developed kinetic model for a multiple reactions that appear in a several differential equations was implemented and solved numerically by ODE45 solver that available in a MATLAB library. The kinetic model results were visualized and displayed graphically by a MATLAB plot functions.

Sweetened acid gas, due to stringent environmental regulations must be processed to reduce further hydrogen sulphide. Nowadays, a new technology to convert acid gas to syngas (AG2STM) is proposed [3,4], but the modified Claus process is still more common than others due to higher sulphur conversion and operation at higher H 2 S concentration and capacity [5]. This implies the importance of declaring a kinetic model for its description and optimization of its operation for better performance and higher conversion, consequently lower sulphur dioxide emissions. There are several researchers are developed a kinetic models of a Clauss process multiple reactions in either a thermal furnace or a catalytic reactors. All researchers validate a kinetic model results with industrial experimental data [1][2][3][4][5]. Among the mentioned researchers Pahlavan et al. [2] simulated a Clauss reactions by a PROMAX V2.0 process simulator.
In this study, it is attempted to perform kinetic analysis of clauss process multiple reactions in a thermal furnace of Clauss process of Mellita gas plant. Kinetic analysis consider a simple reaction mechanism include three multiple reactions.

2
MATLAB PROGRAM OVERVIEW: MATLAB ( Matrix Laboratory) is a generalized mathematical software program with many features such as numerical computations, Simulink and control analysis, visualization of results is available from Mathworks [6] As with any software program, there are few ' rules' and 'codes' which must be followed. The structure of MATLAB consists of M-file or command file type, M-file, there are few general statements to be made about M-file or MATLAB code, is a collection of commands that are executed sequentially, commands can be mathematical operation, function call, flow control statement, and calls to the functions or scripts, the execution of m-file program can be controlled from command file windows or other m-file code. There are two types of m-files, functions and scripts. Function has variables that can passed into and out of the function. Any other variables used inside the function are not saved in memory when the execution of function is finished. Scripts on the other hand, save all their variables in the MATLAB work space. Functions and scripts should be nominated. The first line of a function must contain a function declaration, using a following format:

Built-in ODE Solvers in MATLAB:
As with most differential equation solvers, there a number of different numerical methods that can be used. Essentially, they are split non-stiff and stiff differential equations. A stiff differential equation is one that has a rapidly varying derivative. That is, the derivative changes rapidly and is therefore susceptible to relatively large errors in the numerical analysis. There are also choices within these two broad categories and these choices are based on the type of errors that can tolerated.

Table1. Ode solvers in a MATLAB library
Stiff ode solvers Non-stiff ode solvers

Ode45
Ode23 Ode113 The numbers associated with the names indicates the "order" of the numerical method, the " s" indicates the stiff algorithm and the "t" indicates an implementation of the trapezoidal rule. The best bet is to initially try ode45, which utilizes a 4 th -5 th order Runge-Kutta algorithm. If this takes an unduly long time, then try one of the stiff solvers. If you find that you are having difficulty, the most likely reason is that you have not posed the problem correctly [9].

MATLAB Method of Lines:
Method of lines is based on the concept of converting the partial differential equation (PDE) into a set of ordinary differential equations (ODE) by discretizing only the spatial derivatives using finite difference and leaving the time derivative unchanged [10].
The method of lines (MOL) is generally recognized as a comprehensive and powerful approach to the numerical solution of time-dependent partial differential equations (PDEs). This method proceeds in two separate steps:  spatial derivatives are first replaced with finite difference, finite volume, finite element or other algebraic approximations.  the resulting system of, usually stiff, semi-discrete ( discrete in spacecontinuous in time) ordinary differential equations (ODEs) is integrated in time [11].

Model Validation and visualization of Parameters
It easy to visualize and display results in MATLAB graphically. The plot function is used to create simple plots. The command syntax is: plot(x 1 ,y 1 ,format1,x 2 ,y 2 ,format2,......), x 1 and x 2 are independent variables (usually time), y 1 , y 2 are dependent variables. The format1 and format2 arguments are short combination of characters containing the plot formatted commands [12].

CASE STUDY : CLAUSS RECOVERY UNIT
The modified Claus process is one of the major processes which converts toxic hydrogen sulfide to elemental sulfur from acid gas during sour natural gas processing. The modified Claus process in which hydrogen sulfide is converted to sulfur can be depicted by the following simple overall reaction [1]: The sulfur recovery in the modified Claus process is obtained via the two reaction steps including thermal and catalytic steps. The first step is a thermal step where one third of hydrogen sulfide (H2S) is partially oxidized and sulfur dioxide (SO2) is produced at very high temperatures.
The second step involves the reaction between the un-reacted hydrogen sulfide and sulfur dioxide over a catalytic bed at lower temperatures: Conversion of H 2 S to SO 2 in the modified Claus process is conducted in the reaction furnace (RF) which is a refractory lined cylindrical vessel. Since the combustion process in the RF provides substantial sulfur recovery and affects the entire plant emissions, any modification in the RF affects sulfur conversion and emissions from Claus plant . Therefore any improvement in the reaction furnace has direct impact on reduction of sulfur emissions to the atmosphere [14].

Experimental Data by a PROMAX process simulator:
According to real operation condition in Sulfur recovery unit J 2 , and K 2 at Mellitah complex, the field data have been collected from transmitters and analyzers in distribution control system (DCS). Nowadays, simulation software is one of the most powerful tools in improving the operating conditions of industrial units. it is necessary to improve these software packages or use a better method. Environment has an essential role in Claus unit and its performance affects process efficiency and other equipment performance, so, in the current work, SRU simulation is discussed. Many useful softwares have been improved to study chemical and petroleum industry in terms of the process modelling and simulation, these softwares can be used to simulate the oil and gas processing to enhance the process operation ‫العلوم‬ ‫األسمريت:‬ ‫الجامعت‬ ‫مجلت‬ ‫والتطبيقية‬ ‫األساسية‬ Journal of Alasmarya University: Basic and Applied Sciences and increase the productivity. ProMax and Aspen HYSYS are the most common used software. ProMax has been used as a simulator in this research. It is a powerful engineering simulation software, designed with respect to the program architecture, interface, engineering capability and interactive operation, therefore, ProMax process software is known as one of the best process simulation software. It applies the conception of fluid package that contains all the required data and information to perform chemical and physical properties calculation, therefore, ProMax property is the driver of the process and it must be chosen carefully by allowing the definition of all the information such as property package, components, hypothetical components, interaction parameters, reactions, etc. [15].

Development Of a Reactor Model
Partial differential equations(PDE) play an important role in the mathematical modelling of a biological, physical, or chemical systems. Many engineering problems involve diffusion and reaction processes. These problems include the diffusion of pollutants into the air, water and earth. Each area has its own physics, chemistry, biology and ecology. The modelling of these areas often leads to a partial differential equations (PDE). Such equations may be linear, in which the analytical solutions is often achievable. However, many of diffusion and reaction problems are non-linear and numerical solutions are used to solve and understand these models [16].
There has been a considerable amount of research devoted to developing numerical techniques for solving boundary value problems (PDE) in the chemical engineering literature. The large size of this body of literature reflects the prevalence of (PDE) in chemical process modelling involve diffusion and reaction where their parameter distribution defined on finite domain [17].
Multiple reactions system consisting of a parallel or series reactions that takeplace simultaneously. The development of reactor modelling for multiple reactions need to follow the steps of kinetic analysis that originally published by Fogler [18] and presented below.
Steps of analyzing multiple reactions: • Number each rxn.
. Write the mole balances for each species.
. Determine the rate law for each species in each reaction.
. Relate the rate of reaction of each species to the species for which the rate law is given for each reaction.
. Determine the net rate of formation of each species. . Express rate laws as a function of concentration ,Cj, for the case of no volume change.
. Express rate laws as a function of moles, Nj(batch), or molar flow rates, Fj(flow) when these volume change with reaction.
. Combine all of the above and solve the resulting set of coupled differential equations ( P.F.R, P.B.R, Batch) or (CSTR) The mathematical manupulations translate the multiple reactios by applying basic of chemical reaction kinetics to translate a multiple reactions equations to drdinary differential equations( ODE) that can be solved numerically by a MATLAB ode solver [18]. In this research, the reactor modelling of a reactor furnace was analysed by considering adiabatic plug flow, and the pressure drop is neglected, Consider the following set of reactions [20]: For simplicity the components in the previous set of reactions is renamed as, A,B, C,…., the new set of reactions is written as: A + 0.5B → 0.5C + D A + 1.5B → E + D A + 0.5E ↔ 0.75C + D Write a relative rate of recation for each component in each reaction by applying the general rule that shown below as: The relative of reaction for 1 st reaction is written as: For 2 nd reaction as: For 3 rd reaction as: The general reaction rate expression for selected component of known frequency factor and activation energy are written as:

Kinetic Parameter Estimation
In a non-linear least square analysis, searching for a parameter values that minimize the sum of squre of the diffrences between the measured values and predicted values for all the data points [19]. MATLAB software is used to illustrate this technique for a clauss pocess reactions in sulphur recovery unit at Mellia gas plant.

RESULTS AND DISCUSSION Numerical Solution
MATLAB have a high performance and flexibility for numerical computations and results visualisations. There are several ordinary differential equation solver (ODE) in MATLAB program liberary. The developed kinetic model of differential equations (one per each component) were implemented and solved using ODE45 MATLAB SOLVER that constructed based on four point Ruge-Kuta numerical technique coupled with a method of lines. The MATLAB code was developed by using scripts and built in functions. Figure 1 shows the flow chart of numerical solution algorithm for reactor model  In this work, the influence of air flow rate is studied, by solving a kinetic model for combustion of hydrogen sulphide in a furnace. A typical trend in the time of concentration of residual hydrogen sulphide, intermediates, and carbon dioxide is shown in Figure 3.  As can be seen in Figure 4, the H 2 S concentration decreased since this compound was adsorbed and transformed to intermediate products. In Figure 5, is shown a plot of hydrogen sulphide and sulphate dioxide ratio. it can be noticed from these two figures 4 and 5, their concentrations reaches a maximum at the beginning and then decreased. The time needed to reach a maximum and then decreased is faster during a period of several hours simulation. The air to acid gas flow ratio, initial feed temperature, and steam pressure in the WHB are operating variables that can be changed in operating Claus furnace in order to maintain H 2 S/SO 2 ratio to maximize sulfur conversion efficiency In this paper, kinetic modelling and multiple reaction analysis were applied for Clauss process multiple reactions. Clauss process is used to treatment of acid gases by reducing its concentration before it released to the atmosphere.
The ordinary differential equations of a kinetic model of Clauss process multiple reactions is solved numerically by ode45 solver that available in a MATLAB liberary. The numerical solution was implemented in a MATLAB code by coupling ode45 solver with a method of lines technique. The visualisation of a kinetic model results is shown in graphical representation by a MATLAB plot functions.