Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Reducing variabilities in the mixer inflow concentration Soltanzadeh, Ali 2012

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2012_spring_soltanzadeh_ali.pdf [ 2.2MB ]
Metadata
JSON: 24-1.0072655.json
JSON-LD: 24-1.0072655-ld.json
RDF/XML (Pretty): 24-1.0072655-rdf.xml
RDF/JSON: 24-1.0072655-rdf.json
Turtle: 24-1.0072655-turtle.txt
N-Triples: 24-1.0072655-rdf-ntriples.txt
Original Record: 24-1.0072655-source.json
Full Text
24-1.0072655-fulltext.txt
Citation
24-1.0072655.ris

Full Text

Reducing variabilities in the mixer inflow concentration  by Ali Soltanzadeh Master of Control Engineering, Sharif University of Technology, 2005  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES (Electrical and Computer Engineering)  The University Of British Columbia (Vancouver) March 2012 c Ali Soltanzadeh, 2012 ⃝  Abstract The thesis focuses on mixing fluids, and it has three main parts. The first part investigates mixing in agitated pulp chests. A new model for identifying the performance of mixing in the pulp chest is developed. This newly developed model simplifies the previous pulp chest model and successfully describes non-ideal flow behavior, including channeling and dead-zones which occurs in the agitated pulp chests. The model is verified through experimental data obtained in a laboratoryscale agitated pulp chest. In the second part of the thesis, a novel method for identifying residence time (the average time that it takes for material to exit the mixer) from input and output data is developed. The main benefit of this new method is that it does not require that an individual have prior knowledge of the process to evaluate the mixing inside the flow system. The same idea is used to estimate the higher moments of linear time invariant transfer functions from a batch of input and output data. The ability to estimate higher moments of a transfer function enables one to reconstruct the original transfer function without having knowledge of its structure. This opens new possibilities in non-parametric system identification and residence-time theory. It is also shown how a bound on the delay of all pole transfer functions can be found using the first and second moments. In the third part of the thesis, a new class of mixers is designed. These mixers do not have agitators, yet they are as effective as stirred tanks in reducing fluctuations in their inlet stream. However, they are only effective for reducing the concentration fluctuations in the flow direction (axial mixing). This design splits the fluid into different channels and recombines these streams to achieve mixing. These mixers are especially useful in microfluidic applications where there is a physical ii  limitation in using agitators for mixing; therefore, in the final section of the thesis, a mixer is designed and its geometry is calculated for microfluidic purposes. Moreover, the effectiveness of the mixer in reducing inlet concentration fluctuations is shown through computational fluid dynamic simulations (CFD).  iii  Preface A version of Chapter 2 has been published. Soltanzadeh A., Bhattacharya S., Dumont G.A., and Bennington C.P.J., (2009) Estimation of residence time and channeling in agitated pulp chests, Nordic Pulp and Paper Journal 24(1), p 66-72. The idea for improving the model emerged from a discussion between myself, Professor Dumont and Professor Bennington. I validated the model by using experimental data from the laboratory. The paper was written in part by myself and in part by Dr. Bhattacharya. A version of Chapter 3 has been accepted for publication. Soltanzadeh A., Bennington C.P.J. (posthumous), and Dumont G.A., (Accepted 07/Oct/2010) Direct estimation of residence time from input and output data, Accepted in Canadian Journal of Chemical Engineering, 10 pages. The main contributions of this paper are the two theorems which I was able to prove. Professor Dumont kindly reviewed the proofs and ensured that they were accurate. Professor Bennington assisted us in connecting the research to industrial applications and also provided the funding for the project. The paper was written by myself and then reviewed by my supervisors. A version of Chapter 5 has been submitted for publication. Soltanzadeh A., Dumont G.A., Stoeber B., and Chad P.J. Bennington (posthumous), (2011) Modeling of time-averaging micromixers using parallel channels and a back flow loop. The idea of using parallel channels to accomplish mixing was proposed by Professor Dumont in 1995; however, the ideas of re-circulating part of the flow and emulating the performance of the stirred tank were my contributions to the work. Dr. Stoeber not only helped adapt the design for the purposes of microfluidics, but also proposed the potential application of the mixer in microscale. The first draft iv  of the paper was written by myself. This draft was revised and edited by Dr. Stoeber and myself many times. Professor Dumont gave us useful feedback during the process.  v  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ix  Nomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xii  Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xvi  Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  1.1  Mixing evaluation of a pulp chest . . . . . . . . . . . . . . . . . .  2  1.2  Estimation of residence time without a model . . . . . . . . . . .  8  1.3  Estimation of higher moments of an arbitrary transfer function . .  10  1.4  Design of time averaging micromixer . . . . . . . . . . . . . . .  12  1.5  Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  16  1.6  Summary of contributions . . . . . . . . . . . . . . . . . . . . .  16  2 Estimation of residence time and channeling in agitated pulp chests  18  2.1  Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  18  2.2  Experimental . . . . . . . . . . . . . . . . . . . . . . . . . . . .  19  2.3  Simplified model . . . . . . . . . . . . . . . . . . . . . . . . . .  20  2.4  Results and discussions . . . . . . . . . . . . . . . . . . . . . . .  24  vi  2.5  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  30  3 Direct estimation of residence time from input and output data . . .  32  3.1  Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  32  3.2  Proposed method . . . . . . . . . . . . . . . . . . . . . . . . . .  32  3.3  Results and discussion . . . . . . . . . . . . . . . . . . . . . . .  40  3.4  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  4 Estimation of higher moments of a system from its input and output data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  44  4.1  Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  44  4.2  Proposed method . . . . . . . . . . . . . . . . . . . . . . . . . .  45  4.3  Estimation of an upper and lower limit for the time delay . . . . .  50  4.4  Results and discussion . . . . . . . . . . . . . . . . . . . . . . .  55  4.5  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  5 Modeling of time-averaging micromixers using parallel channels and a back flow loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  5.1  Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  59  5.2  Design of a parallel channel mixer . . . . . . . . . . . . . . . . .  60  5.2.1  Mathematical model for mixing using parallel channels . .  60  5.2.2  Validation of the parallel channel model using CFD . . . .  63  Design of a mixer with recirculation . . . . . . . . . . . . . . . .  66  5.3.1  Mathematical model for mixing with recirculation . . . .  66  5.3.2  Optimization of the mixer transfer function with recirculation 69  5.3.3  Derivation of the mixer dimensions . . . . . . . . . . . .  71  5.3.4  CFD simulation of the mixer with recirculation . . . . . .  77  5.3.5  Validation of CFD results  . . . . . . . . . . . . . . . . .  80  5.4  Results and discussion . . . . . . . . . . . . . . . . . . . . . . .  81  5.5  Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  82  6 Final conclusions and future works . . . . . . . . . . . . . . . . . . .  87  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  92  5.3  vii  A Proof of inequality (4.40) . . . . . . . . . . . . . . . . . . . . . . . .  viii  97  List of Figures Figure 1.1  The frequency response measured for an industrial chest with non-ideal flow compared to the response of an ideally mixed chest ([20]). . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  Figure 1.2  Non-ideal behavior of flow in the chest. . . . . . . . . . . . .  5  Figure 1.3  Continuous-time dynamic model of the stock chest with separate delays. . . . . . . . . . . . . . . . . . . . . . . . . . . .  Figure 1.4  Two different methods of introducing two streams into the mixer, (a) simultaneously (b) sequential. . . . . . . . . . . . . . . .  Figure 2.1  19  Stock chest configurations showing pulp input and exit locations investigated in the laboratory chest [26]. . . . . . . . . .  Figure 2.3  14  Schematic diagram of the scale-model apparatus used for the dynamic tests [26]. . . . . . . . . . . . . . . . . . . . . . . .  Figure 2.2  7  20  Example of scale-model chest input and output signals. Data shown for Configuration 2 with the impeller speed of 1105 rpm, pulp flow of 7.9 L/min and pulp consistency of 3.3%. . .  Figure 2.4  21  Simplified model with one time delay used for identifying the dynamics of pulp stock chests. Mixing in the bypass zone was found to be negligible. τ1 therefore approximates to zero and the transfer function, G′1 (s) to unity. . . . . . . . . . . . . . .  Figure 2.5  23  Comparison of the single- and two-delay models. f is compared at two pulp flow-rates over a range of impeller speeds. The tests were carried out in Configuration 1 with Cm = 3.3%.  ix  24  Figure 2.6  The time constant in the mixing zone is estimated using both models. Note that at N < 1300rpm and Q = 37.1 LPM the one-delay model properly identifies poor mixing in the chest because most of the pulp is short-circuiting. . . . . . . . . . .  Figure 2.7  25  Volume fraction of the chest that is fully mixed estimated for both models under the operating conditions defined above. The single-delay model is able to capture all the dynamic characteristics revealed by the two-delay model. . . . . . . . . . . .  26  Figure 2.8  Low impeller speed can lead to full channeling in the chest. .  28  Figure 2.9  The effect of pulp consistency and chest configuration on parameter estimation using the two models. Bypass fraction, f , is estimated at pulp flow-rate of 37.1 LPM for both chest configurations at two pulp consistencies: (a) Cm = 3.3% and (b) Cm = 2.7%. . . . . . . . . . . . . . . . . . . . . . . . . . . .  29  Figure 2.10 The effect of pulp consistency and chest configuration on the estimation of the fully mixed volume using the one-delay and two-delay models for the operating conditions defined in Figure 2.9. Results for (a) Cm = 3.3% and (b) Cm = 2.7%. . . . . Figure 3.1  Input and output signal (top). Right-hand side of equation (3.1) (bottom)  Figure 3.2  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  35  Input signal (top). Right-hand side of equation (3.1) for first order model and model in example 1. . . . . . . . . . . . . .  Figure 3.3  30  36  Input and output data (top). Right-hand side of equation (3.1) calculated from input and output data and also estimated first order model (bottom). . . . . . . . . . . . . . . . . . . . . .  Figure 3.4  41  Input and output (top). Right-hand side of equation (3.1) calculated from input and output and also estimated first order model (bottom). . . . . . . . . . . . . . . . . . . . . . . . . .  Figure 4.1  Figure 5.1  43  Actual and reconstructed impulse response of the system in Example 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . .  57  Schematic of a mixer based on parallel channels . . . . . . . .  60  x  Figure 5.2  The frequency response of a notch filter designed in the example. 64  Figure 5.3  The step response of a notch filter designed in Section 5.2.1. .  Figure 5.4  The frequency response of a notch filter designed assuming  66  plug flow compared to the revised model with a parabolic profile and diffusion , and estimated from CFD simulations. . . .  67  Figure 5.5  Schematic of a mixer based on parallel channels and recirculation 68  Figure 5.6  The frequency response of a first-order filter representing a typical continuous stirred tank, G f , and the optimized plug flow mixer with recirculation, Gr . . . . . . . . . . . . . . . .  Figure 5.7  72  The frequency response of a first-order filter representing a typical continuous stirred tank, G f , and the plug flow mixer with perturbed parameters, Gr . . . . . . . . . . . . . . . . . .  73  Figure 5.8  The schematic of the open loop mixer in Comsol. . . . . . . .  77  Figure 5.9  Concentration profile across the mixer at t=8 s. . . . . . . . .  79  Figure 5.10 Concentration profile across the mixer at t=18 s. . . . . . . . .  80  Figure 5.11 Outlet concentration profile from plug flow, CFD simulation and its corresponding first order system with the same residence time. . . . . . . . . . . . . . . . . . . . . . . . . . . .  81  Figure 5.12 Snapshot of concentration profile at (a) inlet and (b) outlet of mixer for ω = 15 rad/s and 400 kPa input pressure. . . . . . .  83  Figure 5.13 Ratio of outlet fluorescence to inlet fluorescence in different inlet frequency response for CFD and experimental case.  . .  84  Figure 5.14 The Simulink model for comparing a stirred tank model and plug flow mixer model. . . . . . . . . . . . . . . . . . . . . .  85  Figure 5.15 The time response of a stirred tank and plug flow mixer, (a) input and (b) output concentration. . . . . . . . . . . . . . . .  xi  86  Nomenclature This section explains variables used in equations throughout the thesis.  Chapter 1 Gchest  Chest transfer function between inlet and outlet concentration  f  Bypass fraction  T  Time delay (s)  τ R τdyn t g(t) y(t) u(t) G G′ s µgn νgn Bi  Time constant (s) Recirculation fraction Residence time (s) Time variable (s) Impulse response Output at time t Input at time t Laplace transform of g(t) Derivative of Laplace transform of g(t) Laplace variable The nth moment of the transfer function G(s) The nth normalized moment of the transfer function G(s) ith Laguerre coefficient  xii  Chapter 2 G′  Chest transfer function between inlet and outlet concentration  f  Bypass fraction  τ1 τ2 z T s V f ullymixed VTotal Q LPM  Channeling time constant (s) Mixing time constant (s) Chest transfer function zero Time delay (s) Laplace variable Mixed chest volume (m3 ) Nominal chest volume (m3 ) Flow rate (m3 /s) Liter per minute  Chapter 3 τdyn t g(t) G G′ y(t) u(t) Y Y′ U U′  Residence time (s) Time variable (s) System impulse response Laplace transform of g(t) Derivative of Laplace transform of g(t) Output at time t Input at time t Laplace transform of output data y(t) Derivative of Laplace transform of y(t) Laplace transform of input data u(t) Derivative of Laplace transform of u(t)  Chapter 4 t  Time variable (s)  g(t)  System impulse response  xiii  µgn y(t) u(t) a G H T s G′ G′′ H′ H ′′ µhn α p  The nth moment of the transfer function G(s) Output at time t Input at time t Pole of the first-order model Laplace transform of g(t) Laplace transform of h(t), system without a delay Time delay (s) Laplace variable Derivative of Laplace transform of g(t) Second derivative of Laplace transform of g(t) Derivative of Laplace transform of h(t) Second derivative of Laplace transform of h(t) The nth moment of the transfer function H(s) Pole of the first-order model with delay Pole of the all pole system  Chapter 5 Ti  Time delay in the ith channel (s)  Wi  Width of the ith channel (m)  Li  Length of the ith channel (m)  Qi  Flow rate in the ith channel (m3 /s)  hi  Height of the ith channel (m)  cOi  Output concentration of the ith channel  cI  Mixer inlet concentration  ai  Ratio of the flow rate in the ith channel to the total flow  Q  Total flow rate (m3 /s)  ω G Tr r ωl  Frequency variable (Rad/s) Inlet to outlet transfer function Time delay in the recirculation channel Recirculation ratio Desired mixer cut-off frequency (Rad/s)  xiv  ωh Ri Rtotal µ Vi ∆P Wt D ρ u  The upper range of mixer attenuation (Rad/s) Hydrodynamic resistance of the ith channel (Pa.s/m3 ) Total hydrodynamic resistance of the mixer (Pa.s/m3 ) Viscosity (Pa.s) Horizontal flow velocity of the ith channel (m/s) Pressure drop across the mixer (Pa) Sum of the widths of all channels (m) Diffusion rate (m2 /s) Material density (kg/m3 ) Inlet horizontal flow velocity (m/s)  Chapter 6 D  Cavern diameter (m)  τy cm V  Yield stress (Pa) Pulp consistency Cavern volume (m3 )  xv  Acknowledgments I offer my sincere gratitude to my supervisor, Professor G. A. Dumont, for the freedom he gave me during the course of my PhD. It was an honor to work under his supervision. I am obliged to my departed supervisor, Professor C.P.J Bennington, for introducing me to mixing and providing me with funding for the project. His empty space will be missed forever. I owe my sincere gratitude to Doctor Stoeber for introducing me to microfluidics and all the useful comments he offered while I was writing the fifth chapter of the thesis. I also want to thank my colleagues and officemates, Doctor C. G´omez, Doctor S. Bhattacharya, and Patrick Davis for all the help they gave me. And special thanks to my parents for supporting me morally and financially throughout my life.  xvi  Dedication To the memory of Chad Bennington  xvii  Chapter 1  Introduction Mixing is a ubiquitous process. Despite the frequency of its occurrence in both nature and industry and our intuitive understanding of mixing phenomena, it is hard to scientifically define quantitative indices which can describe different aspects of a mixing process. The aim of this thesis is to investigate mixing from a very narrow but insightful angle. Most of the material in this thesis is based on residence time theory; ‘residence time’ refers to the average time that materials stay in the mixer, and ‘residence time theory’ considers inlet and outlet streams of a mixer in order to study the mixer dynamic. In residence time theory a mathematical model is associated with a mixer dynamic, and the model parameters are estimated from the input and output data. There are several benefits of using this approach. It is usually less complex than other methods of studying a mixer, and it is easy to implement in industry; moreover, accurate measurements of data needed for the analysis are easily obtainable in a typical industrial mixer. On the other hand, results of the analysis are not as detailed as more sophisticated computational or experimental methods e.g. Computational Fluid Dynamics (CFD) or Electrical Resistance Tomography (ERT); however, all or some of these aforementioned methods can be applied simultaneously to the same mixer in order to create a coherent and complete picture of the mixing process. In Chapter 2 the residence time theory has been used to quantify mixing in agitated pulp chests. This chapter presents the efforts of the UBC Dynamic Mixing Group 1  to study these chests. More details about Chapter 2 are given in Section 1.1. Chapter 3, Chapter 4 and Chapter 5 are the main contributions of this thesis. Chapter 3 significantly improves upon the method for identifying residence time originally introduced in Astr¨om and Hagglund [6]. Chapter 4 extends the result of Chapter 3 to higher moments of a transfer function, and Chapter 5 is a continuation on attempts initiated by Nauman et al. [33] for designing static mixers to promote axial mixing (mixing in the flow direction). More detailed descriptions of these chapters are given in Section 1.2, Section 1.3 and Section 1.4.  1.1  Mixing evaluation of a pulp chest  Agitated chests are stirred vessels that are widely used in the pulp and paper industry. Blending different chemicals, providing temporary storage capacity for the pulp stream, removing latency from pulp suspension, and reducing pulp stream variabilities are among the main applications of agitated pulp chests. In pulping operations the primary objective of these units is to suppress variations in suspension properties, such as changes in pulp consistency, Cm , at frequencies higher than the cut-off frequency of the process control loops. The disturbance attenuating abilities of stock chests depends primarily on their mixing efficiency. The machine chest is especially designed for this purpose, and it prepares a well-mixed pulp stock for the paper machine. The paper machine controller has a limited bandwidth for removing high frequency variations in an incoming pulp stream, and failing to remove these fluctuations can diminish paper quality; therefore, it is valuable to develop methods for evaluating pulp chests’ ability for removing unwanted consistency fluctuations. This is illustrated in Fig. 1.1 where a comparison of the frequency responses of an ideal and an industrial chest of equal volume and flow rate show that the upset attenuation of the real chest is significantly worse1 . It is therefore desired to develop methods to properly characterize the dynamics of mixing in pulp chests which can then be used to improve the design of new chests and optimize existing ones. It is important to be aware that there are many devices other than agitated pulp 1 Ideally  mixed chest is the chest that its full volume is agitated by the impeller  2  chests that contribute to mixing in the pulping process (e.g. paper machine headbox and consistency controllers); therefore, it is difficult to isolate the effect of the pulp chest on the quality of the final paper product. Hence, the main goal in Chapter 2 is to compare the performance of the chest during operation with its ideal performance. At present, pulp stock chests are characterized or designed on the basis of previous empirical experience and proprietary information. This is illustrated by the technique summarized by Yackel [54] wherein the impeller speed is estimated based on the momentum flux required to achieve smooth motion over the complete surface of the pulp chest. For process system design, performance is assessed by simply assuming that the chests are ideally mixed [13, 41, 47]. However, recent system design studies by our group [21] show that there can be significant stagnant regions in the chest even when the surface is in complete motion and that deviations from ideal mixing are significant and common in agitated pulp chests. These non-ideal flow behaviors can be categorized as plug flow, channeling, re-circulation and dead zone (Fig. 1.2). These anomalies are in large part due to the complicated rheology of the suspensions. Pulp suspensions form continuous fiber networks with the strong fiber-fiber interactions making them non-Newtonian as evident by their significant yield stress (a minimum amount of shear to make pulp flowing). Moreover, as the suspension mass concentration changes, the interaction between the fibres changes, increasing or decreasing the yield stress [8] and affecting mixing performance. This complex, composition-dependent suspension rheology can cause significant deviations in the flow-field of an ideally mixed vessel. Especially, the presence of yield stress would lead to an agitated volume near the impeller known as a ‘cavern’, and it creates slow moving or even stagnant zones in other parts of the mixer. Stagnant zones, often referred to as ‘dead zones’, are typically found in the chest’s corners. Also, it is possible for the inlet stream to short-circuit the cavern and reach the outlet without passing through the mixing zone; this stream would be referred to as a ‘bypass’ stream. These non-ideal flow behaviors would shrink the mixer-active volume and limit the ability of the chest to remove temporal variations. Different approaches can be taken to study the effect of the aforementioned phenomena on the performance of a chest. Ein-Mozaffari [17] setup an experiment to 3  Figure 1.1: The frequency response measured for an industrial chest with non-ideal flow compared to the response of an ideally mixed chest ([20]). validate the presence of non-ideal behavior in a chest; in this experiment, colored dyes were mixed with pulp suspensions and a video recording was used to visualize the flow in the lab-scale pulp chest. In this recording, the fast-moving region near the impeller and slow-moving regions in other parts of the chest could be recognized; however, this visualization method does not provide detailed information about the flow field inside the chest. To overcome this problem and create a better understanding of the flow field in different parts of the chest, the computational fluid dynamic (CFD) was deployed in G´omez [22]. The complex rheology of pulp suspension makes it very difficult to obtain a reliable flow field from this approach; however, G´omez [22] had some success in describing the flow field in different flow conditions. In other studies Ein-Mozaffari et al. [19–21, 26] illustrated the effect of non-ideal flow for both laboratory- and industrial-scale mixing chests using black-box sys-  4  Figure 1.2: Non-ideal behavior of flow in the chest. tem identification. These studies also showed that not only do portions of the stock chests remain stagnant with very little mixing, but that for higher feed flow rates, slower impeller speeds or higher pulp mass concentrations, a significant increase in the amount of channeling and short-circuiting can occur, adversely affecting mixing performance. In these studies, Ein-Mozaffari et al. developed a parametric model based on a physical interpretation of the flow-field and used this to characterize the nonidealities present. Due to its relevance to the present work, the model is described in some detail with additional information available in the literature [26]. When fluids with a substantial yield stress are agitated, the region immediately surrounding the impeller is well-mixed with high velocity gradients and shear stresses [4, 43, 50] between the active ‘cavern’ and the surrounding fluid. These outer regions have less developed flow-fields and subsequently achieve poorer mixing. The dynamic model ([20]), shown in Fig. 1.3, uses this compartmentalized description of the flow-field. The model visualizes a fraction of the incoming feed-stream, f , bypassing the impeller region and moving directly to the exit while the remaining fraction (1 − f ) reaches the cavern around the impeller. Mixing occurs primarily  5  in the impeller zone, but some secondary mixing could also occur in the bypass stream. Analysis of ideally mixed processes yield first-order differential equations, hence mixing in both the impeller and bypass regions were individually represented by first-order transfer functions: G1 (s) and G2 (s). τ1 and τ2 are the time-constants for the mixing processes in the short-circuiting and the mixing zone, respectively (τ1 ≪ τ2 ). The time delays, T1 and T2 , account for the time taken to transport pulp through these zones, with T1 ≤ T2 . The model also allows for recirculation within the mixing zone, with the recirculating fraction denoted by R. Using Fig. 1.3, the combined transfer function for the chest in continuous-time domain was obtained as:  e−T2 s  Gchest  f e−T1 s (1 − f )(1 − R) 1+τ2 s . = + −T s 1 + τ1 s 1 − Re 2  (1.1)  1+τ2 s  The model reduces to a first-order transfer function, representing an ideally mixed chest, when f = 0 and R = 0 with the time-constant of the mixing zone (τ2 ) apQ . proaching the nominal time constant of the chest, τnom = VTotal  The parameters in equation (1.1) were estimated for a given chest by introducing a set pattern of variation in a suspension property at the inlet of the chest, measuring the chest response to the disturbance at the chest exit, and using the numerical method developed by Kammer et al. [26] to identify the dynamic model. Detailed description of the identification procedure is given in Kammer et al. [26]. While the results from this model allowed identification of the mixing, the form of equation (1.1) introduces a number of difficulties in the identification process: • The presence of the time delay in the denominator in the recirculation loop causes complications while transforming the model from the continuous to its discrete time form. • The presence of two time delays makes the model too complex to consider the use of identification techniques developed for generic systems with delay [39]. Moreover, searching in a two-dimension space for the time delays drastically increases computational demands. • The nonlinear relation between parameters and various constraints on the parameters makes the optimization process very complex. 6  • The procedure does not guarantee the attainment of a global minimum, although a Monte-Carlo simulation showed encouraging results [26]. It was observed that the solution is relatively sensitive to the initial conditions and these initial values had to be re-tuned in some cases. • Finally, being complex and time-consuming the method is only suitable for off-line system identification procedures and cannot be used for real-time identification purposes.  Figure 1.3: Continuous-time dynamic model of the stock chest with separate delays. In Chapter 2 we have revisited the dynamic mixing model in equation (1.1) and have simplified it. The resulting model allows for the use of generic identification methods while removing some of the difficulties associated with the previous model. The results of the simulation from both models are compared to validate the ability of this simplified model to capture the main characteristics of the mixing dynamics identified by the previous model.  7  1.2  Estimation of residence time without a model  The method in Chapter 2 uses a black-box model to identify mixing parameters (e.g. bypass fraction and residence time). As demonstrated in this chapter, the method works well for a chest with simple geometry, but it is very difficult to generalize the method for a flow system with arbitrary geometry; for example, a multiple-impeller industrial chest has very complicated flow dynamics, and it is not easy to model its flow behavior. Therefore, developing a method that describes the chest characteristics and does not rely on the dynamical model is very beneficial. In Chapter 3, such a method is developed and it is shown that it can be used to successfully estimate the residence time of a lab-scaled pulp chest. Residence time is a well known and widely used way to characterize flow systems in various chemical processes [32]. It can describe the degree of fluctuation reduction in mixers or the effectiveness of reaction in reactors. It is, also, an important parameter for designing a controller for a given process [3, 5]. Usually, to estimate the residence time a first-order plus delay model is assumed for the system and then the delay and first-order pole are estimated. In this model, the delay represents the plug flow, and the first-order pole describes mixing in the agitation zone [28]. In order to identify this model most methods use step input. This fits quite well with the reality of industry where step input is the most practical input for exciting a system. This is usually conducted by opening or closing a valve in the input stream of the process. In one method for identifying a first-order plus dead time model from a step input, the graphical characteristics of the step response are used [6]. By drawing a tangent to the step response, the delay is found and the first pole is identified from the settling time. This method is very simple, but it is sensitive to measurement noise. Moreover, in practice it is impossible to have instantaneous changes in the input; there is always a slope corresponding to the limited response time of the valve. Although this ramping time is usually much smaller than the delay and the time constant of the system, it still introduces inaccuracy in the calculation of the parameters. Similar to the previous method, Astr¨om and Hagglund [6] introduced the method  8  of moments. In this method the residence time is calculated from ∫∞  t.g(t) dt τdyn = ∫0 ∞ 0 g(t) dt  (1.2)  where t represents time, τdyn is the residence time, and g(t) is the impulse response of a linear stable casual system. In Laplace domain the above equation can be written as  τdyn = −  G′ (0) G(0)  (1.3)  where G(s) is the Laplace transform of the system. Equation (1.2) has a physical interpretation as follows. If one inserts a significant amount of tracers in the inlet for a very short period then  ∫ ∞g(t) 0  g(t) dt  would be the  probability function of the material coming out of the system at time t; therefore, equation (1.2) implies that residence time is the average exit time of the tracers from the system [32]. To calculate the residence time from equation (1.2) it is necessary to calculate the impulse response of the system; however, by proper rearrangement of equation (1.3) it is possible to calculate the residence time from any arbitrary set of data as ∫∞  ∫  ∞ t.y(t) dt 0 t.u(t) dt ∫ − ∞ 0 y(t) dt 0 u(t) dt  Residence time = ∫0 ∞  (1.4)  where u(t) and y(t) are input and output data respectively. This method is referred to as method of moments [6]. Although the method of moments does not need any prior knowledge about the system and potentially could be very useful, it gives a very poor performance in estimating the residence time. In Bjorklund and Ljung [11] it is suggested that multiplication of input and output data by time is the reason for this poor performance. They argue that the noise is generally constant for input and output data during the time, but the signal is decaying for a stable system; therefore, multiplying data by time deteriorates the signal-to-noise and results in the poor performance of the method. Another possible reason for this method’s inaccuracy which has not previously been mentioned, is the initial time reference. In the lower limit of integrals in equation (1.4) the system should be at rest, but this 9  condition is rarely satisfied at the start of the identification experiment. In Chapter 3 a new approach will be introduced that solves these issues for estimation of residence time. The proposed method can be considered an improvement of the moments method. It can not only be applied to any linear flow system, but the input signals for estimation can also be more general than step or impulse.  1.3  Estimation of higher moments of an arbitrary transfer function  With a slight modification, the method in Chapter 3 can be extended to estimate higher moments of a transfer function. The nth moment of the transfer function G(s) can be calculated as  µgn  ∫∞ n t .g(t) dt = 0∫ ∞ 0  g(t) dt n  = (−1) lim  s→0  d n G(s) dsn  (1.5)  G(s)  where µgn is the nth moment in this equation. Theoretically, it is possible to reconstruct a transfer function using its higher moments. One very basic approach for reconstructing a transfer function from its moments is using the infinite series expansion of G( jω ) based on its higher moments [24, 32]. However, to be able to reconstruct a transfer function from this approach it is necessary to estimate all of its moments which is obviously not feasible; on the other hand, reconstructing a transfer function from a finite number of moments has attracted mathematicians for more than a century because of its practical significance, and wide applications. A complete survey on the subject is given in John et al. [25]. However, in this thesis an old method but familiar for process engineer is introduced [16]. The method is suggested in Nauman and Buffham [32] and it is developed based on a family of orthogonal function known as Laguerre polynomials [1]. Any function can be approximated by weighted summation of these polynomials. There is a close relation between Laguerre weight coefficients and higher moments of a function, but before introducing this relation it is necessary 10  to define normalized moments. The nth normalized moment of a transfer function G(s) , νgn , is defined as  νgn  =  µgn (τdyn )n  .  (1.6)  The first five coefficients of Laguerre expansion of an arbitrary function g(t) can be obtained as B0 = νg0 = 1 B1 = νg0 − νg1 = 0 2νg0 − 4νg1 + νg2 1 νg2 =− + 4 2 4 6νg0 − 18νg1 + 9νg2 − νg3 1 νg2 νg3 B3 = =− + − 36 3 4 36 24νg0 − 96νg1 + 72νg2 − 16νg3 + νg4 νg4 1 νg2 νg3 B4 = =− + − + 576 8 8 36 576 B2 =  (1.7)  and by using these coefficients it is possible to reconstruct the unknown function g(t). It should be emphasized that the optimum reconstruction of a function from its moments is out of the scope of this thesis, and the only aim of this section is to demonstrate that it is possible to reconstruct a function from its moments. The thesis merely concentrates on estimating the higher moments of a transfer function from its input and output data. This new method for estimating higher moments of a transfer function is introduced in Chapter 4. This method can be combined with methods for reconstructing a transfer function from its moments. This combination introduces a new approach in linear system identification. In this approach instead of assuming a model for the system and estimating the model parameters, one would first estimate the higher moments of the system and then reconstruct the system from its moments. The main advantage of this approach is that it does not rely on a specific model structure. In this aspect the method is similar to non-parametric identification methods [30] which do not use a model to extract information about the system from input and output data. The simplicity of non-parametric methods make them appealing for many applications, but they cannot provide detailed information about the system; on the other hand, the improved method of moments can provide any necessary information without assuming a model for the system. 11  1.4  Design of time averaging micromixer  Although industrial pulp chests have been used for many years, performance evaluation of several industrial pulp chests has shown that many of these chests are not effective enough to remove concentration variabilities in the desired frequency range [17]. To improve the mixing performance, it is necessary to increase the amount of energy dissipated in the agitator; this is not easily realizable, and it is essential to develop new mixing concepts for the efficient removal of concentration fluctuations. The purpose of Chapter 5 is to introduce a new mixer design that is capable of removing concentration fluctuations effectively. The idea in this chapter was originally developed for pulp and paper applications [12]; however, the rapid prototyping capability of microelectromechanical systems (MEMS) was the main motivation for designing the mixer at a microscale. On the other hand, this was not the only reason that the proposed mixer was designed at this scale; most of the mixers described in the MEMS literature aim to reduce the fluctuations in the cross direction of the microchannel. The mixer proposed in this chapter introduces a new approach for reducing concentration fluctuations along the microchannel and emulates the performance of stirred tanks. This technology would become crucial if MEMS technology were to grow as an established industry with continuous processes in MEMS scale. The concept of the stirred tank would be very useful for providing storage area in MEMS processes and also for alleviating uncertain changes in the flow of crucial process chemicals. On the other hand, the concept for emulating the behavior of the stirred tank by splitting and recombining the inlet flow is not restricted to MEMS scale. All the derivations and equations in this chapter can be used on a larger scale with the exception of the equations regarding hydrodynamic resistance 2 . These equations are scale-dependent and therefore they should be revisited based on the mixer scale. In microfluidic systems two principal concepts exist to mix two or more fluids at 2 Hydrodynamic resistance of a channel is the constant ratio between pressure drop and flow rate of the channel.  12  desired ratios, in both cases a dosing step that introduces the fluids into a channel is followed by a homogenization step. The fluids can be introduced into a channel side by side at flow rates Qi that correspond to the desired ratios of the components (Fig. 1.4a), mixing then needs to occur across the channel (similar to static mixers). Alternatively, the different fluids can be introduced into the channel in sequence (Fig. 1.4b), and the duty cycle for the individual components then corresponds to the desired ratios [34, 35]; in this case, mixing needs to occur in axial direction along the channel [2]. Most micromixers [23, 37] are designed to achieve mixing across the channel. However, these mixers require an analog control of the flow rates of the fluids, especially if the composition of the mixed fluids needs to vary in time. The second concept requires a much simpler dosing system; here, digital valves can be used that are either fully open or fully closed. The requirements for this mixing configuration are similar to the requirements of traditional applications of stirred tanks. There is one more crucial difference between these two mixing configurations; in the side by side scheme mixing happens instantly. Therefore, changing inlet flow of any of the two streams can be observed almost instantly in the outlet; however, in the sequential configuration mixing takes place after several period of changes between two inlets and therefore it takes time to observe the result of changes in the outlet. Furthermore, these changes are much smoother than changes that happen in the input. This is one of the main benefit of using dynamic mixing concept rather than static mixing concept. Stirred tanks are the oldest and the most widely used dynamic mixers in industry. The idea of stirring to achieve mixing is very intuitive, as it is encountered in many daily life situations. Stirring works fairly well, especially when the material is a low viscosity fluid as there is not much energy consumption associated with stirring to achieve the desirable degree of mixing; on the other hand, mixing through stirring is very challenging for high-viscosity fluids [36]. The typical low Reynolds number associated with microfluidics makes it therefore very difficult to use agitators in micromixers, and it has made it necessary to develop new mixing concepts [46, 51]. Developing new micromixers based on large-scale mixers has been successfully practiced before [9]. The large-scale mixer designed by Nauman et al. [33] suits the purpose of emulating a stirred tank and this design has the potential to be used as a microfluidic mixer. The idea behind this mixer is to 13  Figure 1.4: Two different methods of introducing two streams into the mixer, (a) simultaneously (b) sequential. distribute the fluid in parallel pipes and then re-connect the streams in such a way to assimilate the behavior of a stirred tank. Each pipe in this design emulates part of the stirred tank wash-out function [31], and for achieving that all the pipes must have the same flow rate, but different lengths for introducing different lag times to the stream. These two constraints require that longer pipes should be comprised of several parallel thinner pipes. Nauman et al. [33] implemented these constraints  14  and proposed a design in which the longest pipe is comprised of more than one hundred parallel thinner pipes which makes this implementation unsuitable for microfluidics purposes. The attempt by Nauman et al. [33] is not the first research which suggests using parallel pipes for mixing 3 , but it is the first one that applies appropriate mathematical modeling to design a complete mixer. More recently, Xie et al. [53] used parallel microchannels to separate chemical signals through their frequency domain spectral characteristics. The authors give experimental evidence for how their design is able to act as a notch filter for fluid concentration fluctuations. However, the design methodology is restricted to notch filters and is not applicable to mixers. In the current work the concept of using parallel channels is considered in a more general framework where dimensions of the different parallel flow channels can be determined in order to achieve a desired frequency behavior. In another mixer design by Xie et al. [52] two streams are introduced into the channel in sequence. In this setup mixing is achieved by fast opening and closing the valves for each stream. The main benefit associated with the proposed design in this thesis over the aforementioned design is in the design in [52] if any of the inlet streams is interrupted the mixing performance will drop immediately; however, in the current design the mixer behaves similar to a stirred tank, and therefore the performance of the mixer is not affected by instantaneous changes in the inlet streams’ flow rates. Chapter 5 introduces a design methodology for a mixer that first separates a flow into several parallel channels and then recombines the flow from these channels. The objective is to attenuate high frequency components in concentration fluctuations along the stream that may be introduced by a dosing system as in Fig. 1.4b. There are several similarities and differences between the design in this chapter and the one in [33]. The goal for both designs is to emulate the behavior of a stirred tank; however, the mixer in Chapter 5 is designed in the frequency domain rather than the time domain. The model for each channel is plug flow (piston flow) in the current work while a laminar flow model was used in [33] 4 ; furthermore, the most fundamental difference of our design with [33] is that we recirculate a part of 3 For example in [12] a method for designing such a mixer based on pulp cleaners has been proposed. 4 For more details on modeling of these flow regimes please see [32].  15  the outflow of the mixer and reintroduce it into the mixer inlet as recirculation has been found useful to enhance mixing [14]. In addition, we eliminate the constraint of identical flow rates for each of the parallel pipes, as this did not seem necessary. This makes our design appealing for microfluidics applications; however, as mentioned before the design is not restricted to microfluidics applications and can be implemented in different length scales to reduce fluctuations in concentration.  1.5  Objectives  At the beginning of the project several objectives were defined. The main objective was to develop a model to evaluate pulp chest mixing performance using inlet and outlet consistency data. The purpose of the model was to describe non-ideal physical phenomena occurring in the chest. It was also desired that the model would have a simple structure to make it more useful for industrial applications. The material in Chapter 2, Chapter 3, and Chapter 4 achieves this objective. The next objective was designing a new mixer specifically targeting the reduction of inlet concentration variations. This new type of mixer was meant to substitute the current stirred tank dominating the mixing industry. Chapter 5 is dedicated to designing and simulating this mixer.  1.6  Summary of contributions  Control engineering methods have evolved over the last half century as control engineers have explored time and frequency domains. In different chapters of this thesis methods from both domains are utilized: chapter two and five largely concentrate on the frequency domain approach and chapter three and four use the time domain approach to describe mixing phenomena. In the second chapter a new interpretation of flow behavior in the agitated pulp chest is given. A simple transfer function model is proposed which successfully describes the dynamic behavior of the agitated pulp chest. The new model is validated with the available data from the lab-scale pulp chest. The main advantage of the proposed new model is its simplicity which makes it more desirable for use in 16  the industry. The algorithm for identifying this model is available as a standard tool in MATLAB. The new model is also superior to the previous model [26] as it has the ability to correctly estimate the mixing performance when there is a significant bypass portion in the chest. Further explanation of these topics follows in the next chapter and are also available in Soltanzadeh et al. [44]. Material in the third chapter is one of the main contributions of this thesis. In this chapter a new method for identifying the residence time for an arbitrary system is presented. The method is simple and reliable, and more importantly it does not require prior knowledge of the system structure. It also accepts a wide range of input signals which allows it to easily identify the residence time from a reliable set of data. A version of this chapter is published in Soltanzadeh et al. [45]. Chapter 4 derives a method to identify higher moments of a linear system from a set of input and output data. This method shows reliable results in calculating the first three moments. The method also can be combined with other methods of system identification to verify an estimated model. Moreover, a new relationship between time delay and higher moments is explored in this chapter. In Chapter 5 a new theory for designing mixers is developed and verified. In this new mixing scheme the inlet flow is divided into several parallel channels with different lengths and diameters, and they are then recombined. At the same time, part of the outlet flow is recirculated to the inlet flow. Dimensions of these channels are calculated to enable the mixer to emulate the performance of a stirred tank. This process of splitting and recombining the flow attenuates the concentration variations in the inlet stream, and it also results in a more uniform outlet stream. This theory is then verified through a computational fluid dynamic (CFD) simulation in this chapter.  17  Chapter 2  Estimation of residence time and channeling in agitated pulp chests 2.1  Summary  Agitated pulp chests are widely used in the pulp and paper industry, but their design can be problematic. One main design concern is the existence of flow nonidealities; including channeling, recirculation and the formation of stagnant zones, that adversely impact mixing performance, are difficult to identify and are exacerbated by the non-Newtonian rheology of fibre suspensions. In earlier work in the UBC Dynamic Mixing Group [19–21, 26] Ein-Mozaffari et al. developed a parametric model to identify mixer dynamics and flow non-idealities in these chests. The model was able to evaluate mixing in both laboratory- and industrial-scale vessels but the selection of an appropriate system identification method and the redundancy of model parameters were of concern. Continuing improvement in the understanding of flow characteristics in stock chests has enabled the simplification of this earlier mixing model while also making the parameter identification process more robust. In this chapter the modifications to the model, the physical arguments allowing these simplifications and their impact on the ease and robustness of the solution technique, are discussed. Comparison of the two approaches is presented for a range of chest operating conditions.  18  2.2  Experimental  A laboratory scale model of an industrial chest was designed and used to investigate the influence of operating variables like impeller speed, exit location, chest dimension, pulp consistency and inlet flow on mixing performance of stock chests. The chest was rectangular in shape with dimensions of 0.4 x 0.7 x 0.7m and had a side entry impeller as shown in Fig. 2.1. The two configurations shown in Fig. 2.2 were studied to investigate the influence of the location of the inlet and exit streams relative to the impeller on the performance of the chest. Using this set-up, Ein-Mozaffari et al. [20] carried out extensive tests for a range of impeller speeds (N = 700-1750 rpm), pulp types, pulp mass consistencies (Cm = 2.1- 3.3%) and flow conditions (Q = 7.9-37.1 LPM). The dynamic response of the chest regarding change in inlet concentration was obtained by injecting saline solution into the inlet stream and measuring the variations in the inlet and output conductivity. The flow rate was kept constant during the experiment. For proper identification of the parameters in the mixing model it is necessary  Figure 2.1: Schematic diagram of the scale-model apparatus used for the dynamic tests [26]. to use an input signal that is rich in information and can excite the chest over the range of frequencies in which it is active. In order to do this, instead of employing 19  a single step change for exciting the chest, a binary sequence obtained by digitally switching on and off the saline flow to the inlet stream was used. The input signal was designed to have its energy concentrated at frequencies where the Bode plot is sensitive to parameter variations [30]. This would guarantee that the input signal is rich enough in excitation and identification will lead to valid parameters for the chest. An example of the typical input signal used for the tests is shown in Fig. 2.3. The  Figure 2.2: Stock chest configurations showing pulp input and exit locations investigated in the laboratory chest [26]. low pass characteristic of the chest is apparent in this figure. More details of the setup and experimental procedure are available in [18].  2.3  Simplified model  The model between inlet and outlet concentration of the chest developed by EinMozaffari et al. is described in Section 1.1. This model assumes two flow paths: one path for channeling and one path for mixing. Both path have different time delays and different time constants, and a recirculation flow also occurs inside the mixing path (Fig. 1.3). However, experimental (lab-scale) tests showed that the recirculation parameter, R, was insignificant (∼ = 0) and the time constant in the channeling path was unexpectedly high in some cases. The latter point is inconsistent 20  Figure 2.3: Example of scale-model chest input and output signals. Data shown for Configuration 2 with the impeller speed of 1105 rpm, pulp flow of 7.9 L/min and pulp consistency of 3.3%. with the definition of channeling, which is considered a direct passage of the fluid from input to output without significant mixing. These observations, coupled with the necessity of resolving the issues arising out of the form of equation (1.1) were the main motivation for revisiting the dynamic mixing model. In the new model, the compartmentalized nature of the flow-field in the chest has been retained, but the flow-paths through these two regions have been interpreted differently which allows for substantial simplification of the model. The revised model is shown in Fig. 2.4. Here the time delays in the two paths are assumed to be equal and are merged so that only one delay term is considered. It is also assumed that the time constant in the channeling path is much smaller than the time constant in the mixing path. The physical interpretation of this modification is that the pulp stream flows in plug flow to the entrance to the mixing zone where part of it separates and proceeds to the exit through the bypass zone. Thus both routes have similar delays. In stock chests having side-entry, axial impellers (the common industrial design) a circulation pattern develops where pulp flows out from the impeller discharge to21  wards the opposite wall and then loops back to the impeller suction along the side walls and surface of the suspension. Hence in these chests, especially for Configuration 1 (Fig. 2.2) it can be conjectured that pulp stock added at the surface will be carried towards the impeller by the circulating flow with part of the stream escaping out the chest exit and the remainder entering the impeller suction and thence into the mixing zone. Insignificant mixing happens in the bypass zone which is modeled as a first-order transfer function with time constant, τ1 . Significant mixing is expected to take place in the mixing zone and is modeled as a transfer function with a time constant, τ2 , as seen in Fig. 2.4 (τ1 ≪ τ2 ). On the basis of previous experimental evidence the recirculation loop has been removed (i.e. R= 0). The transfer function for this modified model can be written as G′ (s) = (  1 − f −T s f + )e . 1 + τ1 s 1 + τ2 s  (2.1)  which simplifies to 1 + [ f τ2 + (1 − f )τ1 ]s −T s e (1 + τ1 s)(1 + τ2 s) 1 + zs = e−T s (1 + τ1 s)(1 + τ2 s)  G′ (s) =  (2.2)  with, z = f τ2 + (1 − f )τ1 . This new model has been widely used in chemical engineering and has been shown to simulate flow characteristics reasonably well in a number of industrial processes such as distillation columns and heat exchangers [48]. Unlike the previous (two-delay) model (equation (1.1)) which requires complicated and specialized procedures for identifying the model parameters, for the one-delay model (equation (2.2)) standard routines can be used. In this work, the standard PEM function in MATLAB identification toolbox was used. The results from this new model show a very good fit between the data and the model, capturing the mixing dynamics of the chest predicted by the previous, more rigorous model. These results are described in the next section. There are several significant advantages of this simplified model. Firstly, the model can be identified directly in the continuous time domain. This is beneficial because the delay need not be a multiple of the sample time and no approximation  22  Figure 2.4: Simplified model with one time delay used for identifying the dynamics of pulp stock chests. Mixing in the bypass zone was found to be negligible. τ1 therefore approximates to zero and the transfer function, G′1 (s) to unity. error is introduced during transformation of the model from continuous to discrete time. Secondly, τ1 is almost zero in most cases, and due to the low pass nature of the system the z value in equation (2.2) is smaller than τ2 . As a result the f value will be between 0 and 1 by construction, and there is no necessity of exerting an additional constrain on the magnitude of this parameter. Thirdly, parameter convergence of this model is greatly simplified and the parameters will always converge to a real solution [38]. Moreover, the simplicity and robustness of the model implies that it is very suitable for evaluating online performance of the chest in real-time [40]. This can be very useful for continuous monitoring of mixing performance in sensitive processes. Finally, this model can describe the system with less parameters and thus should be chosen based on the parsimony principle. This principle states that when two or more competing models exist and describe a system equally well, the model with the smallest number of free parameters should be chosen [42]. In the next section, the model in equation (2.1) is identified using the MATLAB PEM function. In this function MATLAB uses a four-stage algorithm to identify delay, zero, and poles in equation (2.1). First it identifies a discrete model using the output error (OE) method, and then uses an autoregressive exogenous (ARX) 23  model to identify the delay. Then it converts this model to continuous time and identifies pole zero and gain with the state space model. The dimension of the state space depends on the integer value of the delay, and the non-integer part of the delay is estimated as a model parameter [29].  2.4  Results and discussions  Q=37.1((LPM), Cm=3.3%, old model Q=37.1(LPM), Cm=3.3%, new model Q=7.9((LPM), Cm=3.3%, old model Q=7.9(LPM), Cm=3.3%, new model  a) 1  Bypass fraction, f(−)  0.8  0.6  0.4  0.2  0 1000  1100  1200  1300  1400 1500 1600 Impeller speed, N (rpm)  1700  1800  1900  2000  Figure 2.5: Comparison of the single- and two-delay models. f is compared at two pulp flow-rates over a range of impeller speeds. The tests were carried out in Configuration 1 with Cm = 3.3%. The dynamics of mixing in the stock chests can be characterized by the fraction of the flow bypassing the mixing zone, f , and/or the fraction of the chest that is fully mixed,  V f ully−mixed VTotal .  While f is obtained directly,  V f ully−mixed VTotal  can be estimated  from f and the other identified parameters. Mixing is negligible in the bypass region hence  V f ully−mixed VTotal  can be obtained simply from the time constant in the mixing 24  1400 Q=37.1((LPM), Cm=3.3%, old model Q=37.1(LPM), Cm=3.3%, new model Q=7.9((LPM), Cm=3.3%, old model Q=7.9(LPM), Cm=3.3%, new model  b)  Time constant in mixing zone, τ2 (s)  1200  1000  800  600  400  200  0 1000  1100  1200  1300  1400 1500 1600 1700 Impeller speed, N (rpm)  1800  1900  2000  Figure 2.6: The time constant in the mixing zone is estimated using both models. Note that at N < 1300rpm and Q = 37.1 LPM the one-delay model properly identifies poor mixing in the chest because most of the pulp is short-circuiting. zone, τ2 , and the flow passing through this region, Q(1 − f ): V f ully−mixed (1 − f )Qτ2 = VTotal VTotal  (2.3)  In the following paragraphs, the values of the two parameters estimated using the simplified one-delay model (equation (2.2) and the two-delay model (equation (1.1)) are compared with each other for selected laboratory tests. A good match between the two models indicates the ability of the simpler model to capture the dynamics identified by the more complex one. As previously mentioned, these values are estimated through a system identification process consisting of estimation and verification steps. Error bars are estimated for the one-delay model based on 25  160  Q=37.1((LPM), Cm=3.3%, old model Q=37.1(LPM), Cm=3.3%, new model Q=7.9((LPM), Cm=3.3%, old model Q=7.9(LPM), Cm=3.3%, new model  c)  140  Vfully mixed/Vtotal%(−)  120 100 80 60 40 20 0 1000  1100  1200  1300  1400 1500 1600 1700 Impeller speed, N (rpm)  1800  1900  2000  Figure 2.7: Volume fraction of the chest that is fully mixed estimated for both models under the operating conditions defined above. The single-delay model is able to capture all the dynamic characteristics revealed by the two-delay model. the accuracy of the estimated model in the system identification process. Figures 2.5, 2.6, and 2.7 show the key parameters identified using the two models for Configuration 1 with a pulp consistency of 3.3% and at flow rates of 37.1 and 7.9 LPM. In Fig. 2.5, the results from the two models are compared on the basis of the bypass fraction, f . It is seen that at the lower flowrate (Q = 7.9 LPM) the amount of channeling decreased continuously with increasing impeller speed and that the trends are the same for both models. When the pulp flowrate is increased to 37.1 LPM, not only does the bypass fraction increase but there is also a sharp change in the chest dynamics with an equally sharp increase in the amount of channeling, to almost 100%, when the impeller speed is reduced only slightly, from 1350 to 1300 rpm. This is an important observation as it suggests that at 26  higher pulp consistencies a small change in operating conditions, such as impeller speed, can cause a dramatic difference in mixing efficiency. The trend identified by the two-delay model ([19]) is reproduced well by the single-delay model, indicating that it capable of capturing the mixing dynamics despite its simplicity. In Figure 2.6 the time constants in the flow zones are directly compared for the two models using the same conditions as in Fig. 2.5. As in the case of f , τ2 also changes continuously with N at the lower flow rate. At the higher flow rate, both models correctly predict the overall reduced levels of mixing (smaller time constants). However, it is only the single-delay model that correctly captures the sudden drop in mixing in the tank. This can be confirmed by inspecting the data, and also from our understanding of chest behavior which complete short-circuiting is expected at this impeller speed (N < 1300 rpm, see Fig. 2.6). The two-delay model, on the other hand, is unable to identify this feature. It should be mentioned that τ1 values for the one-delay model are insignificant (between 0 sec and 5 sec) and the mixing is appropriately represented by τ2 and f values in this figure. To illustrate this matter further, sample data which represents full channeling in the chest is shown in Fig. 2.8. As can be seen in this figure, there is no mixing in the chest and the data shows complete plug flow characteristics. Figure 2.7 presents a comparison of the fraction of the tank volume that is fully-mixed (well-agitated by the impeller), a parameter that can also be used to characterize mixing efficiency. There are some differences in the quantitative values estimated by the two models. The single delay model suggests higher mixing volume for the chest. In many cases the time constant of the bypass zone estimated from the model with two time delays is not negligible; however this time constant has not been considered in the calculation of the mixing volume [21]. As a result the overall trend for the single-delay model is similar to that of the two-delay model. However, the two-delay model predicts inferior mixing performance as the fully mixed volume calculation does not consider the time constant of the channeling region. For these test conditions, it can be concluded that the single-delay model is able to estimate the mixing performance better than the more complex model. One the other hand, for many cases the two-delay model describes data in greater detail as it has an extra parameter. Figures 2.9 and 2.10 present another test of single-delay model. Here the effects 27  3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 input output  1.2 100  200  300  400 time (s)  500  600  700  800  Figure 2.8: Low impeller speed can lead to full channeling in the chest. of pulp consistency and chest configuration on f and  V f ully−mixed VTotal  are investigated.  While industrial chests use Configuration 1 almost exclusively due to operating constraints, investigations with the two-delay model have shown that Configuration 2 provides a better flow path for good mixing [20]. It is apparent from Fig. 2.9 that the bypass fraction in this configuration is close to zero. It is seen from Figs. 2.9a and 2.9b that the one delay model is also able to correctly identify the low levels of short-circuiting in Configuration 2. The effect of pulp consistency on mixing dynamics can be seen by comparing Fig. 2.9a (Cm = 3.3%) and 2.9b (Cm = 2.7%). Mixing is better at lower consistencies. There is good agreement in the fraction of the tank that is fully mixed as identified by both models. Figures 2.10a and 2.10b show the effect of chest configuration and pulp consistency  28  Q=37.1((LPM), Cm=3.3%, old model, config. 1 Q=37.1(LPM), Cm=3.3%, new model, config. 1 Q=37.1((LPM), Cm=3.3%, old model, config. 2 Q=37.1(LPM), Cm=3.3%, new model, config. 2  Bypass fraction, f(−)  1 0.8 0.6 0.4 0.2 0 1000  1100  1200  1300 1400 1500 Impeller speed, N (rpm)  1700  1800  Q=37.1((LPM), Cm=2.7%, old model, config. 1 Q=37.1(LPM), Cm=2.7%, new model, config. 1 Q=37.1((LPM), Cm=2.7%, old model, config. 2 Q=37.1(LPM), Cm=2.7%, new model, config. 2  1 Bypass fraction, f(−)  1600  0.8 0.6 0.4 0.2 0 800  900  1000  1100  1200 1300 1400 Impeller speed, N (rpm)  1500  1600  1700  1800  Figure 2.9: The effect of pulp consistency and chest configuration on parameter estimation using the two models. Bypass fraction, f , is estimated at pulp flow-rate of 37.1 LPM for both chest configurations at two pulp consistencies: (a) Cm = 3.3% and (b) Cm = 2.7%. on the fraction of chest volume that is fully mixed for the conditions described in Fig. 2.9. When Configuration 2 is used or when pulp consistency is decreased, mixing efficiency is improved and a larger volume of the chest remains well-mixed. The improvement in mixing stems from the channeling flow reduction in the second configuration. Also, when pulp suspensions have a lower concentration less energy is needed to agitate the suspension; therefore, the mixing performance is improved while the impeller speed remains the same.  29  120  %(−)  60 40  V  fully mixed  total  80  /V  100  Q=37.1((LPM), Cm=3.3%, old model, config. 1 Q=37.1(LPM), Cm=3.3%, new model, config. 1 Q=37.1((LPM), Cm=3.3%, old model, config. 2 Q=37.1(LPM), Cm=3.3%, new model, config. 2  20 0 1000  1100  1200  1300 1400 1500 Impeller speed, N (rpm)  1600  1700  1800  120 100  60  V  fully mixed  /V  total  %(−)  80  40 Q=37.1((LPM), Cm=2.7%, old model, config. 1 Q=37.1(LPM), Cm=2.7%, new model, config. 1 Q=37.1((LPM), Cm=2.7%, old model, config. 2 Q=37.1(LPM), Cm=2.7%, new model, config. 2  20 0 800  900  1000  1100  1200 1300 1400 Impeller speed, N (rpm)  1500  1600  1700  1800  Figure 2.10: The effect of pulp consistency and chest configuration on the estimation of the fully mixed volume using the one-delay and two-delay models for the operating conditions defined in Figure 2.9. Results for (a) Cm = 3.3% and (b) Cm = 2.7%.  2.5  Conclusions  Issues related to the identification of mixing performance using a two-delay model developed previously [19–21, 26] are addressed in this chapter. The two-delay model visualized a part of the incoming stock flowing into the mixing zone (cavern) around the impeller while part of it bypassed the zone and exited the chest with less mixing with each flow path associated with a different delay. While this model  30  was able to identify dynamic features that agreed with the physical understanding of the system, there were a number of issues in the identification process arising out of the complexity of the model. A new single-delay model has been proposed having a simplified form that leads to improved parameter convergence. Also, the model provides a better estimate of the mixing performance when there is significant channeling in the chest. Moreover, the model is suitable for use as an online tool to continuously monitor chest performance. In the next chapter a novel identification method will be introduced. As will be shown, it is possible to identify residence time without a model.  31  Chapter 3  Direct estimation of residence time from input and output data 3.1  Summary  Residence time is a well known and widely used concept in the process industry, and its estimation is often needed to model or optimize a process. This chapter will demonstrate a new method using suitable measures of input and output data. There are two distinct characteristics of this method: it is not necessary to estimate the residence time from an explicit model, and the constraints on the excitation of the input signal are not as restrictive as common methods.  3.2  Proposed method  To develop the proposed method, it is necessary to prove two theorems: one theorem where the input signal asymptotically decays to zero and one where the input signal does not satisfy this condition. Although the first theorem has been derived in Astr¨om and Hagglund [6], it is useful to prove it here again. The second theorem is the main contribution of this chapter and it will lead to robust identification of the residence time.  Theorem 1 Assume initial conditions are zero, and input and output data have 32  non-zero time averages. Also, assume the input signal u multiplied by time decays to zero asymptotically. Then this relationship holds between residence time and the input and output data of any stable casual linear time invariant system: ∫  ∫  ∞ ∞ t.u(t) dt t.y(t) dt − ∫0 ∞ τdyn = ∫0 ∞ y(t) dt 0 0 u(t) dt  (3.1)  where y and u are the output and input signals, respectively. Proof According to the theorem’s assumption the system is stable and the output multiplied by time decays asymptotically to zero. As a result, the following integrals are well defined  ∫ ∞  ∫ ∞  t.y(t) dt,  t.u(t) dt  0  0  Also, the average of input and output signals is non-zero. Therefore, the following quantities are also well defined ∫∞  ∫  t.y(t) dt 0∞ t.u(t) dt , ∫∞ 0 y(t) dt 0 u(t) dt  ∫0 ∞  To prove the theorem we start from the first term of the right-hand side of (3.1). From the Laplace transform properties ∫∞  t.y(t) dt Y ′ (s) =− Y (s) 0 y(t) dt  ∫0 ∞  .  (3.2)  s=0  From definition of the transfer function G(s) ∫∞  t.y(t) dt (G(s)U(s))′ =− G(s)U(s) s=0 0 y(t) dt G′ (s)U(s) + G(s)U ′ (s) =− G(s)U(s) ( ′ ) G (s) U ′ (s) =− + G(s) U(s) s=0  ∫0 ∞  (3.3) s=0  By rearranging this last equation G′ (s) − G(s)  ∫∞  t.y(t) dt U ′ (s) = ∫0 ∞ + U(s) 0 y(t) dt s=0 33  (3.4) s=0  From the Laplace transform properties ∫∞  t.u(t) dt U ′ (s) =− U(s) 0 u(t) dt  ∫0 ∞  (3.5) s=0  Now by using equations (3.5), residence time definition (equations (1.2) and (1.3)) in (3.4) we obtain  τdyn =  ∫  ∫∞  ∞ t.y(t) dt 0 t.u(t) dt ∫ − ∞ ∞ 0 y(t) dt 0 u(t) dt  ∫0  (3.6)  which proves the theorem. Although the upper limit of integrals in equation (3.1) is infinite, for practical purposes it is sufficient to estimate the residence time from data covering a long but finite period of time. The following example shows this point. Example 3.1. Consider the following system G(s) =  e−700s (1300s + 1) (1057s + 1)(2300s + 1)(1400s + 1)  (3.7)  By using equation (1.3), the residence time is estimated as 4157s for this system. By applying a pulse signal with an amplitude of 1 and a width of 4 hours, the residence time can be estimated from equation (3.1). To calculate the right-hand side of equation (3.1), the data was sampled at 12 sec intervals and then the first order approximation was used to calculate the integrals. The plot of input and output data and the result of estimated residence time is shown in Fig. 3.1. As can be seen in this figure, the residence time has converged to 4159s which is very close to the theoretical value. The difference is caused by the numerical integration of the quantities in the right-hand side of equation (3.1). To emphasize the usefulness of Theorem 1, assume one wants to estimate the residence time by identifying a model for the system. The first difficulty would be finding a model that describes the system well. Ideally this model should have three poles, one zero and one delay; however, this is difficult to realize if there is no prior understanding of the system’s structure. A rich enough signal should then excite the system to identify all five coefficients of the model. After gathering these data, an appropriate method should be deployed to identify the model. This step is not trivial because of the delay in 34  1  input output  0.5  0 0  2  4  6  8  10  6  8  10  (hour)  (sec)  4000  2000  0 0  2  4 (hour)  Figure 3.1: Input and output signal (top). Right-hand side of equation (3.1) (bottom) the model. Finally, the residence time would be calculated from equation (1.3). The result would be prone to error if each of the previously mentioned steps was not carried out properly. By using Theorem 1, one can eliminate these steps and simply identify the residence time from a single bump test. This result illustrates the significance of Theorem 1. Unfortunately, as mentioned in the introduction this method does not give a reliable estimation for a general type of input signals. The right-hand side of equation (3.1) can be calculated for both the model in example 1 and the first-order model with same residence time , 4157s, and the results can then be compared. As can be seen in Fig. 3.2, the expression in the right hand side of equation (3.1) for both systems become close after a transient period. This is expected as the result of Theorem 1 is only valid asymptotically. However, it is not accurate to use the result of Theorem 1 here. As can be seen in Fig. 3.2,  35  the input signal does not satisfy the restrictions needed for validation of Theorem 1 because it does not decay to zero with a rate faster than t. The following theorem extends the result of Theorem 1, and proves that for two systems with the same residence time and the same input excitation the right-hand side of equation (3.1) is asymptotically equal.  1  input output  0.5  0 0  10  20  30  40  50  (hour) Model in example 1 first−order model  (sec)  3000 2000 1000 0 0  10  20  30  40  50  (hour) Figure 3.2: Input signal (top). Right-hand side of equation (3.1) for first order model and model in example 1.  Theorem 2 Assume linear time invariant stable causal systems G1 (s) and G2 (s) which have the same residence time. For any input signal u, the corresponding outputs y1 and y2 of these systems will satisfy the following equation ∫∞  ∫  ∞ t.y1 (t) dt t.y2 (t) dt = ∫0 ∞ ∞ 0 y1 (t) dt 0 y2 (t) dt  ∫0  36  (3.8)  Proof It has been assumed that G1 (s) and G2 (s) are linear time invariant causal systems, then from linear system theory ∫ ∞  y1 (t) =  0  and  ∫ ∞  y2 (t) =  0  u(τ ).g1 (t − τ ) d τ  (3.9)  u(τ ).g2 (t − τ ) d τ  (3.10)  Now substituting (3.9) in the left hand side of (3.8) ∫∞  t.y1 (t) dt = 0 y1 (t) dt  ∫0 ∞  ∫∞ 0  ∫∞  t.(  0  u(τ ).g1 (t − τ ) d τ )dt ∫∞ 0 y1 (t) dt  (3.11)  By bringing t inside the integral and then changing the order of the integration ∫∞  ∫∞∫∞  t.y1 (t) dt = 0 y1 (t) dt  ∫0 ∞  0  t.u(τ ).g1 (t − τ ) dtd τ ∫∞ 0 y1 (t) dt  0  (3.12)  Now u(τ ) is not a function of t; therefore, it can be moved outside the inner integral ∫∞  ∫∞  t.y1 (t) dt = ∞ 0 y1 (t) dt  ∫0  0  ∫∞  u(τ )  ∫0∞ 0  t.g1 (t − τ ) dtd τ y1 (t) dt  (3.13)  Introducing new variable k = t − τ in the inside integral of this last equation we get ∫∞  t.y1 (t) dt = ∞ 0 y1 (t) dt  ∫0  ∫∞ 0  ∫∞  u(τ )  −τ (k + τ ).g1 (k) dkd τ ∫∞ . 0 y1 (t) dt  (3.14)  As τ > 0 and g1 (k) = 0 for k < 0 the lower limit of the last integral can be set to zero, so the last equation can be written as ∫∞  t.y1 (t) dt = ∞ 0 y1 (t) dt  ∫0  ∫∞ 0  ∫∞  u(τ )  37  0∫ (k + τ ).g1 (k) dkd τ ∞ 0 y1 (t) dt  (3.15)  which is equal to ∫∞  t.y1 (t) dt 0 y1 (t) dt ∫∞ ∫∞ ∫∞ 0 u(τ )( 0 k.g1∫(k) dk + 0 τ .g1 (k) dk)d τ = ∞ 0 y1 (t) dt ∫∞ ∫ ∫ ∫∞ τ + 0∞ u(τ ) 0∞ τ .g1 (k) dkd τ 0 u(τ ) 0 k.g1 (k) dkd ∫ = ∞ 0 y1 (t) dt ∫0 ∞  In the last equation  ∫∞ 0  k.g1 (k) dk and  ∫∞ 0  (3.16)  g1 (k) dk are not a function of τ ; therefore,  it is possible to bring them out of the τ integrals which results in ∫∞  t.y1 (t) dt 0 y1 (t) dt ∫∞ ∫∞ ∫ ∫ dk + 0∞ τ .u(τ ) d τ 0∞ g1 (k) dk 0 u(τ ) d τ 0 k.g1 (k) ∫∞ = 0 y1 (t) dt ∫0 ∞  (3.17)  This also can be written as ∫∞  t.y1 (t) dt 0 y1 (t) dt  ∫0 ∞  ∫ ∞  = 0  ∫∞  g1 (k) dk ∗  0  u(τ ) d τ  ∫∞ k.g (k) dk ∫0 ∞ 1  ∫ + 0∞ τ .u(τ ) d τ g (k) dk 1 0 ∫∞ 0 y1 (t) dt  (3.18)  On the other hand, from the theorem’s assumptions the residence time of G1 (s) is equal to the residence time of G2 (s); in other words we have ∫  ∫∞  ∞ t.g1 (t) dt 0 t.g2 (t) dt ∫ = . ∞ ∞ 0 g1 (t) dt 0 g2 (t) dt  ∫0  (3.19)  By substituting (3.19) in (3.18) we get the following result ∫∞  t.y1 (t) dt 0 y1 (t) dt  ∫0 ∞  = =  ∫∞  ∫∞ k.g (k) dk ∫0 ∞ 2  ∫∞ + 0 τ .u(τ ) d τ g (k) dk 0 2 ∫∞ g1 (k) dk ∗ y1 (t) dt 0 ∫∞ 0 ∫ ∫ ∫∞ ∫∞ g1 (k) dk 0 u(τ ) d τ 0 k.g2 (k) dk + 0∞ τ .u(τ ) d τ 0∞ g2 (k) dk ∫0∞ ∫∞ ∗ . 0 g2 (k) dk 0 y1 (t) dt  ∫ ∞  0  u(τ ) d τ  38  (3.20)  Therefore, it has been shown that the left side of (3.8) is equal to (3.20). By starting from the right-hand side of equation (3.8) and by following the previous process, the next equation is achieved ∫∞  t.y2 (t) dt = 0 y2 (t) dt  ∫0 ∞  ∫∞ 0  u(τ ) d τ  ∫∞ 0  ∫  k.g2 (k) dk + 0∞ τ .u(τ ) d τ ∫∞ 0 y2 (t) dt  ∫∞ 0  g2 (k) dk  .  (3.21)  Combining (3.21) and (3.20) ∫∞  ∫  ∫  ∫  ∞ t.y1 (t) dt g1 (t) dt 0∞ y2 (t) dt 0∞ t.y2 (t) dt = ∫0∞ ∗ ∫∞ ∗ ∫∞ ∞ 0 y1 (t) dt 0 g2 (t) dt 0 y1 (t) dt 0 y2 (t) dt  ∫0  (3.22)  , and to prove the theorem it is necessary to show ∫∞  g1 (t) dt ∫0∞ 0 g2 (t) dt  ∫∞  ∗ ∫0∞ 0  y2 (t) dt =1 y1 (t) dt  (3.23)  If we call G1 (s) and G2 (s) the Laplace transform of systems 1 and 2, respectively, and U(s) the Laplace transform of an input signal, the following equations are true for any LTI system ∫ ∞ ∫0 ∞ 0  ∫ ∞  g1 (t) dt = G1 (0),  0  g2 (t) dt = G2 (0) ∫ ∞  y1 (t) dt = G1 (0)U(0),  0  (3.24) y2 (t) dt = G2 (0)U(0)  By simple substitution of these equations in the left hand side of (3.23) the right hand side will be achieved, and this proves the theorem. In the proof of Theorem 2 there is no restriction on the order of systems G1 and G2 . Therefore, it can be concluded that for any arbitrary linear system there is a firstorder system G(s) =  1 τ s+1  which satisfies Theorem 2. τ in this first-order system  would be the residence time of the corresponding system. This is described in the following Corollary. Corollary 1 Consider an arbitrary linear time invariant system G(s) and the first order system H(s) =  1 τ s+1 .  Assume input signal u(t) has been exerted on both  of these systems and corresponding output signals yG (t) and yH (t) have been ob-  39  tained. Then the residence time of G(s) (τdyn ) is equal to  τdyn = arg min τ  (∫ ∞  ∫  ∞ t.yH (t) dt 0 t.yG (t) dt ∫ − ∞ ∞ 0 yH (t) dt 0 yG (t) dt  ∫0  )2 (3.25)  Proof The minimum of this cost function is zero, and in Theorem 2 it was shown that this would occur when τ is equal to the residence time. This minimum is unique as H(s) is first order and only has one parameter. From this Corollary it can be concluded that, in order to identify the residence time, it is necessary to have input and output data from the system. In the next step the same input should be exerted on a first-order model and then the pole of the firstorder model should be estimated in a way that minimizes Equation (3.25). This optimization has only one parameter and it has no constraints. MATLAB function fminsearch can easily solve this kind of problem. The fminsearch function uses the Nelder-Mead simplex method which is a direct search algorithm that does not use gradient implicitly or explicitly [27]. When the optimization search space contains one parameter the simplex is a triangle. At each step of the search, a new point in or near the current simplex is generated using specific predefined rules. The function is then evaluated at the new point, and it is compared with the function’s values at the vertices of the simplex. The vertex to be replaced is determined by the outcome of this comparison. This processes provides a new simplex. This step is then repeated until the diameter of the simplex is less than the specified tolerance. In the next section this method will be used to show how Corollary 1 can be applied to several examples.  3.3  Results and discussion  Two different examples will be considered. First, the residence time of the system in Example 1 will be identified from the result of Corollary 1, and then this method will be applied to real data from a laboratory-scale pulp chest. To find the residence time it is necessary to minimize the cost function in equation (3.25); as previously mentioned the MATLAB function f minsearch has been used for this purpose. Moreover, for faster convergence an initial condition has been introduced for the first-order system. This initial condition has also been found 40  input output  1 0.5 0 0  5  10  15 (hour)  20  25  30  4000 (sec)  3000 2000 1000 0 0  data model 5  10  15 20 (hour)  25  30  Figure 3.3: Input and output data (top). Right-hand side of equation (3.1) calculated from input and output data and also estimated first order model (bottom). during the optimization process. The input and output signals are shown in Fig. 3.3 (top). This input signal consists of periodic signals with a period of 33256s and additive white noise with a power of 0.1. The quantities on the right-hand side of equation (3.1) are shown in Fig. 3.3 (bottom) for the original system and the identified system. As mentioned in Theorem 1, the right-hand side of equation (3.1) is an estimate for residence time, but it is apparent from Fig. 3.3 (bottom) that it does not converge for the type of signal in Fig. 3.3 (top); however, the same quantity calculated from the first-order model and the system in Example 1 are almost tangent to each other in the steady state area. This is the sign that the value for the estimated residence time is correct based on the result of Theorem 2. The value that has been estimated from the optimization is 4122s, which is close to the actual  41  residence time of 4157s. Unlike Fig. 3.1, the input signal used in this example does not decay to zero, and the system is not at rest at time zero. The estimate of Theorem 1 does not converge to any specific value in this case (Fig. 3.3 (bottom) dashed line) but Theorem 2 gives a very good estimate. The initial condition effect was not discussed in Theorem 2, but it is obvious that it is possible to shift time zero to the point where the effect of the initial conditions are insignificant, and then apply the theorem. Therefore, this method does not fail when the system does not start from non-zero initial conditions. This example also shows the method is not sensitive to additive input noise. The method is able to handle additive output noise as well. To decrease the sensitivity against output noise, the filtered version of the output signal should be used rather than the original noisy output. The residence time of the system would then be the estimated residence time from the filtered output minus the residence time of the filter. The last example that concludes this section uses conductivity data from an experiment that was discussed in Chapter 2. The input and output signal and the right-hand side of equation (3.1) are shown in Fig. 3.4. The value calculated for the residence time in this case is 175s. In Chapter 2 the residence time for this data was identified by assuming a first-order model plus delay for the system. The value calculated for these quantities was 13s for delay and 155s for the first-order pole. By adding these values the residence time can be calculated as 168s, which is very close to the value calculated from the proposed method.  3.4  Conclusion  In this chapter a new method for identifying residence time has been proposed. The proposed method requires very few restrictions on the excitation of the input signal. The method only needs several changes in the steady-state part of the input signal, which is the most commonly used input signal in the process industry. The simplicity of this method also makes it appealing for process identification. Users do not need to have a comprehensive understanding of system identification to use this method. The other advantage of this method is having an integral in the cost function which makes it robust against noise. On the other hand, the method in the  42  10  input output  9 8 0  0.05  0.1  0.15 (hour)  0.2  0.25  0.3  8  (sec)  6 4 2 data model  0 −2 0  0.05  0.1  0.15 (hour)  0.2  0.25  0.3  Figure 3.4: Input and output (top). Right-hand side of equation (3.1) calculated from input and output and also estimated first order model (bottom). form presented here cannot readily be used in real-time identification. This can be solved by finding a closed form solution for the optimization problem in Equation (3.25), which could be a topic for future research. In the next chapter the result of Theorem 2 will be extended to higher moments of the transfer function. It will also be shown how a similar method can be used to estimate the higher moments of a system from a set of input and output data.  43  Chapter 4  Estimation of higher moments of a system from its input and output data 4.1  Summary  In this chapter the results of the second theorem of Chapter 3 are extended to the higher moments of a transfer function. Also, an algorithm to estimate the higher moments is proposed. By combining the method introduced in this chapter and available methods for reconstructing a function from its moments, it is possible to identify a system from its input and output data. In other words, this chapter will propose a new approach in system identification that does not depend on a predefined model. In the next section a theorem that extends the result of the second theorem of Chapter 3 is introduced. Also, an algorithm for identifying higher moments of a transfer function is proposed. Moreover, by using this algorithm an upper and lower limit for a time delay in a class of linear systems is derived. In the final section of this chapter the application of these results is demonstrated through an example.  44  4.2  Proposed method  The method is presented through the following theorem and algorithm.  Theorem 1 Consider linear time invariant causal systems G1 (s) and G2 (s). For input signal u, the corresponding outputs y1 and y2 of these systems will satisfy the following equation ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞ 0 y2 (t) dt n−1 ( )  −∑  i=0  µgn2  n i µ i g2  ∫∞ 0  0  y1 (t) dt  ( ) ∫ u(τ ).τ n−i d τ n−1 n i 0∞ u(τ ).τ n−i d τ ∫∞ µg1 ∫ ∞ +∑ = 0 u(t) dt 0 u(t) dt i=0 i  (4.1)  − µgn1  where µgn is the nth moment of a transfer function G(s) as defined in equation (1.5). Proof It has been assumed that G1 (s) and G2 (s) are linear time invariant causal systems. Then from linear system theory ∫ ∞  y1 (t) = and  0  ∫ ∞  y2 (t) =  0  u(τ ).g1 (t − τ ) d τ  (4.2)  u(τ ).g2 (t − τ ) d τ  (4.3)  By substituting (4.2), the second expression in the left hand side of equation (4.1) can be written as ∫∞ n ∫∞ n ∫∞ t .( 0 u(τ ).g1 (t − τ ) d τ )dt 0∫ t .y1 (t) dt ∫∞ = 0 . ∞  y1 (t) dt  0  0  y1 (t) dt  (4.4)  By bringing t n inside this integral and by changing the order of integration the following equation is achieved ∫∞ n ∫∞∫∞ n t .u(τ ).g1 (t − τ ) dtd τ 0∫ t .y1 (t) dt = 0 0 ∫∞ . ∞ 0  y1 (t) dt  0  45  y1 (t) dt  (4.5)  In equation (4.5), u(τ ) is not a function of t; therefore, it can be moved outside the inner integral  ∫∞ n ∫∞ ∫∞ n 0∫ t .y1 (t) dt 0 u(τ ) ∫0 t .g1 (t − τ ) dtd τ = . ∞ ∞  y1 (t) dt  0  0  y1 (t) dt  (4.6)  By introducing a new variable k = t − τ in the inside integral, equation (4.6) can be ∫∞ ∫∞ ∫∞ n n 0 u(τ ) −τ (k + τ ) .g1 (k) dkd τ 0∫ t .y1 (t) dt ∫ = . ∞ ∞  written as  0  y1 (t) dt  0  y1 (t) dt  (4.7)  We also know that τ > 0 and g1 (k) = 0 for k < 0; therefore, the lower limit of the last integral can be set to zero, and the last equation can be written as ∫∞ n ∫∞ ∫∞ n 0 u(τ ) 0 ∫(k + τ ) .g1 (k) dkd τ 0∫ t .y1 (t) dt = ∞ ∞ 0  y1 (t) dt  0  y1 (t) dt  (4.8)  which is the same as ∫ ∞ n (n) n n−i ∫∞ ∫∞ n .g1 (k) dkd τ 0 u(τ ) 0 ∑i=0 i k τ 0∫ t .y1 (t) dt ∫∞ = . ∞ 0  y1 (t) dt  0  y1 (t) dt  (4.9)  Similarly the following equation can be derived for G2 (s) ∫ ∞ n (n) n n−i ∫∞ ∫∞ n .g2 (k) dkd τ 0 u(τ ) 0 ∑i=0 i k τ 0∫ t .y2 (t) dt ∫∞ . = ∞ 0  y2 (t) dt  0  y2 (t) dt  (4.10)  By subtracting the same sides of equations (4.9) and (4.10) from each other, the following equation is achieved ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ = ∞  y2 (t) dt  y1 (t) dt  ∫ ∞ n (n) i n−i ∫∞ ∫ ∞ n (n) i n−i .g2 (k) dkd τ .g1 (k) dkd τ 0 u(τ ) 0 ∑i=0 i k τ 0 u(τ ) 0 ∑i=0 i k τ ∫∞ ∫∞ − 0  ∫∞  0  0  y2 (t) dt  0  y1 (t) dt  (4.11)  46  This equation can be written as ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞ 0  y2 (t) dt  0  y1 (t) dt  ∫ ∞ n−1 (n) i n−i ∫∞ ∫ ∞ n−1 (n) i n−i .g2 (k) dkd τ .g1 (k) dkd τ 0 u(τ ) 0 ∑i=0 i k τ 0 u(τ ) 0 ∑i=0 i k τ ∫∞ ∫∞ + = − ∫∞  y2 (t) dt  0 ∫∞ n ∫ ∫∞ u( u(τ ) 0∞ kn .g1 (k) dkd τ τ ) k .g (k) dkd τ 2 0 0 0 ∫∞ ∫∞ −  ∫∞  0  y2 (t) dt  0  0  y1 (t) dt  y1 (t) dt  (4.12) By changing the order of summation and integration and splitting τ and k integrals, the last equation can be written as ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞  0 y2 (t) dt 0 y1 (t) dt ∫ ∫ ( ) ∫ ∫ ( ) n−1 ∞ u(τ ).τ n−i d τ ∞ n ki .g (k) dk n−1 ∞ u(τ ).τ n−i d τ ∞ n ki .g (k) dk 2 1 0 0 i 0 ∫∞ ∫∞ 0 i − + 0 y2 (t) dt 0 y1 (t) dt i=0 i=0 ∫∞ n ∫∞ n ∫∞ ∫∞ u(τ ) 0 k .g1 (k) dkd τ 0 u(τ )∫ 0 k .g2 (k) dkd τ ∫∞ − 0 ∞ 0 y2 (t) dt 0 y1 (t) dt  ∑  ∑  =  (4.13) From the definition of moments (equation (1.5)) ∫ ∞  t n .g(t) dt =  ∫ ∞  0  0  µgn g(t) dt.  (4.14)  Therefore equation (4.13) can be written as ∫∞ n ∫∞ n 0∫ t .y2 (t) dt 0∫ t .y1 (t) dt − ∞ ∞  0 y2 (t) dt 0 y1 (t) dt ∫ ∫∞ ∫ ( ) ∫ ( ) n−i d τ ∞ n µ i .g (k) dk n−1 n−1 ∞ u(τ ).τ n−i d τ ∞ n µ i .g (k) dk g2 2 g1 1 0 u(τ ).τ 0 0 i ∫∞ ∫∞ 0 i + − 0 y2 (t) dt 0 y1 (t) dt i=0 i=0 ∫∞ ∫∞ ∫∞ ∫∞ u(τ )d τ g2 (k) dk u(τ )d τ 0 g1 (k) dk ∫∞ µgn2 0 ∫ ∞ 0 − µgn1 0 0 y2 (t) dt 0 y1 (t) dt  ∑  ∑  (4.15)  47  =  or ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞ 0  y2 (t) dt  0 y1 (t) dt n−1 ∫ ∞ ( )  ∫∞  ∫∞ g2 (k) dk g1 (k) dk n i n−i 0 − ∫∞ ∗ µg2 u(τ ).τ d τ + ∫0∞ i 0 y2 (t) dt 0 y1 (t) dt i=0 0 ∫∞ ∫∞ ∫∞ ∫∞ u(τ )d τ g2 (k) dk u(τ )d τ 0 g1 (k) dk ∫∞ µgn2 0 ∫ ∞ 0 − µgn1 0 0 y2 (t) dt 0 y1 (t) dt  ∑  n−1 ∫ ∞ ( ) n  ∗∑  i=0 0  i  µgi 1 u(τ ).τ n−i d τ =  (4.16) We also know that for any linear time-invariant system G(s) (see equation (3.24) for proof)  ∫ ∞ 0  ∫∞  g(t)dt = ∫0∞ 0  y(t)dt u(t)dt  (4.17)  By substituting the last equation in equation (4.16) we get ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞ 0 y2 (t) dt n−1 ( )  −∑  i=0  µgn2  n i µ i g2  ∫∞ 0  0  y1 (t) dt  ( ) ∫ u(τ ).τ n−i d τ n−1 n i 0∞ u(τ ).τ n−i d τ +∑ µg1 ∫ ∞ = 0 u(t) dt 0 u(t) dt i=0 i  ∫∞  (4.18)  − µgn1  which proves the theorem. In the proof of Theorem 1 there is no restriction on the order of systems G1 and G2 ; therefore, the first order system G2 (s) =  1 as+1  with a residence time of a will  satisfy Theorem 1. However, for the first order system, G2 , all the higher moments are known and can be calculated as (this can be shown by taking derivatives from G2 (s))  µgn2 = n!an .  48  (4.19)  Therefore, by substituting the higher moments of the first-order system in equation (4.1) the following equation is obtained ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt − 0∫ ∞ ∞ 0 y2 (t) dt n−1 ( )  −∑  i=0  µgn2  n i!ai i  ∫∞ 0  0  y1 (t) dt  ( ) ∫ u(τ ).τ n−i d τ n−1 n i 0∞ u(τ ).τ n−i d τ +∑ = µg1 ∫ ∞ 0 u(t) dt 0 u(t) dt i=0 i  ∫∞  (4.20)  − µgn1  Both sides of equation (4.20) can be squared, and the following equation is obtained as a result ∫∞ n ∫∞ n t .y1 (t) dt 0∫ t .y2 (t) dt ( ∞ − 0∫ ∞ 0 y2 (t) dt n−1 ( )  −∑  i=0  (µgn2  n i!ai i  ∫∞ 0  0  y1 (t) dt  ( ) ∫ u(τ ).τ n−i d τ n−1 n i 0∞ u(τ ).τ n−i d τ 2 ∫∞ +∑ ) = µg1 ∫ ∞ 0 u(t) dt 0 u(t) dt i=0 i  (4.21)  − µgn1 )2  The minimum of the right-hand side of this equation is zero, and it is achieved when µgn2 = µgn1 ; therefore, by minimizing the left-hand side of equation (4.21) with regards to the residence time of the first-order transfer function, the nth moment of the unknown transfer function G1 (s) can be estimated. The following algorithm describes the details. Algorithm 1 1. Start from n=1. 2. Apply the same input signal u(t) to the first order system with residence time a and call the output y2 (t). Call the output data from the real system y1 (t). Then minimize the following cost function and find the residence time ∫∞ (∫ ∞ ) t.y2 (t) dt t.y1 (t) dt 2 0 0 a = arg min ∫ ∞ − ∫∞ a 0 y2 (t) dt 0 y1 (t) dt  3. put µg11 = a. 4. n=n+1.  49  (4.22)  5. Minimize the following equation and find the new a (∫ ∞  ∫  ∞ n t n .y2 (t) dt 0∫ t .y1 (t) dt a =arg min − ∞ ∞ a 0 y2 (t) dt 0 y1 (t) dt ∫∞ ( ) ∫ ) n−1 ( ) u(τ ).τ n−i d τ n−1 n i 0∞ u(τ ).τ n−i d τ 2 n +∑ i!ai 0 ∫ ∞ −∑ µg1 ∫ ∞ 0 u(t) dt 0 u(t) dt i=0 i i=0 i 0∫  (4.23)  6. µgn1 = n!an . 7. Go to 4. It is possible to estimate the higher moments of an arbitrary transfer function from this algorithm and from the system’s input and output data. In the next section another interesting application of Theorem 1 is given. In this section an all-pole linear system with time delay is considered, and a method to calculate an upper limit and lower limit for the delay is proposed. The method uses only first and second moments of the transfer function to estimate the aforementioned upper and lower limit.  4.3  Estimation of an upper and lower limit for the time delay  First it is important to establish the relationship between the time delay and the first and second moments of a transfer function. Assume the system can be modeled as G(s) = e−T s H(s)  (4.24)  where H(s) is an arbitrary linear transfer function without a delay. The first and second moments of G(s) can be written as a function of delay and the  50  first and second moments of H(s). In other words  µg1 = −  G′ (s) G(s)  s=0  Te−T s H(s) − e−T s H ′ (s) = e−T s H(s) H ′ (s) =T− H(s) s=0  (4.25)  s=0  = T + µh1 and  µg2 = =  G′′ (s) G(s)  s=0 2 −T s T e H(s) − 2Te−T s H ′ (s) + e−T s H ′′ (s)  e−T s H(s) H ′′ (s) + = T 2 − 2T H(s) s=0 H(s)  s=0  H ′ (s)  =T  2  (4.26)  s=0  + 2T µh1 + µh2 .  It is possible to estimate the first moment of G1 (s) from the method given in the previous chapter. Therefore, it can be assumed that the estimate of µg1 is available from a set of data. To estimate the second moment of G1 (s) rather using a first e−T2 s α s+1 as G2 (S). µg21 = µg22 ; At the  order model, it is possible to use  In other words, it is desirable to  find α and T2 in such a way that  same time it is possible to find α  and T2 to make  µg11  =  µg12 .  In other words, the first and second moments of G1 (s)  and G2 (s) are going to be equal, and from there it is possible to estimate a unique T2 and α by minimizing the left-hand side of equation (4.1). Now the following equations are true for α and T found from the above process. The first moment of G1 (s) and G2 (s) are equal. In other words  µg11 = µg12 .  51  (4.27)  By substituting equation (4.25) in this equation, the following equation is obtained T + µh1 = T2 + α .  (4.28)  The second moment of G2 (s) can be calculated as  µg22 = T22 + 2T2 α + 2α 2 = (T2 + α )2 + α 2  (4.29)  From equation (4.28) we can write the last equation as  µg22 = (T + µh1 )2 + α 2 = T 2 + 2T µh1 + (µh1 )2 + α 2 .  (4.30)  From the process of finding T2 and α we know that µg21 = µg22 ; therefore, by combing this equation and equation (4.30) the following equation is obtained T 2 + 2T µh1 + (µh1 )2 + α 2 = µg21  (4.31)  and by substituting equation (4.26) in this equation we finally obtain T 2 + 2T µh1 + (µh1 )2 + α 2 = T 2 + 2T µh1 + µh2  α 2 + (µh1 )2 = µh2 .  (4.32)  This equation holds for any arbitrary system H(s). Now consider a special case of a stable all-real-pole transfer function H(s)  µh1 =  ∏n1  1 . pi s + 1  (4.33)  This system can be considered as a series of several first-order systems with the residence time of pi . Therefore, the first moment of H(s) (residence time) is equal to the summation of the residence time of all these first-order systems. In other words  n  µh1 = ∑ pi . 1  52  (4.34)  The second moment of H(s) can be calculated by taking a derivative two times. The first derivative of H(s) is n  H ′ (s) = − ∑ 1  pj . (p j s + 1) ∏n1 pi s + 1  (4.35)  This transfer function is a parallel and series combination of the first-order transfer functions; therefore, its residence time can be calculated without taking a derivative. Using the characteristics of the residence time, µh1′ can be calculated as n  n  1 n  1 n  µh1′ = − ∑ p j (p j + ∑ pi ) (4.36)  = − ∑ (p j ) − (∑ p j ) . 2  2  1  1  However, µh1′ = −µh2 and therefore n  n  1  1  µh2 = ∑ (p j )2 + (∑ p j )2 .  (4.37)  By substituting equation (4.34) into the last equation, we obtain n  µh2 = ∑ (p j )2 + (µh1 )2 .  (4.38)  1  By comparing equation (4.32) and equation (4.38) it is concluded that n  ∑ p2i = α 2 .  (4.39)  1  Also, these inequalities hold for any sequence of positive numbers pi n  n  n  i=1  i=1  i=1  ∑ p2i ≤ ( ∑ pi )2 ≤ n( ∑ p2i )  53  (4.40)  where equality happens for n = 1 (see Appendix A). By using this inequality and equations (4.39) (4.34) we have  α 2 ≤ (µh1 )2 ≤ 2α 2 √ α ≤ µh1 ≤ nα  (4.41)  Using this last inequality and equation (4.28) the following equation is resulted √ T2 + (1 − n)α ≤ T ≤ T2  (4.42)  Also, it is possible to calculate α and T2 from the first and second moments of G1 (s). From equations (4.29), (4.27), and (4.28) the following is obtained  α=  √  µg21 − (µg11 )2 √ T2 = µg11 − µg21 − (µg11 )2  (4.43)  Finally by substituting equation (4.43) in equation (4.42) the following result is achieved  √ √ √ µg11 − n µg21 − (µg11 )2 ≤ T ≤ µg11 − µg21 − (µg11 )2  (4.44)  Therefore, it is possible to find a bound on the delay by estimating the first and second moments of the system and by knowing its order. Example 4.1. Consider following system G(s) =  e−15s . (10s + 1)2 (6s + 1)(16s + 1)  (4.45)  The first and second moments of G(s) can be calculated as µg1 = 57 and µg2 = 3741. The order of this system is four, and from equation (4.44) the upper bound and the lower bound on the delay can be calculated as 35 and 12, respectively. As can be seen, the lower band is very close to the actual value of the delay. In order to calculate the lower band it is necessary to have the system order; however, the system order is not always available. With the exception of the right-hand side of inequality (4.40), all the steps in the 54  above procedure are still valid for the all-pole-stable transfer function with imaginary poles. The inequality (4.40) might or might not be valid for these systems. The validity of the inequality depends on the magnitude of the imaginary part of the imaginary pole and on the relative distance of the real part of the imaginary pole from the other poles. If the value of the imaginary part of this pole is small or the imaginary pole is far from the other poles, there is a good chance that this inequality remains valid. Example 4.2. Consider following system G(s) =  e−20s . ((10 + 4 j)s + 1)((10 − 4 j)s + 1)(12s + 1)(4s + 1)(26s + 1)  (4.46)  The first and second moments of G(s) are µg1 = 82 and µg2 = 7728, respectively. The lower bound and upper bound on T are 11 and 51. If in this example the imaginary part of the pole is increased to 30 (rather than 4), the expression in both parts of the inequality will be imaginary. Therefore, it is not possible to find a bound on the system delay. The process is much more complicated for a system that has zero in its structure, in which equation (4.37) does not have a simple form. Finding a bound for a system with an arbitrary structure can be a topic for future study as it is very useful to have such a bound for a system delay. As mentioned in Chapter 2 identifying the delay is one of the most challenging parts of linear system identification. Having a bound for delay makes the search space smaller and facilitates the system identification process. The focus of the next section is on the implementation of Algorithm 1. In this section, the errors associated with this algorithm and an approach to implement this algorithm are explained .  4.4  Results and discussion  To demonstrate how these methods can be applied in practice, an example similar to that used Chapter 3 is considered. The higher moments are calculated using 55  Algorithm 1 for this example, and the results are compared with the theoretical one. Example 4.3. Consider the system similar to the one used in Example 3.1 (all of the coefficients have been divided by 100 in this new system) G(s) =  e−7s (13s + 1) . (10.57s + 1)(23s + 1)(14s + 1)  (4.47)  By using equation (1.5), the first five moments of this system can be calculated as µg0 = 1, µg1 = 41.57, µg2 = 2395.78, µg3 = 182897.49, and µg4 = 17681000. It is more convenient to use normalized moments rather than the actual ones. The corresponding normalized moments for this transfer function are νg0 = 1, νg1 = 1,  νg2 = 1.38, νg3 = 2.54, and νg4 = 5.921 . It is possible to reconstruct the impulse response from these moments using the Laguerre method describe in section 1.3. The result is shown in Fig. 4.1. The main disadvantage of this method is it does not guarantee a reconstructed impulse response with positive values [32]. To find the higher moments it is necessary to follow Algorithm 1. To minimize the cost function in equation (4.23) the MATLAB function f minsearch has been used. The input signal consists of periodic signals with a period of 600s (about 10 times of residence time). Values that have been estimated from the algorithm for the first five moments are µˆ g0 = 1, µˆ g1 = 41.60, µˆ g2 = 2306.3, µˆ g3 = 196850, and µˆ g4 = 5903300. By comparing these moments with the theoretical ones, it is apparent that the first three moments are very close to the actual ones, but the fourth moment is almost two times higher than the theoretical one. The algorithm encounters very large numbers while calculating the fourth moment, and this results in an error in identification. One way to alleviate this problem is to use more numerically efficient algorithms to calculate the higher moments. Estimating these moments requires estimating large integrals and therefore using efficient algorithms is crucial. The optimization and integrals calculation algorithm used in this thesis are very preliminary, and there is much room for improvements. The main purpose of this example is to show the validity of the algorithm proposed in this chapter. 1 Please  see equation 1.6 for definition of normalized moments.  56  1 reconstructed actual  Impulse Response  0.5  0  −0.5  −1  −1.5 0  1  2  3 time (s)  4  5  6  Figure 4.1: Actual and reconstructed impulse response of the system in Example 4.3.  4.5  Conclusion  The theorem and algorithm proposed in this chapter are opening new possibilities in system identification. By accurately identifying the higher moments and reconstructing the transfer function, it is possible to estimate a linear transfer function. Moreover, the method is very useful because it does not require prior knowledge of the system structure. However, this method, as many other new theories, is far from perfect and can be improved significantly in the future. Also, the reconstruction of the actual function from its moments has not been discussed in this chapter. One main obstacle in this area is that it is necessary to have several higher mo57  ments before successful reconstruction (around 10 moments, see [25] for details). Accurately estimating a large number of moments is very difficult. However, in many mixing processes the model structure is known prior to identification. In these cases the first few moments can be used to estimate the parameters of the mixer. In other words, the value of higher moments is useful and can be used to extract different system parameters [32].  58  Chapter 5  Modeling of time-averaging micromixers using parallel channels and a back flow loop 5.1  Summary  This chapter presents a new method to design continuous mixers for reducing concentration or temperature fluctuations in the flow. Such fluctuations can be caused by the sequential addition of fluids to be mixed into a channel. The mixer is based on separation and reconnection of flow through separate parallel channels. The main contribution of this work which differentiates it from previous ones is that it establishes a mathematical framework in the frequency domain to calculate optimal mixer parameters for efficient mixing. In addition, it uses a backflow loop that permits true low-pass mixing behavior. The method is demonstrated through the design of device parameters and simulation of the mixing performance of a mixer at the microscale. The simulations show that a time-averaging mixer can be designed that has a similar performance as a conventional stirred tank mixer. The material in this chapter is presented as follows. In section 5.2 it is shown that splitting a stream and recombining the parallel flows can only attenuate fluctuations at some specific frequencies; however, in section 5.3.1 it is shown how  59  introducing recirculation enhances the ability to attenuate fluctuations over a wide range of frequencies. In section 5.3.2, the design parameters are optimized in order to attenuate low frequency variations. Finally, a mixer with the optimized design is shown to successfully emulate a stirred tank when the purpose of mixing is to reduce concentration variations along the channel.  5.2  Design of a parallel channel mixer  5.2.1 Mathematical model for mixing using parallel channels  Figure 5.1: Schematic of a mixer based on parallel channels First we establish the mathematical model describing the progression of concentration variation in the parallel plug flows shown in Fig. 5.1. Each channel can be approximately represented by a pure delay which depends on the flow rate, the length and the cross-section of the channel [28]. This assumption ignores the parabolic characteristic of the velocity profile in the channels. Therefore, this approximation represents the worst case because Taylor dispersion caused by the parabolic velocity profile would reduce concentration fluctuations along a channel combined with diffusion [7]. Each individual delay Ti can be calculated as Ti =  Wi hi Li , Qi  (5.1)  where Li is the length, Wi is the width, hi is the height, and Qi is the flow rate 60  in the channel i. The relationship between input concentration, cI (t), and output concentration, cOi (t), of each individual channel cOi (t) = cI (t − Ti )  (5.2)  is characterized by the delay Ti of that channel. Neglecting diffusion in the parallel delay channels, the total output concentration n  cO (t) = ∑ ai cI (t − Ti ),  (5.3)  i=1  depends on the fraction of the total input flow rate ai , entering the ith channel at the input concentration. The total sum of the flow rates in all of these channels should be equal to the total input flow rate; this condition for conservation of mass therefore sets a constraint on the ai as n  ∑ ai = 1 with ai =  i=1  Qi , Q  (5.4)  where Q is the total input flow. Equation (5.3) can also be reproduced in the frequency domain by applying the Fourier transform. Contrary to the time domain, the representation of concentration in the frequency domain needs physical interpretation. The magnitude of the concentration Fourier transform represents the rate at which the concentration changes at each frequency. For example if the concentration does not change at all, the Fourier representation is equal to zero at all frequencies other than zero. If the concentration changes sinusoidally with a constant period the frequency domain representation is also non-zero at that frequency. An arbitrary stream of material can have contribution from all frequencies. By passing the input stream through different processes the frequency characteristics of concentration variations change. For the mixer geometry shown in Fig. 5.1 it is more convenient to represent these changes in the frequency domain than in the time domain. Mixers are designed to introduce uniformity and reduce fluctuations that are present in the input stream. In the time domain this means, the stream coming out of the mixer should have significantly reduced temporal variability compared to the stream going into the 61  mixer. In the frequency domain this corresponds to low-pass filtering. Let us write equation (5.3) in the frequency domain as n  CO ( jω ) = ( ∑ ai e− jω Ti )CI ( jω ) i=1 n  n  i=1  i=1  (5.5)  = ( ∑ ai cos(ω Ti ) − j ∑ ai sin(ω Ti ))CI ( jω ). In order to include the effect of diffusion inside the individual channels, the corresponding low-pass behavior would need to be included with each coefficient ai in equation (5.5). This would complicate the following analysis, and it is not necessary for our discussion as will be shown below. From equation (5.5) the relation between the absolute value of input and output concentrations can be written as ||G( jω )||2 =  n n ||CO ( jω )||2 2 = ( a cos( ω T )) + ( ai sin(ω Ti ))2 i i ∑ ∑ ||CI ( jω )||2 i=1 i=1  (5.6)  where G( jω ) is the system transfer function. For better insight into the significance of equation (5.6) it is useful to find its maximum and minimum values. This equation can be rewritten as  n  n  ||G( jω )||2 = ∑ a2i cos(ω Ti )2 + 2 ∑ i=1 n  ∑  ai a j sin(ω T j )sin(ω Ti )  i=1 j=i+1  n  = ∑ a2i + 2 ∑ i=1  ai a j cos(ω T j )cos(ω Ti )  i=1 j=i+1 n n  + ∑ a2i sin(ω Ti )2 + 2 ∑ i=1 n  n  ∑  n  ∑  ai a j (cos(ω T j )cos(ω Ti )  (5.7)  i=1 j=i+1  + sin(ω T j )sin(ω Ti )) n  n  = ∑ a2i + 2 ∑ i=1  n  ∑  ai a j cos(ω (T j − Ti )).  i=1 j=i+1  For a given set of ai the maximum and minimum of this last equation are determined by its second term. The maximum is achieved if all of the cos(ω (T j − Ti )) terms are equal to 1 (for example at ω = 0). In this case the maximum value can 62  be calculated as n  n  max(||G( jω )||2 ) = ∑ a2i + 2 ∑ i=1 n  n  ∑  ai a j  i=1 j=i+1  = ( ∑ ai )2  (5.8)  i=1  = 1. ||G( jω )||2 is a positive function, and its minimum value is zero. This minimum value can be achieved at a specific frequency by proper design of the Ti for instance in such a way that some of cos( jω (Ti − T j )) are +1 and some −1. Then ||G( jω )||2 is equal to (a1 ± a2 . . . ± an−1 ± an )2 , and if the ai are chosen in such a way that this quantity is equal to zero then ||G( jω )|| will be at its minimum. In a simple example let n = 3 and assume it is desirable to have G( jω ) = 0 at  ω = 100 rad/s. This can be obtained by choosing a1 = a3 = 0.25, a2 = 0.5, T1 = π /100 s, T2 = 2π /100 s, and T3 = 3π /100 s. The frequency response of this plug flow filter is depicted in Fig. 5.2. This example illustrates how the delays in parallel channels can be designed to achieve notch filter behavior for concentration fluctuations. The design given in [53] gives experimental evidence of the feasibility of such a microfluidics filter.  5.2.2 Validation of the parallel channel model using CFD The mixing model derived for the parallel channels is verified through a CFD simulation. A mixer with simple geometry is designed for the above example and the performance of this mixer is investigated using the software Comsol Multiphysics. Unlike the theoretical plug flow model, diffusion has been considered and modeled in the CFD simulation. In Section 5.3.3 it will be shown how the values of the delays Ti and the portions of the flow in each channel ai are related to physical lengths and widths of the parallel channels. For the purpose of this example 100µ m is chosen for the total inlet width, 10µ m for the mixer channel height, 1mPa · s for the dynamic viscosity and 50kPa for the input pressure. The dimensions of the three channels in the example 63  0  magnitude (dB)  −50  −100  −150  −200 0  100  200  300 400 frequency (rad/s)  500  600  Figure 5.2: The frequency response of a notch filter designed in the example. in Section 5.2.1 accordingly are L=[36.2 51.2 62.7] mm W =[18.0 50.9 31.1] µ m To simulate the performance of this mixer a two step procedure was followed in Comsol. In the initial step the steady state Navier-Stokes equations are solved to obtain the flow field in the mixer. In the second step the mixer behavior in response to inlet concentration changes is studied. In this step the transient convection and diffusion equations are solved using the flow field information gathered in the first step. The frequency content of the outlet concentration indicates the frequency range that mixer is not able to remove. In the case of this example, the mixer should remove variations in the range of 100rad/s. The mixer response to a step change in the inlet concentration is shown in Fig. 5.3. The concentration varies at the outlet, and the values shown in this figure are the integral of concentration values across the outlet boundary at each time instant. As it can be seen in this figure the concentration increases smoothly. However, the simple plug flow model in Section 5.2.1 predicts jumps in concentration as the fluid 64  in each channel reaches the outlet, because the parabolic flow profile and diffusion have not been considered. The simplest correction to this model can be achieved by adding first-order behavior to the delay for each channel. This first-order system behavior will guarantee the smooth transient at the outlet. To estimate the time constant of this first-order model we can measure the time that it takes for the concentration to rise from zero and to 63% of its final value, by looking at the step response of each individual channel and calculating its time constant; the new transfer function can be estimated as 0.26 0.49 e− jω 0.024s + e− jω 0.053s 1 + jω 0.023s 1 + jω 0.0472s 0.25 + e− jω 0.1s . 1 + jω 0.0375s  G( jω ) =  (5.9)  To estimate the frequency response of the mixer a sine wave is added to the inlet concentration and the magnitude of the outlet concentration variation is then measured. The result is shown in Fig. 5.4. This figure validates the revised model in equation (5.9) where its prediction is very close to the frequency behavior simulated in CFD. One region that the revised model cannot predict correctly is the notch frequency. The model predicts it would be slightly higher than the value determined through CFD. Also, the attenuation at the notch frequency is much less than the theoretical value (zero) for plug flow because of the presence of diffusion and the parabolic flow profile. Both remove the possibility of exact cancelation of the concentration variations in the three channels. The simple plug flow model that ignores diffusion therefore predicts inferior mixer performance at all frequencies except around the notch frequency. Diffusion is reducing the concentration variations especially at higher frequencies, supporting the low-pass behavior of the designed mixer in the next section.  65  5  outlet concentration  4  3  2  1 CFD simulation plug flow model  0 0  0.05  0.1  0.15  0.2 0.25 time (s)  0.3  0.35  0.4  Figure 5.3: The step response of a notch filter designed in Section 5.2.1.  5.3  Design of a mixer with recirculation  5.3.1 Mathematical model for mixing with recirculation As mentioned previously the combination of flows from different channels alone does not achieve the desired low pass characteristic needed for mixing. To attenuate high frequency fluctuations, part of the output stream can be pumped back to the inlet (Fig. 5.5). Similar to the flow in the forward paths the simplest model assumes the flow in the recirculation path as plug flow, and it assumes that concentration changes across the channel will be negligible when the recirculating fluid combines with the fluid entering the mixer. The new transfer function Gr ( j ω ) =  (1 − r)G( jω ) 1 − re−Tr jω G( jω )  66  (5.10)  0  magnitude (dB)  −50  −100  −150  −200 0  plug flow modified plug flow CFD simulation  50  100 frequency (rad/s)  150  200  Figure 5.4: The frequency response of a notch filter designed assuming plug flow compared to the revised model with a parabolic profile and diffusion , and estimated from CFD simulations. can be found from the transfer function G( jω ) given in equation (5.6), the ratio, r, of recirculated flow rate to the total flow rate, and the time that it takes for the recirculated flow to reach the input Tr . To understand why the recirculation is effective, consider the extreme case where r is close to one. This means the material inside the mixer will be recirculated for a very long time. As a result the exiting material is a very well mixed substance. However, it is impractical to choose r = 1 as both the input flow rate and the output flow rate will need to be zero. To have a better understanding of Gr ( jω ), equation (5.6) is substituted in (5.10) to  67  Figure 5.5: Schematic of a mixer based on parallel channels and recirculation yield ||Gr ( jω )||2 = = =  (1−r)2 ((∑ni=1 ai cos(ω Ti ))2 +(∑ni=1 ai sin(ω Ti ))2 ) (1−r ∑ni=1 ai cos(ω (Tr +Ti )))2 +r2 (∑ni=1 ai sin(ω (Tr +Ti )))2 j (1−r)2 (∑ni=1 a2i +2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j −Ti ))) n n 2 1−2r ∑i=1 ai cos(ω (Tr +Ti ))+r (∑i=1 ai cos(ω (Tr +Ti )))2 +r2 (∑ni=1 ai sin(ω (Tr +Ti )))2 j (1−r)2 (∑ni=1 a2i +2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j −Ti ))) j n 2 1−2r ∑i=1 ai cos(ω (Tr +Ti ))+r (∑ni=1 a2i +2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j −Ti )))  (5.11)  .  For the case where the numerator is non zero (a numerator of zero would correspond to ||Gr ( jω )|| = 0), equation (5.11) becomes ||Gr ( jω )||2 =  (1 − r)2 1−2r ∑ni=1 ai cos(ω (Tr +Ti )) j n 2 ∑i=1 ai +2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j −Ti ))  + r2  .  (5.12)  For ||Gr ( jω )|| to have its maximum attenuation over a wide frequency range it is necessary to solve the following optimization problem for a given r ai , Ti , and Tr = Argmax min ω  1 − 2r ∑ni=1 ai cos(ω (Tr + Ti ))  , j ∑ni=1 a2i + 2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j − Ti ))  68  (5.13)  where the objective is to determine the optimum values for ai , Ti , and Tr that minimize ||Gr ( jω )||; however, the required frequency range over which ||Gr ( jω )|| needs to attenuate fluctuations will depend on the particular application. It is possible to show that for the special case of r = 0.5, equation 5.13 always shows a low pass characteristics independent of the other mixer parameter values. Substituting r = 0.5 in equation (5.12) yields ||Gr ( jω )||2 =  0.52 1−∑ni=1 ai cos(ω (Tr +Ti )) j n 2 ∑i=1 ai +2 ∑i̸= j ∑nj=1 ai a j cos(ω (T j −Ti ))  + 0.52  .  (5.14)  This transfer function is one at ω = 0, and it is smaller than one at other frequencies since the first expression in the denominator is always positive. It is easy to show that 1 − ∑ni=1 ai cos(ω (Tr + Ti )) is always positive; ∑ni=1 ai = 1 and cos(ω (Tr + Ti )) cannot be equal to +1 for all Ti′ s . This shows that at least for the special case of r = 0.5 the transfer function Gr ( jω ) always has low-pass filter behavior.  5.3.2 Optimization of the mixer transfer function with recirculation As mentioned in section 5.3.1 the optimal parameters can be determined according to the expression (5.13); however, in this section a different approach is taken to calculate the design parameters. The basic idea is to calculate the mixer parameters in such a way that the transfer function Gr ( jω ) shows a first-order low-pass characteristic, i.e. the plug flow mixer will behave like an ideal stirred tank in spite of its different geometry. This requirement can be formulated as the following optimization problem ai , Ti , and Tr = Arg min ||G f ( jω ) − Gr ( jω )||, 0<ω <ωh  (5.15)  where Gr ( jω ) is the transfer function in equation (5.10), ωh is the maximum frequency that the mixer is designed for, and G f ( jω ) is a first-order low-pass filter. The transfer function G f ( jω ) =  69  1 1 + j ωωl  (5.16)  of the first-order filter includes the cut-off frequency ωl . The following requirements were chosen without any loss of generality. The maximum frequency ωh is assumed to be thirty times larger than ωl , and it is assumed fluctuations above this frequency can be attenuated by diffusion. Also, the value of r = 0.5 is chosen as it has been shown in the previous section that there is a guaranteed low-pass behavior for this value of r; moreover, the mixer was designed in such a way to approximate a normalized low-pass filter with cut-off frequency, ωl , equal to 1 rad/s; this normalization facilitates adjusting design parameters through simple scaling for different desired criteria. The actual optimization problem then becomes ai , Ti , and Tr = Arg min  0<ω <30  1 0.5 ∑ni=1 ai e− jω Ti , ω − 1 + j ωl 1 − 0.5 ∑ni=1 ai e− jω (Ti +Tr )  sub ject to 0 < T1 < . . . < Ti < . . . < Tn ,  (5.17)  0 < Tr , 0 < ai , n  ∑ ai = 1.  i=1  This optimization problem can be solved using the MATLAB fmincon function. The optimization problem is not convex, and the solution depends on the chosen initial condition. To find the optimum values for the parameters more than one hundred random initial values are generated. It is obvious that many of these initial values will result in the same extremum point. The initial values for the flow portions (ai ) and the delays (Ti ) are chosen between zero and one and zero and ten respectively. The sets that result in the smallest value on the right hand side of equation (5.17) are retained; in the next step the sets that result in feasible physical implementation are chosen from these candidates (details of the process is given in the next section). There is also a trade-off between the accuracy of approximating a first-order transfer function and the complexity of the mixer when choosing the value for n. The number of channels n = 6 was found to be a good compromise 70  and the following coefficients were calculated based on this choice: a = [0.1352 0.1121 0.2445 0.2829 0.1741 0.0511] T = [0.1815 0.3614 0.5331 0.6849 0.8348 0.983] s Tr = 0.1052 s The sixth component of the flow ratio is 0.0511 which is very small. This shows that there is not much benefit in using more than six channels as in the current design one of the pipes already only carries 5% of the total flow. The magnitudes of Gr ( jω ) and G f ( jω ) are compared in Fig. 5.6. The figure shows that the two transfer functions are very similar to each other in the desired frequency range (0 < ω < 30 rad/s). The value of the right hand side of equation (5.17) for these coefficients is 210. This value is much smaller than the typical 105 achieved by substituting arbitrary ai , Ti , and Tr in this equation. To investigate the sensitivity of the design, delay values are changed by less than 5% and the frequency response of the new Gr is calculated. The result is shown in Fig. 5.7. The new transfer function behaves similarly at very low frequencies, but shows much worse performance at higher frequencies. Therefore, it is crucial to make delays and flow ratios as close as possible to the values calculated in this section while implementing the design.  5.3.3 Derivation of the mixer dimensions In this section a method is presented to convert the delays, Ti , and the flow portions, ai , resulting from the optimization in section 5.3.2 to the mixer geometry. The hydrodynamic resistance of the channels is used to calculate their dimensions. The channels in Fig. 5.1 are all parallel to each other; therefore, the total resistance of the forward path can be written as Rtotal =  1 ∑ni=1 R1i  71  ,  (5.18)  0  Gr  −5  Magnitude (dB)  −10 −15 −20 −25 −30  Gf  −35 −40 −45 0  20  40 60 Frequency response (rad/s)  80  100  Figure 5.6: The frequency response of a first-order filter representing a typical continuous stirred tank, G f , and the optimized plug flow mixer with recirculation, Gr . where Ri is the hydrodynamic resistance of the ith channel. In order to relate the resistance of each channel to its geometry it is assumed that the height of each channel is much smaller than its width. Then from White [49] the resistance of the straight microfluidics channel i is given by Ri =  12µ Li , Wi h3i  (5.19)  where µ is the viscosity of the fluid. Equation (5.18) now can be written as Rtotal =  1 Wi h3  ∑ni=1 12µ Li i  72  .  (5.20)  0 −5  Magnitude (dB)  −10 −15 −20 −25 −30 −35 −40 −45 0  Gr with perturbed coefficients Gf 20  40 60 Frequency response (rad/s)  80  100  Figure 5.7: The frequency response of a first-order filter representing a typical continuous stirred tank, G f , and the plug flow mixer with perturbed parameters, Gr . Also, Ti is given as Li Vi Li Ai = Qi  Ti =  73  (5.21)  where Vi is the flow velocity, and Ai is the area of the ith channel. Combining equations (5.21) and (5.4) yields Li Ai Q j Ti = Tj Qi L j A j LiWi hi Q j = Qi L j W j h j LiWi hi a j = ai L j W j h j a j LiWi = . ai L j W j  (5.22)  In deriving this last equation it is assumed that all channels have the same height hi =h j =h; this assumption is typically valid for microchannels due to their fabrication process. The pressure drop along a channel ∆P = Ri Qi  (5.23)  is given by its flow rate Qi and hydraulic resistance Ri . Combining equations (5.23) and (5.4) yields ai =  Qi Q  =  ∆P Ri ∆P Rtotal  =  Rtotal . Ri  74  (5.24)  Substituting equation (5.20) in (5.24) gives  ai = = =  Wi h3 12µ Li W h3 ∑nj=1 12µj L j  Wi Li ∑nj=1  (5.25)  Wj Lj  1 LW  ∑nj=1 LijWji  .  Combining equations (5.25) and (5.22) yields ai =  1 ai TiW 2 ∑nj=1 a j TjWj2 i  ,  (5.26)  From there the following n − 1 non linear equations can be extracted: n W2 Wi2 1 j ( − 1) = ∑ , ai Ti ai a T j̸=i j j  (5.27)  where the only unknowns are the channel widths, Wi ; solving for these unknowns requires one additional equation n  Wt = ∑ Wi ,  (5.28)  i=1  where Wt is the total width of all channels combined. Wt is chosen by the designer as it determines the scale of the mixer. The Li are then found by rewriting equation  75  (5.21) as LiWi h Qi LiWi h = Qai LiWi hRtotal = ∆Pai LiWi h = . W h3 ∑nj=1 12µj L j ∆Pai  Ti =  (5.29)  This will give the required n equations to determine the Li . These non-linear equations were treated as an optimization problem. The value of Wi and Li are found using fmincon MATLAB function to minimize the difference between two sides of equations (5.27), and (5.29), respectively. For the purpose of this chapter the mixer was designed for microfluidics applications, and values ∆P = 50 kPa, h = 10 µ m, µ = 10−3 Pa·s, Wt = 1 mm were chosen as design parameters. Using the calculated values for Ti and ai , the following values for Li and Wi were found L=[8700 12300 14900 16900 18700 20200] µ m W =[80 93 247 323 220 70] µ m to yield a cut-off frequency of 1 rad/s. In these calculations the effect of diffusion and the parabolic profile have not been considered. As a reference we also evaluated mixing in a single straight channel including both the parabolic velocity profile and diffusion in order to compare this to the performance of the mixer. To investigate this problem a CFD simulation is carried out and the length of the channel is chosen as a weighted average of the above six lengths. In other words, L=  ∑ni=1 LiWi ∑ni=1 Wi  = 1579µ m.  76  (5.30)  The width of the channel is chosen as the inlet width, the diffusion coefficient as 10−8 m2 /s and all other parameters were similar as the ones used to derive the dimension of the closed loop mixer. The value of the inlet concentration attenuation that has measured for this channel at ωh = 30 rad/s is −29.5dB which is much less than the expected attenuation of the mixer without diffusion 0.8 (or −1.9dB).  5.3.4 CFD simulation of the mixer with recirculation  Figure 5.8: The schematic of the open loop mixer in Comsol. The mixer designed in section 5.3.3 is simulated using Comsol 3.4. The same two step process mentioned in Section 5.2.2 also is followed for this simulation. However, the presence of the recirculation adds an extra degree of complexity to the simulation; to simplify the simulation, the flow field is studied in two different steps. In the first step the flow field in the parallel six channels with the dimensions as determined in Section 5.3.3 is found without considering any recirculation. These channels are implemented as shown in Fig. 5.8. The following values have been used as boundary conditions and simulation parameters: material density ρ =1000 kg/m3 ; material dynamic viscosity µ =0.001 Pa·s; inlet horizontal velocity u=0.015 m/s; diffusion coefficient D=10−8 m2 /s. It is obvious that the only way to connect several channels with different lengths and widths in parallel is to wind the longer ones. In order to achieve the same flow resistance ratios or delay ratios for all channels as determined in Section 5.3.1 and Section 5.3.2 the channel lengths were adjusted by trial and error, until the observed delay ratios in the CFD simulation were close to the ones in section 5.3.2. 77  The delays were determined by observing the concentration changes of the channel outlets and they may be associated with small inaccuracies. The trial an error process is necessary as the hydrodynamic resistance calculated in Section 5.3.3 is for straight channel. Therefore, the length and width values calculated in this section can only be used as initial design values for implemented non-straight channels. The determined lengths and widths are L=[17300 25627 37933 43248 45180 41993] µ m W =[80 93 247 323 220 70] µ m. The corresponding delays that we measured are T≈[0.6 1.1 1.7 2.1 2.5 3.7]s. In the next step the recirculation path is added. The recirculation path connects the ceiling of the outlet channel to the inlet channel and has the same length as the shortest forward channel. In reality a pump will be implemented in the design, and will provide the back pressure needed for the recirculation; however, this is very challenging in simulation as it is almost impossible to incorporate a real pump as it adds an immense degree of complexity to the simulation. Instead, we used the Slip velocity1 boundary condition provided in the Comsol 3.4 software. By applying this boundary condition the flow can be pushed back to the inlet by the slip velocity on the walls of the recirculation channel. This velocity was found through another trial and error process. As mentioned in section 5.3.1 in this design half of the flow is supposed to be recirculated to the inlet; therefore, the slip velocity was adjusted to Uwall =0.075 m/s along the entire length of the recirculation channel in order to achieve an identical flow rate in the recirculation channel as at the outlet. The delay that has been measured in this channel is 1 s. In the final step the inlet concentration is changed from zero to five and concentration variations are observed across the mixer. This step is very computationally intensive as the convection diffusion equations are solved in the transient mode. To 1 The  velocity difference between the wall and boundary fluid layer.  78  Figure 5.9: Concentration profile across the mixer at t=8 s. evaluate the performance of the mixer the time at which the outlet concentration reaches 86% of its final value is determined. This value is equivalent to two times the residence time for a stirred tank (or twice the time constant of any first order system). The concentration profile at different moments is shown in Figs. 5.9 and 5.10. Also, the outlet concentration is shown in Fig. 5.11. At t = 11.62, the outlet concentration is almost 86% of its final value. From this the mixer residence time is estimated to be 5.8 s. A first-order response for a system with this residence time is also depicted in Fig 5.11. The first-order model and the mixer show very similar behavior in response to a step change in their inlet concentration. However, the plug flow model predicts a faster rise time as it does not consider the parabolic profile and diffusion. The mixer was originally designed for a residence time around 3.5 s. This number is still close to the residence time of 5.8 s that is observed in CFD. The main reason for the difference may be the presence of the parabolic profile and diffusion in the simulation. As it is shown in section 5.2.2 these phenomena increase the amount of mixing and therefore the effective residence time is higher than expected.  79  Figure 5.10: Concentration profile across the mixer at t=18 s.  5.3.5 Validation of CFD results To validate CFD results described in the previous section two different measures are taken: mesh refinement and experimental validation. To be certain that the solution is mesh-independent a coarse mesh with 190,000 elements is generated. The flow field inside the mixer is then found using this preliminary mesh. In the next step the mesh is refined to 500,000 elements. The flow field is again calculated using this mesh. There are no differences in the flow field values calculated using these meshes. Indeed, the value of the flow rate in different channels shows less than a one percent. Therefore, it is concluded that the solution is not dependent on the mesh size. To validate the CFD results the mixer is fabricated without the recirculation channel. The design was fabricated using Polydimethylsiloxane (PDMS) channels. In the next step clear water and water labeled with a florescent dye were added to the inlet, and the fluorescence of the dye at the inlet and at the outlet was measured. The concentration of the dye was periodically changed at the inlet by alternating the type of fluid entering the system and the magnitude of inlet and outlet changes were measured at different frequencies. A snapshot of the inlet and outlet of the design is shown in Fig. 5.12. The same process was carried out in CFD and the magnitude ratio of outlet fluorescence to inlet fluorescence is compared in Fig. 5.13 at different inlet frequencies. As can be seen in Fig. 5.13 CFD results and 80  4.5  outlet concentration  4 3.5 3 2.5 2 1.5 1  plug flow model mixer first−order model  0.5 0 0  5  10 time (s)  15  Figure 5.11: Outlet concentration profile from plug flow, CFD simulation and its corresponding first order system with the same residence time. experimental results agree well with one another in particular at low frequencies. This confirms and validates the CFD simulation conducted in the previous section.  5.4  Results and discussion  Fig. 5.14 shows the block diagram used for a Simulink simulation of the mixer performance in response to a periodic square signal with the coefficients of the transfer function calculated in section 5.3.3 with the following coefficients a = [0.1352 0.1121 0.2445 0.2829 0.1741 0.0511] T = [0.1815 0.3614 0.5331 0.6849 0.8348 0.983] s Tr = 0.1052 s This simulation also includes a first-order model for an ideal stirred tank according to equation (5.16) for comparison. The input signal has a period of 10s, and the 81  first-order transfer function has a corner frequency of 1rad/s. Fig. 5.15 shows the input concentration (a) and the output concentrations from the first-order model and from the plug flow mixer (b). This example shows that applying a stream of material with changing concentration to these systems will result in a similar output for both mixers. According to Fig. 5.15, both systems are expected to have nearly identical transient behavior, and a plug flow mixer and a stirred tank can show comparable mixing performance.  5.5  Conclusions  In this work we have developed an approach to designing plug flow mixers to have near first-order behavior. We have simulated the mixing performance of a mixer with the designed parameters, consisting of six parallel streams. The mixing performance was very close to that of a stirred tank mixer, which was modeled as a first-order system. This shows that time-averaging mixers can substitute for stirred tanks whenever the aim of mixing is to decrease variations in input concentration such as the dosing system shown in Fig. 1.4b. In our model we ignored Taylor dispersion and diffusion, which would have contributed to enhance mixing further. However, in a preliminary step we developed a model to incorporate these phenomena when the step response of the channels are available. Although the optimized design in this chapter was demonstrated on a micro-mixer, the design approach presented in section 5.3.3 can be applied to other length scales. This will enable a new design space for mixers.  82  Figure 5.12: Snapshot of concentration profile at (a) inlet and (b) outlet of mixer for ω = 15 rad/s and 400 kPa input pressure.  83  Figure 5.13: Ratio of outlet fluorescence to inlet fluorescence in different inlet frequency response for CFD and experimental case.  84  Figure 5.14: The Simulink model for comparing a stirred tank model and plug flow mixer model.  85  1  0.5  0 0  10  20  30 40 time (s)  50  60  70  1 plug flow stirred tank 0.5  0 0  10  20  30 40 time (s)  50  60  70  Figure 5.15: The time response of a stirred tank and plug flow mixer, (a) input and (b) output concentration.  86  Chapter 6  Final conclusions and future works Although mixing phenomena may look deceivingly simple, the scientific description and qualification of their impact on a given process is fairly complicated; therefore, this thesis focused only on one aspect of mixing: reducing variabilities in the mixer inflow concentration. In Chapter 2 the problem of finding an appropriate model for a pulp chest is answered. The proposed model is able to quantify the overall performance of a laboratory-scaled pulp chest. Despite these encouraging results, the final goal of the project is to evaluate the mixing performance of an industrial chest, preferably in real time. Part of our progress to achieve this goal is published in [10]. In this paper the mixing performance of an industrial latency chest is evaluated and compared using computational fluid dynamics (CFD) and the black-box system identification method (the method described in Chapter 2). Although this paper demonstrates acceptable agreement between the two methods applied to evaluate this industrial chest, it is not accurate to generalize the results. This chest was not the best study case for mixing. As described in [10], the geometry of the chest is complicated and it was originally designed for latency removal from pulp suspension. It is therefore necessary to conduct additional experiments on different industrial pulp chests to validate the model described in Chapter 2; however, before conducting more experiments on industrial pulp chests it is necessary to describe 87  the differences between the lab and industrial experimental setup. Chapter 2 is built upon the previously developed lab-scaled experimental setup developed in our group [21]. The purpose of designing this experiment was to imitate real industrial chest conditions, and to observe the same flow behavior in the lab; unfortunately, it is not accurate to compare the results of this experiment with an industrial pulp chest without further considerations. In Section 2.2 the input signal in the lab-scaled experiment is produced by changing the concentration of the saline solution injected to the inlet stream. This changes the conductivity of the inlet stream, and by measuring the input and output conductivity signals it is possible to identify the model for the chest; on the other hand, in the real industrial chest it is very difficult to change and measure conductivity in the inlet and outlet streams. However, most industrial chests are equipped with an inlet dilution valve, and by manipulating the amount of water going into the chest it is possible to change the inlet consistency. At a glance, it seems the results of Chapter 2 are applicable to an industrial chest by using consistency data rather than conductivity data; however, this would only be true if the dynamic of the chest were not affected by consistency variations during the experiment. Therefore, it is necessary to determine the extent to which it is possible to change the inlet consistency without significantly changing the chest dynamic. As mentioned in section 1.1 the residence time and channeling in the agitated pulp chest are closely related to the ‘cavern’ volume, the well-mixed and well-agitated zone around the impeller [50]. The size of the cavern depends on several factors, one of which is yield stress of the fluid (τy ). The correlation derived by Amanullah et al. [4] predicts that cavern diameter (D) is a function of the inverse of yield stress; in other words, D2 ∝  1 . τy  (6.1)  This equation predicts that in the same chest the cavern for a fluid with higher yield stress is smaller than the fluid with lower yield stress; furthermore, the following power law between pulp consistency (Cm ) and yield stress has been derived in literature [15]  τy = aCm b  88  (6.2)  where a is very close to one and b varies between 2.4 to 4 depending on the specific pulp type. Combining equations (6.1) and (6.2) will give the following equation D2 ∝  1 . Cm b  (6.3)  Assuming the average value of three for b and having in mind that cavern volume (V ) is a function of cavern diameter cube, the following correlation is achieved V∝  1 Cm 4.5  .  (6.4)  The minimum amount of detectable changes in consistency is about 5%. Assuming a 5% decrease in consistency, the ratio of the new cavern volume (Vnew ) to the old cavern volume (Vold ) would be Vnew = (1.05)4.5Vold ≃ 1.25Vold .  (6.5)  As a result, a 5% change in the pulp consistency may lead to a 25% change in the cavern volume; therefore, the experiment itself can disrupt the chest dynamic and introduce error into the estimation of the mixing zone. In other words, the model introduced in Chapter 2 has more than a 20% error while estimating the residence time for a typical industrial chest; however, this error would be less significant if there were a well-developed cavern in the whole chest, in which case equation 6.1 is not accurate anymore and changes in the consistency would not affect the flow behavior. Therefore, the mixed volume of the chest is not only a function of the volume and the flow but also it is a strong function of consistency. In the future, this speculation must be verified with more data from industrial chests. Another aspect of this project that has not been answered is the mixer effect on the quality of the final paper product. In response to this question, a trial experiment was designed for the PAPRICAN paper machine. The goal of this experiment was to study the quality of the final paper product in relation to the amount of energy expended in the pulp mixers. Unfortunately, it was not possible to finish this experiment during the course of this project, but it is important to run similar experiments in the future to quantify the effects of mixing on the quality of the final  89  paper product. Chapter 3 introduces a new method for identifying residence time. The method does not require prior information about the structure of the system and only requires system input and output data. In this sense the method is very appealing, as one of the most difficult steps in system identification is finding the appropriate structure for the system; however, the method in its current form is not appropriate for online system identification. Developing a recursive form of the method is a potential topic for future research. It Chapter 4 the results of Chapter 3 is generalized in order to estimate the higher moments of any transfer function from input and output data (Chapter 4). The method is fairly fast and very easy to use. It has many potentials and can be considered a serious contender for current available identification methods. However, the optimization algorithm used this thesis is a standard MATLAB optimization algorithm and may not be the best choice for solving the optimization problem to find higher moments. Also, the integration method used in this thesis to calculate different integrals is very preliminary and more sophisticated algorithm can be used. Using more efficient numerical algorithms are crucial especially for calculating the moments higher than three which involves multiplication, summation and integration of very big numbers. Also, it is important to study the effect of input signal on accuracy of estimation and comparing the results by using different input signal type. As many other identification method the accuracy of the estimation can be affected by the type and frequency components of input signal; however, in general this method is very less sensitive to the input signal type than many other methods. Also, it works better with signals having long period of steady states. These sort of signals are encountered regularly in the process industry and the method can offer many benefits in this field. Chapter 5 employs a fresh perspective for solving the problem of reducing variabilities in the inlet stream. In this chapter a new mixer is introduced which achieves mixing by time averaging of the inlet stream. Time averaging in this mixer is accomplished by guiding the inlet stream into separate channels and recombining the outlet flow of the channels. Also, it is necessary that part of the outflow be recirculated to the inlet by using a pump. The main advantage of the mixer is that it does not have an agitator. The energy expenditure in this mixer is through the pressure 90  drop across the mixer to push the fluid through different channels, and also through the recirculating part of the outlet fluid to the inlet through a pump. The lump sum of energy consumption through these two processes could be potentially less than a stirred tank to achieve the same amount of mixing; on the other hand, the mixer has not yet been fabricated completely, but the initial results of the fabrication of the device are very promising. Complete fabrication of the mixer in different scales and evaluations of their performance are topics that should be pursued in the future.  91  Bibliography [1] M. Abramowitz and I. A. Stegun. Ch. 22 in Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. 9th edition, 1972. → pages 10 [2] A. Ainla, I. Zen, O. Orwar, and A. Jesorka. A microfluidic diluter based on pulse width flow modulation. Analytical chemistry, 81:5549–5556, 2009. → pages 13 [3] J. Alvarez and P. Gonzalez. Constructive control of continous polymer reactors. J. Proc. Cont., 17:463–476, 2007. → pages 8 [4] A. Amanullah, S. Hjorth, and A. Nienow. A new mathematical model to predict cavern diameters in highly shear thinning power law liquids using axial flow impellers. Chem. Eng. Sci., 53:455, 1998. → pages 5, 88 ˚ om and T. H˝agglund. Revisiting the ziegler-nichols step response [5] K. Astr˝ method for pid control. J. Proc. Cont., 14:635–650, 2004. → pages 8 [6] K. Astr¨om and T. Hagglund. PID controllers: theory, design, and tuning. Instrument Society of America, 2nd edition, 1995. → pages 2, 8, 9, 32 [7] D. A. Beard. Residence time theory. J. Appl. Phys., 89:4667–4669, 2001. → pages 60 [8] C. Bennington. Mixers and mixing in Pulp bleaching: principles and practice. Atlanta: TAPPI Press, 1996. → pages 3 [9] A. Bertsch, S. Heimgartner, P. Cousseau, and P. Renaud. Static micromixers based on large-scale industrial mixer geometry. Lab Chip, 1:56–60, 2001. → pages 13 [10] S. Bhattacharya, C. Gomez, A. Soltanzadeh, F. Taghipour, C. P. J. Bennington, and G. A. Dumont. Computational modeling of industrial pulp 92  stock chests. Canadian Journal of Chemical Engineering, 88:295–305, 2010. → pages 87 [11] S. Bjorklund and L. Ljung. A review of time-delay estimation techniques. Technical Report LiTH-ISY-R-2554, Department of Automatic Control, Lund Institute of Technology, 2003. → pages 9 [12] D. Bonnewijn. Study of disturbance propagation and damping in a paper machine approach system via omsim dynamic modeling. Master thesis, Department of Mathematcs, and Pulp and Paper Center, Universite Catholique De Louvain, and University of British Columbia, 1997. → pages 12, 15 [13] J. Brown. The dynamic behavior of paper mill stock chests. Southern Pulp Pap. Manf., pages 103–112, 1968. → pages 3 [14] Y. C. Chung, Y. L. Hsu, C. P. Jen, M. C. Lud, and Y. C. Lin. Design of passive mixers utilizing microfluidic self-circulation in the mixing chamber. Lab Chip, 4:70–77, 2004. → pages 16 [15] B. Dalpke and R. J. Kerekes. The influence of fibre properties on the apparent yield stress of flocculated pulp suspensions. Journal of pulp and paper science, 31:39–43, 2005. → pages 88 [16] G. A. Dumont, C. C. Zervos, and G. L. Pageau. Laguerre-based adaptive control of ph in an industrial bleach plant extraction stage. Automatica, 26: 781–787, 1990. → pages 10 [17] F. Ein-Mozaffari. Macroscale mixing and dynamic behavior of agitated pulp stock chests. Phd thesis, Department of Chemical and Bioligical Engineering, University of British Columbia, 2002. → pages 3, 12 [18] F. Ein-Mozaffari, G. Dumont, and C. Bennington. Performance and design of agitated stock chests. In In: Proc. PACWEST Conference, page 3B6, Whistler, Canada, 2001. → pages 20 [19] F. Ein-Mozaffari, L. Kammer, G. Dumont, and C. Bennington. Dynamic modelling of agitated pulp stock chests. TAPPI J., 2:13, 2003. → pages 4, 18, 27, 30 [20] F. Ein-Mozaffari, L. Kammer, G. Dumont, and C. Bennington. The effect of operating conditions and design parameters on the dynamic behavior of agitated pulp stock chests. Can. J. Chem. Eng., 82:154, 2004. → pages ix, 4, 5, 19, 28 93  [21] F. Ein-Mozaffari, C. Bennington, and G. Dumont. Suspension yield stress and the dynamic response of agitated pulp chests. Chem. Eng. Sci., 60:2399, 2005. → pages 3, 4, 18, 27, 30, 88 [22] C. G´omez. Numerical and experimental investigation of macro-scale mixing applied to pulp fibre suspensions. Phd thesis, Department of Chemical and Bioligical Engineering, University of British Columbia, 2009. → pages 4 [23] V. Hessel, H. L˝owe, and F. Sch˝onfeld. Micromixersa review on passive and active mixing principles. Chem. Eng. Sci., 60:2479–2501, 2005. → pages 13 [24] C. Jeffreson. Dynamic testing- a unification. Chemical Engineering Science, 25:1319–1329, 1970. → pages 10 ¨ ul, and D. Th´evenin. Techniques for the [25] V. John, I. Angelov, A. Onc¨ reconstruction of a distribution from a finite number of its moments. Chemical Engineering Science, 2007:2890–2904, 2007. → pages 10, 58 [26] L. Kammer, F. Ein-Mozaffari, G. Dumont, and C. Bennington. Identification of channeling and recirculation parameters of agitated pulp stock chests. J. Proc. Con., 15:31, 2005. → pages ix, 4, 5, 6, 7, 17, 18, 19, 20, 30 [27] J. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright. Convergence properties of the nelder-mead simplex method in low dimensions. SIAM Journal of Optimization, 9:112–147, 1998. → pages 40 [28] O. Levenspiel. The chemical reactor engineering omnibook. Oregon St Univ Bookstore, 1st edition, 1993. → pages 8, 60 [29] L. Ljung. System Identification Toolbox-For Use with MATLAB. 3th edition, 1995. → pages 24 [30] L. Ljung. System Identification - Theory For the User. PTR Prentice Hall, 2nd edition, 1999. → pages 11, 20 [31] E. B. Nauman. residence time thoery (100th anniversary). Ind. Eng. Chem. Res., 47:3752–3766, 2008. → pages 14 [32] E. B. Nauman and B. A. Buffham. Mixing in continuous flow systems. John Wiley & Sons, 1st edition, 1983. → pages 8, 9, 10, 15, 56, 58 [33] E. B. Nauman, D. Kothari, and K. Nigam. Static mixer to promote axial mixing. Transaction on Instituation of Chemical Engineering, 80:1–5, 2002. → pages 2, 13, 14, 15 94  [34] N.-T. Nguyen and X. Huang. Mixing in microchannels based on hydrodynamic focusing and time-interleaved segmentation: modelling and experiment. Lab Chip, 5:1320–1326, 2005. → pages 13 [35] N.-T. Nguyen and X. Huang. Modelling, fabrication and characterization of a polymeric micromixer based on sequential segmentation. Biomedical Microdevices, 8:133–139, 2006. → pages 13 [36] J. M. Ottino. The kinematics of mixing: stretching, chaose, and transport. Cambridge University Press, 1st edition, 1989. → pages 13 [37] J. M. Ottino and S. Wiggins. Introduction: Mixing in microfluidics. Phil. Trans. A, 362:923–935, 2004. → pages 13 [38] A. Pearson and C. Wuu. Decoupled delay estimation in the identification of differential delay systems. Automatica, 20:761, 1984. → pages 23 [39] A. Rad, W. Lo, and K. Tsang. Simultaneous online identification of rational dynamics and time delays: a correlation-based approach. IEEE Trans. on Control System Tech., 11:957, 2003. → pages 6 [40] X. Ren, A. Rad, P. Chan, and W. Lo. Online identification of continuous-time systems with unknown time delay. IEEE Trans. on Automatic control, 50:1418, 2005. → pages 23 [41] E. Reynolds, J. Gibbon, and D. Attwood. Smoothing quality variations in storage chests holding paper stock. Trans. Inst. Chem. Engrs 42a, page 13, 1964. → pages 3 [42] T. Soderstrom and P. Stoica. System identification. Prentice Hall, 1989. → pages 23 [43] J. Solomon, T. Elson, A. Nienow, and G. Pace. Cavern size in agitated fluids with a yield stress. Chem. Eng. Comm., 11:143, 1981. → pages 5 [44] A. Soltanzadeh, S. Bhattacharya, G. A. Dumont, and C. P. Bennington. Estimation of residence time and channeling in agitated pulp chests. Nord. Pulp Pap. Res. J., 24:66–71, 2009. → pages 17 [45] A. Soltanzadeh, C. P. J. Bennington, and G. A. Dumont. Direct estimation of residence time from input and output data. The Canadian journal of Chemical Engineering, 2011. doi:10.1002/cjce.20525. → pages 17 [46] B. Stoeber, D. Liepmann, and S. J. Muller. Strategy for active mixing in microdevices. Phy. Rev. E, 75:066314, 2007. → pages 13 95  [47] O. Walker and A. Cholette. Determination of optimum size and efficiency of stock chests, part i: The ideal chest. Pulp Pap. Mag. Canada, page 113, 1958. → pages 3 [48] D. White. Parameter estimation in time delay systems. Phd thesis, Princeton University, NJ, 1976. → pages 22 [49] F. White. Viscous fluid flow. McGraw-Hill Book Company, 1st edition, 1974. → pages 72 [50] K. Wichterle and O. Wein. Threshold of mixing of non-newtonian liquids. Int. Chem. Eng., 21:116, 1981. → pages 5, 88 [51] H. Xia, Z. Wang, Y. Koh, and K. May. A microfluidic mixer with self-excited turbulent fluid motion for wide viscosity ratio applications. Lab Chip, 10:1712–1716, 2010. → pages 13 [52] Y. Xie, F. Azizi, and C. H. Mastrangelo. Compact high-frequency mixing module for microfluidic chips. Proceedings of the 3rd IEEE Int. Conf. on Nano/Micro Engineered and Molecular Systems, 1:916–920, 2008. → pages 15 [53] Y. Xie, Y. Wang, L. Chen, and C. H. Mastrangelo. Fourier microfluids. Lab Chip, 8:779–785, 2008. → pages 15, 63 [54] D. Yackel. Pulp and Paper Agitation: the history, mechanics and process. Atlanta: TAPPI Press, 1990. → pages 3  96  Appendix A  Proof of inequality (4.40) In this appendix it is shown that the following inequality is true for any sequence of positive numbers n  n  n  i=1  i=1  i=1  ∑ p2i ≤ ( ∑ pi )2 ≤ n( ∑ p2i ).  (A.1)  Proof First we start from the inequality on the left. n  n  i=1 n  i=1 n  i=1  i=1 n n  ∑ p2i ≤ ( ∑ pi )2 n  ∑ p2i ≤ ∑ p2i + 2 ∑  0 ≤ 2∑  ∑  2  ∑  i=1 j=1,i̸= j  2pi p j  i=1 j=1,i̸= j  97  pi p j  (A.2)  which is true as pi are positive numbers. This proves the left-hand side of the inequality. Now for the right-hand side n  n  ( ∑ pi )2 ≤ n( ∑ p2i ) i=1  0≤  i=1 n n 2 (n − 1) p2i − 2 pi p j i=1 i=1 j=1,i̸= j n n (pi − p j )2 i=1 j=1,i̸= j  0≤∑  ∑  ∑ ∑  ∑  which is positive as it is the summation of several positive numbers.  98  (A.3)  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0072655/manifest

Comment

Related Items