Download Mathematical Modeling - Advanced Engineering Analysis | CHEE 502 and more Study notes Chemistry in PDF only on Docsity!
CHEE 502
ADVANCED CHEMICAL ENGINEERING ANALYSIS
CLASS NOTES – FALL 2008
Eduardo Sáez
INTRODUCTION
MATHEMATICAL MODELING
"A mathematical model is a representation, in mathematical terms, of certain aspects of a nonmathematical system" (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 (e.g. 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 (e.g. 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 "model" 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 and/or energy conservation equations. If the equations are differential equations, find the appropriate boundary and/or initial conditions to be applied. (3) Determine the physical properties and constitutive equations to be used (e.g, 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 verified, either a priori or a posteriori. (5) Simplify the equations to form a coherent, well-posed system. (6) Solve the equations, either analytically or numerically.
(7) Before the solution is used, the model must be validated. Typically, this involves finding 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 (i.e., 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.
dissipation in the reactor; i.e., 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):
bychemical reaction
massofspeciesis
The rateat which
ofspecies thereactor
The rateat whichmass
speciesin thereactor
thetotalmassofachemical
Thetimerateofchangeof
created
net
enters
net
The terms on the right-hand 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
bychemical reaction
molesofAare
The rateat which
ofA thereactor
The rateat whichmoles totalmolesofAin thereactor
Thetimerateofchangeofthe
created
net
enter
net
To express equation (2) in mathematical form, we first 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 1 for notation):
Qc Qc R V dt
d (Vc ) Ai A A
A = − + (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 → Pr oducts (4)
then the reaction is first order in A and the rate of reaction is
R (^) A = −kc A (5)
where k is the reaction rate constant, which can be expressed as a function of temperature using Arrhenius equation,
RT
Ea k k 0 e
− = (6)
where k 0 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,
A Ai
A (^) c V
Q
c V
Q
k dt
dc =
+ ^ + (7)
which will be subject to an initial condition,
cA = cA 0 , t=0 (8)
The quantity tR=V/Q is the residence time in the reactor and its use in equation (7) leads to
Ai R
A R
A (^) c t
c t
k dt
dc =
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
y f(t ) dt
dy
The integrating factor technique is based on the use of the following identity
e y dt
dy e dt
d (ye t ) t t = α + α α
α (11)
Note that, if we multiply equation (10) by eαt^ , then the left-hand side of the resulting equation is precisely the right-hand side of equation (11). Therefore, equation (10) becomes
e f(t ) dt
d (yeαt^ ) αt = (12)
If y(t) is subject to the initial condition
y=y 0 , t=0 (13)
t
cA
cA
cA
cAss
cA0>cAss
cA0 Ai R
A R
RT
Ea 0
A (^) c t
c t
k e dt
dc
− (23)
RT A r
Ea p (^) dt QCp(Ti T) Vk 0 e c H
dT ρVC =ρ − − ∆
− (24)
This is a system of two first-order, coupled, nonlinear ODEs. They are subject to the initial conditions
cA = cA 0 , t=0 (25)
T = T 0 , t=0 (26)
Before we explore the solution of this problem, we will formulate it in dimensionless form. We define the following dimensionless variables:
Ai
A c
c x = (27)
T i
T
θ = (note the need to use absolute temperatures) (28)
t R
t τ = (29)
Equations (23) to (26) become
θ
γ − = − − α τ
1 x xe d
dx (30)
θ
−γ = −θ− αβ τ
θ 1 xe d
d (31)
x = x 0 , τ=0 (32)
θ =θ 0 , τ=0 (33)
where
α =tR k 0 (34)
explore the solution of this equation using the Newton-Raphson 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 (∆Hr>0 ⇒ β>0) The graphical solution will be found by plotting the functions on both sides of equation (42)
vs. θ, and finding the point at which both functions intersect. Note that, in this case, both functions are positive since β>0. The two sides of equation (42) can be physically interpreted as follows. The left-hand side,
Ef (θ )= 1 − θ (43)
represents, in dimensionless form, the energy lost by the fluid that goes through the reactor, while the right-hand side,
−γ θ
−γ θ
αβ θ = /
/ r (^1) e
e E ( ) (44)
represents the thermal energy consumed by the reaction progress. The steady state solution is the value of θ for which Er=Ef. To find the solution graphically, we plot Ef vs. θ, which will be a
straight line of slope -1, and Er vs. θ, which will be a monotonically increasing function that at large θ tends asymptotically to αβ/(1+α). 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 (∆Hr<0 ⇒ β<0) In this case it is convenient to change the signs in the definition of the functions that intersect in equation (42). Let
Ef (θ) =θ− 1 (45)
Now this function represents the energy gained by the fluid that goes through the reactor, while
−γ θ
−γ θ
−αβ θ = /
/ r 1 e
e E ( ) (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 α, β, and γ. 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
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.
Ef(θ)
θ
Er(θ)
solution
Figure 3. Qualitative illustration of the graphical solution of equation (42) for an endothermic reaction.
(a) (b)
Figure 4. Qualitative illustration of the graphical solution of equation (42) for an exothermic reaction: (a) single solution, (b) multiple solutions.
Ef(θ)
θ
Er(θ)
1
Ef(θ)
θ
Er(θ) SS
SS
SS
and x.
θ
P
P
P
Ef Er
Ef>Er
Er>Ef
Ef>Er
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 x=xi, where i denotes the number of the iteration being performed, the function is evaluated at x=xi+1, and the second term is truncated to produce the approximation
f (xi + 1 )≈ f(xi)+f'(xi)(xi+ 1 −xi ) (51)
The value xi is a known value of x for which f(xi)≠0. Equation (51) is applied to find a new value, xi+1, aiming at f(xi+1)=0. Using this in equation (51) and rearranging leads to
f'(x )
f(x) x x i
i i + 1 = i− (52)
Because equation (51) is an approximation, this equation does not guarantee that f(xi+1)=0. However, xi+1 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 y-x plane, equation (52) is equivalent to extrapolating to y=0 a straight line y=y(x) that passes through the point (xi,f(xi)) and whose slope is f'(xi) (Figure 6). To estimate error propagation in the NR method, let xs be an exact solution,
f(xs)=0 (53)
Figure 6. Graphical representation of the Newton-Raphson method. From a point x=xi, a better approximation to the solution (xi+1) is found by extrapolating (for this case; in other cases this step may be an interpolation) to y=0 a straight line that coincides with the function f(x) at x=xi 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 x=xs, with x 0 =xi, and use equation (53) to get
i i s i f"()(xs xi)^2 2!
0 = f(x)+f'(x)(x −x)+ ξ − (54)
where ξ is between xi and xs. We define the absolute error of the method at a given iteration by
Ei = xs−x i (55)
Subtracting equation (51) from equation (54) (recalling that we are letting f(xi+1)=0) leads to
i s i 1 f"( )(xs xi)^2 2!
0 = f'(x)(x −x+ )+ ξ − (56)
from which we get
s i^2 i
s i 1 (x x) 2 f'(x)
f"( ) x x −
ξ − (^) + = − (57)
(36) to get
α=1.101× 1012 , β=-0.200, γ=30.
To solve this problem, we apply the NR method to solve equation (42). First of all, we rewrite the equation as
−γ θ
−γ θ
αβ θ =θ− + /
/ 1 e
e f ( ) 1 (60)
so that equation (42) is equivalent to f(θ)=0. Differentiation of equation (60) with respect to θ yields, after manipulations,
/ 2 2
(e )
e f '( ) 1 +α θ
αβγ θ = + γ θ
γ θ (61)
Application of the NR method (equation 52) for different initial values of θ 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 10-9^ 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 (^) θ x SS1 1.0209 0. SS2 1.1000 0. SS3 1.1647 0.
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 first-order ODEs: initial value problems Consider the following system of two ODEs with two dependent variables, y 1 and y 2 ,
f(t,y,y ) dt
dy 1 1 2
f (t,y,y ) dt
dy 2 1 2
Number of iterations
θ
Figure 7. Convergence to the steady-state 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, 1.08 and 1.2.
subject to the initial conditions
y 1 =y1,0, t=0 (64)
y 2 =y2,0, t=0 (65)
Since the boundary conditions are specified at the same boundary (t=0), 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 finite difference methods. The simplest of these is Euler's method.
y"( ) O(h ) 2
h (^22) ξ = (69)
(note that y"(ξ) will tend to y"(ti) as h→0). We then can write equation (68) as follows,
y (^) i + 1 =yi+hy'(ti)+O(h^2 ) (70)
Now consider that the function y satisfies the initial value problem
f(t,y ) dt
dy = (71)
y=y 0 , t=0 (72)
Substituting y' from equation (71) into equation (70) leads to
y (^) i + 1 =yi+hf(ti,yi)+O(h^2 ) (73)
Euler's method is based on truncating the second-order term in this equation:
yˆ^ i + 1 =yˆi+hf(ti,yˆi ) (74)
where the circumflex will be used to denote approximate value of the dependent variable. The
sequential application of equation (74) from i=0, 1, 2… (starting from yˆ 0 = y 0 , the initial
condition) leads to calculation of approximate values of the solution at the discretization points. Note that the values of yˆ (^) i+ 1 depend only on the value of the dependent variable in the previous
discretization point ( yˆ ). This means that the calculation can be performed explicitly fromi
equation (74) regardless of the complexity of the function f. Because of this characteristic, Euler's 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 Euler's method (equation 73) to a system of two ODEs (equations 62 and 63). The appropriate discretized equations are:
yˆ^1 , i+ 1 =yˆ 1 ,i+hf 1 (ti,yˆ 1 ,i,yˆ 2 ,i ), i=0, 1, 2… (75)
yˆ^2 , i+ 1 =yˆ 2 ,i+hf 2 (ti,yˆ 1 ,i,yˆ 2 ,i ), i=0, 1, 2… (76)
Direct application of these equations to equations (30) and (31) leads to
= + − − α θ
−γ
ˆ xˆ (^) i 1 xˆi h 1 xˆi xˆie , i=0, 1, 2… (77)
θ =θ + −θ − αβ θ
−γ
ˆ ˆ (^) i 1 ˆi h 1 ˆi xˆie , i=0, 1, 2… (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 Euler's method, consider the production of propylene glycol with initially no propylene oxide in the reactor and a reactor temperature equal to the feed temperature:
x 0 =0 (79)
θ 0 =1 (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 h=0.02 is exact for the scale of the plot. Note that the solution for h=1 is quite inaccurate, but h=0.25 leads to a solution that is relatively close to the exact solution.
h= h=0. h=0. h=0.
x
τ
Figure 9. Dimensionless concentration vs. time curves generated by Euler's method for various values of h for the production of ethylene glycol in an adiabatic CSTR with initial conditions (79) and (80).
