# Solution of Fick's Second Law of Diffusion Equation

152 views (last 30 days)
Soumyadeep on 9 Jan 2023
Edited: Torsten on 11 Jan 2023
I require the code to solve the equation ∂C/∂t=D.[(∂^2 C)/(∂x^2 )+(∂^2 C)/(∂x^2 )+(∂^2 C)/(∂x^2 )]
I have value for the diffusion coefficient D, I also have value for the rate of variation of concentration ∂C/∂t, all I need is the distribution of concentration C along x, y and z axes respectively.
Torsten on 9 Jan 2023
There is no internal source, only inward diffusion of a gas takes place.
So the boundary condition is C = C_gas on all six faces of the cuboid and initial concentration is C = 0 within the cuboid ?
Soumyadeep on 9 Jan 2023
Exactly. Thats correct.

Torsten on 9 Jan 2023
Edited: Torsten on 9 Jan 2023
Then discretize d^2C/dx^2, d^2C/dy^2 and d^2C/dz^2 and use the method of lines to integrate
dC(i,j,k)/dt = D*((C(i+1,j,k)-2*C(i,j,k)+C(i-1,j,k))/dx^2 + (C(i,j+1,k)-2*C(i,j,k)+C(i,j-1,k))/dy^2 + (C(i,j,k+1)-2*C(i,j,k)+C(i,j,k-1))/dz^2)
(2<=i<=nx-1 , 2<=j<=ny-1, 2<=k<=nz-1)
with MATLAB's ODE15S.
Soumyadeep on 11 Jan 2023
Is there a way to write this code like this?
for j>=2, i<=nx-1
for j>=2, j<=ny-1
for k>=2, k<=ny-1
dC(i,j,k)/dt = D*((C(i+1,j,k)-2*C(i,j,k)+C(i-1,j,k))/dx^2 + (C(i,j+1,k)-2*C(i,j,k)+C(i,j-1,k))/dy^2 + (C(i,j,k+1)-2*C(i,j,k)+C(i,j,k-1))/dz^2)
end
end
end
If I want to write this code like this then what do I need to initialise?
Soumyadeep on 11 Jan 2023
I wanted to obtain the 3d discretisation of x, y, z axis together to make the 3d distribution of concentration. Is there a way to write this code like this? For that how to I proceed? I intend to solve it in form of loop, such that i provide the initial c values, the loop runs and calculates all points and I plot the graph at the end. Is there a way to make it that way?

### More Answers (1)

J. Alex Lee on 10 Jan 2023
Edited: J. Alex Lee on 10 Jan 2023
I agree with the general idea in Torsten's answer using FD and method of lines. Specifically what ode solver to use on the discretized system I think depends. Granted explicit Euler is not generally practical for reasons of accuracy in addition to stability in general, but this answer is just to close the loop against Torsten's categorical objection.
With this simple problem you can find stable sets of discretization for specific values of diffusion coefficient and box dimensions, and you can then decide if those discretizations are acceptable for you w.r.t. practicality of number of steps you need to take and accuracy.
Explicit Euler should be stable as long as the condition is met:
Here is an example in 1D to illustrate: vanilla diffusion on unit domain with both boundaries fixed to unity, and initial condition of zero everywhere else.
In 1D with space indexed by subscript i and time indexed by superscript, to keep notation simple, for the generic internal node i, in a way to emphasize that the finite difference coefficients can be pulled outside, explicit Euler is, after solving for the "next" time step:
Let's put our units in terms of mm and hr.
For discretization of the 1D domain 100mm long (the first dimension of the original problem) into N nodes
N = 50;
x = linspace(0,100,N)'; % mm
So is
dx2 = (x(2)-x(1))^2
dx2 = 4.1649
The longest time I'm interested in is
tEnd = 96; % hours
Then as long as I choose the pair of D and such that the above condition is met, we should be good to go. Physically, we can expect the concentration to approach some significant fraction of the steady state by 96 hours if the diffusion coefficient is
D = 30; % mm^2/hr
and to satisfy stability the min timestep we need is
dt = dx2/D * 0.5
dt = 0.0694
and the number of time steps required to get to it is
nSteps = round(tEnd/dt)
nSteps = 1383
This is trivial, but of course the actual numbers depend on the actual value of D, and the desired spatial resolution of the box, and some tolerance on accuracy.
So to solve the problem, construct the finite difference coefficient matrix
F = spdiags(-2*ones(N,1),0,N,N) ...
+ spdiags(ones(N,1),+1,N,N) ...
+ spdiags(ones(N,1),-1,N,N);
The simple fixed (Dirichlet) BCs just mean zero out the coefficients for the end nodes
F(1,1:2) = 0;
F(N,end-1:end) = 0;
Initial condition:
c = zeros(N,1);
c(1) = 1;
c(N) = 1;
execute and plot 10 intermediate steps
figure(1); hold on;
A = F*dt*D/dx2;
for k = 1:nSteps
c = c + A*c;
if mod(k,round(nSteps/10))==0
plot(x,c,'.-')
end
end
Is this accurate enough? Well, you'd have to check against better integrators, and decide on the trade-off in runtimes.
Soumyadeep on 11 Jan 2023
Yes it is definitely helpful Alex.
Torsten on 11 Jan 2023
Edited: Torsten on 11 Jan 2023
I will assert that it isn't about "conduction or diffusion" problems though.
Each system of ODEs that contains a discretization of 2nd order derivatives is stiff. The reference examples are diffusion and heat conduction problems.
Concerning the method to use in the above case:
The fastest solution will be to supply the Jacobian for ODE15S directly in sparse form as you did it for the 1d-case, tell the solver that it is constant over the complete time span and take advantage of the adaptive time stepping of the integrator.
But of course one will first have to read a little about how to supply a constant Jacobian in sparse form for ODE15S in the options structure.

### Categories

Find more on Mathematics and Optimization in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!