# Simulate Fick's 1st and 2nd laws of mass diffusion.

12 Ansichten (letzte 30 Tage)
JAEYEON am 19 Mai 2024
Beantwortet: Torsten am 19 Mai 2024
Dear a person who is looking this.
I am trying to simulate the mass diffusion based on Fick's laws on 2D plane.
The initial condition is that the concentration on the left wall is 1.
But the result seems to be odd as below.
I guess the problem is from the code for the 2nd law. So I have tried to modify that like
C(i,j) = C(i,j) ...
- dt * ((Jx(i, j) - Jx(i-1,j))/(dx) ...
+ (Jy(i, j) - Jy(i,j-1))/(dy));
But the result was the same.
How can I solve this problem. I really want some help..
Thank you for reading this.
======Here is my code (you can operate by just copy and paste)=======
% Parameters
Nx = 20; % Number of grid points in x-direction
Ny = 20; % Number of grid points in y-direction
D = 0.1; % Diffusion coefficient
dx = 1.0; % Spatial step in x-direction
dy = 1.0; % Spatial step in y-direction
dt = 0.1; % Time step
T = 200.0; % Total simulation time
C = zeros(Nx, Ny); % Initial concentration
Jx = zeros(Nx, Ny);
Jy = zeros(Nx, Ny);
C(:, 1) = 1.0; % Initial condition: concentrated mass at the center
% Time loop
for t = 0:dt:T
% Create a copy of the concentration array to hold the new values
% C_new = C;
% Iterate over each interior point
% flux
% Fick's 1st law: J = -D * del c
for i = 2:Nx-1
for j = 2:Ny-1
Jx(i,j) = - D * (C(i,j) - C(i-1,j))/dx - D * (C(i+1,j) - C(i,j))/dx ;
Jy(i,j) = - D * (C(i,j) - C(i,j-1))/dy - D * (C(i,j+1) - C(i,j))/dy ;
end
end
% concentration
% Fick's 2nd law: partial dc/dt = - div J
for i = 2:Nx-1
for j = 2:Ny-1
C(i,j) = C(i,j) ...
- dt * ((Jx(i+1, j) - Jx(i-1,j))/(2*dx) ...
+ (Jy(i, j+1) - Jy(i,j-1))/(2*dy));
end
end
% Optionally plot the concentration at each time step
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
colormap('jet');
title(sprintf('Time = %.2f', t));
drawnow;
end
% Final plot
slice_view = squeeze(C(:,:)); % View a slice in the middle of the z-plane
imagesc(slice_view);
colorbar;
title('Final Concentration');
xlabel('X');
ylabel('Y');
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Akzeptierte Antwort

Torsten am 19 Mai 2024
The correct update for inner grid points is
C_new(i,j) = C_old(i,j) + dt * (D_x*(C_old(i-1,j)-2*C_old(i,j)+C_old(i+1,j))/dx^2 + D_y*(C_old(i,j-1)-2*C_old(i,j)+C_old(i,j+1))/dy^2 )
And since you don't change the values for C for i=1,i=Nx,j=1,j=Ny, they always have a value of 0. Thus you boundary condition is C = 0 at all the four edges.
##### 0 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Physics finden Sie in Help Center und File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by