Plug Flow Reactor with axial dispersion ODE (error)

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)

Answer 1

[Sam Chak]

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

[Cai Li Song]

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.