- Nam Sun Wang
- Department of Chemical Engineering
- University of Maryland
- College Park, MD 20742-2111
- 301-405-1910 (voice)
- 301-314-9126 (FAX)
e-mail: nsw@eng.umd.edu
©2002 by Nam Sun Wang
- Department of Chemical Engineering
| Due | Homework | Solutions (Computer Files) |
|---|---|---|
| 9/09 | Homework #1 | |
| 9/16 | Homework #2 | |
| 9/25 | Homework #3 | |
| 10/02 | Homework #4 | |
| 10/09 | Homework #5 | |
| 10/16 | Homework #6 | |
| 10/23 | Homework #7 | |
| 10/30 | Homework #8 | |
| 11/06 | Homework #9 | |
| 11/13 | Homework #10 | |
| 11/20 | Homework #11 | |
| 11/27 | Homework #12 | |
| 12/04 | Homework #13 | |
| 12/11 | Homework #14 | |
axi
yi = ---------
1+(a-1)xi
yi = axi/(1+(a-1)xi) ... if your browser fails to line up the above equation.
where xi and yi are the mole fractions of
the more volatile component in the liquid phase and vapor phase,
respectively, at the ith stage, and a
is the relative volatility between the two components.
Write a program to input x1 (which is called the bottom composition), yn (which is called the top composition), and a (which is called the partition coefficient). Your program will then proceed to find n, the number of stages needed to distill a binary mixture from x1 to yn. Finally, print the answer and plot the above equation and connect x and y in a staircase manner (which is called the McCabe-Thiele diagram).
Background. Start with a liquid feed to the 1st stage with x1; then from the above equation find y1, the composition of the vapor phase in equilibrium with the liquid. Condense the vapor from the first stage, and that becomes the liquid feed to the second stage (i.e., x2=y1). Now that you know x2, find from the above equation y2... Repeat this process for i=1,2,...,n until yn exceeds the specified value. Moonshine example: you start with a fermentation broth containing 5% alcohol and you wish to find how many stages are needed to distill it to 95% under total reflux. Note that the top composition is a function of the number of stages, the partition coefficient, the bottom composition... We do not have to be able to express it nicely in terms of sine, cosine, etc. in order to define a function.
Solution:
The solution appears in one of the links in my ENCH250 web page.
x2*y" + x*y' + x2*y = 0Or we can define it in terms of Taylor's series:
J0(x) = 1 - (x/2)2/(1!)2 + (x/2)4/(2!)2 - (x/2)6/(3!)2 +...Write a function that returns J0(x) accurate to at least 5 significant digits via Taylor's series expansion.
Solution:
In the Mathcad version, I test it with x=0.0, 0.1, ..., 10.0 and plot
the Taylor's series expansion and J0(x) on the same graph.
Notice how clean the Mathcad code is -- thus my reason for relying on it
for my lecture notes and future homework solutions, even though Mathcad is
about 100 times slower than Fortran and it sometimes give wrong answers.
1 óp
J0(x) = --- ô cos(x*sinq) dq
p õ0
J0(x) = 1/p * integral (from 0 to p) cos(x*sinq) dq ... In case your internet browser fails to display the above equation
Numerically calculate the integral to at least 5
significant digits. Compare your answer to the value you get from a table
or from your calculator or mathematical package.
Solution:
See the last problem.
Solution:
Intelligent selection of data spacing helps minimize error.
Never extrapolate.
Solution:
This is identical to the last problem, except for the function
and interval. A polynomial interpolation with 21 Chebychev
interpolation points seems to do the job. Missing from problem
statement is the number of points and where to place these
points. You need to base your judgement. Note that the Matlab
file here uses the formula based on matrix inverse, where the
matrix to be inverted is almost singular for large n, say, n=26;
thus, the results are not to be trusted. The lecture note warns,
"In practice, we do not find the polynomial coefficients from the
above formula because the columns of matrix X are somewhat
correlated. Thus, the matrix is often ill-conditioned and nearly
singular, especially when n is large." Many of you did not heed
this warning. Note that at interpolation points, the polynomial
passes through the given function exactly. If this is not the
case, you are seeing only computational artifacts.
x y --------- 0.0 0.00 0.1 0.95 0.2 0.90 0.3 0.95 0.4 0.10 0.5 0.05 0.6 0.05 0.8 0.20 1.0 1.00If one single continuous polynomial is not to your liking, try piecewise polynomial (e.g., quadratic or cubic spline).
Solution:
This problem tries to convey that we should not fit data blindly.
When we connect data points with a curve, we sometimes create trends
that do not exist.
Solution:
Zeros of J0 are important in eigenvalue-eigenvector problems.
The following file finds zeros of other Bessel's functions as well.
R*T a(T)
P = --- - -----------------
V-b V*(V+b) + b*(V-b)
where
R2*Tc2
a(T) = 0.45724*------*alpha(T)
Pc
R*Tc
b = 0.07780*----
Pc
alpha(T) = [ 1 + k*(1-sqrt(T/Tc)) ]2
k = 0.37464 + 1.54226*w - 0.26992*w2
The critical parameters for various compounds have been tabulated, and
the following are the values for water.
Tc = 647.3 K (Critical temperature) Pc = 22.048 MPa (Critical pressure) Vc = 0.056 m3/kmol (Critical volume -- not used in this problem) w = 0.344 (acentric factor) R = 8.314x10-3 MPa m3/kmol K (Gas Constant)
fugacityliquid=fugacitygaswhere fugacity is defined as:
fugacity óP p*V/(R*T) - 1
ln(--------) = ô ------------- dp
P õ0 p
ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp (if you don't see the integral symbol above)
Applying the Peng-Robinson equation of state to the definition of
fugacity results in:
fugacity A Z+(1+sqrt(2))*B
ln(--------) = (Z-1) - ln(Z-B) - -----------*ln[---------------]
P 2*sqrt(2)*B Z+(1-sqrt(2))*B
where
Z = compressibility factor = P*V/(R*T)
A = a(T)*P/(R2*T2)
B = b*P/(R*T)
With this additional
fugacity equation, the pressure we find is the saturation
vapor pressure in equilibrium with the liquid, and that pressure
depends on only one variable, not two variables. (Each equation
or phase reduces the degree of freedom by one.) Now, P(V,T)=Psat(T).
Compute and plot the vapor pressure of water versus temperature
from 0°C up to the critical
temperature, Tc. Finally, the boiling point is the
temperature at which the vapor pressure equals the barometric
pressure, which is normally 1atm at sea level but lower
at high elevations. Calculate the boiling point of water at
1atm.
Solution:
The basic concept is the following.
The following worksheet is a bit messy, because Mathcad does not
return polynomial roots in a systematic manner (e.g., clean
ordering of real and complex roots). Otherwise, it is quite
straight forward.
Diffusion = Reaction d2y --- = f(y) Boundary Conditions: at x=0 dy/dx=0; at x=1 y=1 dx2Use the following saturation (Michaelis-Menten) kinetics:
f(y)=a*y/(b+y)Find a solution via the finite difference method, where the first derivative dy/dx and the second derivative d2y/dx2 are approximated as:
dy/dx ==> 0.5(yi+1-yi-1)/Dx d2y/dx2 ==> (yi+1-2yi+yi-1)/(Dx)2where Dx is the step size. You may choose your own values for a and b. In my solution, I set a=1 and b=1. How small does Dx have to be? In other words, how many points yi do you need?
Solution:
Express the above ODE as a set of algebraic equations of the form
g(y)=0 and solve for y numerically with Newton's method.
We can handle practically any ODEs this way, no matter how complicated they are.
Mathcad program shows very fast convergence (within a few iterations).
The Matlab program below gives y(0)=0.778 for 100 intervals and y(0)=0.777 for
200 intervals. Thus, 100 subintervals shall suffice.
Solution:
In my solution, I first tried powers of T and powers of wt%.
I was not satisfied with the results. I then added the
cross-product terms, such as T*wt%, T2*wt%,
T*wt%2, T2*wt%2. The results
got a bit bettter, but I still was not completely satisfied.
Then I remembered from thermodynamics and every day experience
that viscosity decreases exponentially with temperature...
Finally, I am happy.
Example: At T=50°C & wt%=40, predicted viscosity=2.503 centipoise (measured viscosity=2.506 centipoise).
At T=70°C & wt%=40, predicted viscosity=1.620 centipoise (measured viscosity=1.614 centipoise).
Because the axes are not linear, you should regress the axes scale versus the coordinate at several points to obtain a relationship between them. I lifted the coordinates given below from p2-321 graph by viewing the .gif image file of the graph in an image program (e.g., Photoshop or Paint Shop) and by reading the coordinates under the cursor. The (x,y) coordinates and the viscosity from Table 2-326 has been transcribed into an ASCII file here.
-------------------------
Temperature Coordinate
(Deg C) x y
-------------------------
-100 129 143
0 129 386
100 129 559
200 129 687
300 129 795
400 130 881
500 130 958
600 130 1025
700 130 1085
800 131 1138
900 131 1188
1000 132 1232
-------------------------
Likewise, we regress viscosity versus coordinates
-------------------------
Viscosity Coordinate
(10^-7 poise) x y
-------------------------
10000 814 141
8000 814 222
6000 814 325
5000 815 393
4000 815 473
3000 816 580
2000 816 723
1500 817 827
1000 817 974
800 817 1058
600 817 1162
500 817 1230
-------------------------
Coordinates at the four corners of the grid
-------------------------
Grid Coordinate
X Y x y
-------------------------
0 0 343 1130
18 0 773 1131
0 30 341 422
18 30 771 422
-------------------------
Solution:
For this problem, I modified the
specific heat example.
I slightly re-arranged the X-Y grid data file into such a format that Mathcad
can read (and changed names of compounds that start with a number, e.g.,
3-Methylene-1-butene and 2,2,3-Trimethylbutane).
Diffusion = Reaction d2y --- = f(y) Boundary Conditions: at x=0 dy/dx=0; at x=1 y=1 dx2Use the following saturation (Michaelis-Menten) kinetics:
f(y)=a*y/(b+y)Add 10% random noise to your simulated results of y; the random noise represents measurement errors on concentration y at various positions. From the noisy data (and now pretend you do not know what the parameters were), estimate the reaction rate parameters. For example, start simulation with a=1 and b=1, and estimate a and b from noisy y via regression.
Solution:
I modified the solution from the previous homework and added the noise
portion and performed nonlinear regression.
fi(x)=J0(lix) where li=i for i=1,2,..,5
Solution:
Notice the number of times the functions cross
the x-axis (i.e., number of zeros) within x=[0, 1].
fi(x)=J0(lix) for i=1,2,..,5where li are the first five roots of J0(x)=0 (which you had already calculated in a previous homework assignment).
Solution:
The given functions are already orthogonal, if we define the
scalar product with a weighting function w(x)=x. We simply
normalize the length. Notice the number of times the functions cross
the x-axis (i.e., number of zeros) within x=[0, 1].
fi(x)=J0(lix) for i=1,2,..,5where li are the first five roots of J0(x)=0
Solution:
Note that
J0(lix)=0 at x=1. (After all,
that is how the zeros li are defined.)
Thus, as we increase the number of terms, the approximation becomes better,
but the approximated function is always 0 at x=1. This is because the
sum of many zero terms is zero, no matter how many terms are present.
x2y" + x*y' + x2*y = a*y for x=[0,1] B.C.: y(0)=0 y(1)=1Start with, say, a=1 (but I can easily change it to be anything by entering a different value in my worksheet). An analytical solution exists for the above ODE. You may want to compare your approximation to the analytical solution. (This problem is similar to the heat conduction example on the 10/09's class web page.)
Solution:
The purpose of this problem is to demonstrate the concept of
expressing a vector as a linear combination of limited number of
vectors. I am not saying that we should go through this kind of
trouble when we know what the solution is. If the given ODE is
much more complicated and an analytical solution does not exist,
we can continue to apply the same method. Approximating a
function with a linear combination of basis functions reduces
expressing a given function with a small set of numbers. It
reduces solution of an ODE into solution of an algebraic
equation. It reduces a PDE into an ODE... Many researchers are
occupied with the method of weighted residuals (MWR).
a0=-9.2026E+00 intercept a1= 4.2360E-03 total enrollment a2= 9.9039E-02 research $ a3= 1.9728E+00 student/faculty ratio a4=-1.2351E-01 acceptance rate a5= 1.1434E-01 quantitative GRE score
Solution:
Mole Total
Fraction Concentration -----Spectral Intensity-----
(M)
y<1> y<2> x<1> x<2> x<30>
------------------------------------------------------
0.75 3.5230 1.6600 1.5830 ... 2.5220
: : : : :
: : : : :
Here is one such set of experimental data (nadh.dat) for you to work with.
The first column of the data file contains the mole fraction, the
second column the total concentration, and columns three and
beyond the spectral intensity at a series of wavelengths. Thus,
the dependent variables are the sample composition, i.e., mole
fraction and total concentration;
Y = [y<1> y<2>]and the independent variables are the sample color, i.e., spectral intensities at thirty wavelengths.
X = [x<1> x<2> ... x<30>]Before you start, you may want to center each variable around the mean value and rescale with the standard deviation so that the new scaled variables are roughly of order 1, i.e., y<1>~[-1, 1]. Because the independent variables are correlated, the scalar product matrix of the various vectors x<j>, i.e., the matrix XTX, is singular (or nearly singular), and naive regression based the following normal equation does not work. You will have trouble evaluating (XTX)-1.
Naive regression: mole fraction y<1> = a11*x<1> + a21*x<2> + a31*x<3> + ... + a30,1*x<30> + error<1> Naive regression: total conc. y<2> = a12*x<1> + a22*x<2> + a32*x<3> + ... + a30,2*x<30> + error<2> Regress Y against X: Y=X*a+error Normal equation provides the solution: a=(XTX)-1*XT*YFind the eigenvalues for the square matrix XTX and list them in decreasing order. Also find the associated normalized eigenvectors, v<1>, v<2>, v<3>.... How many independent eigenvectors are there? Show that all of these eigenvectors are mutually orthogonal (i.e., v<i>T*v<i>=1 for i=j, vT*v=I). It is better to describe each sample in terms of values (scores) along these mutually orthogonal eigenvectors (loadings) rather than x<i>.
score<i>=X*v<i>In other words, we employ a new coordinate system constructed out of eigenvectors.
mole fraction y<1> = a11*score<1> + a21*score<2> + a31*score<3> + ... total conc. y<2> = a12*score<1> + a22*score<2> + a32*score<3> + ...Find the coefficients. How many terms to you need to describe adequately mole fraction and total concentration? Note that ai1*score<i> is the projection of the vector y<1> onto the vector score<i>.
Projection of y<1> onto score<i> = ai1*score<i> = (y<1>,score<i>)/(score<i>,score<i>)*score<i>Thus, the coefficient ai1 is
ai1 = (y<1>,score<i>)/(score<i>,score<i>) ai1 = (y<1>,score<i>) if score<i> is normalized, i.e., (score<i>,score<i>)=1Finally, provide the regression equation y(spectral intensities) to predict composition (mole fraction and total concentration) from color.
Solution:
dx1/dt = 1.5*x1 + x4 dx2/dt = 0.75*x1 + 2*x2 + 0.5*x4 dx3/dt = 0.875*x1 - 0.5*x2 + 3*x3 - 0.25*x4 dx4/dt = -0.25*x1 + 2.5*x4Note that we can write the above set of first ODEs compactly as:
dx/dt=A*x x(0)=x0Find the eigenvalues and eigenvectors of A and exp(At). Solve the ODE via the eigenvalue-eigenvector approach. Try a few initial conditions.
Solution:
L*y = [ (x2-1)d2/dx2 + x*(d/dx) ]*y y(-1)=+/-1 y(1)=1 (or any symmetrical boundary condition)Is this linear transformation L symmetrical? Note that the definition of a symmetrical transformation involves scalar product; thus, you need to come up with a suitable definition. How many eigenvalues are there? How many eigenfunctions are there? Plot the first few eigenfunctions. Can you find a systematic trend with regards to the number of zeros within the interval of interest? One reason for solving this type of problem is, of course, when we encounter just such an ODE in our work. Another reason is simply to generate a set of orthogonal basis functions, because these eigenfunctions provide a set of basis to represent any other continuous functions as a linear combination of the eigenfunctions.
f(x) = S aj*vjFind the coefficients aj. Try f(x)=J0(x) in x=[0,10] as an example. Another example: f(x)=0 in x=[-1 1]. (Since our eigenfunctions are nonzero at x=-1 or x=1, are we even allowed to represent f(x)=0?
Solution:
drag ~ v2*sin2y*(cross-section area)We have set up all the relevant equations in class. Solve them and plot the side profile of the bullet. Common sense tells us that a long slender bullet like the one used by the recent sniper travels farther than a short fat one. Compare the drag on an optimized bullet of 5 length unit and 10 length unit. (That is, provide the value of the quantity we try to minimize.)
Solution:
The problem statement does not specify the final position xf.
We either specify it or impose a condition to determine it.
Min J=w*s(1)-b(1)
where w=the weighting factor (Try w=1.1 and w=2, and compare the paths)
s(x)=the distance traveled between x=0 and x
b(x)=the amount of berries picked between s=0 and s(x).
Find the path y(x) you should take if the density of berries (or
equivalently the probability of encountering them) is given by
the following function:
g(x,y)=exp[-(x-0.25)2-(y-0.5)2]
Solution:
Similar problems include: a bible salesman traveling to sell
most bibles, my kids collecting candies on Halloween night, a
politican campaigning through his/her district to collect votes,
light traveling through space, sunlight being deflected by the
atmosphere, spacecraft traveling through gravitational field, an
electron or ions traveling trough a magnetic field...
e*y" + (1+x)*y' + y = 0 B.C. y(0)=y(1)=1
Solution:
e*y" + sqrt(x)*y' - y = 0 B.C. y(0)=0 y(1)=e2
Solution:
e*y" - x2*y' - y = 0 B.C. y(0)=y(1)=1
Solution:
|