Solution of Fick's Second Law of Diffusion Equation
126 次查看(过去 30 天)
显示 更早的评论
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.
8 个评论
Torsten
2023-1-9
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 ?
采纳的回答
Torsten
2023-1-9
编辑:Torsten
2023-1-9
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.
更多回答(1 个)
J. Alex Lee
2023-1-10
编辑:J. Alex Lee
2023-1-10
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1257627/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1257632/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1257637/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1257642/image.png)
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) ...
+ 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.
6 个评论
Torsten
2023-1-11
编辑:Torsten
2023-1-11
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!