I am trying to solve PDE using FDM, with initial and boundary conditions. I have written a code for that, can you please check whether I have correctly discretize it or not?

D = 1.5e-6; %m^2/s ; diffusion coefficient

C0= 50; % molecules/ml;

T = 60; % Total time to simulate

R1 = 5;

% Spatial grid

num_points_R = 60;

dr_R = R1 / num_points_R;

r_R1 = 0:dr_R:R1;

% Temporal grid

Nt = 60;

dt = T /Nt;

% Create concentration arrays

Ce_R = zeros(num_points_R+1, Nt+1);

% Initialize concentration profiles at t=0, Ce_R=C0 i.e

Ce_R(:, 1) = C0;

% Set boundary condition at r=0 for Ce_R

Ce_R(1, :) = 50;

for t = 2:Nt

% Update the concentration in region Re using FDM (Only radial loop)

for r = 2:num_points_R

d2Ce_dr2 = (Ce_R(r + 1, t-1) - 2 * Ce_R(r, t-1) + Ce_Re(r - 1, t-1)) / dr_R^2;

Ce_R(r, t) = Ce_R(r, t-1) + D * dt * (d2Ce_dr2);

end

end

Answer 1

[SANKALP DEV]

Hi Aiman,

I understand that you are seeking assistance in verifying if your code is correctly discretized.

On investigating the attached code, it appears that you have correctly discretized the PDE using the Finite Difference Method (FDM) for the given diffusion problem.

The code you provided correctly implements the diffusion equation in cylindrical coordinates, '∂C/∂t = D * (∂²C/∂r²)'.

The line of code,

Ce_R(r, t) = Ce_R(r, t-1) + D * dt * d2Ce_dr2; 

corresponds to the finite difference approximation of the diffusion equation. It calculates the concentration at the current time step ‘t’ based on the concentration at the previous time step ‘t-1’, the diffusion coefficient ‘D’, the time step size ‘dt’, and the second derivative of the concentration with respect to ‘r’.

Hence, the code is exhibiting the intended behaviour.

Hope this answers your question.

[Aiman]

Thanks very much SANKALP DEV for for confirming the correctness of the code's discretization.