I tried to do an advection-diffusion equation but I'm not sure if it's a good idea.

I tried to solve the advection-diffusion equation from the equation dc/dt - v(dc/dy) = D(d^2c/dx^2), which is the equation for methylene blue degradation in microchannel reactor, using the tridiagonal matrix method. However, the graph of ct/c0 values that I got (I started by guessing the values of k,K) does not look very similar to the experiment. Can you please check where I went wrong, or should I not have used the tridiagonal matrix method in the first place?

This is the code I wrote.

clear; close all; clc;
%% === Parameters ===
Lx = 1e-3;   % x-direction length (m)
Ly = 0.04;   % y-direction length (m)
Nx = 40;     % number of x points
Ny = 80;    % number of y points
dx = Lx / (Nx-1);
dy = Ly / (Ny-1);
D = 1e-9;    % diffusion coefficient (m^2/s)
C0 = 6.253e-3;  % initial concentration (mol/m3)
k = 2.25e-2;  % reaction rate (m/s)
K_ads = 70;    % adsorption equilibrium constant (m^3/mol)
rho = 1000;    % density (kg/m3)
g = 9.81;      % gravity (m/s2)
mu = 1e-3;     % dynamic viscosity (Pa.s)
H = 1e-3;      % characteristic length in velocity equation (m)
dt = 5;      % time step (s)
Tmax = 180*60; % total time (s)
Nt = round(Tmax/dt);
x = linspace(0, Lx, Nx);
y = linspace(0, Ly, Ny);
[X, Y] = meshgrid(x, y);
%% === Experimental data ===
t_exp = [0 30 60 90 120 150 180]*60; % sec
ca_ratio_exp = [1.00 0.6768 0.5753 0.4924 0.4142 0.2742 0.2264];
%% === Velocity Profile ===
% v(x) = (rho * g)/(2 * mu) * (H * x - x^2)
v = (rho * g)/(2*mu) * (H*X - X.^2); % size Ny x Nx
%% === Initial Condition ===
C = C0 * ones(Ny, Nx); % C(y,x)
%% === Preallocate ===
idx = @(i,j) (j-1)*Nx + i;  % helper for linear indexing
N = Nx*Ny;
A = sparse(N,N);
b = zeros(N,1);
rx = D*dt/dx^2;
Ct_out = zeros(Nt+1,1);
time = (0:Nt)*dt;
% Store initial outlet concentration
Ct_out(1) = mean(C(end,:));
%% === Time Loop ===
for tstep = 1:Nt
    C_old = C;
    Cvec = reshape(C_old,[],1);
    A(:,:) = 0; b(:) = 0;
    for j = 1:Ny
        for i = 1:Nx
            n = idx(i,j);
            vx = v(j,i); % v(x) depends on x only
            if i==1  % x=0 Dirichlet
                A(n,n) = 1;
                b(n) = C0;
            elseif j==1 % y=0 Dirichlet
                A(n,n) = 1;
                b(n) = C0;
            elseif i==Nx % x=Lx Robin condition (reaction)
                nm1 = idx(i-1,j);
                A(n,n) = 1 + D*dt/dx + dt * k*K_ads/(1 + K_ads*C_old(j,i)); % ใช้ approx r_new
                A(n,nm1) = -D*dt/dx;
                b(n) = C_old(j,i);
            elseif j==Ny % y=Ly (outlet)
                jm1 = idx(i,j-1);
                im1 = idx(i-1,j);
                ip1 = idx(i+1,j);
                % Upwind scheme at outlet
                A(n,n)  = 1 + 2*rx + vx*dt/dy;
                A(n,im1) = -rx;
                A(n,ip1) = -rx;
                A(n,jm1) = -vx*dt/dy;
                b(n) = C_old(j,i);
            else % Interior points
                jm1 = idx(i,j-1);
                jp1 = idx(i,j+1);
                im1 = idx(i-1,j);
                ip1 = idx(i+1,j);
                A(n,n)   = 1 + 2*rx + vx*dt/dy;
                A(n,im1) = -rx;
                A(n,ip1) = -rx;
                A(n,jm1) = -vx*dt/dy;
                b(n) = C_old(j,i);
            end
        end
    end
    % Solve
    Cvec_new = A\b;
    C = reshape(Cvec_new,Ny,Nx);
    % Store outlet
    Ct_out(tstep+1) = mean(C(end,:));
end
%% === Plot Results ===
% Concentration profile at final time
figure;
contourf(Y,X,C,30,'LineColor','none');
colormap(jet);
colorbar;
xlabel('y (m)');
ylabel('x (m)');
title('Concentration Profile at 180 min');
set(gca,'FontSize',14);

% Velocity profile
figure;
contourf(Y,X,v,30,'LineColor','none');
colormap(jet);
colorbar;
xlabel('y (m)');
ylabel('x (m)');
title('Velocity Profile');
set(gca,'FontSize',14);

% Outlet Concentration
figure;
plot(time/60,Ct_out/C0,'b','LineWidth',2, 'DisplayName','Simulation');
hold on;
plot(t_exp/60,ca_ratio_exp,'ro','MarkerSize',8,'DisplayName','Experimental');
xlabel('Time (min)');
ylabel('C_t / C_0');
title('Outlet Concentration vs Time');
grid on;
set(gca,'FontSize',14);

Answer 1

[Rahul]

Hi Sasiwimol,

I understand that you are trying to simulate a degradation iprocess in microchannel reactor using the advection-diffusion equation in MATLAB, and are trying to match your results with experimental data. From the provided MATLAB code, it seems that you are trying to implement a 2D transport model where advection occurs along the y-axis and diffusion along the x-axis.

The discretized form of the governing equation has been approached using finite differences i.e., an upwind scheme for advection and an implicit treatment for diffusion, which is an ideal approach for problems involving directional transport and stiff reactions.

However, here are a few changes that might help improve the numerical behavior and physical representation:

1. Incorporating the Local Velocity Profile

Since the velocity is position-dependent in the x-direction, the advection term can benefit from being evaluated using the local velocity value rather than a uniform one. For example, the following code helps capture the effect of varying flow more accurately across the domain.:

vx = v(j, i);  % use local velocity at this point
A(n, jm1) = -vx * dt / dy;
A(n, j)   = 1 + vx * dt / dy;

2. Refining the Boundary Condition at the Catalyst Wall

The implementation at the wall where the reaction occurs (typically the boundary at 'x=Lx​') may be improved by clearly combining both diffusion and reaction kinetics. A representative form of the same can be implemented as shown below:

A(n, j)   = 1 + D * dt / dx + dt * k * K_ads / (1 + K_ads * C_old(j, i));
A(n, jm1) = -D * dt / dx;

This ensures that the wall flux incorporates both the transport into the surface and the adsorption-based reaction.

3. Ensuring Numerical Stability (CFL Condition)

It might be helpful to verify whether the chosen time step satisfies the following CFL condition:

v⋅(dt/dy) < 1

If the CFL number exceeds 1, reducing the time step (e.g., to 1 second or smaller) could improve the accuracy and reduce artificial diffusion in the results.

4. Matching Simulation to Experimental Trends

To better align with the experimental 'Ct/C0'​ values, it may be beneficial to explore a parameter sweep for the reaction rate constant 'k' and the adsorption coefficient 'K_ads'. Here is how you can perform a simple parameter scan:

k_values = [0.01, 0.05, 0.1];
K_ads_values = [0.1, 0.5, 1.0];
% Loop through combinations and compare Ct/C0 results

Alternatively, parameter fitting techniques using MATLAB functions like 'lsqcurvefit' or 'fmincon' can help in estimating optimal values by minimizing the deviation from experimental data.

5. Using 'pdepe' for Simpler Setup

If the 2D domain can be approximated as a 1D problem along the flow direction (e.g., using an effective diffusion coefficient), MATLAB’s 'pdepe' solver might offer a more concise and stable framework for solving such parabolic PDEs with boundary conditions representing surface reactions. This could also make it easier to implement and test multiple parameter sets.

For more information regarding the usage of 'pdepe' function, you can refer to the following documentation link:

• Solve 1-D parabolic and elliptic PDEs using 'pdepe' function: https://www.mathworks.com/help/matlab/ref/pdepe.html