How to implemented boundary condition in solving convection diffusion equation in cylindical coordinate

For C_c at r = R_c, you can consider that node as C(length(r1),j)

For C_a at r = R_c, you can consider that node as C(length(r1)+1,j)

Then, via central difference,

dC_c/dr (at r = R_c) can be written as (C(length(r1)+1,j)-C(length(r1)-1,j))/(del_r)

Similarly, dC_a/dr (at r = R_c) can be written as (C(length(r1)+2,j)-C(length(r1),j))/(del_r)

Try to suitably implement the above in the condition mass flux continuity condition at r = R_c, and check what happens. Also, ensure that the other conditions are correctly implemented, and then run the code.

I want to know how to write this in matlab so that I will not get error. I have aslo attached my matlab file here

clc;
close all;
clear all;

% Constants for Diffusion and Convection
D = 1e-7;  % Diffusion coefficient in m^2/s
dz = 0.0002;  % Axial step size (z-direction)
dr = 2e-6;  % Radial step size (r-direction)
delta = 1e-4;  % Maximum radial distance (domain size)

% Radial and Axial Grids
r = 0:dr:delta;  % Radial nodal points
z = 0:dz:0.02;   % Axial nodal points

% Velocity in the z-direction
% Initializing the velocity profile in cylindrical coordinates
Rc = 1e-4; % Core radius (in meters)
Ra = 2e-4; % Annular radius (in meters)
mu_c = 0.001; % Viscosity of the core fluid (Pa·s)
mu_a = 0.003; % Viscosity of the annular fluid (Pa·s)
dP=707.71;% Axial pressure gradient (Pa/m)
L= .05;  
Dhc=(2*Rc)

for i= 1:length(r)
Vc_z(i) = (1/4) *( -dP/L) .* ((Rc^2 - r(i)^2) ./ mu_c + (Ra^2 - Rc^2) / mu_a);
end
Vc_avg =(1 / 4) * (-dP/L) * ((1 / mu_c) * (2 / 3) *(Rc^2 )+ ((Ra^2 - Rc^2) / mu_a));
Vc= (Vc_z(i))/(Vc_avg);
rho=1000;
Re_c=((Vc_avg*rho*Dhc)/ mu_c);
Re_c = 10.62;
Sc_c= (mu_c/(rho*D));
Sc_c=500;
A=(mu_c/(rho*D))*((Vc_avg*rho*Dhc)/ mu_c)*(Ra/Dhc)*(Ra/L)

% Initializing concentration matrix
%C(2:length(r)-1) = 0
%C = zeros(length(r), length(z));

% Boundary Conditions
%C(:, 1) = 0;       % C = 0 at z = 0 for all r
%C(1, :) = 1;       % C = 1 at r = 0 (centerline)
%C(end, :) = 0;     % C = 0 at r = delta (outer boundary)

%C (i,j) is the conc. at ith value of r and jth value of z
for i =1:length(r)
    C(i,1)=0;
end

for j=1:length(z)
    C(1,j)=1;
    C(length(r),j)=0;
end 

% Precomputing coefficients for the finite difference scheme
for i = 2:(length(r)-1)
    %update concentration at each interior node
    a(i) = (((2 * dz) / (A * Vc * dr^2)) + ((dz) / (A * Vc * r(i) * dr))); % For C(i+1,j)
    b(i) = ((1-(4 * dz ))/ (A * Vc * dr^2));                              % For C(i,j)
    c(i) = (2 * dz) / (A * Vc * dr^2) - (dz) / (A * Vc * r(i) * dr); % For C(i-1,j)
end

% Time-stepping loop for axial propagation
 for i = 2:(length(r)-1)
    for j = 1:(length(z))
        % Applying the finite difference equation in cylindrical coordinates
        %C(i, j+1) = b(i) * C(i, j) + a(i) * C(i+1, j) + c(i) * C(i-1, j);
        C(i, j+1) = a(i) * C(i+1, j) + b(i) * C(i, j)+ c(i) * C(i-1, j);
    end
end

% Compute CFL number (optional for stability check)
CFL = (D * dz) / (dr*dr);

% Plotting the concentration profile at the last axial point
%plot(r, C(:, end), 'r');
plot(r, C(:,length(z)), 'r')
xlabel('Radial Position (r)');
ylabel('Concentration (C)');
title('1-D Convection-Diffusion in Cylindrical Coordinates');
legend('Concentration Profile');

[Torsten]

What is the PDE you are trying to solve and what are your boundary conditions ?

What is C_c and what is C_a ?