— PART TWO TRANSIENT BEHAVIOR OF PROCESSES— CHAPTER 3 Laplace Transforms In Chapter 2 we developed a number of mathematical models that describe the dynamic operation of selected processes. Solving such models, that is, finding the output variables as functions of time for some change in the input variable(s), requires either analytical or numerical integration of the differential equations. Sometimes considerable effort is involved in obtaining the solutions. One important class of models includes systems described by linear differential equations. Such linear systems represent the starting point for many analytical techniques in process control, In this chapter we introduce a mathematical tool, the Laplace transform, that can significantly reduce the effort required to solve linear differential equation models analytically. A major benefit is that this transformation converts differential equations to algebraic equations, which can simplify the mathematical manipula- tions required to obtain a solution First we define the Laplace transform and show how it can be used to convert the elements of any linear differential equation—both the time derivatives and the inputs, which are functions of time—into a standard set of transforms that can be placed in a compact table. How a simple differential equation can be Laplace transformed by referring to this table is then illustrated. Since many practical problems involve functions that are not found in the standard tables, we next demonstrate how any transform can be expanded into functions that are found in the table. This technique, called partial fraction expansion, is the key part of a general technique that can be used to solve virtually any linear ordinary differential equation problem. Some important general properties of Laplace transforms are then presented, and, finally, we illustrate the use of these techniques for a practical modeling situation. 3.1 THE LAPLACE TRANSFORM OF REPRESENTATIVE FUNCTIONS The Laplace transform of a function f(s) is defined as Fos) = S900] = [ster (1) where F(s) is the symbol for the Laplace transform, f(1) is some function of time, and y is an operator, defined by the integral. The function f(z) must satisfy mild 4344 LAPLACE TRANSFORMS: conditions [1] which include being piecewise continuous for 0 < t < %; this re- quirement almost always holds for functions that are useful in modeling and control. When the integration is performed, the transform becomes a function of the Laplace transform variable s which is a complex variable. The inverse Laplace transform (2°) operates on the function F(s) and converts it to f(t). Notice that F(s) contains no information about f(t) for ¢ < 0. Hence f(t) = £-'{F(s)} is not defined for <0. One of the important properties of the Laplace transform and the inverse Laplace transform is that they are linear operators; a linear operator satisfies the general superposition principle: Hax(t) + by(t)) = a5(x(t)) + bS(y(t)) (3-2) where ‘¥ denotes a particular operation to be performed, such as differentiation or integration with respect to time. If ‘* = 2, then Eq. 3-2 would be 2(ax(t) + by(t)) = aX(s) + BY(s) (3-3) Therefore, the Laplace transform of a sum of functions is the sum of the individual Laplace transforms; in addition, multiplicative constants can be factored out of the operator, as shown in (3-3). In this book we are more concerned with operational aspects of Laplace trans- forms, that is, using them to obtain solutions of /inear differential equations. For more details on mathematical aspects of the Laplace transform, the text by Churchill [1] is recommended. Before we proceed to solution techniques, the application of Eq. 3-1 should be discussed. The Laplace transform can be calculated easily for most simple functions, as shown below. Constant Function. For f(t) = a (a constant), 2(a) GB-4) Step Function. The unit step function, defined as _fo <0 S@ = {? 120 (3-5) is an important input that is used frequently in process dynamics and control applications. The Laplace transform of the unit step function is similar to that obtained for the constant a above. ets(o) = + 6-6) If the step magnitude is a, the Laplace transform is a/s. The step function incor- porates the idea of initial time, zero time, or time zero; all of these terms are ways to refer to the time at which this function changes from 0 to 1. To avoid any ambiguity concerning the value of the step function at ¢ = 0 (it has two values), we hereafter consider f(t = 0) to be the value of the function approached from the positive side, t = 0*.3.1 The Laplace Transform of Representative Functions 45 Derivatives. The transform of a first derivative is important because such deriv- atives appear in linear differential equations. This transform is (dfdt) = f (dfidt)e-* dt (3-7) 0 Integrating by parts, e(dfide) [ floes dt + feo] (3-8) 0 0 sx(f) — (0) = sF(s) — f(0) G-9) where F(s) is the symbol for the Laplace transform of (1). Generally, the point at which we start keeping time for a solution is arbitrary. Model solutions are ordinarily obtained assuming that time starts (i.e., t = 0) at the moment the process model is first perturbed. For example, if the process initially is assumed to be at steady state and an input changes by a step function, time zero is taken to be the moment at which the input changes in magnitude. In many process modeling applications, functions are often defined so that they are zero at initial time, that is, att = 0, f(t) = 0, or f(0) = 0. In these cases, (3-9) reduces to ¢(dfidt) = sF(s). The Laplace transform for higher-order derivatives can be found using Eq. 3-9. To calculate 8[f"(A)], we define a new variable ( = df/df) such that af) _ . : : (2) (# ) = s(s) — (0) (3-10) From the definition of (= df/dt), (s) = sF(s) — f(0) (3-11) Substituting into Eq. 3-10, s[sF(s) — f(O)] — 6(0) (3-12) = s°F(s) — sf(0) — f') (3-13) where f’(0) denotes the value of dffdtat ¢ = 0, The Laplace transform for derivatives higher than second order can be found by the same procedure. An nth-order derivative, when transformed, yields a series of (n + 1) terms: (2) = s*F(s) = s"-¥9(0)= 5-290) = « = sf-DQ0) = f-D) where f‘(0) is the ith derivative evaluated at ¢ = 0. If n = 2, Eq. 3-13 is obtained (3-14) Exponential Functions. The Laplace transform of an exponential function is important because exponential functions appear in the solution to all linear differ- ential equations. For a negative exponential, e~*', with b > 0 le e-Me-t dt = li oO dt (3-15) fi B 1 ‘ 1 See - pele Ol = a5 G16)46 LAPLACE TRANSFORMS: The Laplace transform for b < 0 is unbounded if s < b; therefore, the real part of s must be restricted to be larger than ~6 for the integral to be finite. This condition is satisfied for all problems we consider in this book. Trigonometric Functions. In modeling processes and in studying control sys- tems, there are many other important time functions, such as the trigonometric functions, cos wt and sin w/, where @ is the frequency in radians per unit time. The Laplace transform of cos wf or sin w/ can be calculated using integration by parts. An alternative method is to use the Euler identity V=1 (3-17) : Ghul ate a cos wf = ——>——., = 2 2 and to apply (3-1). Since the Laplace transform of a sum of two functions is the sum of the Laplace transforms, L(cos wt) = eel") + Ae(e-Fe) (3-18) Using Eqs. 3-15 and 3-16 gives, X(cos wt) (3-19) Note that the use of imaginary variables above was merely a device to avoid integration by parts; imaginary numbers do not appear in the final result. To find (sin wt), we can use a similar approach. Table 3.1 lists some important Laplace transform pairs that occur in the solu tion of linear differential equations. For a more extensive list of transforms, see Ref. 2. Note that in all the transform cases derived above, F(s) is a ratio of polynomials in s, a so-called rational form. There are some important cases when nonpolynomial (nonrational) forms appear. One such case is discussed next. The Rectangular Pulse Function. An illustration of the rectangular pulse is shown in Fig. 3.1. The pulse has height h and width ¢,. This type of signal might be used to depict the opening and closing of a valve regulating flow into a tank The flow rate would be held at A for a duration of ¢,, units of time. The area under the curve in Fig. 3.1 could be interpreted as the amount of material delivered to the tank (=Az,,). Mathematically, the function f(t) is defined as follows: oe 0 f(Q=th O 0 only. Note that an exponential term in F(s) results. A special case of (3-22) is the unit rectangular pulse, where A = 1/ty. In this case, the area under the pulse is unity. Impulse Function. A limiting case of the unit rectangular pulse is the impulse, or Dirac delta function, which has the symbol 8(¢). This function is obtained when t, > 0 while keeping the area under the pulse constant. A pulse of infinite height and infinitesimal width but with a finite area of unity is obtained. Mathematically, this can be accomplished by substituting A = 1/t,, into (3-22); the Laplace transform of 8(t) is 9((0) = tim La -e (3-23) to by i) 0 tw Time, ¢ Figure 3.1. The rectangular pulse function.82 Solution of Diflerential Equations by Laplace Transform Techniques 49 Equation 3-23 is an indeterminate form that can be evaluated by application of L’'Hospital’s rule. (Also spelled L’Hopital) Taking derivatives with respect to 4, of both numerator and denominator, £(8(t) = Jim 1 (3-24) If the impulse magnitude (i.e., area t,,/) is taken to be equal to a rather than unity, then 9(a8(t)) = a (3-25) as given in Table 3.1. The unit impulse function may also be interpreted as the time derivative of the unit step function S(t). ‘A physical example of an impulse function might be given by rapid injection of dye or tracer into a fluid stream, where f(t) would correspond to the concen- tration or the flow rate of the tracer. This type of signal is sometimes used in process testing, for example, to obtain the residence time distribution of a piece of equip- ment, as illustrated in Section 3.5. The response of a process to a unit impulse is called its impulse response. 3.2 SOLUTION OF DIFFERENTIAL EQUATIONS BY LAPLACE TRANSFORM TECHNIQUES In the previous section we developed the techniques required to obtain the Laplace transform of each element in a linear ordinary differential equation. Table 3.1 lists important functions of time, including derivatives, and their Laplace transform equivalents. Since the Laplace transform converts any function f(t) to F(s) and the inverse Laplace transform converts F(s)back to f(t), the table provides a simple way to carry out these transformations. The procedure that we will use to solve a differential equation is quite simple: Laplace transform both sides of the differential equation, substituting values for the initial conditions in the derivative transforms. Rearrange the resulting algebraic equation, finding the transform of the dependent (output) variable. Find the inverse of the transformed output variable. The method is best illustrated by means of several examples. EXAMPLE 3.1 Solve the differential equation, 5a +4y=2 yQ)=1 (3-26) using Laplace transforms. Solution First take the Laplace transform of both sides of Eq. 3-26: dy v(s2+4y) = 92 3-27 ( dt ») (2) B-27) Using the principle of superposition, each term can be transformed individually: d ’ (s 2) + 9(4y) = £2) (3-28)50 LAPLACE TRANSFORMS «(s a) =5¢ (2) = 5(6¥(s) = 1) = 58¥(s)— 5 (3-29 2(4y) = 4e(y) = 4¥(s) (3-30) #@) =2 (331) Substitute the individual terms: S8¥(s) = 5 + 4¥(5) = 2 (332) Rearrange (3-32) and factor out ¥(s): Y(s)(5s + 4) = 54 2 (3-33) or ¥(s) = es (3-34) Having obtained an explicit relation for ¥(s), now take the inverse Laplace trans- form of Eq. 3-34: (3-35) 2£'[Y)] = £ | 5s +2 | ss + 4) The inverse Laplace transform of the right side of (3-35) can be found by using Table 3.1. First divide the numerator and denominator by 5 to put all factors in the s + b form corresponding to the table entries: s+ 0.4 lz + = 0, and b; = 0.4, the solution can be written immediately, y(t) = 0.5 + 0.5e-°8 (3-37) Note that in solving (3-26) both the forcing function (the constant 2 on the right side) and the initial condition have been incorporated easily and directly. As with any differential equation solution, (3-37) should be checked to make sure it satisfies the initial condition and that it satisfies the original differential equation." Next we apply the Laplace transform solution to a higher-order differential equation. ‘In some books, Eq. 3-37 would be written as y(t) = 0.5 + 0.5e"*8(0) to emphasize the fact that the time function e~* is defined only for = 0, that is, when S(f) = 1. However, we have already noted that Y(s) is defined only for 1 = 0; hence the solution y(0) = £- [¥(s)] clearly has meaning only for ¢ > 0. Inclusion of the unit step function is redundant,33 Partial Fraction Expansion 51 EXAMPLE 3.2 Solve the ordinary differential equation ay y dy 7p t get UT + by =1 (3-38) with initial conditions y(0) = y‘(0) = y"(0) = Solution Take Laplace transforms, term by term, using Table 3.1: ¢ (2) = s¥(s) de 2 (5 #2) = 6s°¥(s) e(11%) = £ (u “| = ls¥(s) £6y) = 6Y(s) : s) = 5 Rearranging and factoring out Y(s), we obtain ¥(s)(s? + 65? + 11s + 6) = t (3-39) ul ¥(s) = (3-40) s(6 + 6s? + Ty + 6) To invert (3-40) to find y(t), we must find a similar expression in Table 3.1. Unfortunately no formula in the table has a fourth-order polynomial in the de- nominator. This example will be continued later, after we develop the techniques necessary to generalize the method. ‘We conclude from this example that a general transform may not exactly match any of the entries in Table 3.1. In the case of higher-order differential equations, this problem will always arise because the order of the denominator polynomial (characteristic polynomial) of the transform is equal to the order of the original differential equation and no table entries are higher than third order in the de- nominator. Hence, the method used in Examples 3.1 and 3.2 must be made more general. It is simply not practical to expand the number of entries in the table ad infinitum. Rather, we choose to develop a method based on elementary transform building blocks. This procedure, called partial fraction expansion, is discussed in the next section. 3.3. PARTIAL FRACTION EXPANSION When a high-order denominator polynomial arises in a Laplace transform solution, we first factor it. Note that the denominator polynomial consists of terms arising from the differential equation (its characteristic polynomial) plus terms contributed by the inputs. The factors of the characteristic polynomial may have to be obtained52 LAPLACE TRANSFORMS: by finding the roots of the characteristic equation which is obtained by setting the characteristic polynomial equal to zero. The input factors usually are quite simple. Once the factors are obtained, the Laplace transform is then expanded into partial fractions. As an example, consider s+5 Y(s) = G+ D6rd (3-41) which can be expanded into the sum of two partial fractions: sts oO & GtD@O+) stl s+4 where a, and a> are unspecified coefficients that must satisfy Eq. 3-42, The ex pansion in (3-42) indicates that the original denominator polynomial has been factored into a product of first-order terms. In general, for every partial fraction expansion, there will be a unique set of the a; There are several methods for calculating the appropriate values of a and aa: Method 1. Multiply both sides of (3-42) by (s + 1)(s + 4) (3-42) s+5=aj(s +4) + als + 1) (3-43) Equating coefficients of each power of s, we have slay ta, =1 (3-4da) 5%: doy + a = 5 (3-44b) Solving for a, and a, simultaneously yields @, = 4.02 = 4. Method 2. Since Eq. 3-42 should hold for all values of s, we can specify two values of s and solve for the two constants: 0 = -lu, - a) (3-45a) Jo, + a (3-45b) fag = 4 Solving, Method 3. The fastest and most popular method is called the Heaviside expansion. In this method we multiply both sides of the equation by one of the denom- inator terms (s + 5,) and then set s —b,, which causes all terms except ‘one to be multiplied by zero. Multiplying Eq. 3-42 by s + 1 and then letting s = -1, gives s+5 of se 4 |e 3 Similarly, after multiplying by (s + 4) and letting s = —4, the expansion gives ao = a = os 3 Thus, the coefficients essentially can be found by simple calculations For a more general transform, where the factors are real and distinct (no complex or repeated factors appear), the following expansion formula can be used: Ms) _ __N(s) oi YS) = DG) To+5) meth (3-46)39 Partial Fraction Expansion 53 where D(s), an nth-order polynomial, is the denominator of the transform. D(s) is the characteristic polynomial. N(s), the numerator, has a maximum order of n — 1, The ith coefficient can be calculated using the Heaviside expansion Ns) a, = (s +b; 3-47: oO oon Alternatively, an expansion for real, distinct factors may be written as N' N‘ yo) = MO. VO (48) © Tas +1) In this case the coefficients are calculated by (3-49) Note that several entries in Table 3.1 follow the 75 + 1 format. We now are in position to use the methods of partial fraction expansion to complete the solution of Example 3.2. EXAMPLE 3.2 (CONTINUED) First factor the denominator of Eq. 3-40 into a product of first-order terms (n = 4 in Eq. 3-46). Simple factors, as in this case, rarely occur: s(s3 + 6s? + Ils + 6) = s(s + 1)(s + 2)(s + 3) (3-50) This result determines the four terms that will appear in the partial fraction ex- pansion, namely 1 a1 a ay : s) = a . 3-51 EG) 8 nena) 6G son se es The Heaviside expansion method gives a, = 4,0; = —}, a; = tay, = —4. Once the transform has been expanded into a sum of first-order terms, merely invert each term individually using Table 3.1: 2 s+3 o+ 1 1 1 1 Se et ae, 21 el — (') (44) +4: (4) 8 (4) he} hee denn (3-52) y(t) = &[Y(} = 8 (# 42, 12 1/6 ) gs ostl Equation 3-52 is thus the solution y(s) to the differential equation (3-38). The a,'s are simply the coefficients of the solution. Equation 3-52 also satisfies the54 LAPLACE TRANSFORMS: three initial conditions of the differential equation. The reader should verify this result. General Solution Procedure Having solved Example 3.2, we now can state a general procedure to solve ordinary differential equations using Laplace transforms. The procedure consists of four steps: Step 1. Take the Laplace transform of both sides of the differential equation. Step 2. Solve for ¥(s). If the expression for ¥(s) does not appear in Table 3.1, Step 3a. Factor the characteristic equation polynomial. Step 3b. Perform the partial fraction expansion. Step 4. Use the inverse Laplace transform relations to find y(1). Figure 3.2 illustrates the four steps for solving an ordinary differential equation; note that solution of the differential equation involves use of Laplace transforms as an intermediate step. Step 3 can be bypassed if the transform found in Step 2 matches an entry in Table 3.1. In carrying out Step 3, other cases that we have not covered can arise. Both repeated factors and complex factors require a modified partial fraction expansion procedure. Time Laplace ‘domain domain o.d.e. Step 1 Take Laplace Initial transform conditions (able 3.0) Y Step 2 ro Ss _ Nie) ; Yu) = Me | | I y Step 3 | Factor D(a), I perform partial i fraction expansion | Step 4 Solution Take inverse y(t) Laplace transform (Table 3.1) Figure 3.2. The general procedure for solving an ordinary differential equation using Laplace transforms.3.3 Partial Fraction Expansion = 55. Repeated Factors If a denominator term s + 6 occurs r times in the denominator, r terms must be included in the expansion that incorporate the s + 6 factor a s+ 1 oe a, +233 b 5 a (3-53) in addition to the other factors. Repeated factors arise very infrequently in models of real systems. They are most often encountered in hypothetical models, for example, in the model of a cascade of identical stages. Because real system models exhibit slight physical differences from one stage to the next, the factors in real system cascade models will be only approximately equal, not truly repeated in the mathematical sense. EXAMPLE 3.3 For s+ ay mw oo ge 3-54) 0) Era) er GR oO G54) set up the partial fractions and evaluate their coefficients. Solution To find a in (3-54), the Heaviside rule cannot be used for multiplication by (s + 2), since 5 = —2 causes the second term on the right side to be unbounded, rather than 0 as desired. We therefore employ the Heaviside expansion method for the other two coefficients («, and a;) that can be evaluated normally, then solve for a, by arbitrarily selecting some other value of s. Multiplying (3-54) by (s + 2)° and letting s = —2 yields 1 a (3-55) Multiplying (3-54) by s and letting s = 0 yields stl 1 ae 3. see 4) 4 _. Substituting the value 5 = —1 in (3-54) yields (3-57) 1 a, = -= (3-58) An alternative approach to find a, is to use differentiation of the transform. Equa- tion 3-54 is multiplied by s(s + 2), SH1= ay(s + 2)s + as + a(s + 2)? (3-59) Then (3-59) is differentiated twice with respect to s, 0 = 2a, + 205; so thata; = —a; = —} (3-60)56 = LAPLACE TRANSFORMS: Differentiation in this case is tantamount to equating powers of s, as demonstrated earlier. The differentiation approach illustrated above can be used as the basis of a more general method to evaluate the coefficients of repeated factors. If the denominator polynomial D(s) contains the repeated factor (s + b)’, first form the quantity 265) = NO (6 + by = (5 + bymtay + (6 + BY to D(s) (3-61) + a, + (s + b)[other partial fractions) Setting s = ~b will generate a, directly. Differentiating Q(s) with respect to 8 and letting s = ~b generates a,_;. Successive differentiation a total of r — 1 times will generate all a, i = 1,2,..., r from which we obtain the general expression _ 1 d%QCs a, = aa ; i=0,...,7-1 (3-62) For i = 0 in (3-62), 0! is defined to be 1 and the zeroth derivative of Q(s) is defined to be simply Q(s) itself. Returning to the problem in Example 3.3, (3-63) from which (3-64a) (3-64b) Complex Factors Another important case arises when the factored characteristic polynomial yields terms of the form 2 +c st owhere Deca, Ss? + dis + dy 4 : In this case the denominator cannot be written as the product of two real factors. However, we can put it into factorable form by completing the square st + dis + dy =3.3 Partial Fraction Expansion 57 Hence, it is possible to rewrite the denominator in the form of two complex factors st t dys + dy = (s + b + jo) (s + b — jo) (3-66) where b = 4 and w= Note that the quantity under the radical is positive, according to the earlier as- sumption. The complex factors in Eq. 3-66 lead to a pair of complex terms in the partial fraction equation, as follows: a + JB), a + iB 2_ + other terms (3-67) ¥(s) = 0) = bt jet s+b— jo Note that one factor is the complex conjugate of the other. Appearance of these terms implies an oscillatory behavior: Terms of the form e~" sin wf and e~"' cos wt arise after combining the inverse transforms e~(*/*)' and e~(*~/»", The sine and cosine terms yield periodic, oscillatory behavior which ultimately damps to zero (e-* — 0) if b is positive. Although dealing with complex factors is more tedious than dealing with real factors, the Heaviside expansion (3-47) can still be employed. First we examine inherent requirements on the a,, B,, b, and w. If the two complex terms in (3-67) are combined into a single fraction, _ Ucn + a2) + f(Bi + Ba)Js , orb + Biw + arb — Bow 7 (s+ bY + @ * (s+ bP +o (3-68) (Bib — aw + ow + Brd) y G+ by +e . Since Y(s) is a transform of real quantities, the right side must also contain only real quantities, Hence By = — Bz and a, = a:. Equation 3-68 then simplifies to — 1 + JBi oi — JBi 7 u YS) = b+ ja ts hb jo* 6-69) The two complex numbers now combine to yield _ als +b) + 2Brw oD (5 + bY + w (3-70) When the individual terms in (3-69) are inverted, y(t) = aye 4 jBye-P-Fee on + ayel-Ptie — jBel-btiow eo. Rearranging gives y(t) = aye-M(e-Mt + elt) + jBye-M(e-™! = eM) +++ (3-72) With the use of the identities _ _ cos wt = —s and sine EMO) ayy Eq. 3-72 becomes y(t) = 2aye-' cos wt + 2B, sin wt + =~ (3-74)58 LAPLACE TRANSFORMS. Therefore, once a, and B; are evaluated using partial fraction decomposition, the inverse transform can be written out immediately. An alternative partial fraction form that avoids complex algebra is as follows: Let _a(s +b) ta . ¥(s) = (s + bY + w? _ Using Table 3.1, this fraction form inverts to y(t) = ae" cos wt + Be sin ot + +++ (3-76) However, the coefficients a, and a, must be found by simultaneous equations rather than by the Heaviside expansion, as is shown in Example 3.4. EXAMPLE 3.4 Find the inverse Laplace transform of stl YG) = SG aa 45) (3-77) Solution Using the quadratic formula, the factors of s? + 4s + 5S are found to be (s + 2+ j) and (s + 2 — j), so that Fe stl YO) - Sy aas Ses 2+ pe FP oD (3-78) The partial fraction expansion is, therefore, stl Cr S14 2 YS) = Sea ess tS First, use the Heaviside expansion to evaluate a3 and 83; multiply by s + 2 + j and lets = —2 — j: stl SF 2— Dinas _o22-j+1 _ -1-j | (2-72) B= 6 Rationalize the complex number, multiply (3-80) by (8 + 6)/(8 + 6j), to obtain -2- 14 100 a3 + JBs = (3-80) «3 + jBs = = 0.02 = 0.14) (3-81) —0.02 and By = —0.14, Now evaluate the repeated root using the formula (3-62): stl O0 = ya 45 ) stl a = O(s) ea Pras (3-83) 0.23.4 Other Laplace Transform Properties 59 2 - a oa en @ 2 a + ae = 5-0 = 0.04 Use Eq. 3-79 and Table 3.1 to obtain the corresponding time-domain expression: y(O) = 0.04 + 0.21 — 0.04e-* cos ¢ — 0.28% sin t (3-85) The alternative partial fraction form for (3-77), the form that completely avoids using complex factors, is stl s*(s? + ds + 5) Multiply both sides of Eq. 3-86 by s%(s? + 4s + 5) and collect terms: 5+ 1 = (a + 045)5? + (oy + ce + as? + (Sa, + 4oq)s + Se, (3-87) ass + a6 7 S +4545 (3-86) Equate coefficients of like powers of s: Say tas =0 (3-88a) 4oy + a2 + a, = 0 (3-88b) st: Sa; + 4a) = 1 (3-88c) ood (3-88d) Solving simultaneously, we obtain a, = 0.04, 0; = 0.2, a5 = —0.04, a, = —0.36. The inverse Laplace transform of Y(s) is 0.04 0.2 0.045 — 0.36 oe (es » 1 (02 e-1 . : y= 2+ (288) 5 e+ (22) 4 2 (EO) ay Before using Table 3.1, the denominator term (s? + 4s + 5) must be converted by completing the square to (s + 2)? + 1°; the numerator is —0.04(s + 9). Note that, to match the expressions in Table 3.1, the argument of the last term in (3-89) must be written as =0.04s = 0.36 _ =0.04(s + 2), __ = 0.28 ++i @+3+1 G+A+1 This procedure yields the same results as given in Eq. 3-85. ‘ s (3-90) 3.4 OTHER LAPLACE TRANSFORM PROPERTIES Final Value Theorem The asymptotic value of y(s) for large values of time y(%) can be found from (3-91) as long as lim y(t) exists, that is providing the limit exists for all Re(s) = 0: y() = lim[s¥)] (1) Equation 3-91 can be proved using the relation for the Laplace transform of a derivative (Eq. 3-9): iF 2 edt = s¥(s) ~ y(0) (3-92) Taking the limit as s > 0 and assuming that dy/dr is continuous and s¥(s) has a60 LAPLACE TRANSFORMS limit for all Re(s) = 0, we find [oa = lim[s¥(s)] — y(0) (3-93) a Integrating the left side and simplifying yields y(*) = lim [s¥(s)] (3-94) =o If the transform corresponds to a function of time that is unbounded for 1 », Eq. 3-94 will give erroneous results. For example, if ¥(s) = 1/(s — 5), Eq. 3-94 predicts y(~) = 0. Note that Eq. 3-9, which is the basis of (3-92), requires that lim y(¢— ~) exist. In this case y(t) = e*, which is unbounded for t > =. However, we can see that (3-92) does not apply here without having to evaluate (0), by noting that s¥(s) = s/(s — 5) does not have a limit for some real value of s = 0, in particular, for s = 5. Initial Value Theorem In analogy to the final value theorem, the initial value theorem can be stated as y(0) = lim[sY(s)] (3-95) The proof of this theorem is similar to the development in (3-92) through (3-94). It also requires continuous functions; however, the limit is taken as s > ~. The proof is left to the reader as an exercise. EXAMPLE 3.5 Apply the initial and final value theorems to the transform derived in Example 31: = y(s) = _ (3-34) Solution Initial value y(0) = lim [s¥(s)] = lim Bee =1 (3-96) Final value | 4s (3-97) a — The initial value of 1 corresponds to the initial condition given in Eq. 3-26. The final value of 0.5 agrees with the time-domain solution in Eq. 3-37. Both theorems are useful for checking mathematical errors that may occur in the course of ob- taining Laplace transform solutions.3.4 Other Laplace Transform Properties 61 Transform of an Integral Occasionally, it is necessary to find the Laplace transform of a function that is integrated with respect to time. In general, by application of the definition (Eq. #1), e{fneyae} [ { {near} ew a (3-98) 2 fe freee |, + ala f(t) dt (3-99) after integration by parts. The first term in (3-99) yields 0 when evaluated at both the upper and lower limits, as long as f(¢*) possesses a transform (is bounded). The integral in the second term obviously represents the definition of the Laplace transform of f(¢). Hence, {f rear} = 1 FG) (3-100) Note that Laplace transformation of an integral function of time leads to division of the transformed function by s. We have already seen that transformation of time derivatives leads to the reverse procedure, that is, multiplication of the transform bys. Time Delay (Translation in Time) Functions that exhibit time delay represent an important case in process modeling and control. Time delays are commonly encountered in process control problems because of the transport time required for a fluid to flow through piping. Consider the stirred-tank heater example presented in Chapter 1. Suppose one thermocouple is located at the outflow point of the stirred tank, and a second thermocouple is fo ° @) ae om 0 7 @) Figure 3.3. A time function with and without time delay. (The initial value has been subtracted from both functions.) (a) Original function (no delay); (b) function with delay.62 —_LAPLACE TRANSFORMS immersed in the fluid a short distance (10 m) downstream. The heater is off initially and, at time zero, it is turned on. If there is no fluid mixing in the pipe (the fluid is in plug flow) and if no heat losses occur from the pipe between the two ther- mocouples, the shape of the temperature response of the two sensors should be identical. However, the second sensor response will be translated in time, that is, it will exhibit a time delay. If the fluid velocity is 1 m/s, the time delay (t = L/v) is 10 s. If we denote f(t) as the transient temperature response at the first sensor and g(t) as the temperature response at the second sensor, Fig. 3.3 shows typical behavior of the functions f() (no delay) and g(t) (added delay of fa units). In Fig. 3.3 we have plotted only the response; in other words, we have subtracted the initial sensor output value. The function g(¢) is 0 for ¢< f. Therefore, we can state that BO) = Ft — &) Se — &) (3-101) which is the function f(t) delayed by ( time units. Here we include the unit step function S(t — f) to denote explicitly that g(t) = 0 for all values of ¢ < . If £ (f() = F(s). then © (800) = £ (FE ~ W)8(¢— 19) = [7 Fle ~ SUE ~ whe de [ise = (Oe dt + e fle = gen de 0 A 7 J “f(t = eM-We-Med(t — fp) (3-102) Since (t — t) is now the artificial variable of integration, which can be replaced by * £(g@) =e ie (Pye dt (3-103) lo yielding the Real Translation Theorem G(s) = £ (F(t ~ G)S(e — )) = e-*% F(s) (3-104) Consequently, the time delay introduces a nonrational (nonpolynomial) element into the transform In inverting a transform that contains an e~ element (time-delay term), the following procedure will yield results easily and also avoid the pitfalls of dealing with translated (shifted) time arguments. Starting with the Laplace transform ¥(s) = e“F(s) (3-105) 1. Invert F(s) in the usual manner, that is, perform partial fraction expansion, and so forth, to find f(t). 2. Find y(t) = f( — t) S(t — &) by replacing the t argument, wherever it appears in f(¢), by (¢ — ); then multiply the entire function by the shifted unit step function. EXAMPLE 3.6 Find the inverse transform of lie a. (3-106)3.5 A Transient Response Example 63 Solution Equation 3-106 first can be split into two terms ¥(s) = Yi(s) + ¥xls) (3-107) -— te © (4s + DGs +)” Gs + NGs +1) From Table 3.1, the inverse transform of Y;(s) can be obtained directly: DiC) ee ce (3-109) Since Y,(s) = e-?* ¥\(s), its inverse transform can be written out immediately by replacing ¢ by (¢ — 2) in (3-109), then multiplying by the shifted step function: yo(0) = fe 9 — e25]8(4 — 2) (3-110) (3-108) Hence, the complete inverse transform is y(t) = et = et a feta -P3]8(¢ — 2) (111) Equation 3-111 can be numerically evaluated without difficulty for particular values of, noting that the bracketed terms are multiplied by 0 (the value of the unit step function) when t < 2, and by 1 when f= 2. An equivalent and simpler method is to evaluate the contributions from the bracketed time functions only when the time arguments are nonnegative. An alternative way of writing Eq. 3-111 is as two relations, each one applicable over a particular interval of time: w= ee Osr<2 (3-112) and y(t) = et — et $ fet-24 — 9-203] = ed ee) = UL + 8) = 2,6487e-""4 — 2.9477 t= 2 (3-113) Note that (3-112) and (3-113) give equivalent results for ¢ = 2 since, in this case, y(t) is continuous at ¢ 3.5 A TRANSIENT RESPONSE EXAMPLE In Chapter 4 we will develop a standardized approach for the use of Laplace transforms in solving transient problems. That approach will unify the way process models are manipulated after transforming them, and it will further simplify the way initial conditions and inputs (forcing functions) are handled. However, we already have the tools to analyze an example of a transient response situation in some detail. The example below illustrates many features of Laplace transform methods in investigating the dynamic characteristics of a physical process. EXAMPLE 3.7 The Ideal Gas Company has two fixed-volume stirred-tank reactors connected in series as shown in Fig. 3.4. The three IGC engineers who are responsible for reactor operations—Larry, Moe, and Curly—are concerned about the adequacy of mixing in the two tanks and want to run a tracer test on the system to determine if dead zones and/or channeling exist in the reactors.64 = —_LAPLACE TRANSFORMS Stage 1 Stage 2 Figure 3.4. Two-stage stirred-tank reactor system, Their idea is to operate the reactors at a temperature low enough that reaction will not occur and to apply a rectangular pulse in the reactant concentration to the first stage for test purposes. In this way, available instrumentation on the second- stage outflow line can be used without modification to measure reactant (tracer) concentration. Before performing the test, the engineers would like to have a good idea of the results that should be expected if perfect mixing actually is accomplished in the reactors. A rectangular pulse input for the change in reactant concentration will be used with the restriction that the resulting output concentration changes must be large enough to be measured precisely. The process data and operating conditions required to model the reactor tracer test are given in Table 3.2. Figure 3.5 shows the proposed pulse change of 0.25 min duration that can be made while maintaining the total reactor input flow rate constant. As part of the theoretical solution, Larry, Moe, and Curly would like to know how closely the rectangular pulse response might be approximated by the system response to an impulse of equivalent magnitude. Based on all of these considerations, they need to obtain the following information: (a) The magnitude of an impulse input equivalent to the rectangular pulse of Fig 3.5. (b) The impulse and pulse responses of the reactant concentration leaving the first stage. (c) The impulse and pulse responses of the reactant concentration leaving the second stage. Solution The reactor model for a single-stage CSTR was given in Chapter 2 as v « = q(c; — 0) ~ Vke (2-34) Table 3.2 Two-Stage Stirred-Tank Reactor Process and Operating Data Volume of Stage 1 = 4m Volume of Stage 2 = 3m Total flow rate . 2 m/min Nominal reactant concentration in input . 1 kg mol/m>3.5 A Transient Response Example 65 0 0.25 Time (rin) Figure 3.5. Proposed input pulse in reactant concentration. where c is the reactant concentration. Since the reaction term can be neglected in this example (k = 0), the stages serve simply as continuous-flow mixers. Two simple material balance equations are required to model the two stages: de 4 + 2c, = 2¢, (3-114) dey 3S + 2a = 2e (3-115) If the system initially is at steady state (0) = (0) = ¢(0) = 1 kg mol/m$ (3-116) (a) The pulse input is described by 1 1<0 6 0s1<0.25 (3-117) 1 120.25 A convenient way to interpret (3-117) is as a constant input of 1 added to a rectangular pulse of height = 5: cf = 1 + S(rectangular pulse with height = 1 and width = 0.25) (3-118) The magnitude of an impulse input that is equivalent to the time-varying portion of (3-118) is kg mol M = 5——— x 0.25 min = 1.25 kg mol + min a 3 Note that M is simply the integral of the rectangular pulse. Therefore, the impulse input in concentration ¢? would be it) = 1 + 1.258(¢) (3-119) Although the units of M may have little physical meaning, we note that Eq. 3-114 is written with units of kg mol/min. In evaluating its right side, we see that 5m kg mol - min aM = 27 x 125 SE = 2.5 kg mol can be interpreted as the amount of additional reactant fed into the reactor in cither the rectangular pulse or the impulse.66 LAPLACE TRANSFORMS: (b) The impulse response of Stage 1 is obtained by Laplace transforming (3-114), using c,(0) = 4sC\(s) — 4(1) + 2C\(s) = 2C(s) (3-120) Rearranging we obtain C,(s): 2 a) = ia eet err G21) The transform of the impulse input in (3-119) is CXs) = 2 + 1.25 (3-122) Substituting (3-122) into (3-121), we have 2 6.5 2, _65_ s(4s +2) 4542 Equation 3-123 does not correspond exactly to any entries in Table 3.1. However, putting the denominator in 7s + 1 form yields CHS) = (3-123) 1 3.25 |. es OO) = Sasa * 3541 G-124) which can be directly inverted using the table to yield GX) = 1 =e"? + 1.625e-H = 1 + 0.6250 "2 (3-125) The rectangular pulse response is obtained in the same way. The transform of the input pulse (3-117) is given by (3-22) so that 1 5(1 — eo cr) = 4+ ’ @-126) Substituting (3-126) into (3-121) and solving for C¥(s) yields 4 R 1e-28 + _ He _ CW) = 49 + as +d ~ Sas +) @-127) Again, we have to put (3-127) into a form suitable for inversion oo : gen CUS) = F544 * Ss + ~ ss _ Before inverting (3-128), note that the term containing e~°2* will involve a translation in time. Utilizing the procedure discussed above, the following inverse transform is obtained h(t) = e+ OL = e-8) = SIL — eP2974S(¢ — 0.25) (3-129) Note that there are two solutions; for ¢ < 0.25 min (or f,) the right-most term, including the time delay, is zero in the time solution. Hence cA(t) =e"? + 61 — et?) = 6 — Se"? < 0.25 min (3-130) eh) = #24 (1 — e-#2) — SCL — e025) 1 = Sen? + Sete +0252 = 1 + 0.6657e-"? = 0.25 min (3-131),3.5 A Transient Response Example 67 Plots of (3-125), (3-130), and (3-131) are given in Fig. 3.6. Note that the rectangular pulse response approximates the impulse response fairly well for ¢ > 0.25 min. Obviously, the approximation cannot be very good before 1 = 0.25 because the full effect of the rectangular pulse is not felt until that time, whereas the full effect of the hypothetical impulse is felt immediately at ¢ = 0. (c) For the impulse response of Stage 2, Laplace transform (3-115), using (0) = 1 3sCx(s) — 3(1) + 2C,(s) = 2C,(s) (3-132) Rearrange to obtain C,(s) 3 3s + 2 5 For the input to (3-133), substitute the Laplace transform of the output from Stage 1, namely (3-124) CAs) = a6) 133) 3 a 5 OO) = 9 td l= +) *is+ | G34) which can be rearranged to a form suitable for inversion 3.25 cy) = = : (3-135) Tas +1 7° s5s + )Qs41) * (5s + DQs4+) Since each term in (3-135) appears as an entry in Table 3.1, partial fraction ex- pansion is not required. 3.25 le) 1 25e2) | 2c (3-136) For the rectangular pulse response of Stage 2, substitute the Laplace transform of =e at (1 5e-t5 — 20-12 AO = eS + [1 + gig se 2e | + = — = Impulse input Rectangular pulse input 1.50 1 (kg mol?) 1.25 Time (rin) Figure 3.6. Reactor Stage 1 response.68 LAPLACE TRANSFORMS, the appropriate stage output, Eq. 3-128, into Eq. 3-133 to obtain LS 2 (oe —=—Os Ds + 1) 6 5e-028 SS ee : s(.3s+DQs +1) s(.3s+D@s +1) C2 Again, the right-most term in (3-137) must be excluded from the inverted result or included, depending on whether ¢ < 0.25 min or not. The calculation of the inverse transform of (3-137) giv cht) = 6 + 1Se-"!S ~ 20" 1 < 0.25 (3-138) S(t) = 1 = 2.7204e-"5 + 2.663e-1? se (3-139) Plots of Eqs. 3-136, 3-138, and 3-139 are given in Fig. 3.7. The rectangular pulse response is virtually indistinguishable from the impulse response for ¢ larger than about 1 min. Hence, Larry, Moe, and Curly can use the simpler impulse response solution to compare with real data obtained from the reactor when forced by a rectangular pulse, as long as ¢= 1 min. The maximum expected value of c(t) is approximately 1.25 kg mol/m’, This value should be compared with the nominal concentration before and after the test ( = 1.0 kg mol/m?) to determine if the instrumentation is precise enough to record the change in concentration. If the change in concentration is too small, then the pulse amplitude or the pulse width or both must be increased. Since this system is linear, multiplying the pulse magnitude () by a factor of four would yield a maximum concentration of reactant in the second stage of about 2.0 (since the difference between initial and maximum concentration will be four times as large). On the other hand, the solutions obtained above strictly apply only for ¢, = 0.25 min. Hence, the effect of a fourfold increase in 1, can be predicted only by resolving the model response for ,, = 1 min, Qualitatively, we know that the maximum value of c, will increase with an increase in t,. Since the impulse response model is a reasonably good approximation with f, = 0.25 min, we might expect that small changes in the pulse width will yield an approximately proportional 13 = — = Impulse input Rectangular pulse input 12 | 2 (kg movi?) | | | 1 7 Tirme (ruin) Figure 3.7. Reactor Stage 2 response.Exercises 69 effect on the maximum concentration change. This argument is based on a pro- portional increase in the approximately equivalent impulse input. A quantitative verification is left as an exercise. SUMMARY In this chapter we have considered the application of Laplace transform techniques to solve linear differential equations. Although this material may be a review for some readers, an attempt has been made to concentrate on the important properties of the Laplace transform and its inverse, and to point out the techniques that make manipulation of transforms easier and less prone to error. The final example il- lustrates the use of a number of these techniques in a practical setting. The use of Laplace transforms to obtain solutions can be extended to models. consisting of multiple simultaneous differential and algebraic equations. However, before discussing further complications, it makes sense to discuss the definition and use of transfer functions. The conversion of differential equation models into transfer function models, covered in the next chapter, represents an important simplification in the methodology, one that can be exploited considerably in process modeling and control system design. REFERENCES 1, Churchill, R. V., Operational Mathematies, McGraw-Hill, New York, 1972 2. Nixon, F. E., Handbook of Laplace Transformation—Fundamental Applications, Tables and Ex amples, Prentice-Hall, Englewood Clifls, NJ, 1965. EXERCISES 3.1. Use Eq. 3-1 to show that the Laplace transform of (a) eM sinwt is Gupte b b) eb : ine -..«. 3.2. Find the Laplace transform of x(r) = cosh(at) sin(wt) by substituting for the hyperbolic function and then using Table 3.1. 3.3. An input function that is sometimes used to test the dynamics of processes is the so- called half-sine-wave pulse shown in the drawing. Here 0 1<0 x(t) = | Asin or _ 0 Fer @ Find X().70 LAPLACE TRANSFORMS: 3.4, Calculate the Laplace transforms of the graphical input signals in the accompanying figure 5 qa fo) fo, | | | a = 7 | oO 0 ° 2 6 0 2 6 erin) (rin) @ ®) 3.5. The start-up procedure for a batch reactor includes a heating step where the reactor temperature is gradually heated to the nominal operating temperature of 75°C. The desired temperature profile 7(1) is shown in the drawing, What is T(s)? Teo tkmin) 3.6. Using partial fraction expansion where required, find x(¢) for ' / s(s + 1) ooa4 OXO- Tae aery O*%O- GSP (b) Xs) = 4 @ x6) == (FQ + 3 + 4) +541 3.7. Expand each of the following s-domain functions into partial fractions: (s + 2y(5 + 3) FAG + 5G + 6) 1 [Gs + 1? + IP ( + 2) (©) ¥(s) = (d) ¥(s) = 3.8. Solve the following equation for (0) using Laplace transforms: [var = 22 yay = 1Exercises 71 3.9. For each of the following functions X(s), what can you say about x(t) (0 = 15 =) without solving for x(¢)? In other words, what are x(0) and x(%)? Is x(t) converging or diverging? Is x(t) smooth or oscillatory? _ tee @) X0) = Byes + EF 4H _ los? - 3 b) XO) = Bes $1 FD 165 + 5 0X60) - Fy Find the mathematical form of the solutions to the following equations, that is, de- termine the form of the time solution but do not numerically evaluate the coefficients. Be sure to obtain correct arguments of all exponential and periodic functions. Use the Final Value Theorem to obtain the solution for large values of time. In all cases x(0) = (0) = O and the dots denote differentiation with respect to time (a) €+4e+8r= 10 (©) H+ 4x = 10 (b) + 4+ 4x = 10 (ff) + 4 10 (c) ¥+4i + 3x = 10 (g) ¥- 4 + 3x = 10 @d) e+ 4e+ x= 10 Ch) &~ 4k + 4x = 10 3.11. Which solutions of the following equations will exhibit convergent behavior? Which oscillatory? dx i # dx a | 4=2=—CiCcidCzdiCzsCNW ax dx @ ataer4 Note: All of the above differential equations have one common factor in their char- acteristic equations. 3.12, For a system described by the following equations a oo Gtety=e* x0) ay 7 = Geytde=0 (0) =0 (a) Find X(s). Be sure to eliminate ¥(3) from this expression. (b) Determine the final value of x(t) as > ». (©) Find y() 3.13. Find the complete time-domain solutions for the following differential equations using Laplace transforms: ax dx(0) (@) Gat ae =e with x(0) = 0, 0, BO-9 dt de dx dt a@x | dx dx(0) dt (b) - 12x =sin3r x(0) = 0 | iéééié_—_ a + 2x, — (d) x0) = 0, (0) = 0 dx, Be x, + Day = sine72 LAPLACE TRANSFORMS 3.14, Find the solution to the following differential equation: fbx bx = sine x(0) = £(0) = 0 Use the final value theorem to determine x(1) as ¢ +» from X(s). What can you say about this result? Note: Dots denote differentiation with respect to time 3.15. Find the solution of dx — + 4x = f(t) di 0 1<0 where f()= jh O5¢<1/h 0 t= Wh x0) = 0 Plot the solution for values of ft = 1, 10, 100, and the limiting solution (/t > %) from 1 = 0 tor = 2. Put all plots on the same graph. 3.16. Three stirred tanks in series are used in a reactor train (see drawing). The flow rate into the system of some inert species is maintained constant while tracer tests are conducted. Assuming that mixing in each tank is perfect and volumes are constant: (a) Derive model expressions for the concentration of tracer leaving each tank. ¢; is the concentration of tracer entering the first tank. (b) If c;has been constant and equal to zero for a long time and an operator suddenly injects a large amount of tracer material in the inlet to tank 1, what will be the form of c3(t) (i.e., what kind of time functions will be involved) if 1.V=¥,=Vs3 2.V# V2 # Vs. (c) If the amount of tracer injected is unknown, is it possible to back-calculate the amount from experimental data? How? ‘ u | [= ta = oly 2 4 4 s vi Ha Vo ve 3.17. The system pictured is used to dilute a concentrated caustic solution. It is started up with pure water in the mix tank. (a) Derive a mathematical model to describe the concentration of caustic in the mix tank effluent c. You may assume that the caustic solution has approximately the same density as the water and that volume in the tank is constant. (b) Simplify the model for conditions which give rise to c, >> c. When will this occur? (©) Find the solution c(t) for both q,(t) and q-(t) constant with the assumption ee >>e Volumetric rate of caustic: 4.(0) Volumetric rate of water: q(t) Volume of stirred tank: VExercises: 73 Concentrated | caustic solution M Pure water : =| aw 3.18. An electrical heating element is being tested in a large mass of water at temperature T,. Determine how the element temperature 7, varies in time if the power dissipation rate in the element Q changes suddenly from 0 to ©. The element initially is at the water temperature 7,,. You may assume that the entire element is at the same tem- perature at any instant in time and that the convective heat transfer coefficient is constant. What is the relation between the unsteady-state model for this case and that of the stirred-tank heater described by Eqs. 2-10 and 2-11? 3.19. I. M. Appelpolscher has just returned home from a hard day at work and is looking forward to a cool, refreshing drink of his favorite beer, Old Froth and Slosh. Unfor- tunately, the only can available has been left sitting on the kitchen counter. Since sunlight from the kitchen window has been warming the beer alll afternoon, it is hardly fit to drink. He is interested in cooling the beer to the desired temperature of 60 °F as quickly as possible (a) Appelpolscher immediately springs into action by placing the can in the freezer section of the refrigerator, which has a temperature of —6 °F. If the initial tem- perature of the beer is 85 °F, estimate how long it will take to cool it to 60 °F. (You may assume that a can of beer weighs 0.8 Ib and that the heat capacity is ~1 Btu/lb °F). (b) Isthere a faster method of cooling the beer, using only the resources that normally would be found in Appelpolscher’s house? Justify your answer. 3.20, A stirred-tank reactor is operated with a feed mixture containing reactant A at a mass concentration cj, The feed flow rate is w,, as shown below in the drawing. Under certain conditions the system operates according to the model °74 LAPLACE TRANSFORMS = wea — pV ke, (a) For cases where the feed flow rate and feed concentration may vary and the volume is not fixed, simplify the model to one or more equations that do not contain product derivatives. The density may be assumed to be constant. Is the model in a satisfactory form for Laplace transform operations? Why or why not? (b) For the case where the feed flow rate has been steady at W, for some time, determine how c, changes with time if a step change in ca; is made from 4, to aa. List all assumptions necessary to solve the problem using Laplace transform techniques.