Numerical solution technique

The present chapter provides an overview of several numerical techniques that can be used to solve model equations of ordinary and partial differential type, both of which are frequently encountered in multiphase catalytic reactor analysis and design. Brief theories of the ordinary differential equation solution methods are provided. The techniques and software involved in the numerical solution of partial differential equation sets, which allow accurate prediction of nonreactive and reactive transport phenomena in conventional and nonconventional geometries, are explained briefly. The chapter is concluded with two case studies that demonstrate the application of numerical solution techniques in modeling and simulation of hydrocar-bon-to-hydrogen conversions in catalytic packed-bed and heat-exchange integrated microchannel reactors. 

1 Techniques for the numerical solution of ordinary differential equations 

Consider the following ordinary differential equation (ODE) and the initial condition  

Numerical solution of the ODEs describing initial value problems is possible by explicit and implicit techniques, which are described in the Sections 11.1.1 and 11.1.2, respectively. It is worth noting that the techniques are formulated for the solution of a single equation (Equation 11.1), but they can be used for solving multiple ODEs as well. Theoretical background of these methods as well as their stabilities are described elsewhere [1,2] and will not be discussed here. 

---

The use of PB modeling by practitioners has been hmited for two reasons. First, in many cases the kinetic parameters for the models have been difficult to predict and are veiy sensitive to operating conditions. Second, the PB equations are complex and difficult to solve. However, recent advances in understanding of granulation micromechanics, as well as better numerical solution techniques and faster computers, means that the use of PB models by practitioners should expand. 

The three modes of numerical solution techniques are finite difference, finite element, and spectral methods. These methods perform the following steps ... 

Problems of inclusions in solids are also treated by exact elasticity approaches such as Muskhelishvili s complex-variable-mapping techniques [3-9]. In addition, numerical solution techniques such as finite elements and finite differences have been used extensively. 

This equation is coupled to the component balances in Equation (3.9) and with an equation for the pressure e.g., one of Equations (3.14), (3.15), (3.17). There are A +2 equations and some auxiliary algebraic equations to be solved simultaneously. Numerical solution techniques are similar to those used in Section 3.1 for variable-density PFRs. The dependent variables are the component fluxes , the enthalpy H, and the pressure P. A necessary auxiliary equation is the thermodynamic relationship that gives enthalpy as a function of temperature, pressure, and composition. Equation (5.16) with Tref=0 is the simplest example of this relationship and is usually adequate for preliminary calculations. 

As noted earlier, a,>, will usually be independent of r, but the numerical solution techniques that follow can easily accommodate the more general case. The radial boundary conditions are... 

Analytical solutions of the self-preserving distribution do exist for some coalescence kernels, and such behavior is sometimes seen in practice (see Fig. 40). For most practical applications, numerical solutions to the population balance are necessary. Several numerical solution techniques have been proposed. It is usual to break the size range into discrete intervals and then solve the series of ordinary differential equations that result. A geometric discretization reduces the number of size intervals (and equations) that are required. Litster, Smit and Hounslow (1995) give a general discretized population balance for nucleation, growth and coalescence. Figure 41 illustrates the evolution of the size distribution for coalescence alone, based on the kernel of Ennis Adetayo (1994). 

A comparison of the benefits and drawbacks of common numerical solution techniques for complex, nonlinear partial differential equation models is given in Table II. Note that it is common and in some cases necessary to use a combination of the techniques in the different dimensions of the model. 

Numerical Solution Techniques for Partial Differential Equations Arising in Packed Bed Reactor Modeling... 

The resultant equations are non-linear and in this general case numerical solution techniques must be used. However, there exists a special case where an analytical solution may be obtained. If the increase in biomass concentration during flow through the reactor is small then an average value for the biomass concentration, independent of the distance Z along the fermenter, may be used. The material balance for the substrate over the reactor element may then be written ... 

The classification of numerical solution techniques lends itself excellently to the system theory classification as well. 

Typically, the numerical solutions techniques used are very specific to the problem. Particularly challenging problems include moving front problems where concentration profiles, for example, may vary widely over a short distance but may not change much at other spatial locations. The spatial discretization must be small close to the front for accuracy and numerical stability, but must be larger at other locations to reduce computation time. Various adaptive grid techniques to change the spatial step sizes have been developed for these problems. One of the more common codes to solve fluid-flow-related problems is FLUENT. 

However because the rate law expressions for iron sulfide formation (equations 8 and lO) are non-linear the differential equations for H2S and the iron sulfides are not amenable to explicit solution Thus it is important to develop an equation for that can be incorporated in a numerical solution technique such as that of Runge-Kutta (15) Fortunately an appiropriate differential equation for can be developed firom charge balance considerations Here it is assumed that dissolved substances other than those listed in Table I are not affected by diagenesis If this is true, then a charge balance difference equation can be written (16) ... 

Numerical Solutions For many practical applications, numerical solutions to the population balance are necessary. Several numerical solution techniques have been proposed. It is usual to break the... 

Several physicochemical models of ion exchange that link diffuse-layer theory and various models of surface adsorption exist (9, 10, 14, 15). The difficulty in calculating the diffuse-layer sorption in the presence of mixed electrolytes by using analytical methods, and the sometimes over simplified representation of surface sorption have hindered the development and application of these models. The advances in numerical solution techniques and representations of surface chemical reactions embodied in modem surface complexation mod-... 

We saw above that the concentration gradient at an electrode will be linear with respect to the spatial coordinate perpendicular to the electrode surface if the anode/cathode cell were operated at a constant current density and if the fluid velocity were zero. In actuality, there will always be some bulk liquid electrolyte stirring during current flow, either an imposed forced convection velocity or a natural convection fluid motion due to changes in the reacting species concentration and fluid density near the electrode surface. In electrochemical systems with fluid flow, the mass transfer and hydrodynamic fluid flow equations are coupled and the solution of the relevant differential equations is often a formidable task, involving complex mathematical and/or numerical solution techniques. The concept of a stagnant diffusion layer or Nemst layer parallel and adjacent to the electrode surface is often used to simplify the analysis of convective mass transfer in... 

The need to use specific numerical values for the rate constants and initial conditions is a weakness of numerical solutions. If the specific values change, then the numerical solution must be repeated. Analytical solutions usually apply to all values of the input parameters, but special cases are sometimes needed. Recall the special case needed for Uo = bo in Equation 1.32. Numerical solution techniques do not have this problem, and the problem of specificity with respect to numerical values can be minimized or overcome through the judicious use of dimensionless variables, as... 

This equation is coupled to the component balances of Equation 3.9 and with an equation for the pressure, for example, one of Equations 3.18, 3.19, or 3.21. There are iV + 2 equations and some auxiliary algebraic equations that need to be solved simultaneously. Numerical solution techniques are similar to those used in Section... 

Analytical expressions similar to those for spherical particles have been derived for infinite-length cylinders in perpendicular incidence as well as in oblique incidence, for elliptic cylinders, and for spheroids (see Refs. 168 and 169). With increasing complexity of the shape of the particle, even with as little change as from sphere to spheroid, the analytical solution to the problem becomes formidable. Then, the use of numerical solution techniques may be preferable to analytical techniques. 

For systems with transfer functions that are very difficult to factor and consequently very hard to complete the frequency response analysis, Luyben [Ref. 7] discusses various numerical solution techniques. He has also included a computer program in FORTRAN which uses the stepping technique to develop the Bode and Nyquist plots for a distillation column. More details on the philosophy of the Ziegler-Nichols tuning method can be found in the original work ... 

---