Solution of Fick's Second Law of Diffusion Equation

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]

So my guess is you mean

dC/dt = D*(d^2C/dx^2 + d^2C/dy^2 + d^2C/dz^2)

instead of

dC/dt = D*(d^2C/dx^2 + d^2C/dx^2 + d^2C/dx^2)

?

Use the PDE toolbox.

[Soumyadeep]

In PDE Toolbox 2D is available not 3D

[Torsten]

Could you specify your model in more detail ?

Volume on which you want to solve the equation, boundary conditions, possible internal sources, time span over which you want to integrate ...

[Soumyadeep]

The model is a 3D Fick's second law of diffusion model to calculate material concentration gradient locally. For this an experiment is performed for inward diffusion of a material into the sample in gaseous state. The rate of change of concentration is a value available experimentally. So with time, a new concentration is obtained from the experiment. Diffusion coefficient is available in literature. For the purpose of accuracy of results, the diffusion coefficient is evaluated by heating the material to elevated temperature. The scientific result sought is the distribution of concentration along x, y and z axes in the sample, which is why solution of the equation dC/dt = D*(d^2C/dx^2 + d^2C/dy^2 + d^2C/dz^2) is sought. With dC/dt value changing with time, and D remaining constant, the distribution of concentration is the target of the investigation.

[Torsten]

To solve the diffusion equation, you need the dimensions of the volume of integration (lengths in x-, y- and z-direction), D and the boundary conditions for C on the 6 faces of the cuboid. That's what I've been asking for.

[Soumyadeep]

Volume on which I am solving the equation is 100mm*20mm*10mm. Boundary conditions being that material is fixed on all sides, so it is static. There is no internal source, only inward diffusion of a gas takes place. I would like to integrate over 24h, 48h, 72h, and 96h respectively.

[Torsten]

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]

Exactly. Thats correct.

Answer 1

[Torsten]

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.

[J. Alex Lee]

just a comment that, while straightforward in concept, implementation can be very important in 3D for reasonably resolved meshes. once you have C, you could also entertain the possibility of Euler steps instead of ode15s, which could save a lot of time and/or effort.

[Torsten]

@J. Alex Lee

Explicit Euler for a stiff system of ODEs over a time span of 96 h ? I think it's not a good idea.

Of course, this equation in 3d can cause memory problems with ODE15S. Thus an iterative method with implicit Euler could be an alternative.

[J. Alex Lee]

@Torsten, I am not seeing the source of stiffness - the only timescale I see is through D since boundary conditions are fixed. Certainly CFL condition (sorry if I've got the name wrong) about time steps and mesh sizes will apply...

[Torsten]

@J. Alex Lee

Heat conduction or diffusion problems are always stiff.

[J. Alex Lee]

Perhaps we are not agreed on what "stiff" means...

But anyway the important part of my comment was about how to construct the coefficient matrix efficiently for 3D cases with proper use of sparse commands.

[Soumyadeep]

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]

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?

Answer 2

[J. Alex Lee]

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]

I am thankful to you for this answer.

[Torsten]

@J. Alex Lee

That's what I meant by "stiff":

N = 50;
F = spdiags(-2*ones(N,1),0,N,N) ...
    + spdiags(ones(N,1),+1,N,N) ...
    + spdiags(ones(N,1),-1,N,N);
S = eig(F);
Smax = max(abs(S))

Smax = 3.9962

Smin = min(abs(S))

Smin = 0.0038

Smax/Smin

ans = 1.0535e+03

[Soumyadeep]

I thank you heartily for this answer.

[J. Alex Lee]

@Torsten, Ok I see where you're coming from now. I take your point that "methods of line are stiff" in the sense that the system of equations resulting from the space discretization is stiff, as measured by ratio of eigenvalues in the resulting matrix F (or probably more to the point A). I will assert that it isn't about "conduction or diffusion" problems though.

With my background, I tend to think more of stiffness resulting from the physics, e.g., reaction diffusion problems where reaction rate and diffusion rate scales are disparate. You can also introduce stiffness in forcing terms like boundary conditions and source terms. Stiffness in these senses seem to me less straightforward than that introduced by discretization in this simple case, and present challenges beyond what we are talking about here, which is essentially about tradeoff between accuracy, time step size, and spatial resolution.

To the point of this discussion, every choice in numerical methods represents a trade-off, so I think it is important to lay out first "what the investigator cares about". What's missing, and I'm trying to supply, is a question about runtime. Unquestionably, more sophisticated integrators can be much more efficient with time step size given a fixed spatial discretization and diffusion coefficient for some target accuracy. What's unclear to me is how the trade-off plays out with run time. To illustrate (and disclaimer, I have no idea how the numbers actually come out, and surely it depends on a lot of details): you may be able to take 10 times larger time steps, but if each time step takes 10 times longer to calculate, what have you gained?

But I guess my point is look how trivial it is to set up an explicit Euler, is it good enough for your needs?

[Soumyadeep]

Yes it is definitely helpful Alex.

[Torsten]

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.