필터 지우기
필터 지우기

Plug Flow Reactor with axial dispersion ODE (error)

조회 수: 8 (최근 30일)
Cai Li Song
Cai Li Song 2024년 3월 27일
댓글: Cai Li Song 2024년 3월 27일
I am trying to solve mass balance with axial dispersion (governing equations in the file attached) related to methane conversion to hydrogen in a plug flow reactor (or more precisely a bubble column in axial dispersion).
Here's my matlab code:
function dX_CH4dz = SBC_ADM(z, X_CH4)
% Define constants
T = 1040 + 273.15;
Tm = 1020 + 273.15;
P = 101325;
n_CH40 = 10;
n_H20 = 0;
n_I0 = 0;
MW_ni = 0.05869;
MW_bi = 0.20898;
MW_CH4 = 0.01604;
MW_H2 = 0.002016;
MW_I = 0.02802; % Nitrogen
R = 8.314;
g = 9.81;
D = 0.0127;
% Define parameters (as a function of conversion)
sigma = 0.1; % surface tension of liquid
rho_ni = 9207 - 0.7818 * T;
rho_bi = 10698 - 1.2064 * T;
n_T0 = n_CH40 + n_H20 + n_I0;
MW = ((n_CH40 * (1 - X_CH4)) / (n_T0 + (n_CH40 * X_CH4))) * MW_CH4 + ((n_H20 + 2 * n_H20 * X_CH4) / (n_T0 + (n_CH40 * X_CH4))) * MW_H2 + (n_I0 / (n_T0 + (n_CH40 * X_CH4))) * MW_I;
rho_l = (0.27 * MW_ni + 0.73 * MW_bi) / ((0.27 * MW_ni / rho_ni) + (0.73 * MW_bi / rho_bi)); % density of liquid
mu_l = 1.7e7 * (rho_l^(2/3)) * (Tm^0.5) * (MW^(-1/6)) * exp(2.65 * (Tm^1.27) / R * ((1 / T) - (1 / Tm))); % dynamic viscosity of liquid
rho_g = P * MW / (R * T); % density of gas
N_mu = mu_l ./ ((rho_l * sigma * sqrt(sigma ./ (g * (rho_l - rho_g))))^0.5); % dimensionless viscosity number range
D_H_star = D ./ (sqrt(sigma ./ (g * (rho_l - rho_g))));; % dimensionless hydraulic equivalent of diameter of the column
% Define correlation
V_g = (n_T0 + (n_CH40 * X_CH4)) * R * T / P; % volumetric flow rate
U_g = 4 * V_g / (pi * D^2); % gas velocity
U_g_plus = U_g ./ ((sigma * g * (rho_l - rho_g) ./ (rho_l)^2))^0.25; % dimensionless gas velocity
% Calculate V_d_plus based on the provided conditions
V_d1_plus = 0; % Initialize V_d1_plus to a default value
if U_g_plus >= 0.5
V_d2_plus = 0;
if N_mu <= 2.25e-3
if D_H_star <= 30
V_d1_plus = 0.0019 * D_H_star^0.809 * (rho_g / rho_l)^(-0.157) * N_mu^(-0.562);
elseif D_H_star > 30
V_d1_plus = 0.030 * (rho_g / rho_l)^(-0.157) * N_mu^(-0.562);
end
elseif N_mu > 2.25e-3
V_d1_plus = 0.92 * (rho_g / rho_l)^(-0.157);
end
elseif U_g_plus < 0.5
% Initialize variables
alpha0 = 0.1;
V_d2_plus = sqrt(2) * (1 - alpha0)^(1.75);
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d2_plus);
% Iteratively calculate alpha and Vgj_plus until tolerance is met
while abs(alpha - alpha0) > 0.001
% Update previous values
alpha0 = alpha;
V_d2_plus = sqrt(2) * (1 - alpha0)^(1.75);
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d2_plus);
end
V_d2_plus = sqrt(2) * (1 - alpha)^(1.75);
end
V_d_plus = (V_d2_plus * exp(-1.39 * U_g_plus) + V_d1_plus * (1 - exp(-1.39 * U_g_plus))); % dimensionless bubble drift velocity
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d_plus); % cross-sectional area
D_ax = 50 * ((U_g ./ alpha)^3) * D^1.5; % axial dispersion coefficient
R_c = 4.73e4 * exp(-208000 / (R * T)) * ((n_CH40 * (1 - X_CH4)) / V_g); % reaction rate
% Methane mass balance equation
% Second order ODE (X_CH4 as a function of Z)
dX_CH4dz = [
X_CH4(2); % First element of dX_CH4dz
(U_g ./ (D_ax * alpha)) .* X_CH4(2) - (alpha * R_c * MW * V_g) ./ (D_ax * alpha * n_CH40)
];
end
L_c = 1; % length of catalyst
[z, X_CH4] = ode45(@SBC_ADM, [0 L_c], [0; 0]);
When I run the code, i get this error message:
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in SBC_ADM (line 73)
dX_CH4dz = [
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in SBC_ADM_ode (line 4)
[z, X_CH4] = ode45(@SBC_ADM, [0 L_c], [0; 0]);
I suspect the error is in the writing of the second order ODE, please assist and give suggestion how I can rectify this, thank you (Note: I don't have symbolic matlab toolbox, so nothing of SYMS please)

채택된 답변

Sam Chak
Sam Chak 2024년 3월 27일
There doesn't seem to be any coding issue with the 2nd-order ODE 'dX_CH4dz' when I substituted all parameters with unity scalars:
dX_CH4dz = [X_CH4(2); % First element of dX_CH4dz
(U_g./(D_ax*alpha)).*X_CH4(2) - (alpha*R_c*MW*V_g)./(D_ax*alpha*n_CH40)];
Therefore, the matrix dimension mismatch error must be triggered by some parameters {U_g, D_ax, alpha, R_c, MW, V_g, n_CH40} that are in matrix form. If you trace them individually, you'll notice that the issue arises from the state vector term, X_CH4, which appears in these three parameters: MW, V_g, and R_c. The vector contains two state variables, X_CH4(1) and X_CH4(2).
I have made a 'temporary fix' in the code below by modifying all instances of X_CH4 to X_CH4(1), but it is important for you to review and ensure that the equations for MW, V_g, and R_c are corrected as described in the PDF.
L_c = 1; % length of catalyst
[z, X_CH4] = ode45(@SBC_ADM, [0 L_c], [0; 0]);
plot(z, X_CH4), grid on
xlabel('z'), ylabel('X_{CH4}')
legend({'$X_{\mathrm{CH4}}(1)$', '$X_{\mathrm{CH4}}(2)$'}, 'interpreter', 'LaTeX')
function dX_CH4dz = SBC_ADM(z, X_CH4)
% Define constants
T = 1040 + 273.15;
Tm = 1020 + 273.15;
P = 101325;
n_CH40 = 10;
n_H20 = 0;
n_I0 = 0;
MW_ni = 0.05869;
MW_bi = 0.20898;
MW_CH4 = 0.01604;
MW_H2 = 0.002016;
MW_I = 0.02802; % Nitrogen
R = 8.314;
g = 9.81;
D = 0.0127;
% Define parameters (as a function of conversion)
sigma = 0.1; % surface tension of liquid
rho_ni = 9207 - 0.7818 * T;
rho_bi = 10698 - 1.2064 * T;
n_T0 = n_CH40 + n_H20 + n_I0;
MW = ( (n_CH40*(1 - X_CH4(1))) / (n_T0 + (n_CH40*X_CH4(1)))) * MW_CH4 + ((n_H20 + 2 * n_H20 * X_CH4(1)) / (n_T0 + (n_CH40 * X_CH4(1)))) * MW_H2 + (n_I0 / (n_T0 + (n_CH40 * X_CH4(1)))) * MW_I;
rho_l = (0.27 * MW_ni + 0.73 * MW_bi) / ((0.27 * MW_ni / rho_ni) + (0.73 * MW_bi / rho_bi)); % density of liquid
mu_l = 1.7e7 * (rho_l^(2/3)) * (Tm^0.5) * (MW^(-1/6)) * exp(2.65 * (Tm^1.27) / R * ((1 / T) - (1 / Tm))); % dynamic viscosity of liquid
rho_g = P*MW/(R*T); % density of gas
N_mu = mu_l ./ ((rho_l * sigma * sqrt(sigma ./ (g * (rho_l - rho_g))))^0.5); % dimensionless viscosity number range
D_H_star= D ./ (sqrt(sigma ./ (g * (rho_l - rho_g))));; % dimensionless hydraulic equivalent of diameter of the column
% Define correlation
V_g = (n_T0 + n_CH40*X_CH4(1))*R*T/P; % volumetric flow rate
U_g = 4*V_g/(pi*D^2); % gas velocity
U_g_plus= U_g/((sigma*g*(rho_l - rho_g) ./ (rho_l)^2))^0.25; % dimensionless gas velocity
% Calculate V_d_plus based on the provided conditions
V_d1_plus = 0; % Initialize V_d1_plus to a default value
if U_g_plus >= 0.5
V_d2_plus = 0;
if N_mu <= 2.25e-3
if D_H_star <= 30
V_d1_plus = 0.0019 * D_H_star^0.809 * (rho_g / rho_l)^(-0.157) * N_mu^(-0.562);
elseif D_H_star > 30
V_d1_plus = 0.030 * (rho_g / rho_l)^(-0.157) * N_mu^(-0.562);
end
elseif N_mu > 2.25e-3
V_d1_plus = 0.92 * (rho_g / rho_l)^(-0.157);
end
elseif U_g_plus < 0.5
% Initialize variables
alpha0 = 0.1;
V_d2_plus = sqrt(2) * (1 - alpha0)^(1.75);
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d2_plus);
% Iteratively calculate alpha and Vgj_plus until tolerance is met
while abs(alpha - alpha0) > 0.001
% Update previous values
alpha0 = alpha;
V_d2_plus = sqrt(2) * (1 - alpha0)^(1.75);
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d2_plus);
end
V_d2_plus = sqrt(2) * (1 - alpha)^(1.75);
end
V_d_plus = (V_d2_plus * exp(-1.39 * U_g_plus) + V_d1_plus * (1 - exp(-1.39 * U_g_plus))); % dimensionless bubble drift velocity
alpha = U_g_plus / ((1.2 - 0.2 * sqrt(rho_g / rho_l)) * U_g_plus + V_d_plus); % cross-sectional area
D_ax = 50*((U_g/alpha)^3)*D^1.5; % axial dispersion coefficient
R_c = 4.73e4 * exp(-208000 / (R * T)) * ((n_CH40 * (1 - X_CH4(1))) / V_g); % reaction rate
% Methane mass balance equation
% Second order ODE (X_CH4 as a function of Z)
dX_CH4dz = [X_CH4(2); % First element of dX_CH4dz
(U_g./(D_ax*alpha)).*X_CH4(2) - (alpha*R_c*MW*V_g)./(D_ax*alpha*n_CH40)];
end
  댓글 수: 1
Cai Li Song
Cai Li Song 2024년 3월 27일
Dear @Sam Chak,
Thanks for the explanation, now I understand I should have defined X_CH4 as X_CH4(1) as corrected. This is very helpful for me to understand my error, appreciate your effort in debugging. Cheers.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

제품

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by