Solution of Fick's Second Law of Diffusion Equation
186 views (last 30 days)
Show older comments
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
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.
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
dx2 = (x(2)-x(1))^2
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
and the number of time steps required to get to it is
nSteps = round(tEnd/dt)
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) ...
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;
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;
Is this accurate enough? Well, you'd have to check against better integrators, and decide on the trade-off in runtimes.
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.
Find more on General PDEs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!