CHEE 390 Computational Methods in Chemical Engineering - Project #2 Transient Diffusion in Hygroscopic Cylinders
Project #2
Transient Diffusion in Hygroscopic Cylinders
Instructor: Professor Phillip Servio (398-1026) Instructions
This computer project must be completed using MATLAB®. You must work individually and submit a report electronically as a PDF along with all the m-files required to execute your code. All the files must be contained within 1 zip file. Instructions necessary to execute your code must be imbedded as comments in your m-files. In order to submit your assignment, you will click on the reply to all in the project submission email that I will send out approximately 1 week before the assignment is due.
Problem Statement
You work for a company that specializes in designing equipment in order to dehumidify air. In the past, the company used silica beads but you and your team have developed a novel material with superior performance in the shape of a cylinder. The process can be seen in figure 1, utilizing a square cuboid chamber design. When the cylindrical desiccant becomes saturated inside the chamber, they are replaced with new ones. The saturated cylinders, which can not absorb anymore water, are then dried to their original state and ready for reuse.
Air + Water Air
Hygroscopic Cylinders of
length, L
Figure 1: Schematic of the square cuboid chamber used to dehumidify air with novel hygroscopic cylinders.
You are now asked to write a program using MATLAB® to predict the spatial concentration dependence of the hygroscopic cylinders in the square cuboid chamber, from r = [0 , r0] as a function of time during the absorption process. The equation for the dimensionless concentration of water in a hygroscopic cylinder during
the absorption process,θ*w, as a function of the dimensionless time parameter, Fom, and the dimensionless position, r*, is given by the following relations:
* Cw,i−Cw ∞ 2 *
θw =C −C =1−∑Gnexp(−ζnFo)Jo(ζnr) (1)
r0
w ,i
w,∞ n=1
o
Gn = 2 J1(ζn) (4)
Fom = DABt r 2o
(2) r* = r (3)
ζn J 20(ζn)+J 21(ζn)
The values of ζn are the positive roots (0 < ζn < ∞) of the following transcendental equation:
0 = Bi−ζ J1(ζn) (5) n J0(ζn)
where
Bim = hmr
DAB
(6)
Equation 5 must be solved using a numerical root finding technique such as Brent’s Method. Brent’s method is a popular technique used for finding roots of equations of the form F(x)=0. Although it is more algebraically complicated than standard root finding methods that we have discussed in class, it has convergence characteristics (superlinear) that are comparable to the Newton-Raphson method (quadratic). Brent’s method uses interpolation and will reliably locate a root once an interval bracketing a sign change has been found. Unfortunately, just like the Newton-Raphson Method, Brent’s method also behaves badly in the vicinity of a singularity!
After bracketing a sign change, the method then uses the bisection method or linear interpolation to locate a point, x2, which lies in-between the bracketed interval x1 & x3 such that x1 < x2 < x3. The method then uses inverse quadratic interpolation of the 3 points on F(x) associated with these positions to iteratively improve the value of the root that has been bracketed. The inverse quadratic function which passes through the 3 points is directly given by the Lagrange interpolation polynomial and is shown in the algorithm below. It is your responsibility to modify the below algorithm to deal with singularities that cause instability in Brent’s Method (Point III). NOTE: Please explain in the report section why you would use the bisection method in this scenario over incremental search or linear interpolation.
Brent’s Method
-
Choose an increment ∆x which is small compared to the interval being searched [x0 , xf].
-
Locate 2 positions using incremental search and label them x1 & x3 (x3 – x1 = ∆x) such that y1 ∙ y3 < 0 or
until x > xf .
-
Find the position x2 (x1 < x2 < x3) where the straight line connecting (x1, y1) & (x3 ,y3) intersects the x-
axis.
x2 =(x1 y3−x3 y1) (y3 − y1)
IV. Find a better estimate of the root, x4, where x4=x2+qp
-
Evaluate y4=f(x4)
-
Set x2=x4 & y2=y4
p=s[t(r−t)(x3−x2)−(1−r)(x2−x1)] q=(r−1)(s−1)(t−1)
r=y2 s=y2 & t=y1 y3 y1 y3
IF x1<x4<x2 x3=x2 & y3=y2
ELSE
x1=x2 & y1=y2
VII. Repeat steps IV, V & VI until |p/q|⩽Tol x4
Report
Each person is required to submit a report. Some of the results that you can expect are given below in
Figure 2, which shows the dimensionless concentration,θ*w, vs dimensionless position, r*, for varying Fourier
numbers, Fom. You must vary the Biot number, Bim = {0.01, 0.1, 1, 10, 1000}, as well as the Fourier number, Fom = [0 , 2] and comment on the validity of the first term approximation assumption, Fom > 0.2, and the lumped capacitance assumption, Bim < 0.1, in your discussion.
The first term approximation assumption implies that if the dimensionless time parameter is greater than 0.2, Fom > 0.2, only 1 term of the infinite series is required to give an adequate estimate of the concentration profile. In other words, you only need 1 eigenvalue! The lumped capacitance assumption states that if the Biot number is less that 0.1, Bim < 0.1, the concentration profile is independent of spatial effects.
Your program must also have the capability of telling the user how long it takes for the centerline dimensionless concentration,θ*w(0,Fom), to reach a specific value. This is important because the centerline concentration has the slowest response with respect to time compared to any other position. In other words, all other concentrations at any other location will be higher in this process than the centerline value, see figure 2.
The format of the report should be as follows:
-
Title page
-
Objective
-
Flowchart
-
Results & Discussion
-
Conclusion
-
Nomenclature
-
