COMMUNICATIONS
An Uncoupling Technique for the Study of M uItivariabIe NonIinea r Ti me-De Iay Systems An uncoupling technique i s presented b y which relatively complex nonlinear multivariable systems with time delays can be decomposed into more manageable problems, of linear and nonlinear character, respectively. The linear part can be treated b y standard transform techniques and convenient approximations, while the nonlinear part of reduced dimensionality i s amenable to a separate analysis. A reactor-separator-recycle process i s offered as an illustrative example, and regional stability results are shown.
T h i s communication presents and illustrates a technique whereby the stability of multivariable systems described by linear and nonlinear differential-diff erence equations can be studied by uncoupling the linear and nonlinear equations and performing separate stability analyses on the two equation sets that result. It treats the dynamics of physical systems that may be represented by equations of the form
where .If
G(s)
=
Is - A
- C A j exp(-e,s)
(6)
j= 1
and H(sj
=
B
+
.M j=1
B, exp(-ejs)
(7)
The variables z(t) can be obtained from eq 5 by the inverse transformation where y, y(t - e,), and f are vectors of dimension S Jand the e, are the values of 111 time delaj-s. Equation 1 is composed of J nonlinear equations and K ( = S - J ) linear differentialdifference equatioiis. If those elements of y(fj which correspond to the linear f functions are renamed as z ( t ) and the remaining elements of y ( t j are renamed as x(tj, eq 1 can be rewritten as
and
z ( t ) = $-'[G(s)-'H(s)X(s)]
+ &-'[G(s)-
Z(011
(8)
whereupon substitution of eq 8 into eq 2 produces eq 4. The straightforward substitution ordinarily calls for a specification of X(sj ; however, this can be avoided by approximating [G(s)-'H(s)] = C(s)
(9)
in such a manner bhatan indirect method for inversion of C (s) . X(s) exists. When t'he linear system is stable, this transfer function matrix can be approximated by the summation of exponentials L
C(s)
(10)
Ci exp(-Pis)
= i=1
(3) This transformation has the advantage that eq 3 , being linear, can be solved in implicit form. Such a solution for z ( t ) can be substituted into eq 2 to reduce the order of the system from A ' to J , yielding the reduced nonlinear set dX dt
=
f[x(t), x(t
- ej)]
(4)
There are nom only J equations to be considered in the nonlinear analysis. This simplification is particularly attractive when it is used in conjunction with the approach of Razumikhin (195S, 1959, 1960) for direct stability determination, because the bounding technique suggested in Razumikhin's theorem (see Landis and Perlmutter, 1972) avoids the need for an explicit determination of the relationship between x(t) and ~ ( -t ej). To detail this manipulation, one considers that in eq 3 , x(t) and x(t - e,) are treated as forcing vectors for thez(tj response variables. Laplace transformation will then give Z(sj = G(sj-'H(sjX(s)
+ G(s)-'z(O)
476 Ind. Eng. Chem. Fundam., Vol. 12, No. 4, 1973
(5)
If in addition the initial condition vector z(0) proximation for eq 8 becomes
= 0,
the ap-
L
z(t)
=
$-1[C(s)X(s)]
=
CiX(t
- Pi)
(11)
i=l
If z(0j is nonzero and the G(s)-' is in simple form, then standard direct or indirect inversion procedures can be used to evaluat'e the term G(s)-'z(O). Otherwise, a function of the form of C(s) can also be used to approximate G (s) - I . The accuracy of the approximation eq 11 in describing the relationship between z(t) and x(t) depends on the range of frequencies contained in the forcing vector x(t) and on the number of terms in the series approximation (Landis, 1970). When the transfer function is a decreasing function of s over the whole range of s, the computational effort may be reduced by setting
c1 = cz
= ... =
CL
(12)
and (13)
Then
Heat Exchanger
[G-'(s)H(s)],=o = CiL
4 91
(14)
T,
which provides a means of directly evaluating the coefficients in the esponential aplrosimation in terms of the constants of the transfer function. If the transfer function for the linear equations lends itself to the simplified approximation in eq 12 and the C, can be calculated from eq 14, stability criteria can be formulated without the necessity of evaluating the p,. The practical application of this uncoupling technique is not restricted to these transfer functions since methods of evaluation of the coefficients and exponents in the approximation are available (Hamming, 1962) ; however, many chemical engineering process units, such as those considered in the following esample, are described by linear energy and material balances which are inherently stable and do in fact have transfer functions that decrease nioiiotonically with increasing s.
Ci
CSTR
\
Ta
Ci
-
v
T
Application to a Process Analysis
For illustration purposes, consider a hypothetical chemical engineering process consisting of a continuous flow stirredtank reactor (CSTR),a separator, and a heat exchanger in the recycle line (see Figure 1). Studies of the local stability of similar processes have been carried out by Bilous and Amundson (1955) and by Gilliland, et al. (1964); however, the effect of transport lag on system behavior was not considered. .Issuming with the usual restrictions that the system components can be modeled by simple mass and energy conservation statements yields for the reactor
A H V R - C ( T - T,)
dr
which is in the form of eq 2 and 3
with
(16)
for the separator H dCi - = (Pi dt
+ q)C(t - Si) - piCi - K(C1 - QCz) h dCz - = K ( C1 dt
(17)
and
- QC2) - PCZ
and for the recycle line with heat exchanger
ct
=
(A-) Cl(t - e,) + ( 1 Ca ) (20) Y
+ 91
Tl(t -
P
ez) exp(- Cd 1
+ 41
+ (-L) P + 91
The z(t) are the state variables for the separator, and the x(t) correspond to those for the reactor. Following the prescription of eq 6 and 7
["'
G-'H(s)
TO (21)
These equations call be written in ternis of dimensionless deviation variables and dimensionless parameters by using the definitions: xl = ( C ' , / A H ~ O ) (T T , ) ;xz = (c - cs)/CO; 71 = (Cp'AHCo)(Ti - Tis);zz = (Ci - Ci,)/Co; 23 = (C, C2s)/Co;r = ( R - Rs)V,qCo; and T = qt/ 8. By these substitutions, eq l%21 reduce t o
=
exp(jhs)/(s
-
ad
0 b d s - a331 exp(-Sls)lD bzza23 exp(-els)/ D
1
(28)
where D = (s - azz)(s - a33) - a32a23. The characteristic roots of 24-26 have negative real parts for all realistic values of the system parameters; hence with the initial conditions zero, the approximation equations (9) and (10) may be used to reduce the system of equations to the pair
Ind. Eng. Chem. Fundam., Vol. 12, No.
4, 1973 477
-
103.3
Ij-
*
x2
t0.1
0
-. 1-
t
+-O.' 0.134
for which stability is ensured is about & 6 O C while the maximum concentration disturbance is 0.1 mole/l. This elliptical boundary for the region of stability is also the boundary for the allowable initial curves for the system. The figure shows only a region of asymptotic stability for the two variables of the reactor, but since there exist five dependent variables in the process model, the region shown should be related to the three remaining state variables for the separator. If a disturbance were to occur within the reactor such that its extent were contained within the contour shown in Figure 2, the separator Fould have a bounded input disturbance and the entire system would be asymptotically stable. The uncoupling technique presented here is of particular importance in the regional stability analysis of nonlinear equations where a finite region of stability is often determined by search methods, and a reduction of the number of variables over which the search must be carried out is a great advantage. I n the fifth-order reactor-separator process the variable space searched was reduced from five variables to two. The uncoupling can be useful, however, in other applications, and is not limited to stability analysis alone.
0.354
Concentrotion, C ( m o l e s / l i t e r )
Figure 2. Region of stability for numerical example
Since the transfer functions in eq 28 are decreasing functions in s, the C, can be evaluated directly from eq 14. If one elects to follow Razumikhin, a region of stability may be established by applying the symmetric quadratic Liapunov function
Nomenclature
(30)
v = XTPX
The criterion for stability which results is a nonlinear algebraic inequality which is independent of the time delays. The details of the formulation and the numerical evaluation of the inequality constraint can be found in the work of Landis (1970). Numerical Example
Suppose the system parameters for a reaction-separationrecycle system are C, = 0.48 cal/cm3"K, V = 100 l., q1 = 1.67 l./min, Co = 2.2 g-moles/l., T o = 273"K, U = 1920 ca1j"K min, 9 = 4.0, p = 25 l./min, C1 = 0, K = 25 l./min, H = 50 l., and h = 50 1. The kinetics follow a first-order Xrrhenius-form rate equation
R
=
koC exp(- Q / T )
with ko = 2.5 X 10" min-l, Q = 8.25 X 103"K, and AH = 2.4 X l o 4cal/g-mole. For this set of constants there exist the following three equilibrium states State 1: T
=
T 1 = 92.2"C; C
=
0.134;
CI
=
0.537;
Cz = 0.106 mole/l. State 2: T
=
TI
=
50.2"C; C = 1.359;
C1
=
5.434;
Cz = 1.087 moles/l. State 3 : T
=
T I = 3.3"C; C
=
2.713; C1 = 10.852; Cz = 2.176 moles/l.
I t is desired to determine a region of stability about the first equilibrium state by applying the Razumikhin technique detailed by Landis and Perlmutter (1972). Figure 2 shows the region of asymptotic stability which results when the two eq 29 are analyzed via the Liapunov function of eq 30 for which the matrix P = [ 0.0250 0.0428
0.04281 0.0847
was obtained from local stability considerations. The maximum temperature disturbance from the equilibrium state 478
Ind. Eng. Chem. Fundam., Vol. 12, No.
4, 1973
A
coefficient matrix for linear equations -V(p ql)/q(H h )
a22
= = = = =
C
= concentration, moles/l.
Ci
=
all
+
+
-iVqdqH) - ( K V / p H ) a23 +KV/qH a32 KV/qh a33 = - i+KV/qh) - ,(V/h) B = coefficient matrix for linear equations bll = -all b22 = V(q qi)/qH
+
coefficient matrix for the exponential approximation
C,
= specific heat G(s) = transfer function defined by eq 5
H(s) = transfer function defined by eq 6 H , h = volumes of the two phases in the separator AH = heat of reaction I = identity matrix = mass transfer coefficient between phases (includes K area of contact) = frequency factor for first-order kinetics ko 2 = Laplace transform operator P = Liapunov coefficient matrix Q = activation energy divided by gas constant for Arrhenius-form kinetics q = inlet flow rate = recycle flow rate = reaction rate = deviation of reaction rate from the steady-stat,e value r T = temperature t = time = heat transfer coefficient for reactor U Lrl = heat transfer coefficient for heat exchanger V = reactor volume v = Liapunov function = nonlinear state vector of dimension J x = general state vector of dimension N y = linear state vector of dimension K z
%
GREEKLETTERS p i = exponent matrix for the exponential approximation o1 = time delay between reactor and separator o2 = time delay in recycle line = (pt/V), dimensionless time variable T 6 = equilibrium separator factor SUBSCRIPTS ASD SUPERSCRIPTS 0 = fresh feed condition i = reactor inlet conditions s = steady state 1 = recycle conditions
c
=
w
=
wall temperature of heat exchanger wall temperature of reactor
Literature Cited
Bilous, O., Amundson, N. R., A.I.Ch.E. J . 1, 513 (1955). Gilliland, E. R., Gould, L. A . , Boyle, T. J., JACC Proc. 5 , 140 ( 1964). Hamming, R. W., “Xumerical Methods for Scientists and Engineers,” RZcGraw-Hill, New York, N. Y., 1962. LandiL, J. G . , Ph.D. Thesis, University of Pennsylvania, 1970. Landis, J. G., Perlmutter, D. D., A.I.Ch.E. J . 18, 380 (1972): Razuniikhin. B. S..Ph.D. Thesis, Institut Medhaniki Akademiva Nauk USgR, 1958.
Razumikhin, B. S., A p p l . Math. Mech. 22, 215 (1959). Razumikhin, B. S., Automat. Remote Contr. 21, 515 (1960). J . G. LANDIS DANIEL D. P E R L N U T T E R * Department of Chemical and Biochemical Engineering L‘niversity of Pennsylvania Philadelphia, Pa. 192 74 RECIXVLD for review XTay 1.5, 1972 ACCEPTLD J d y 31, 1973 This research was supported by a grant from the Sational Science Foundation.
Effect on Heat Transfer Due to a Particle in Motion through Thermal Boundary layer over a Flat Plate A mechanism i s outlined to explain the heat transfer due to movement of a particle through the thermal boundary layer on a flat plate. It consists of the experimental determination of the improvement in heat transfer due to film scraping and b y particle convection. The heat picked up b y the particle was determined experimentally; a quantitative explanation of its variation i s given. The experimental values of the heat transfer coefficient for a moving spherical particle have been found to support literature relationships.
T h e study of the effect on heat transfer due t’o movement of a particle through the thermal boundary layer on a flat plate has a significant practical application in many fields of engineering. The most important of these is the field of fluidized bed heat transfer. The heat transfer coefficient in a fluidized bed is many times greater than the coefficient from wall to fluid alone. There are many reports in the literature explaining this fact (Botterill and Williams, 1963; Levenspiel and Kalton, 1954; Nickley and Fairbanks, 1955; Van Heerden, et al., 1951; TVicke and Fetting, 1966). According to Zabrodsky (1966), the thickness of the laminar sublayer is reduced anywhere from one-sixth of a particle diameter to a n arbitrary diameter due to particle movement; this increases the heat transfer coefficient. The particles themselves gain heat from their surroundings as they approach the heated surface, and when they leave the thermal region they give up their heat to the bulk. This is termed “particle convection heat transfer.” The effect of such transfer as a plausible mechanism was undertaken by Kashmund and Smith (1965), Ziegler and Brazelton (1964), and others. Though these authors have demonstrated the simultaneous presence of the two effects experimentally, their relative importance has not yet been investigated. It is the intention of this paper to determine experimentally the improvement in heat transfer due to film scraping and the heat picked up b y a particle when it moves through the thermal sublayer. Many experimental studies have been undertaken with regard to the motion of a particle in a fluidized bed. Some research workers have made qualitative statenients about the motion being raiidom and others state that it is not random; both views are supported by experimental evidence. Handley’s (1957) photographic study with glass beads in a liquid fluidized bed shows that the particles move smoothly throughout the bed and as such they do not have a finite residence time a t the bed wall. .Is suggested b y Mori aiid Wakhumara (19681,
when the particles are in downward motion along the wall, it is thought that’ they may pass through the fluid film adjacent to the wall. I n doing so, the particles pick up heat aiid also decrease the thickness of the thermal boundary layer. To isolate the effects of the particle heat conv’ection and the scraping of the film adjacent to the wall, the experimental setup we have used involves a single spherical particle moving in a circular path. During the course of its motion it passes through the boundary layer, adjacent to a heated platinum foil surface. Experimental Section
The details of the experimental setup are sketched in Figure 1. The particles, made of copper, were 0.4975, 0.3925, 0.2999,
0.29, and 0.25 cni iri diameter. A copper-constantan thermocouple was spot-welded to each of the spheres. The spheres themselves formed the ends of the thermocouple hot junction. The thermocouple was drawn through a groove in a silver steel rod of lja-iii. diameter and the copper-coristantaii wire was attached to two copper rings fitted in the rod. The rod rested in a bearing mounted in a pipe a t one end, and the other end was coupled to a worm and wheel arrangement. The worm could be rotated a t various speeds via a motor and set of reduction gears, thus achieving various velocities for the spherical particle. The platinum foil used as a heater was activated by a de electrical source. The position of the foil across which the particle moved measured 1.485 x 0.8275 x 0.00254 em. All the thermocouples were standardized a t the start of the experiment. The platinum foil was also calibrated to calculate its temperature from the resistance measurement. The voltage drop across the foil when there was no niovement of the particle was measured by a de microvoltmeter. Whenever the particle moved across the foil, its temperature fluctuations were recorded by measuring its voltage drop by an X-Y recorder. The particle temperature was measured by the thermoInd. Eng. Chem. Fundam., Vol. 12, No. 4, 1973
479
