Extending explicit and linearly implicit ODE solvers for index-1 DAEs

Nonlinear differential-algebraic equations (DAE) are typically solved using implicit stiff solvers based on backward difference formula or RADAU formula, requiring a Newton–Raphson approach for the nonlinear equations or using Rosenbrock methods specifically designed for DAEs. Consistent initial conditions are essential for determining numeric solutions for systems of DAEs. Very few systems of DAEs can be solved using explicit ODE solvers. This paper applies a single-step approach to system initialization and simulation allowing for systems of DAEs to be solved using explicit (and linearly implicit) ODE solvers without a priori knowledge of the exact initial conditions for the algebraic variables. Along with using a combined process for initialization and simulation, many physical systems represented through large systems of DAEs can be solved in a more robust and efficient manner without the need for nonlinear solvers. The proposed approach extends the usability of explicit and linearly implicit ODE solvers and removes the requirement of Newton–Raphson type iteration.     

Modified bounded homotopies to enable a narrow bounding zone

The aim of this paper is to introduce modifications that enhance the usability of the bounded homotopy methods proposed by Paloschi [1995. Bounded homotopies to solve systems of algebraic nonlinear equations. Computers and Chemical Engineering 19, 1243-1254; 1997. Bounded homotopies to solve systems of sparse algebraic nonlinear equations. Computers and Chemical Engineering 21, 531-541], especially in the area of chemical engineering. In modified bounded homotopies, the homotopy path is tracked by exploiting mapped variables instead of unmapped ones. Path tracking based on mapped variables makes it significantly easier to track the bounded homotopy path even though the bounding zone has to be narrow. Mapping also improves the bounding effect of bounded homotopies and makes it possible to avoid unreasonable variable values in homotopy path tracking. The performance of the modifications is illustrated with test cases. These examples clearly show that the modifications enlarge the capability and accuracy of bounded homotopies when solving both small- and large-scale sets of nonlinear equations describing chemical engineering problems.     

Towards a quasilinear process simulator—I. fundamental ideas

Process simulation in the steady state is approached by setting up a global set of nonlinear algebraic equations to represent the plant. These equations are solved by a quasi-linear method, whereby a series of linearised systems are solved whose limiting solution is identical to that of the original set of nonlinear equations. With proper care in the formulation of the linearised approximations, second-order convergence can be achieved in a region sufficiently close to the solution. The Newton-Raphson method falls into this category. However, in general, precautions have to be taken in order to get within the convergent region. Precautions necessary to assure satisfactory solution of the equations are described, but the possibility of multiple solutions is not excluded.     

LEGENDRE FUNCTION APPROXIMATIONS OF ORDINARY DIFFERENTIAL EQUATION AND APPLICATION TO CONTINUOUS CRYSTALLIZATION PROCESSES

The technique of the shifted Legendre integral transformation associated with its series expansion is applied to solve population balance equations describing the continuous crystallization processes. The ordinary differential equation of population balance is transformed into a series of algebraic equations of the expansion coefficients which can be easily solved. Illustrative examples for realistic cases are given to simulate the systems. Satisfactory computational results are obtained when are compared with those of the exact solutions. The method is straightforward and the computational algorithm is effective.     

Steady state multiplicity of chemically reacting systems. The method of computation

To analyse the steady state multiplicity of chemically reacting systems of different types it is necessary to find all roots of systems of non-linear equations describing steady states of the systems. This can be solved rather simply in the case of a reducible system to a single equation. We present a numerical method of nonlocal solution for systems of non-linear equations which does not require computing with the Jacobian matrix. The method is illustrated by joint equilibrium absorption of two species for a semi-empirical model of induced inhomogeneity of the catalyst surface and by computation of the steady states of the stirred tank reactor in which the reaction A → B → C occurs.     

LEGENDRE FUNCTION APPROXIMATIONS OF ORDINARY DIFFERENTIAL EQUATION AND APPLICATION TO CONTINUOUS CRYSTALLIZATION PROCESSES

The technique of the shifted Legendre integral transformation associated with its series expansion is applied to solve population balance equations describing the continuous crystallization processes. The ordinary differential equation of population balance is transformed into a series of algebraic equations of the expansion coefficients which can be easily solved. Illustrative examples for realistic cases are given to simulate the systems. Satisfactory computational results are obtained when are compared with those of the exact solutions. The method is straightforward and the computational algorithm is effective.     

Experimental methods in chemical engineering: Computational fluid dynamics/finite volume method—CFD/FVM

Computational fluid dynamics (CFD) applies numerical methods to solve transport phenomena problems. These include, for example, problems related to fluid flow comprising the Navier–Stokes transport equations for either compressible or incompressible fluids, together with turbulence models and continuity equations for single and multi-component (reacting and inert) systems. The design space is first segmented into discrete volume elements (meshing). The finite volume method, the subject of this article, discretizes the equations in time and space to produce a set of non-linear algebraic expressions that are assigned to each volume element—cell. The system of equations is solved iteratively with algorithms like the semi-implicit method for pressure-linked equations (SIMPLE) and the pressure implicit splitting of operators (PISO). CFD is especially useful for testing multiple design elements because it is often faster and cheaper than experiments. The downside is that this numerical method is based on models that require validation to check their accuracy. According to a bibliometric analysis, the broad research domains in chemical engineering include: (1) dynamics and CFD-DEM, (2) fluid flow, heat transfer, and turbulence, (3) mass transfer and combustion, (4) ventilation and the environment, and (5) design and optimization. Here, we review the basic theoretical concepts of CFD and illustrate how to set up a problem in the open-source software OpenFOAM to isomerize n-butane to i-butane in a notched reactor under turbulent conditions. We simulated the problem with 1000, 4000, and 16 000 cells. According to the Richardson extrapolation, the simulation underestimates the adiabatic temperature rise by 7% with 16 000 cells.     

Optimal operation of multicomponent batch distillation—multiperiod formulation and solution

A method is proposed to determine optimal multiperiod operation policies for binary and general multicomponent batch distillation of a given feed mixture, with several main products and intermediate off-specification cuts. A two-level optimal control formulation is presented so as to maximize a general profit function for the multiperiod operation, subject to general constraints. The solution of this problem determines the optimal amount of each main and off cut, the optimal duration of each distillation step and the optimal reflux ratio profiles during each production period. The outer level optimization maximizes the profit function by manipulating carefully selected decision variables. These are chosen in such a manner that the need of specifying the mole fractions of all the components in the products, as required by previous methods is avoided. For values of the decision variables fixed by the outer loop, the multiperiod operation is decomposed into a sequence of independent optimal control problems, one for each production step. In the inner loop, a minimum time problem is then solved for each step to generate the optimal reflux ratio values, reflux switching times and duration of the step. The procedure permits the use of very general distillation models described by differential and algebraic equations, including rigorous thermodynamics if desired. The model equations are integrated by using an efficient Gear's type method, the inner loop optimal control problems are solved using a variational method, and all optimisations are solved using a robust and efficient successive quadratic programming code (Chen, Ph.D. Thesis, Imperial College, 1988).

Several example problems (involving binary and multicomponent mixtures) are used to demonstrate the idea and to show the effect of the cost functions used (in particular the value of the main products) on the optimal solutions.     

A collocation solver for systems of boundary-value differential/algebraic equations

A spline collocation procedure has been implemented to solve systems of multi-point boundary-value problems for mixed order differential/algebraic equations (BVP-DAEs). This implementation is based on the modification of an existing code “COLSYS”, which solves boundary-value ordinary differential equations (BVP-ODEs) alone. With minor changes to the original COLSYS code, the resulting procedure offers robust solutions to systems of DAEs, in addition to ODEs. A numerical example with application in chemical engineering is shown to demonstrate the superior ability of this procedure over other software.     

Simultaneous optimization and solution of systems described by differential/algebraic equations

Engineering approaches to the solution of constrained variational problems often involve converting the problem into a nonlinear programming (NLP) problem and solving it using current NLP methods. These methods usually use a sequential optimization and solution strategy. We propose a method, using piecewise constant functions for the independent variables, that combines the technologies of quasi-Newton optimization algorithms and global spline collocation to simultaneously optimize and integrate systems described by differential/algebraic equations. A computer implementable algorithm is discussed and three test problems are solved. The algorithm allows the solution of a more general class of optimization problems than previous methods employing this strategy.     

A model-based control framework for industrial batch crystallization processes

Dynamic optimization is applied for throughput maximization of a semi-industrial batch crystallization process. The control strategy is based on a non-linear moment model. The dynamic model, consisting of a set of differential and algebraic equations, is optimized using the simultaneous optimization approach in which all the state and input trajectories are parameterized. The resulting problem is subsequently solved by a non-linear programming algorithm.     

Concentration and temperature profiles for complex reactions in porous catalysts

A method of computation for stationary concentration and temperature profiles in porous catalysts was developed. These calculations can be used with highly complex reactions and with any type of rate equation or transport model. The differential equations are transformed by the collocation method to non-linear algebraic equations. These are solved by a method which uses the time-derivatives of the variables. The method finds stationary points for all problems which can be formulated as ? = f(y) and is demonstrated on two examples.     

Fluid dynamic stability of fluidised suspensions: the particle bed model

The equations of change for a fluidised suspension are closed by application of a simple model for the interaction of a single particle with the fluid; this relationship delivers the total force on elements of the particle phase including a fluid dynamic formulation for the 'particle phase pressure gradient'; numerical solutions of the non-linear equations, for small imposed perturbations of voidage, show clearly the development of shocks (bubbles) for unstable (aggregate) beds and the return to the homogeneous state for stable (particulate) systems; the linearised equations yield the general criterion for bed stability as a simple algebraic expression and quantitative predictions of disturbance propagation velocities in good agreement with published data.     

A solution method for solving coupled nonlinear boundary-value problems

A Strum-Liouville integral transform technique is novelly applied to solve system of coupled nonlinear boundary-value problems approximately. The systems of differential equations consist of a linear differential operator and a nonlinear function of the dependent variables. To illustrate the potential of this technique we consider an example which comes from the modeling of diffusion and nonlinear chemical reaction systems in chemical engineering. The approximate solutions obtained by our technique agree surprising well with the numerically exact solutions obtained by the orthogonal collocation technique. To improve the approximation an iteration scheme in transform space is also defined.

Scope—Today, mathematical modeling of physical phenomena often produces (single or coupled) nonlinear differential equations. The true physical situation can, in many cases, be more closely described if the differential equations are allowed to be nonlinear. However, nonlinear differential equations are generally too difficult to be solved analytically apart from a few “tricks” or substitutions which apply only to a handful of equations [1]. An alternative approach is to look for a method which will reduce the problem, via analytical techniques, to a point where a “simple” computer program can solve the rest of the problem. The method introduced in this paper belongs to this class of solution techniques.

The method, which in this paper is applied to solving coupled nonlinear boundary-value problems, is a generalization of an idea in a paper by Do and Bailey [2] who apply it to a single nonlinear differential equation of boundary-value type. The equations, to which the technique is applied, arise from Fick's law diffusion into a porous solid and nonlinear reaction within the solid.

The solution method employs a Strum-Loiuville integral transform and to account for the nonlinear part an approximation is introduced. An iteration scheme is defined to improved the accuracy of the solution. The system of coupled nonlinear differential equations is reduced to a system of coupled nonlinear algebraic equations which is solved using a Newton-Raphson process. Finally, the solution is expressed as an infinite series, which is summed using a computer.

In response to papers by Do and Bailey [3] and Do and Weiland [4], Jerri [5] has tried to put this method on a more mathematical footing, and he shows that this method is a special case of a more general technique he has devised. Jerri uses the idea of Fourier transforms and convolution products to justify his method. The results for the example he considered are good, but he did not state how many iterations he required to obtain the solutions reported.

Conclusions and Significance—This paper has presented a very powerful method of solving boundary-value problems with linear operators and a nonlinear function of the dependent variable. The method works well for a single equation or coupled equations and can handle any kind of nonlinear function. We have shown through extensive numerical calculation the accuracy of this solution method, where the accuracy is measured in terms of a ratio of norms. In most cases an error of 4% can be achieved with just one iteration (Tables 2 and 3). Even though the present method has been applied to problems which have arisen from the modeling of chemical engineering problems, it would also be applicable to differential equations arising in other areas, provided they are of the same form.     

Multi-scale analysis of exotic dynamics in surface catalyzed reactions—II: Quantitative parameter space analysis of an extended Langmuir-Hinshelwood reaction scheme

The multi-time-scale approach to modelling heterogeneously catalysed reaction systems developed in a companion paper [1] is applied to a simple Langmuir-Hinshelwood mechanistic scheme for an overall reaction 2A(g) + B(g)→products supplemented by a slow buffering step involving a chemisorbed inert species. Based on leading-order approximations to the full governing equations and the use of non-linear algebraic equations in the definition of system manifolds and prediction lociin appropriate phase planes, regions of different dynamic behaviors (stationary states, single and multipeak oscillations) are qualitatively and quantitatively delineated as a function of residence time. The bifurcation predictions from the leading-order approximations achieve remarkable agreement both with numerical integration results as well as at local bifurcation points obtained via Hopf analysis of the full equations. Numerical simulations reveal some very complex (“chaotic”?) behavior in the multipeak region, with an indication that the onset of chaos may be via intermittency. It is expected that the application of multi-time-scale analysis to specific well-characterized reaction systems will provide fuller understanding of past and adequate theoretical guidance for future work.     

Automatic differentiation‐based quadrature method of moments for solving population balance equations

The quadrature method of moments (QMOM) is a promising tool for the solution of population balance equations. QMOM requires solving differential algebraic equations (DAEs) consisting of ordinary differential equations related to the evolution of moments and nonlinear algebraic equations resulting from the quadrature approximation of moments. The available techniques for QMOM are computationally expensive and are able to solve for only a few moments due to numerical robustness deficiencies. In this article, the use of automatic differentiation (AD) is proposed for solution of DAEs arising in QMOM. In the proposed method, the variables of interest are approximated using high‐order Taylor series. The use of AD and Taylor series gives rise to algebraic equations, which can be solved sequentially to obtain high‐fidelity solution of the DAEs. Benchmark examples involving different mechanisms are used to demonstrate the superior accuracy, computational advantage, and robustness of AD‐QMOM over the existing state‐of‐the‐art technique, that is, DAE‐QMOM. © 2011 American Institute of Chemical Engineers AIChE J, 2012     

The application of the parametric continuation method for determining steady state diagrams in chemical engineering

The paper discusses a simple method of using the parametric continuation method to designate complex diagrams of steady states generated by systems and apparatuses used in chemical engineering. The main advantage of the discussed approach is the fact that it does not require the installation of huge professional IT systems. The deliberations are illustrated by examples of a tank reactor model described by the algebraic equations and of a tubular reactor model with longitudinal dispersion, described by the differential equations of the second order. Both models render solutions in the form of the so called multiple steady states.     

Direct methods for consistent initialization of DAE systems

Rigorous techniques for consistent initialization of differential–algebraic systems (DAEs) are in general difficult to apply. Even worse, since initialization algorithms generally require the solution of a set of non-linear algebraic equations, they may not even converge. Hence, simpler approaches for initialization are of considerable interest in the numerical implementation of DAE models. In this contribution, it is proposed a simple and effective initialization technique, which is based on the assumption that the dynamic response of a mathematical model may be approximated by the response to discontinuous perturbations of a similar system, originally in steady state. Determining the analogous system and smoothening adequately those perturbations are the key aspects for the consistent initialization of the DAEs involved. The technique presented has been tested in several examples, always yielding encouraging results.     

Sliding motion of discontinuous dynamical systems described by semi-implicit index one differential algebraic equations

We have proposed an approach to derive a continuous system of differential algebraic equations (DAE) of index one that is dynamically equivalent to a discontinuous index one DAE system. This involves augmenting the convex combination of the ordinary differential equations with the algebraic equations from individual models. This result is proved by using the implicit function theorem. This procedure is illustrated with the help of an ideal gas-liquid system in which the algebraic variables can be expressed as explicit functions of differential variables. It is also demonstrated with an example from a soft-drink manufacturing process, in which, it is difficult to express the algebraic variables as explicit functions of differential variables. Through computer simulation, it is shown that the equivalent dynamic DAE system and the discontinuous DAE system have identical solutions. The proposed method is several orders of magnitude more efficient than the procedure that works with the discontinuous system of DAEs.