Tag: batch reactor

用数值方法解决常微分或者偏微分方程要比分析解要更加容易理解和处理，而matlab的优势在于拥有强大的数值处理能力，这个学期有reaction engineering的课程，老师要求使用matlab来解决一些反应器的问题，由于之前一直用matlab做一些简单的计算，从来没有在里面写过脚本或者方程，所以第一次写起来还是有些难度，但好在matlab有比较详细的说明文档，最终还是能够写出来一些东西。下面的脚本就是处理non-isothermal Batch Reactor的温度和产率随温度变化的趋势。

%% Script for Problem 6.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Constants : k_m, T_m, Ea_R, Qloss, HR, VR
%% Units: cal, gmol, min, K, L
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k_m = 0.01725;   %% L/mol.min
T_m = 300;       %% K
Ea_R = 2660;     %% K
Qloss = -41700;  %% cal/min
HR = -10000;     %% cal/mol A
VR = 1200;       %% L

time(1) = 0;
r = 0;
Ca(1) = 2;
dT_dt = 0;
T(1) = 300;
k = 0;
x(1) = 0;
i = 1;

while (T(i) > 299)
    k = k_m * exp(-Ea_R * (1/T(i) - 1/T_m));
    r = k * Ca(i);
    Ca(i+1) = Ca(i) + 0.1 * (-r);
    dT_dt = (-HR * r * VR + Qloss)/(VR * 2 * Ca(1) * 20);
    T(i+1) = T(i) + 0.1 * dT_dt;
    time(i+1) = time(i) + 0.1;
    x(i+1) = 100 * (1 - Ca(i) / Ca(1));
    i = i + 1;
end
subplot(2,1,1); plot(time, T)
title('T vs time')
xlabel('time/min')
ylabel('T/Kelvin')
subplot(2,1,2); plot(time, x)
title('Conversion vs time')
xlabel('time/min')
ylabel('x')
Tmax = max(T)

处理之后的趋势图