Class Note for CHEE 502 with Professor Saez at UA
Class Note for CHEE 502 with Professor Saez at UA
Popular in Course
Popular in Department
This 39 page Class Notes was uploaded by an elite notetaker on Friday February 6, 2015. The Class Notes belongs to a course at University of Arizona taught by a professor in Fall. Since its upload, it has received 19 views.
Reviews for Class Note for CHEE 502 with Professor Saez at UA
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 02/06/15
CHEE 502 ADVANCED CHEMICAL ENGINEERING ANALYSIS CLASS NOTES FALL 2008 Eduardo saez INTRODUCTION MATHEMATICAL MODELING quotA mathematical model is a representation in mathematical terms of certain aspects of a nonmathematical systemquot Aris 1999 In the context of Chemical Engineering the nonmathematical system to which Aris refers would be a physicochemical process that might involve basic concepts from Thermodynamics Transport Phenomena Reaction Engineering and Process Design A mathematical model is developed to accomplish one or more of the following functions 1 Validate established principles or hypotheses 2 Quantify physicochemical parameters of the materials and process eg kinetic rate constants transport properties thermodynamic properties that may subsequently be used for other purposes 3 Predict how the process will behave under normal operating conditions 4 Predict the sensitivity of the process to changes in physical properties or operating parameters 5 Design and evaluate similar processes eg scale up 6 Design control systems for a specific process To develop a model the following steps are undertaken 1 Conceptual description of the process Usually this step involves the formulation of a quotmodelquot problem that we hope that can represent the real physical problem that we are trying to model This step usually involves formulating simplifying assumptions or hypothesis about the system s behavior 2 Establish the conservation principles to be used Find the appropriate version of the mass andor energy conservation equations If the equations are differential equations find the appropriate boundary andor initial conditions to be applied 3 Determine the physical properties and constitutive equations to be used eg ideal gas equation of state Fourier s law of heat conduction 4 Make plausible assumptions to simplify these equations as much as possible The assumptions must be veri ed either a priori or a posterior 5 Simplify the equations to form a coherent wellposed system 6 Solve the equations either analytically or numerically Advanced Chemical Engineering Analysis 7 Fall 2008 l 7 Before the solution is used the model must be validated Typically this involves nding a similar but simpler problem for which a solution is known and then simplifying the model to verify that the simpler problem is adequately described 8 Since the model never can be validated exactly for the problem considered at this point there are no assurances that the mathematical solution found is correct To increase the level of confidence in the model study the response of the model to variation in conditions and parameters ie a sensitivity analysis to be used for validation purposes 9 Apply the model The main focus of this course will be point 6 above but the rest of the steps in model development will be presented when relevant The course is organized in terms of specific models Each chapter will start with the formulation of a model and then the mathematical tools necessary to solve the model will be explored The following mathematical techniques will be explored Analytical methods 1 Separation of variables eigenfunction expansion 2 Combination of variables similarity transform 3 Laplace transformation 4 Moment technique 5 Perturbation methods and asymptotic analysis Numerical methods 1 Solution of nonlinear equations 2 Solution of initial value problems 3 Finite difference methods for partial differential equations Advanced Chemical Engineering Analysis Fall 2008 2 MODEL 1 The Adiabatic ContinuousFlow Stirred Tank Reactor CSTR Mathematical Aspects 1 Solution of nonlinear algebraic equations 2 Solution of simple ordinary differential equations 3 Numerical solution of initial value problems The mole balance on a CSTR Consider a chemical reaction taking place in a perfectly mixed tank with perfect thermal insulation so that there is no heat transfer between the reaction mixture and the surroundings A continuous ow rate of reactant liquid mixture is fed into the reactor and a continuous ow exits the reactor Figure 1 We will consider that the system under study is subjected to the following assumptions and conditions 1 The physical properties of the liquid phase density specific heat etc are assumed to be constant and uniform ie independent of temperature and liquid composition We expect this assumption to hold only when there are moderate temperature changes and for systems in which reactants and products are present at low concentrations 2 The reactor is perfectly mixed so that intensive properties of the exit ow are identical to intensive properties of the reacting mixture everywhere in the reactor 3 The volumetric flow rates of material entering and leaving the reactor are the same This will ensure a steady state with respect to total mass in the reactor since the densities are the same at inlet and outlet equal volumetric ow rate implies equal mass flow rate and the principle of mass conservation states that the total mass in the reactor must remain constant The pressure in the reactor will be considered constant also Q7 Tia cAi V T P CA Q volumetric ow rate cAi inlet species A molar concentration Ti inlet temperature CAI species A molar concentration in reactor V reacting mixture volume Q T P CA Figure 1 Continuouslystirred tank reactor CSTR and notation used in the mole balance of species A 4 The reactor is adiabatic all energy interactions between the reactor and its surroundings occur through inlet and outlet flows In addition we will consider that there is no viscous Advanced Chemical Engineering Analysis 7 Fall 2008 3 dissipation in the reactor ie no heating as a consequence of viscous friction To build a model for this process we will apply two conservation principles conservation of mass and conservation of energy The principle of mass conservation applied to a chemical species in the reactor states this is sometimes called a macroscopic mass balance The time rate of change of The net rate at which mass the total mass of a chem1cal of species enters the reactor spec1es 1n the reactor 1 The net rate at wh1ch mass of species is created by chemical reaction The terms on the righthand side of this equation will be positive whenever mass of species enters or is created by chemical reaction and negative whenever mass of species leaves or is destroyed by chemical reaction In reactive systems the species balance is expressed more commonly in terms of moles Let A be the chemical species for which we write the mass conservation principle To convert equation 1 to moles we divide both sides by the molecular weight of A and the balance becomes The time rate of change of the total moles of A in the reactor The net rate at which moles of A enter the reactor The net rate at which 2 moles of A are created by chemical reaction To express equation 2 in mathematical form we rst will define an intensive property to quantify the reaction rate Let RA be the rate of reaction of A representing the net number of moles of A that are created by chemical reaction per unit time and per unit volume of reaction mixture Equation 2 can be written as refer to Figure l for notation dVcA dt QcAi QGA RAV 3 In general the rate of reaction will depend on the concentrations of the various species that participate in the chemical reactions that occur in the mixture For example if species A undergoes an irreversible elementary decomposition A gt Products 4 then the reaction is rst order in A and the rate of reaction is RA k A 5 Advanced Chemical Engineering Analysis 7 Fall 2008 4 where k is the reaction rate constant which can be expressed as a function of temperature using Arrhenius equation k koe RT 6 where k0 is a constant Ea is the activation energy of the reaction R is the gas constant and T is the absolute temperature in the reactor First let us consider the case in which the temperature is constant then k will be constant and equation 3 can be solved analytically Knowing that this is a constant volume system equation 3 can be rearranged as follows d if kcAcAi 7 which will be subject to an initial condition 0A 0A0 F0 8 The quantity tRVQ is the residence time in the reactor and its use in equation 7 leads to dciAkicAicAi 9 dt tR tR Equation 9 can be integrated by direct separation and integration Alternatively the equation can be solved using the integrating factor technique which is applicable to first order ordinary differential equations of the form Ly dt ocy ft 10 The integrating factor technique is based on the use of the following identity It Lye zewdle0 tocy 11 dt dt Note that if we multiply equation 10 by e0 then the lefthand side of the resulting equation is precisely the righthand side of equation 1 1 Therefore equation 10 becomes dtw ewm 12 If yt is subject to the initial condition yy0 t0 l3 Advanced Chemical Engineering Analysis 7 Fall 2008 5 then equation 12 can be separated and integrated to get t yew y0 Je tftdt 14 0 which yields the solution t yt e o t yo Je0 tftdt 15 0 Applying this method to equation 9 leads to after manipulations t 7lktR 7 cA A1 c 0 cm e tR 16 1 ktR 1 ktR This equation predicts that a steady state will be reached at long times lim cA t 2 em cASS steady state 17 t gtltgto lktR In terms of this steady state concentration equation 16 can be rewritten as t 7lktR 7 0A mm 0A0 cAss 6 R 18 Figure 2 shows the behavior of this solution the concentration will increase or decrease monotonically with time depending on whether the initial concentration in the reactor is lower or higher than the steady state value respectively The steady state value will be approached asymptotically The energy balance on a CSTR Next we will consider that the reaction takes place in an adiabatic reactor In general this implies that the temperature will change with time Equation 9 still will be valid but it will contain 2 dependent variables cA and T from equation 6 To complete the formulation of the problem we need to apply the principle of energy conservation An energy balance in the reactor can be expressed as follows The time rate of change of 19 The net rate at which the total energy in the reactor energy enters the reactor The application of this energy balance to the reactor is presented in detail in Appendix A Here we will state the final equation obtained equation A29 Advanced Chemical Engineering Analysis 7 Fall 2008 6 0A 0A0 0 Ass cA0 t Figure 2 Qualitative behavior of the concentration of A with time when A undergoes a rst order consumption in an isothermal CSTR dT R pvcpaszcpm T VV7AAHI q 20 A where p and Cp are the density and heat capacity of the mixture respectively Cp has been assumed constant q is the net heat transfer rate to the reactor energy per unit time AHr is the enthalpy of reaction at T and P enthalpy change in the reaction per unit mole of A reacted defined as the difference between enthalpy of products and enthalpy of reactants and 1 if A is a reactant A 21 1 1f A 15 a product For the case of reaction 4 using the reaction rate 5 in equation 20 leads to dT pVCpaz pQCpTi T chAAHr q 22 The formulation of the problem of an adiabatic reactor with a firstorder reactions is obtained by substituting equation 6 into equations 9 and 22 letting q0 to get Advanced Chemical Engineering Analysis 7 Fall 2008 7 EJ M koe RT ti cA 2cAi 23 R R dT 7E4 pvcp a chpTi T Vk0e RT c AAHr 24 This is a system of two rstorder coupled nonlinear ODEs They are subject to the initial conditions 0A 0A0 F0 25 T To t0 26 Before we explore the solution of this problem we will formulate it in dimensionless form We define the following dimensionless variables x ci 27 Ai T 6 7 note the need to use absolute temperatures 28 i 17 L 29 tR 1 dlZl X OCXC 9 30 d1 1 d e ocsxe e 31 d1 XX0T0 32 eeor0 33 where 0LtRk0 34 Advanced Chemical Engineering Analysis 7 Fall 2008 8 AHrcAi 5 pCpTi 35 note that Blt0 for exothermic reaction and Bgt0 for endothermic reaction here we will assume that the enthalpy of reaction can be considered constant Ea 36 Y ITi on 6 0A1 To 6 2 38 0 TE The solution of the system of ODEs will yield the evolution of temperature and concentration of A in the reactor There is no known analytical solution to this problem although a relation between dimensionless concentration and temperature can be found analytically as shown in Appendix B Later in this chapter we will explore a numerical solution to this problem but for now let us focus on the steady state solution which will have practical significance as defining a continuous operating condition for the reactor If we assume that a steady state solution exists then this solution can be found by setting the time derivatives equal to zero in equations 30 and 31 39 40 This is a system with two nonlinear equations and two unknowns x and 6 that will represent the steadystate dimensionless concentration and temperature This system cannot be solved explicitly but it can be reduced to a single equation by solving equation 39 for x 1 X1oce l e 41 and substituting this into equation 40 to get a single equation with a single unknown 6 779 1 XBe 7 42 Hwy9 The solution of this equation can be found by a numerical method Any numerical method will require an iterative procedure that must start with an initial guess for the solution We will Advanced Chemical Engineering Analysis 7 Fall 2008 9 explore the solution of this equation using the NewtonRaphson method later in this chapter Another method to solve nonlinear equations like this one is through a graphical construction Even though this type of approach might be inaccurate and inefficient it allows one to analyze the solution from a physical perspective We will explore this option first The type of result obtained will depend on whether the reaction is exothermic or endothermic so we will treat the two cases separately 1 Endothermic reaction AHrgt0 3 Bgt0 The graphical solution will be found by plotting the functions on both sides of equation 42 vs 6 and finding the point at which both functions intersect Note that in this case both functions are positive since Bgt0 The two sides of equation 42 can be physically interpreted as follows The lefthand side Efe1 e 43 represents in dimensionless form the energy lost by the uid that goes through the reactor while the righthand side OCBe Y e E e r 1oce Y9 44 represents the thermal energy consumed by the reaction progress The steady state solution is the value of 6 for which ErEf To find the solution graphically we plot Ef vs 6 which will be a straight line of slope 1 and Er vs 6 which will be a monotonically increasing function that at large 6 tends asymptotically to XBIHX The solution is the intersection of the two functions This is illustrated in Figure 3 Clearly the shapes of the two functions assure that there will be a unique solution 2 Exothermic reaction AHrlt0 3 Blt0 In this case it is convenient to change the signs in the de nition of the functions that intersect in equation 42 Let Ef991 45 Now this function represents the energy gained by the uid that goes through the reactor while OCBe Y e E e r 1oce Y9 46 represents the thermal energy liberated by the reaction progress The type of graphical solution in this case depends on the values of the parameters X B and 7 Two different cases are possible as illustrated in Figure 4 The fact that three steady state solutions may be possible Figure 4b for an exothermic reaction complicates the modeling effort considerably Two important questions can be asked at this point 1 which steady state if any will be achieved for a particular operation of the reactor and 2 are all three steady states physically realizable The answer to question 1 will be based on the temporal evolution of the system which will be treated later The main factor to Advanced Chemical Engineering Analysis 7 Fall 2008 10 consider in answering question 2 is the stability of the steady state A steady state condition is stable if when the system experiences a perturbation that takes it away from the steady state it tends to evolve in time back to the original steady state Since all systems experience perturbations to their operating conditions unstable steady states will not be achieved in practice 1 Er9 1319 solution 9 Figure 3 Qualitative illustration of the graphical solution of equation 42 for an endothermic reaction Er9 l Ede 9 9 a b Figure 4 Qualitative illustration of the graphical solution of equation 42 for an exothermic reaction a single solution b multiple solutions Advanced Chemical Engineering Analysis 7 Fall 2008 l 1 To analyze the stability of the steady states in Figure 4b let us consider what happens to the system when the reactor temperature experiences a perturbation Ifthe system is displaced from a steady state then the energy balance must consider transient terms so we must go back to equation 31 In terms of the definitions in equations 45 and 46 equation 31 can be written as de iZEreEfe 47 d1 Take the case in which a perturbation increases slightly the temperature of the reactor Figure 5 A slight increase in the steady state temperature of 881 Figure 4b will take the system to a point at which EfgtEr P1 in Figure 5 According to equation 47 the energy balance dictates that the temperature should decrease with time from that perturbed point d6d17lt0 This will bring the system back to 81 We say that 81 is a stable steady state For SSZ on the other hand a slight increase in temperature will take the system to a point at which ErgtEf According to equation 47 this means that d6d17gt0 which implies that the temperature will increase with time and therefore instead of going back to SSZ the system will evolve towards SS3 We say that SSZ is an unstable steady state The same arguments show that SS3 is stable Also perturbations that decrease the temperature will have similar results In the case of SSZ a slight decrease in temperature will induce an evolution of the system towards 881 Hence even though SSZ is a solution of the steady state energy balance this point cannot be achieved in practice We will return to this point later in the chapter The graphical solution method is not the best way to solve a nonlinear equation such as equation 42 because of its relatively low accuracy Here we will consider a numerical method to solve nonlinear equations the NewtonRaphson method Solution of a single nonlinear equation the Newton Raphson NR method Consider that we are looking for the value of X that satisfies the equation fX0 48 where f is a nonlinear function of X that is continuous and that has a continuous derivative with respect to X The NewtonRaphson method is based on a Taylor series eXpansion of the function fX around a base point XX0 1 fXfXof39XoX Xo5fquotXoX Xo2 quot39 49 where the primes denote derivatives with respect to X The series can be truncated at any point by using Taylor39s theorem The NR method is based on a truncation of the series at the second order term According to Taylor39s theorem fX1 Xof39XoXXofquot XXo2 50 is an eXact representation of the function for a specific unknown value of E that lies between X0 Advanced Chemical Engineering Analysis 7 Fall 2008 12 and x 55 11 quotU 6 Figure 5 Perturbations to the steady state of an adiabatic reactor Figure 4 that increase slightly the reactor temperature P1 P2 and P3 are the perturbed states The arrows indicate the trend of the system evolution with time in each case The NR method is an iterative procedure in which equation 50 is applied at a base point xxi where i denotes the number of the iteration being performed the function is evaluated at xxi1 and the second term is truncated to produce the approximation fXi1 fXif39XiXi1Xi 51 The value xi is a known value of x for which fxi 0 Equation 51 is applied to nd a new value xi1 aiming at fxi10 Using this in equation 51 and rearranging leads to fXi f39Xi 52 Xi1Xi Because equation 51 is an approximation this equation does not guarantee that fxi10 However xi1 should be a better approximation to the solution than xi which means that successive application of equation 52 might eventually converge to the solution Graphically in a yx plane equation 52 is equivalent to extrapolating to y0 a straight line yyx that passes through the point xifxi and whose slope is fxi Figure 6 To estimate error propagation in the NR method let xs be an exact solution Advanced Chemical Engineering Analysis 7 Fall 2008 13 fxs0 53 y A YfXif XiX39Xi i solution Xli Xi1 gtX yfX Figure 6 Graphical representation of the NewtonRaphson method From a point xxi a better approximation to the solution xiH is found by extrapolating for this case in other cases this step may be an interpolation to y0 a straight line that coincides with the function fx at xxi and has the same derivative as the function at that point For a given step in the iterative procedure an estimate of the solution is the value xi Since equation 50 is exact we can evaluate it at xxs with x0xi and use equation 53 to get 0fXif39XiXs Xifquot xs Xi2 54 where E is between xi and xs We de ne the absolute error of the method at a given iteration by E lxs xi 55 Subtracting equation 51 from equation 54 recalling that we are letting fxi10 leads to 0f39ltxigtltxs xi1gtfquotlt gtltxs xigt2 56 from which we get fquot5 2fXi Xs Xi2 57 Xs Xi1 Advanced Chemical Engineering Analysis 7 Fall 2008 14 Taking the absolute value of this equation and using the definition of error in equation 55 leads to fquot5 2 2f xiE 58 i Ei1 2 This equation states that the error at a given iteration is directly proportional to the square of the error at the previous iteration This behavior is called quadratic convergence The success of the NR method in terms of the number of iterations required to find a solution and on the fact that a solution can be found at all will depend on how close the initial guess is to the solution For cases in which there are multiple solutions the NR method should be able to predict all the solutions from appropriate different choices of the initial guess Example Production of propylene glycol from propylene oxide In the presence of sulfuric acid which acts as a catalyst propylene oxide reacts with water to produce propylene glycol according to the reaction CH2 7CH7CH3 H20 gtCH2 7CH7CH3 59 O OH OH In excess water the kinetics of this reaction is rst order in propylene oxide A concentration Furthermore the reaction is exothermic so that multiple steady states are possible in an adiabatic CSTR Kinetic and thermodynamic data for this reaction are E 49058K k0321gtlt109 s39l R k k pCp 239507 AHr 102000 m3K kmol A Consider that the CSTR is operated at the following conditions 3 mfg Q35gtlt10 3n III S cAi 23 T5297 K v12 m3 The residence time in the reactor is 3 tR XL343s Q 35X10393m3s With the information given we can calculate the dimensionless constants from equations 34 to Advanced Chemical Engineering Analysis 7 Fall 2008 15 36 to get a1101gtlt10 520200 F3050 To solve this problem we apply the NR method to solve equation 42 First of all we rewrite the equation as OCBe Ye f e 9 1 106e 7e 60 so that equation 42 is equivalent to f60 Differentiation of equation 60 with respect to 6 yields after manipulations ocBer 9 1 f e 1 7 eY9oc2 62 61 Application of the NR method equation 52 for different initial values of 6 yields three different solutions Values obtained for each iteration are plotted in Figure 7 for the three different solutions Convergence is very fast in this case after 4 iterations the error is below 10399 for each case The dimensionless temperature for each solution and the dimensionless concentration calculated from equation 41 are shown in Table 1 Table 1 Three steady state solutions for the adiabatic conversion of propylene oxide to propylene glycol Steady state 6 x SS1 10209 0896 SS2 11000 0500 SS3 11647 0177 In this case as explained before SS2 is unstable but both SS1 and SS3 are physically achievable In principle SS3 would be the most attractive solution because of the high conversion achieved Note that there is a large difference in conversion between SS1 and SS3 It is evident that deciding which steady state will be achieved in a particular operation is of great practical significance In general the final steady state of processes with multiple steady states can be determined by performing a simulation of the temporal evolution of the process from a specific initial condition In this case this would involve the solution of equations 30 and 31 subject to initial conditions 32 and 33 Solution of this system of ODEs must be done numerically In what follows we explore numerical methods for the solution of this type of problem Numerical methods for rst order ODEs initial value problems Consider the following system of two ODEs with two dependent variables y1 and yz Advanced Chemical Engineering Analysis 7 Fall 2008 16 df1tay1y2 62 3T2f2ty1y2 63 12 f 118 7 116 7 3 114 7 112 7 11 7 391 E E 1 108 E 106 7 104 7 102 7 E U E 1 L 0 1 2 3 4 Number of iterations Figure 7 Convergence to the steadystate dimensionless temperature in the application of the NR method for the adiabatic conversion of propylene oxide to propylene glycol The three different steady state solutions are generated from three different initial estimates 1 108 and 12 subject to the initial conditions y1y170 t0 64 YZy20 F0 65 Since the boundary conditions are specified at the same boundary t0 this is an initial value problem IVP which defines a class of problems whose solution can be found by a certain type of numerical algorithms The numerical methods that we will explore for IVPs are based on the use of Taylor series expansions to represent the dependent variables and on the discretization of the solution domain into specific subdomains The methods fall into the category of what are called nite difference methods The simplest of these is Euler39s method Advanced Chemical Engineering Analysis 7 Fall 2008 17 1 Euler39s method Consider a Taylor series expansion of a function yt about a point tt0 We will use Taylor39s theorem to nd a rstorder expansion plus a residual term as we did in equation 50 2 YtYt0t0tt0 tt02 66 where E is between to and t Finite difference methods seek a solution of a differential equation not as a function of the independent variable but rather as discreet values of the function at speci c points in the solution domain For a solution yyt tgt0 we discretize the domain into an array of points For example these points could be equally spaced along the t axis as shown in Figure 8 The objective of the nite difference method is to nd approximate values of the dependent variable at the discretization points t t00 tlzh tzzzh t39 Figure 8 Discretization of the independent variable t into equallyspaced intervals nite differences of length h de ned by the discretization points to t1 Consider the Taylor series expansion 66 using as base point t0ti and evaluating the series at tzti1 d 1 12 Wm yltt1gtd lttigtltt 11 ti t 11 tigt2 67 Since we have hZ y11 y1 hy39ti 7yquot 68 where we are using yiyti The idea of the method is to use relatively small values for h We can see that as the value of h is decreased if yt is a function with bounded second derivatives the last term on the right hand side of equation 68 will be less and less important We say that this term is quotof order h2quot to represent that it will be proportional to h2 as haO This is denoted as Advanced Chemical Engineering Analysis 7 Fall 2008 18 h 2 7 5 0012 69 note that yquotE will tend to yquotti as h gt0 We then can write equation 68 as follows Yi1 yihy39ti0h2 70 Now consider that the function y satisfies the initial value problem dy 72f t 71 dt y yy0t0 72 Substituting y39 from equation 71 into equation 70 leads to Yi1 yihftiyi0h2 73 Euler39s method is based on truncating the secondorder term in this equation 7i1 7ihfti i 74 where the circum ex will be used to denote approximate value of the dependent variable The sequential application of equation 74 from i0 l 2 starting from 90 y0 the initial condition leads to calculation of approximate values of the solution at the discretization points Note that the values of 7i1 depend only on the value of the dependent variable in the previous discretization point yi This means that the calculation can be performed explicitly from equation 74 regardless of the complexity of the function f Because of this characteristic Euler39s method is an explicit method At this point we are ready to solve the transient mole and energy balances for the adiabatic CSTR First we generalize Euler39s method equation 73 to a system of two ODEs equations 62 and 63 The appropriate discretized equations are hm m hf1 01311323 i0 1 2m 75 71m 72 hf2ti 1i 2ia i0 1 2w 76 Direct application of these equations to equations 30 and 31 leads to i1 ih l xi chie ei i0l2 77 Advanced Chemical Engineering Analysis 7 Fall 2008 19 l i1 ih1 i oc3xie 9i i012 78 The evolution of the dimensionless concentration and temperature and eventually the steady state that is reached will depend on the initial conditions As an illustration of the use of Euler39s method consider the production of propylene glycol with initially no propylene oxide in the reactor and a reactor temperature equal to the feed temperature x00 79 601 80 It is expected that the solution found by the application of equations 77 and 78 from the initial conditions 79 and 80 converges to the exact solution at low enough values of h However an appropriate value of h cannot be defined a priori The usual procedure employed is to generate solutions for various decreasing values of h until the solution becomes independent of h For example Figure 9 shows the solution found for the dimensionless concentration as a function of dimensionless times for various values of h The solution for h002 is exact for the scale of the plot Note that the solution for h1 is quite inaccurate but h025 leads to a solution that is relatively close to the exact solution T Figure 9 Dimensionless concentration vs time curves generated by Euler39s method for various values of h for the production of ethylene glycol in an adiabatic CSTR with initial conditions 79 and 80 Advanced Chemical Engineering Analysis 7 Fall 2008 09 08 07 06 05 04 03 02 01 90090 90108 90120 a 12 11 09 08 602120 602108 602090 Figure 10 Evolution of dimensionless concentration a and temperature b for the b 20 production of ethylene glycol in an adiabatic CSTR with X00 and various initial temperatures Advanced Chemical Engineering Analysis 7 Fall 2008 21 One of the advantages of being able to generate the transient solution is to find out which steady state the system will approach Figure 10 shows the exact solution for various choices of the initial temperature showing how choice of initial condition dictates whether 81 or SS3 Table 1 will be achieved Error analysis The numerical solution found by Euler39s method differs from the exact solution due to the truncation of highorder terms in equation 73 In this section we will analyze the magnitude of the error of the approximate solution First we rewrite equation 68 as Yi1 yihftiyidi1 81 112 di1 yquot 82 where tiltEltti1 The parameter di1 is called local truncation error because it is the term of the expansion that is dropped to generate Euler39s method equation 74 The best way to characterize the accuracy of the method is by the error which is defined as the difference between the exact and approximate values of the function 61 Yi 1 83 Since the method starts with the initial condition it is obvious that e00 After the calculation at the first discretization point the error will be equal to the local truncation error since this will represent the only difference between the exact and approximate solutions equation 81 see also Figure 11 However for the calculation at the second discretization point and beyond the error will be increased by the local truncation error but errors made in previous points will propagate making ei greater in magnitude than di for iZZ Figure 11 To analyze quantitatively the propagation of errors we start by subtracting equation 74 from equation 81 Using the de nition of the error 83 this leads to ei1ZeihftiaYi fti ildi1 84 The error is best characterized in terms of the absolute value of the difference between the exact and approximate solutions Taking the absolute value of both sides of this equation and recalling the identity labclSlallbllcl 85 leads to lei1 S leil hfti Yi fti il ldi1 86 Using the mean value theorem we can state that Advanced Chemical Engineering Analysis 7 Fall 2008 22 ftiaYiftia i 3f 87 yr 39i 3y am y A yX exact solution t Figure 11 Propagation of errors in sequential numerical methods for IVPs the local truncation error adds to the overall error at each step but since the starting values for each step are approximate values already the total error will include accumulation of successive truncation errors where m is between yi and 371 Here we will assume that the partial derivative in equation 87 is bounded Let L be an upper bound for the magnitude of the derivative so that Of S L for all t 88 3y From equations 83 and 87 we can then state ftia Yr fti ml S Lleil 89 To nd a bound for the local truncation error we use equation 82 Let M be an upper bound for the magnitude of yquott lyquott s M for all t 90 We then have Advanced Chemical Engineering Analysis 7 Fall 2008 23 112 ldi1 SM 91 The combination of inequalities 86 88 and 91 leads to 112 leillslei1hL7M 92 We can evaluate this inequality sequentially starting at i0 recall that e00 to get 112 M s 7M 93 2 2 3 e2 3 e1 1hLhMSh2Mh ML 94 2 2 3 2 4 e3 3 e2 lhLLMS3h M3h MLL h M 95 2 2 2 2 which leads to 39h2M ei S 1 higher order terms 96 But when we evaluate the solution at a given ti the number of steps required to get there i is i2 97 so that ei S tihM higher order terms 98 This implies that ei 0h 99 That is even though the local truncation error is of Oh2 the propagation of errors in Euler39s method leads to an actual error of 0h This reduction in the order of the error by one unit always happens in sequential finite difference methods Because of this we say that Euler39s method is of rst order To illustrate the variation of the error at a point in time consider the example treated before of Advanced Chemical Engineering Analysis 7 Fall 2008 24 production of ethylene glycol in an adiabatic CSTR Figure 12 shows how the error at 171 changes with h in a loglog scale Equation 99 implies a linear relation in this scale with slope 1 As we can see the results follow this trend closely especially at low values ofh Euler39s method is very convenient due to its simplicity However the accuracy of the method can be improved by improving the level of approximation over that in equation 73 In what follows we will explore other methods to solve IVPs that have lower error orders 6 10E00 10E 01 r l 10E 02 I 10E 03 10E 04 I 10E 05 10E 06 7 D 10E 07 1 0001 001 01 h Figure 12 Euler39s method error at 171 as a function of h for the production of ethylene glycol in an adiabatic CSTR with initial conditions x00 and 601 The error is calculated from equation 83 using as exact solution the Euler solution at h10394 Filled symbols error in x open symbols error in 6 The solid line is a straight line of slope 1 2 Taylor39s second order method The simplest way to improve the level of approximation over that of Euler39s method is to increase the number of terms in the Taylor series used to represent the solution equation 67 For example including an additional term in the series leads to dy 1 d2y 1 d3y Yti1 Ytiatit H1 ti5dt72tit H1 ti2 F 0 H1 ti3 100 This equation can be written as 112 Yi1ZYihy39ti7yquotti0h3 101 Advanced Chemical Engineering Analysis 7 Fall 2008 25 To solve the basic IVP given by equations 71 and 72 we first notice that Y39 KEY 102 which implies that df t y quot 103 y dt It is important to recognize that the derivative on the righthand side of this equation is still a total derivative with respect to t since yyt To evaluate this derivative we start by expressing the total differential of f by Bf Bf df 7dtid 104 at 3y y Dividing by dt leads to g E 105 dt at By dt Using equations 102 and 103 we obtain y aiaif 106 Bt By which can be substituted in equation 101 to find h2 Bf Bf Yi1 Yi hftiYii tiyitiYiftiy1 Oh3 107 2 Bt By Truncation of this series leads to an approximation with local truncation error of Oh3 A A A h2 Bf A Bf A A Yi1 Yi hfaiYiEtiaYitiYiftiYi 108 This equation is used sequentially from the initial condition The method is called Taylor39s second order method The error is of Oh2 so it performs better than Euler39s method as h is decreased However the increase in accuracy for a given h comes with a price computational effort has increased considerably for each step In fact if the function f and its derivatives are complex enough most of the computational effort of this type of method will be associated with function and derivative evaluations Each step in the calculation required one function evaluation for Euler39s method 74 whereas now it requires three function and derivative evaluations equation 108 This means that the computational effort may have increased by a factor of 3 Advanced Chemical Engineering Analysis 7 Fall 2008 26 that is assuming that the derivatives of f in equation 108 are both nonzero Even so Taylor39s second order method is better than Euler39s method in terms of the accuracy of the solution for a given level of computational effort especially for fine discretization grids low h Nevertheless it would be interesting to find out if it is possible to develop a method with local truncation error of Oh3 but that requires less computational effort In what follows we explore a class of methods that accomplishes this 3 RungeKutta methods Methods based on Taylor series expansions such as Euler39s and Taylor39s lead to an approximation that has the form m a hSntiampi 109 where Sn is a series here n denotes the order of the method whose algebraic complexity increases with n For example for Taylor39s second order method we have equation 108 A h 3f A 3f A A s2 faiYi tiaYiitiYiftiYi 110 2 3t By This particular series contains the function f and its two firstorder partial derivatives From this point of view it looks like a Taylor series expansion for f In fact for a function of two variables the firstorder Taylor series expansion around a point t0y0 is given by 1 2 at axon to2 fty ftoyoftoYOt t0toYOy YO 111 e 1 lyyo2 i TIttoyyo 2 3y2 Btay where E is between to and t and 1 is between yo and y If we use as base point a generic point ty and evaluate the series at a point tocyB equation 111 can be written as Bf Bf ft0 YB f0YELY0 LYBR TI 112 where the residual R is given by 1 32f 1 32f 32f R if 2 if 2 7 113 Em 2 at iron 2 3y arms may 511th where E is between t and t0 and 1 is between y and yB If we compare equation 112 with equation 110 we can see that if we neglect the residual we can make the approximation 32tia9i fti09il5 114 Advanced Chemical Engineering Analysis 7 Fall 2008 27 by letting 0c h E 115 Bft191 116 Since XBOh then ROh2 which means that neglecting R will not incur in an error of higher order than already made in the truncation that leads to equation 110 Using equation 114 we can now propose a scheme based on equation 109 i1 ihft10lt il5 117 This new method preserves the order of Taylor39s second order method local truncation error of Oh3 and error of Oh2 while requiring only two function evaluations equations 116 and 117 one less than Taylor39s method Even though this might not seem like an exceptional advantage the savings in computational effort of this approach with respect to Taylor39s methods increase substantially as the order of the method is increased The explicit finitedifference method represented by equation 1 17 is a secondorder RungeKutta method One of the most widely used methods to solve IVPs is the fourthorder RungeKutta method RK4 Its derivation follows the ideas developed above and will not be presented here The basic equation is yi yik12k22k3k4 118 where k1 hft1amp1 119 k2 hftih2yik12 120 k3 hfti h2yi k22 121 k4 hfti hyi k3 122 This method has a local truncation error of Oh5 and an error of Oh4 Application of higher order methods to systems of ODEs The higherorder methods presented in the previous section can be extended to solve systems of ODEs such as that represented by equations 62 and 63 In the application of Taylor39s secondorder method to systems of equations the derivation of the appropriate expression for the second derivative of the dependent variables ie the steps equivalent to equations 103 to 106 Advanced Chemical Engineering Analysis 7 Fall 2008 28 gets more complicated due to the possible dependence of the functions on the rest of the dependent variables For example for equations 62 and 63 we can show that Bf YI 1 311 7f 7f 123 at 3y1 1 3y2 2 and 312 312 312 277fl7f2 124 at 3Y1 3Y2 so that the equations for Taylor39s secondorder method for a system of two ODEs are A A A A h2 Bf A A Bf A A A A y111 y11hf1t1y11y21i4t1y11y21401y11y21f1t1y11y21 2 at Byl 125 Bfl A A A A it1y11y21f2t1y11y21 3Y2 A A A A h2 sz A A sz A A A A y211 y21 hf2t1y11y217it1y115y21Tt1y11y21f1t1y11y21 Bf yl 126 i2t1 11 21f2t1 11 21 3Y2 The extension of the RK4 to systems simply consists in calculating the k factors 119 to 122 for each ODE separately For example for the system of two ODEs considered here a set of k39s is calculated in each step for each ODE as follows k1j hfjti 1i 2iaj12 127 k2j hfjti h2 1i k11 2amp2i k12 2 j12 128 kij hfjtih2y1ik212y2ik222j12 129 k4j hfjt1h 11k31a 21k32j12 130 and the values of the dependent variables for the next step are calculated by A A 1 Yji1 Yji gk1j 2k2j 2k3j k4j112 131 Figure 13 shows a comparison of the error at a point obtained in the solution of the problem stated by equations 30 and 31 with data for production of ethylene glycol These results show Advanced Chemical Engineering Analysis 7 Fall 2008 the clear advantage of higherorder methods at low h39s 10E01 10E 01 o 10E 03 r We 539 10E 05 10E 07 10E 09 10E 11 10E 13 A O Euler El Taylor 2nd order A RungeKutta 4th order 10E 15 0001 001 01 h 29 Figure 13 Errors at 171 as a function of h for the production of ethylene glycol in an adiabatic CSTR with initial conditions X00 and 601 The straight lines are arbitrary lines placed on the results for each method with slopes 1 Euler 2 Taylor 2nd order 4 RK4 Advanced Chemical Engineering Analysis 7 Fall 2008 30 Appendix A Derivation of the thermal energy balance on a CSTR In this appendix we will show how an energy balance in the reactor equation 19 leads to the nal energy balance used equation 20 First of all we make the following assumptions 1 kinetic energy changes are negligible 2 viscous dissipation is negligible and 3 there is no mechanical work performed in the reactor this would mean that the shaft work done by stirring devices has a negligible contribution to the energy balance These assumptions decouple thermal energy from mechanical energy so that all effects due to mechanical energy interactions are neglected Let n be the total moles in the reactor at a given time and U be the molar internal energy of the mixture inside the reactor internal energy per unit mol of mixture The energy accumulation term can be expressed as The time rate of change of dnU A 1 the total energy in the reactor dt Input and output of energy will be due to input and output of mass and in general heat exchange with the surroundings so that The net rate at which FiHi FeHq A2 energy enters the reactor where F represents total molar ow rate molestime H is molar enthalpy enthalpymole the subindices i and e denote inlet and exit respectively and q is the heat transfer rate energytime from the surroundings to the reactor Note that we have used the fact that HeH the molar enthalpy of the mixture inside the reactor since perfect mixing has been assumed The molar internal energy of the mixture in the reactor is related to the enthalpy by nU nH PV A3 so that dnU dnH dPV A4 dt dt dt Constant pressure and volume makes the second term on the righthand side of this equation vanish Collecting terms equations Al and A2 into the energy balance leads to dnH dt FiHi FeH q A5 To quantify the enthalpies we will consider that the mixture is an ideal solution for simplicity That is the enthalpy of the mixture is the sum of the enthalpies of its N components ie the enthalpy of mixing is negligible Advanced Chemical Engineering Analysis 7 Fall 2008 31 N nH Z nJ39HJ39 A6 jl Similarly N FiHi FDqu A7 Jl where iji represents the molar ow rate of component in the inlet ow Also N FeH Z FjeHj A8 j1 Substituting equations A6 to A8 into equation A5 leads to N dHJ39 dnj N N 21 nj dt HJ39 dt ZleiHji Zle Hj q A9 J J J The mole balance for species j can be written as follows equivalent to equation 3 dnj dt Fji FjeRjV A10 Substituting equation A10 into equation A9 and rearranging leads to N dHJ N N an zZFjiHji Hj VZRjHjq All F1 dt 397 j1 Since components are converted by the chemical reactions reference states for enthalpy calculations cannot be component speci c It is customary to choose the chemical elements that compose the species in the mixture as a basis to establish a reference state since elements will be conserved neglecting element transmutation by nuclear reactions Therefore the enthalpies of all compounds are calculated with respect to an elementspeci c reference state the enthalpy of all pure elements at their natural state at 25 C and 1 atm standard conditions is zero Because of this the enthalpy of each component inside the reactor and at the inlet ow can be expressed as follows 0 HJAHfJAHJ A12 Hji AHJ AHji A13 Advanced Chemical Engineering Analysis 7 Fall 2008 32 where AH J is the molar enthalpy of formation of component at standard conditions 25 C and 1 atm and AHJ39 and AHji are the molar enthalpy changes of component between the standard state and the reactor and inlet conditions respectively Substitution of equations A12 and A13 into equation A11 yields after manipulations note that the enthalpies of formation are constant N dAH N N J 0 Zn dt ZFJ1AHJ1 AHJ VZRJAHMAHJq A14 J21 11 11 To calculate AH we use the molar specific heat at constant pressure of componentj ij T AHJ IijdT A15 Tref where TreFZSOC Note that if there is a change of phase between the standard state and the state of the component in the reactor an enthalpy difference for the phase change would have to be added to this equation but this term will eventually cancel out Similarly AHJJ AHJ is the enthalpy change of component between the inlet and outlet conditions in the reactor so that T AHD AHJ jcpde A16 Ti The reaction term in equation A14 represents the change in enthalpy per unit time that occurs as the reactions proceed Here for simplicity we will consider that there is only one chemical reaction occurring and that this reaction involves component A The reaction will be written so that the stoichiometric coefficient for each compound j Vj is positive if the compound is a product and negative if it is a reactant and the stoichiometric coefficient for A is unity i1 For example in a multicomponent mixture where the following reaction occurs A2B gtC A17 we have VA1 VB2 Vc1 and Vj0 for all other components in the mixture 072A B C Note that for all components a stoichiometric balance yields VJ RJTRA A18 A Using this equation we can state that Advanced Chemical Engineering Analysis 7 Fall 2008 33 N N R ZRJAH AHJ7AZ vJAH AHJ A19 gtJ VA 1 J F1 J The summation term on the righthand side of this equation represents the change in enthalpy that occurs because of the reaction per mole of A reacted at the conditions at which the reaction is taking place This is the enthalpy of reaction N AHr ZvJAHgJ AHj A20 j1 Note that since the stoichiometric coefficients are positive for reaction products AHr represents an enthalpy difference of products minus reactants This means that AHr lt0 implies that thermal energy is quotliberatedquot by the reaction exothermic reaction whereas when AHr gt0 implies that the reaction quotconsumesquot thermal energy endothermic reaction Using equations A15 A16 A 19 and A20 equation A14 becomes N N T A dT A R anijgz ZFJ39J IijdT VviAAHrq A21 F1 j1 Ti A Since the mixture is an ideal solution we can define an average molar heat capacity for the mixture as A N A Cp ZXJij A22 an pj n p A23 j1 N T A T A EFF IijdT Fi IdeT A24 F1 Ti Ti where Fi is the total inlet molar ow rate Equation A21 can be written as T ndelz Fi JdeT VRiAAHrq A25 dt T VA 1 Finally this equation is usually expressed in terms of heat capacity Cp instead of molar heat Advanced Chemical Engineering Analysis 7 Fall 2008 34 capacity For the mixture in the reactor we can write n p mcp pVCp A26 Similarly Fi p chp A27 Equation A25 then becomes dT T R A pVCp E pQ IdeT VTAHI q A28 A Ti If the heat capacity can be approximated to be constant in the range of temperatures considered we have dT R pvcpaszcpm T VV7AAHI q A29 A which is equation 20 Advanced Chemical Engineering Analysis 7 Fall 2008 35 Appendix B An analytical relation between dimensionless concentration and temperature in the adiabatic CSTR with first order reaction Even though there is no known analytical solution for the ODE system represented by equations 30 to 33 an analytical relation between X and 6 can be obtained To find this relation multiply equation 30 by B and then subtract equation 31 This eliminates the nonlinear term yielding dX d6 7 1 1 6 131 5dr dr 5 X Now let uBX 6 B2 Using this transformation equation Bl can be written as dLHZIH 133 d1 The solution of this equation is uAe TB l B4 where A is an integration constant which can be evaluated from the initial conditions 32 and 33 Doing this and substituting back equation B2 into equation B4 leads to the final result after manipulations el5Xl3Xo90e Tl5116 T 35 This relates X to 6 for each 17 Advanced Chemical Engineering Analysis 7 Fall 2008 36 Appendix C A note on the solution of multiple linear non homogeneous first order ODEs and linear non homogeneous second order ODEs Multiple firstorder ODEs are common in reactor design problems If reaction kinetics is nonlinear or the energy balance is one of the equations the solution of the system of equations must be obtained numerically For problems that result in systems of linear firstorder ODEs analytical solutions can be found This appendix sketches how to solve this type of system for the case of two coupled ODEs Consider the system of ODEs d fzax yf1t C1 dy dt yx 5yf2t 02 where the coef cients will be considered constants for simplicity This system of two rstorder ODEs can be transformed into a single secondorder ODE as follows First take the derivative of equation C2 with respect to t 2 d y d15iy 7 C3 dt2 dt dt dt Substitution of dxdt from equation Cl into this equation leads to d2y dy df2 7 57 Xx 7 C4 dtz dt 73y Y Yfl dt Now to eliminate x from this equation we substitute it from equation C2 This leads to after rearrangements d2y dy d 7 5 7 5 f 7 C5 dtz 06 dt 06 73W 39Yfl 0 2 dt This shows how two firstorder ODEs can be transformed into a secondorder ODE The new ODE is linear and nonhomogeneous with constant coefficients We will explore now the solution of this type of equation Equation C5 is a particular case of the ODE d2 d a1tda2tyft C6 This differential equation is nonhomogeneous because yEO does not satisfy it as long as ft 0 Advanced Chemical Engineering Analysis 7 Fall 2008 37 for at least one point The associated homogeneous equation is dZYh dyh 7a t ia t 0 C7 dtz 1 dt 2 yh It can be shown that this equation has two independent solutions y1t and y2t consequence of it being second order Since the equation is linear any function that is a linear combination of the two solutions is also a solution Hence the general solution of equation C7 is YhtC1Y1tC2y2t CS where C1 and C2 are arbitrary constants This can be proven by substituting equation C8 into equation C7 and remembering that y1t and y2t are basic solutions The solution C8 is undetermined until two additional boundary conditions are used to evaluate the constants C1 and C2 Going back to the original equation C6 it is sometimes possible to find a specific function that satisfies the differential equation by inspection or by the application of specific techniques This solution of the nonhomogeneous equation is called particular solution ypt This solution does not necessarily satisfy the boundary conditions The general solution of equation C6 will be yt yht ypt C9 since substitution of this equation into equation C6 leads to the identity dZYh dyh dzyp dyp a tia t 7a tia t ft ClO m2 1 dt 2yh dt 10m 2yp 0 ft Example Consider the ODE dzy dy 475 t2 C11 dt2 dt y The homogeneous equation is dZYh dyh 4 5 0 C12 dtz dt yh To solve this equation we find the roots of the characteristic polynomial r2 4r50 c13 Advanced Chemical Engineering Analysis 7 Fall 2008 38 which are ri2iVl6 202ii Cl4 where i J l Independent solutions of equation C 12 can be shown to be yhtC1eTt C2er7t C15 To deal with the imaginary part of the eXponentials we make use of Euler39s formula eix cosxisinx C16 Using the roots from equation C14 and this formula equation C15 can be rewritten as yh t C1e2tcos t isin t C2e2tcos t i sin t C17 which can be expressed in terms of new constants as follows yhtC3e2t costC4e2t sint C18 A particular solution to equation Cll has polynomial form When the non homogeneous part of the equation is a polynomial it can be shown that the particular solution of a linear ODE with constant coef cients is a polynomial of the same degree as the non homogeneous term We then propose yttt p AZBC C19 Substituting this equation into equation CH and equating terms of the same power in the result the coef cients A B and C can be found In this case this leads to t2 8t 22 t777 C20 ypo 5 25 125 The solution of equation C l l is tC e2tcostc e2tsint 2 C 21 y 3 4 5 25 125 39 The only remaining thing to do is to find the constants C3 and C4 from the application of two boundary conditions
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'