Project ReportProblem Solving Using Computers (25CHME21I)
Prepared from the lecture notes and the project brief provided in class.
This report follows the lecture methods for matrices, MATLAB plotting, Excel goal-seek / Solver, material balance,
and PFR modeling. The project brief requires a pipe-network solution by the Cross Hardy method and a PFR
simulation in MATLAB and Excel.
Q1. Pipe Network Analysis
The network has 19 pipes and 8 loops. The project brief asks for a Hardy Cross solution using the Darcy–Weisbach
equation, with a summary table containing Q, L, D, ε/D, velocity, Reynolds number, friction factor, K, head loss, and
n·hL/Q. It also specifies that the diagonal entries of the loop matrix are positive sums of n·hL/Q and the off-diagonal
entries are negative shared-pipe terms.
Pipe network diagram from the project brief:
Recommended Excel setup for the first iteration:
Pipe No. Q (m3/s) L (m) D (m) ε/D V (m/s) Re f K hL,
n·hL/Q
1 0.04814 457 0.305 0.00087
4444444
44444
2 0.04166 305 0.203 0.00104
6666666
666664
3 0.03611 366 0.203 0.0013
1111111
11111
4 0.00183 610 0.203 0.0013
3333333
3333333
Use the lecture formulas exactly as given in the project: V = Q / (πD²/4), Re = ρDV/μ, f = 64/Re for laminar flow, f =
0.25 / [log10((ε/D)/3.7 + 5.74/Re^0.9)]² for turbulent flow, K = 8fL / (π²gD⁵), hL = KQ², and n = 2 for the Hardy Cross
term n·hL/Q. The loop corrections are obtained from ΔQ = A⁻¹F.
The loop-by-loop matrix is then repeated for about six iterations until the change in flowrates becomes negligible,
which is exactly the workflow described in the project brief.
Q2. Plug Flow Reactor (PFR) Analysis
The project gives the reaction A + C → B with rate law rB = k·CA·CC, and the PFR model u dCi/dz = ri. The resulting
ODE system is: dCA/dz = -kCA·CC/u, dCC/dz = -kCA·CC/u, and dCB/dz = +kCA·CC/u. The lecture material on
batch/PFR modeling and MATLAB ODE solving is the same style used here.
Given conditions for part 2.1 / 2.2:
Variable Value
CA(0) 1 kmol/m³
CC(0) 0.5 kmol/m³
CB(0) 0 kmol/m³
u 0.5 m/s
k 0.3 m³/kmol·s
L 10 m
Because CA − CC stays constant, CA = CC + 0.5 for part 2.1. That makes the exact isothermal solution one-
dimensional: dCC/dz = -(k/u) CC(CC+0.5). With k/u = 0.6, the closed-form result is CC(z) = 0.25 e^(-0.3z) / [1 − 0.5
e^(-0.3z)], CA(z) = CC(z) + 0.5, and CB(z) = 1 − CA(z).
At z = 10 m, the exact isothermal outlet concentrations are CA = 0.5128, CC = 0.0128, and CB = 0.4872 kmol/m³. A
forward-Euler step size of 0.1 m gives CA ≈ 0.5119, CC ≈ 0.0119, and CB ≈ 0.4881 kmol/m³, which is close to the
exact profile.
Concentration profile for the isothermal case:
Euler error relative to the exact solution:
Selected comparison points for part 2.2:
z (m) CA exact CA Euler CC exact CC Euler CB exact CB Euler
0 1.000000 1.000000 0.500000 0.500000 0.000000 0.000000
2 0.689090 0.683699 0.189090 0.183699 0.310910 0.316301
4 0.588649 0.585072 0.088649 0.085072 0.411351 0.414928
6 0.545048 0.542821 0.045048 0.042821 0.454952 0.457179
8 0.523757 0.522383 0.023757 0.022383 0.476243 0.477617
10 0.512765 0.511920 0.012765 0.011920 0.487235 0.488080
For part 2.3, the project changes CC(0) to 1 kmol/m³ and uses the Arrhenius equation with A = 10⁵ m³/kmol·s, Ea =
40,000 J/mol, and R = 8.314 J/mol·K. At the safety limit T = 400 K, the rate constant is k = 0.5975 m³/kmol·s. With
CA(0) = CC(0) = 1, the solution becomes CA = CC = 1 / [1 + (k/u)z].
To reach 90% conversion, CA must fall from 1 to 0.1 kmol/m³, so the minimum reactor length is L = 7.53 m. This is
the shortest length because the temperature is capped at 400 K, so the maximum allowable rate constant is used.
Non-isothermal design at the 400 K safety limit:
Appendix A. MATLAB code for the PFR plots
% Q2.1 and Q2.2: Isothermal PFR
u = 0.5;
k = 0.3;
h = 0.1;
z = 0:h:10; CA = zeros(size(z)); CC = zeros(size(z)); CB = zeros(size(z));
CA(1) = 1; CC(1) = 0.5; CB(1) = 0; for i = 1:length(z)-1 r = k*CA(i)*CC(i)/u; CA(i+1) =
CA(i) - h*r; CC(i+1) = CC(i) - h*r; CB(i+1) = CB(i) + h*r;
end CC_exact = 0.25*exp(-0.3*z) ./ (1 - 0.5*exp(-0.3*z));
CA_exact = CC_exact + 0.5;
CB_exact = 1 - CA_exact; plot(z,CA_exact,'-o',z,CA,'--',z,CC_exact,'-
o',z,CC,'--',z,CB_exact,'-o',z,CB,'--');
xlabel('Reactor length, z (m)');
ylabel('Concentration (kmol/m^3)');
title('Isothermal PFR: exact vs Euler');
grid on
legend('CA exact','CA Euler','CC exact','CC Euler','CB exact','CB Euler','Location','best')
Appendix B. MATLAB code for the non-isothermal design
A = 1e5;
Ea = 40000;
R = 8.314;
T = 400;
u = 0.5;
k = A*exp(-Ea/(R*T));
z90 = (1/0.1 - 1)/(k/u); z = 0:0.01:10;
CA = 1 ./ (1 + (k/u)*z);
CC = CA;
CB = 1 - CA; plot(z,CA,z,CB,z,CC);
xlabel('Reactor length, z (m)');
ylabel('Concentration (kmol/m^3)');
title('Non-isothermal design at T = 400 K');
grid on
legend('CA','CB','CC','Location','best')
Appendix C. Excel formulas
Use these formulas in Excel and fill downward by row: V = Q/(PI()/4*D^2), Re = rho*D*V/mu, f =
IF(Re<=2100,64/Re,0.25/(LOG10((eD/3.7)+(5.74/Re^0.9))^2)), K = 8*f*L/(PI()^2*g*D^5), hL = K*Q^2, and [Link]/Q =
2*hL/Q.