


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 holdup profiles for threedimensional laminar
bubbly twophase 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 pressurecorrection 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 twolevel
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 subproblems 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 threedimensional vertical cylindrical pipe.
Comparison with other authors is included. Theoretically, there exists
a possibility that the initially onedimensional wave generates, via
the nonlinear interactiom, threedimensional 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 wallforce and
the laterallift 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 fixedbed 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 welldefined
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 nonlinear 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 steadystate and dynamic behaviors of a pilot plant fixedbed
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 nonlinear 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 nonlinear
algebraic equations and were solved with the GaussianSeidel 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 steadystate operation. It also reasonably captured the
dynamic behaviors of reactor startup 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 facetoface 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,
19.
 FENG DING, Department of Electrical and Computer Engineering, University
of Alberta, Edmonton, Alberta T6G 2V4
Parameter estimation for dualrate 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 dualrate systems in
which the output sampling period is an integer multiple of the input
updating period, we obtain frequencydomain 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 leastsquares (LS) algorithm,
instrumental variable least squares (IVLS) algorithm, and
forgettingfactor least squares (FFLS) algorithm, using directly the
finite dualrate inputoutput 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 waterlevel system.
 MICHAEL G. FORBES AND FRASER FORBES, University of Alberta
Regulating discretetime 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 nonuniform 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
linearGaussianquadratic 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 closedloop 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 discretetime
stochastic processes is approached with a focus on the PDF of the
closedloop process. The discretetime process model is given by
equation (1),
x(t+1)=f 
æ è

x(t),u(t) 
ö ø

+w(t) 
 (1) 
where x(·) Î Â^{n} 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 (· , ·) Î Â^{n} is a smooth mapping; and t is
an integervalued time index. Note that equation (1) can
represent inputoutput 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, highorder
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
closedloop 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 closedloop dynamics.
In this work, equation (3), that gives the stationary PDF,
p(·), for the closedloop process (2), is examined
in detail.
p(z)= 
ó õ

Â^{n}

p_{w} 
æ è

z 
~
f

(x) 
ö ø

p(x) dx 
 (3) 
While this integral equation gives an exact implicit relationship
between the closedloop dynamics, [(f)\tilde] (·), and the
PDF of the closedloop process, in general no explicit solution exists.
An approximate solution to this equation is developed by
parameterization of both the process PDF and the closedloop 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 closedloop process is selected and
the closedloop dynamics are calculated using the approximation. The
closedloop dynamics may then be inverted to give the feedback control
law. This approach to feedback control design is referred to as
PDFshaping. The second application is in optimal control design with
nonquadratic, nonsymmetric objectives. Here, the parameters in the
closedloop 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 closedloop dynamics and the parameters in the
process PDF. This suboptimal 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 lightoff
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 convectionfree electrolytic systems occurs via
diffusion and electromigration. Primary electromigration is driven by
an encompassing electrical potential gradient often produced by
separation of halfcell 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
plantmodel 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 trialanderror
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 TS Fuzzy
descriptor systems via LMI approach

In this paper, the problem of H_{¥}Control based on the state
feedback for TS 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
TS 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 TS fuzzy descriptor
systems have been given.
 NEZAM MAHDAVIAMIRI, 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 (AbaffyBroydenSpedicato) 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_{2}S and S0_{2} 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 TiO2coated 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 areaspecific 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 airphase 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 UVA lamp.
However, the efficiency of photon capture is not significantly lower,
and remains substantially higher than TiO2 films on flatplate
geometries.

