Hi @Cai Li Song
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
