---

Industrial Mathematics / Mathématiques industrielles
(Org: Biao Huang, Yanping Lin and/et Shijie Liu)

---

CLAUDIA CAIA, Department of Mathematical and Statistical Sciences, University of Alberta, Edmonton, Alberta  T6G 2G1
Finite element method for 3D multiphase flow modelling: validation and numerical results

A fundamental problem in multiphase flows is the prediction of the velocity and the gas hold-up profiles for three-dimensional laminar bubbly two-phase flow in a pipe.

The derivation of the multiphase flow model assumes an interpenetrating continua hypothesis. The set of continuity and momentum equations is derived through statistical averaging of each phase. The expressions for the coefficients are chosen so that the numerical simulations to be in good agreement with some available experimental data. The simple Lagrangian tracking models, even though more accurate, are limited to small void fractions.

The time discretization based on a splitting scheme. The convective terms in both the continuity and momentum equations are discretized by a characteristic method. The remaining problem id discretized by a scheme that stems from a pressure-correction projection scheme for incompressible flows. The incompressibility in the present case is imposed on an averaged velocity field. The spacial discretization is accomplished wirh a stabilized version of the finite element method. The standard element consists of 10 velocity nodes and for 4 pressure nodes. The stabilized finite element is based on a two-level hierarchical decomposition of the approximation space. It introduces two more terms in our formulation: a stabilization term and a shock capturing term. Our approach yields fully symmetric matrices for both the velocity and the pressure sub-problems that are solved by a conjugate gradient method.

To validate the efficiency of the stabilization and shock capturing terms we used the so called sedimentation problem, which is a difficult benchmark problem and we compared our results with other authors. Since the final goal of this study is to perform numerical simulations of bubble columns, and since the model is stronlgy nonlinear, we tested its performance on a problem involving a nonlinear void fraction wave propagation in a three-dimensional vertical cylindrical pipe. Comparison with other authors is included. Theoretically, there exists a possibility that the initially one-dimensional wave generates, via the nonlinear interactiom, three-dimensional waves that could dissipate enough energy to yield a stable (3D) wave solution. This hypothesis has been verified through humerical simulations. In real bubble columns, the bubbles expand due to the loss of hydrostatic head, and therefore we incorporate an additional term in the continuity equation for the dispersed phase, that accounts for this increase in the void fraction. We determined the magnitude of this term so that the results to match the available experimental data. The influence of the wall-force and the lateral-lift force also used by other authors is investigated. We also determined the threshold of instabilty and compare it with the results provided by other authors. The experiments show that close to the transition to agitated flow, the flow can no longer be regarded as uniform, since the bubbles become clustered towards the center of the tube. That is the point where our numerical model is no longer valid.

JINWEN CHEN AND ZBIGNIEW RING, National Center for Upgrading Technology
Modeling fixed-bed chemical reactors with high heat release

Fixed bed reactors, in which gases flow through cylindrical beds packed with catalyst particles to undergo chemical reactions, are widely used in chemical, petrochemical and petroleum refining processes. The design and performance optimizations of these commercial processes rely on successful reactor simulation using mathematical models based on mass, heat and momentum balances. Such reactor models consist of a set of partial/ordinary differential equations (PDEs/ODEs) with well-defined initial and boundary conditions. However if a large amount of heat is released due to chemical reactions, the numerical computation of these PDEs/ODEs could become difficult due to the stiffness of the equations, which is caused by the highly non-linear reaction term coupled with the strong heat release. Consequently, serious challenges could be expected in process design, control and operation.

In this study, a mathematical reactor model was developed to simulate the steady-state and dynamic behaviors of a pilot plant fixed-bed reactor used for hydrotreating of cracked naphtha, a highly exothermic process. An efficient numerical method was used to overcome the stiffness of the model equations. The model, consisting of three PDEs, describes the profiles of temperature and reactant concentration in the reactor. The three PDEs, together with the corresponding initial and boundary conditions, were converted into dimensionless forms which were further transformed into a set of ODEs by discretion in axial and radial directions using a finite difference method. A second order, backward finite difference method was used for the discretion to overcome the stiffness of the PDEs in the axial direction caused by the highly non-linear reaction term, coupled with the high heat release due to olefin hydrogenation. The ODEs were solved with an integration package LSODA. Under steady state, these ODEs became non-linear algebraic equations and were solved with the Gaussian-Seidel iteration method. These computation methods proved to be suitable and efficient in dealing with the problem in question. The model successfully simulated the temperature and reaction concentration profiles in the reactor under steady-state operation. It also reasonably captured the dynamic behaviors of reactor start-up and the reactor temperature response to a step change in gas flow.

CHUWANG CHENG, University of Alberta
Fault tolerant control systems with stochastic integral quadratic constraints

This paper concerns the fault tolerant control problem, in which the system component failure and the fault detection and isolation (FDI) scheme are characterized by two Markov random processes. The model uncertainties and noise/disturbance are also treated in the same framework. A more general description of system uncertainties is proposed in this paper and a less conservative approach for dealing with such uncertainties is presented. In addition to the stochastic stability of the fault tolerant control system, the system performance using a Stochastic Integral Quadratic Constraint (SIQC) is also achieved.

KEVIN DAVIES, Lakehead
On symmetric pyramidal finite elements

(joint work with Liping Liu and Michal Krízek)

Finite element method is considered to be one of the most effective numerical methods for solving problems in mathematical physics and engineering. Pyramidal elements are used for a face-to-face connection of tetrahedral finite elements with block elements. This is a very useful tool for joining tetrahedral meshes with 3D rectangular meshes, which has a lot of practical applications in conforming finite element discretizations of domains in which some parts can be decomposed into blocks.

In [2], the author presents a special family of mortar pyramidal finite elements. They are composed of two tetrahedra which cause an artificial anisotropy in solving isotropic problems (compare with [1, p. 38]). In this talk we introduce a new family of pyramidal elements, which are composed of four tetrahedra. They have more symmetries and the same number of degrees of freedom (function values at vertices) as pyramidal elements in [2], but aviod the artificial anisotropy. We also show the results by computing the numerical solutions for Poisson equation by decomposing a unit cube into N×N×N subcubes and each subcube into 6 pyramids. Finally, we compare the results using our highly symmetric pyramidal elements and nonsymmetric elements.

References

[1]

M. Krízek and P. Neittaanmäki, Finite Element Approximation of Variational Problems and Applications. Pitman Monographs and Surveys in Pure and Applied Math. 50, Longman Scientific & Technical, Harlow, 1990.

[2]

C. Wieners, Conforming discretizations on tetrahedrons, pyramids, prisms and hexahedrons. Univ. Stuttgart., Bericht 97/15, 1997, 1-9.

FENG DING, Department of Electrical and Computer Engineering, University of Alberta, Edmonton, Alberta  T6G 2V4
Parameter estimation for dual-rate systems with finite measurement data

This paper on parameter estimation is motivated by practical consideration that the output sampling rate is often limited and that the data length is finite. In particular, for dual-rate systems in which the output sampling period is an integer multiple of the input updating period, we obtain frequency-domain models, study the properties of the least squares type algorithms in detail in the stochastic framework, and derive convergence rates and upper bounds of parameter estimation errors (PEE) for the least-squares (LS) algorithm, instrumental variable least squares (IVLS) algorithm, and forgetting-factor least squares (FFLS) algorithm, using directly the finite dual-rate input-output data. The analysis indicates that the mean square PEE upper bounds of LS and IVLS algorithms are proportional to 1/k and converge to zero as the data length k increases, and the PEE upper bound of the FFLS algorithm approaches a finite constant. Finally, we illustrate and verify the theoretical findings with example systems, including an experimental water-level system.

MICHAEL G. FORBES AND FRASER FORBES, University of Alberta
Regulating discrete-time stochastic processes: Focusing on the probability density function

In the chemical and manufacturing industries, a common challenge is the regulation of processes subject to unmeasurable, random disturbances. These disturbances may by present by a number of different causes including non-uniform feedstocks, sensor noise, fluctuating ambient conditions, etc. The majority of work in the chemical process control literature focuses on regulation of linear processes subject to disturbances with a Gaussian distribution. Standard approaches to controller design involve quadratic objective functions and result in linear feedback control laws. An important advantage of the linear-Gaussian-quadratic framework is that in the derivation of the control law design techniques, it is not necessary to make explicit use of the probability density function (PDF) of the closed-loop process. In alternate frameworks (nonlinear processes, nonquadratic objectives, etc.), explicit knowledge of process PDF becomes necessary and nonlinear feedback control laws result.

Here, the topic of regulatory feedback control design for discrete-time stochastic processes is approached with a focus on the PDF of the closed-loop process. The discrete-time process model is given by equation (1),

x(t+1)=f æ è x(t),u(t) ö ø +w(t)

(1)

where x(·) Î Âⁿ is the process `state'; u(·) Î Â^m is the input; w(·) Î Â^w is an independent, identically distributed random noise sequence with a PDF p_w(·); f (· , ·) Î Âⁿ is a smooth mapping; and t is an integer-valued time index. Note that equation (1) can represent input-output models of processes as well as those that recognize a true internal process state. Further, equation (1) can represent a wide range of linear and nonlinear processes including those with delay, coloured noise, high-order dynamics, multiple inputs and multiple outputs.

The process given by equation (1), when controlled by a `state' feedback, u (t)=k(x(t)), has the closed-loop form:

x(t+1)=f æ è x(t),k æ è x (t) ö ø ö ø +w(t)= ~ f   æ è x(t) ö ø +wt(t)

(2)

and the mapping [(f)\tilde] (·) is referred to as the closed-loop dynamics.

In this work, equation (3), that gives the stationary PDF, p(·), for the closed-loop process (2), is examined in detail.

p(z)= ó õ Ân  pw æ è z- ~ f   (x) ö ø p(x) dx

(3)

While this integral equation gives an exact implicit relationship between the closed-loop dynamics, [(f)\tilde] (·), and the PDF of the closed-loop process, in general no explicit solution exists. An approximate solution to this equation is developed by parameterization of both the process PDF and the closed-loop dynamics. These parameterizations allow for reduction of the integral equation to a set of algebraic equations in the parameters. By choosing values for either set of parameters, the other set can then by obtained. Thus, the approximation technique has a number of different applications in feedback control design.

Two applications will be presented. In the first, a desired PDF for the closed-loop process is selected and the closed-loop dynamics are calculated using the approximation. The closed-loop dynamics may then be inverted to give the feedback control law. This approach to feedback control design is referred to as PDF-shaping. The second application is in optimal control design with nonquadratic, nonsymmetric objectives. Here, the parameters in the closed-loop dynamics may be explicitly included in the objective function as well as influencing the objective through their relationship with the probability density function. Thus in optimizing the objective it is necessary to use the relationship between the parameters in the closed-loop dynamics and the parameters in the process PDF. This sub-optimal but effective control design technique will be illustrated with a small case study.

In addition to the technique derivations and simulation examples, comments on the stability of the resulting controlled processes and accuracy of the approximation technique are given.

ROBERT HAYES, University of Alberta
Modelling and optimisation of the automotive catalytic converter

The automotive catalytic converter has become ubiquitous in Europe and North America as a pollution abatement device for gasoline and diesel engines. The driving force behind this reactor development is environmental legislation, which is becoming increasingly stringent. Much of the early work on the development of the catalytic converter was based on the design of catalysts and reactor systems capable of meeting the demanding requirements placed upon them. For gasoline engines, it is necessary to oxidize hydrocarbons and carbon monoxide and to reduce nitric oxide. To achieve both the oxidation and reduction requirements it is necessary to operate the engine at the exact stoichiomnetric point, and thus the oxygen required for the oxidation reaction is supplied by the NOX, which is correspondingly reduced. In practice, the engine oscillates about the stoichiomentric point, alternating between rich and lean operation. Excess oxygen in the exhaust gas in the lean cycle is stored in the catalyst and released during the rich cycle, thus preserving a time averaged stoichiomentric balance. This type of catalyst is called the three way catalyst (TWC), and forms the basis of most gasoline engine converter catalysts in use today. For diesel engines, which operate in the lean combustion environment, the emphasis is on oxidation only, although NOX reduction is desired, and remains an area of active research.

Reactor technology in catalytic converters is almost entirely based on the honeycomb monolith reactor, although some early designs used packed beds. Monolith reactors provide a high surface area to volume ratio and give a small pressure drop, which is essential in the operation of the engine. The catalyst is typically contained in a thin washcoat layer deposited on the walls of the monolith channels.

The increasing stringincy of government regulation, combined with a desire on the part of catalyst manufacturers to reduce devolment costs and to shorten design cycles, has led to an increasing interest in the use of computer aided design to assist in catalytic converter development and understanding. The simulation of the catalytic converter presents a number of modelling challenges. Methods of responding to some of these challenges is the topic of this presentation.

It is a prerequisite of reactor modelling to have a reasonably representative set of kinetic equations. For the TWC, there are essentially five or six overall or global reactions, which are comprised of a plethora of elementary mechanistic reaction steps. The current commonly used reaction scheme for modelling reactors is based on a set of LHHW inspired kinetic models, with adaptations for such features as catalyst oxygen storage. Such a set of equations contains of the order twenty to thirty adjustable parametrs, or more, depending on the level of complexity of the model. These parameters are empirical, and must be determined from experimental data. The determination of these parameters represents a challenging optimisation problem, which classical gradient based methods can have difficulty in dealing with. Current practice is to fix most of the variables and adjust the remaining few manually. Pattern search algorithms offer the advantage of a gradient free method for determining optimal parameter values. These methods are useful, for example, when the derivatives are not available, or expensive to evaluate with sufficient precision. This is often the case when the objective function is given by a "black box", such as, for example, the output of a simulation, a PDE solver, or the results of an experiment. Although they potentially be expensive in terms of CPU time, they also offer the advantage of the application of parallel computing, which can significantly reduce total wall clock turn around time. The application of a global pattern search algorithm will be demonstrated for the solution of several converter problems.

A key operating variable in the operation of a catalytic converter is the resistance to diffusion offered by the washcoat. Because of the complex reaction mechanisms and kinetic inhibition terms, the role of diffusion is not inherently obvious, especially on the light-off performance under cold start conditions. Furthermore, the development of effectiveness factor relationships for the rate expressions is challenging. The use of effectiveness factors to model the reactor performance is desirable to reduce the execution time for full converter reactor models, which otherwise would require massive amounts of computer power to run. The development of some effectiveness factor models for the LHHW type kinetics will be presented. This development uses the GPS optimization algorithm used earlier.

More complex mechanistic kinetic models are the next level of development in the modelling of catalytic converters. Such models contain a large number of degrees of freedom, which require special computational methods to deal with. Furthermore, many of the elementary reaction steps lead to stiff systems of equations, which require a robust solution method to solve. Some results obtained for detailed mechanistic models will be presented, and the computation algorithms used will be described.

KEVIN L. HEPPNER, Department of Chemical Engineering, University of Saskatchewan, Saskatoon, Saskatchewan  S7N 5A9
Computation of mass transport in electrolytic systems

Mass transport in convection-free electrolytic systems occurs via diffusion and electromigration. Primary electromigration is driven by an encompassing electrical potential gradient often produced by separation of half-cell electrodes. However, variations in the mobility and charge of dissolved ionic species induce the formation of a net solution charge density that drives diffusion potential electromigration. This paper mathematically quantifies the effect of diffusion potential electromigration on mass transport in electrolytic systems.

BIAO HUANG, Department of Chemical and Materials Engineering, University of Alberta
Control of chaos in a convective loop system

Dynamics of Newtonian fluids that involve chemical reaction and heat transfer offer a highly nonlinear model of flow systems with no uncertainties in the model. The distributed and highly nonlinear nature of the model however offers challenges even in predicting the behavior of such systems. Identification and developing control strategies of such mass and heat transportation is a bigger challenge that is attracting lot of attention currently. The system becomes more complicated in terms of prediction and regulating it to a fixed operating condition when it shows chaotic behavior. Such type of chaotic system can have multiple equilibrium points and may jump suddenly from one state to another without any external influence and may cause the system to show unstable behavior even leading to the shutting down of plants. A convective loop is one such system in which a fluid circulates freely inside a closed circular pipe. The circulating fluid works as a transport media of heat from a source to a sink. First order lumped parameter modeling of this system leads to a set of nonlinear ordinary differential equations (Lorenz type). Depending on Rayleigh number this system can show chaotic behavior. Using nonlinear model predictive control algorithm, control action (perturbation of heating rate from its nominal value) can be calculated to stabilize the system to its equilibrium point. The main target is to study the performance of different nonlinear and linear control. The idea is to develop control strategies based on the simplified lumped parameter models and apply it to distributed models of the system to investigate its robustness for different uncertain conditions like plant-model mismatch.

SHIJIE LIU, University of Alberta, Edmonton, Alberta
Parametric estimation and error structure

In this paper, the basic regression approaches have been reviewed. The parametric estimation has been discussed in association with the error structure of the data series. Different regression models have been discussed. In particular, one common type of regression model: differential regression model has been discussed in detail and an example has been used to further the discussions. In the Physical Chemistry and Reaction Engineering literature, the integral methods and the differential methods have been "standardized" for solving the differential regression model problems. Discussions on these two "standard methods" are made with a real kinetic estimation problem. Finally, a viable approach to parametric estimation for the differential regression model problem is presented. The proposed method consists of a general numerical integrator and the utilization of the solver program in Excel. With the proposed method, the trial-and-error nature and the tedious calculations of the traditional methods can be eliminated. Above all, the parameters estimated with the proposed method directly reflect the quality of the data series and thus can be trusted over those obtained from the traditional methods.

XIAODONG LIU, Dalian Maritime University
New approaches to H_¥-control for T-S Fuzzy descriptor systems via LMI approach

In this paper, the problem of H_¥-Control based on the state feedback for T-S fuzzy descriptor systems has been studied. First a new sufficient condition in the form of a single negative definite matrix which guarantees the existence of the state feedback H_¥ control for the T-S fuzzy descriptor systems has been proposed. The conditions in this paper are not only simple but also consider the interactions among the subsystems. Secondly based on the LMIs, the controller designing methods for the stabilization and the H_¥-Controller for the T-S fuzzy descriptor systems have been given.

NEZAM MAHDAVI-AMIRI, Sharif University of Technology, Tehran, Iran and York University, Toronto, Ontario  M3J 1P3
ABS solution of full row rank linear inequality systems over real and integer spaces

ABS (Abaffy-Broyden-Spedicato) methods have been used extensively to solve linear systems of equations in the past two decades. We have recently used these methods to solve linear Diophantine systems of equations as well as systems of linear inequalities having full row rank. We explain how the ABS approach is used to compute the general solution of these systems over both the real and integer spaces.

EDGAR TAMAYO, Syncrude Canada Ltd.

ANDREW WILLIS, Syncrude Canada Ltd., Mildred Lake Plant, Fort McMurray, Alberta  T9H 3L1
Condition monitoring of flow instrumentation through inferential modeling

In many processes found in industry, input streams of reagents are combined and to form a product. Inference of the input reagent statistics based on those of the process output can thus be viewed as an inverse problem. The results of this inversion may then be used to check the operation of instrumentation on the input streams.

In this paper we take the modified Claus process as an example, in which H₂S and S0₂ are combined to form Sulphur and water. The inversion previously mentioned is approached from a Bayesian perspective and a numerical approach is discussed that allows for this. In the special case of Gaussian statistics the Bayesian solution is developed in closed form and shown to reduce to a simple least squares estimator on a weighted metric space (i.e. (L2,W)). By iterating the prior an empirical Bayes scheme can be developed to converge the solution. In the linear case this is shown to correspond to iterating the metric weight, W, until convergence.

Results on simulated and real data are given in which the algorithm generally converges within a few iterations. The results are shown to be equivalent to a `regularised' form of Partial Least Squares inversion.

JASON ZHANG, University of New Brunswick, New Brunswick
Radiation modeling of an air phase corrugated plate photocatalytic reactor

The radiative energy absorption profiles on catalyst surfaces in an air phase corrugated plate photocatalytic reactor are mathematically modeled and numerically simulated based on the emission characteristics of a typical UV fluorescent lamp and solar radiation, the refraction indexes of the reactor materials, the spectral absorption characteristics of the catalyst, geometric optics, and analytical geometry. The effects of dimensions of the TiO2-coated corrugated plates, multiple reflections, radiation sources, and refraction by the reaction medium were examined.

Recapture of the reflected photons was found to result in higher photon absorption efficiency as well as more uniform local area-specific rate of energy absorption, which were shown to be dependent strongly on the dimensions of the corrugated plates. Compared to water phase systems, the radiation fields in an air-phase corrugated plate reactors were predicted to be less uniform due to the lack of refraction by the reaction medium, especially when the light source is a UV-A lamp. However, the efficiency of photon capture is not significantly lower, and remains substantially higher than TiO2 films on flat-plate geometries.