フィルターのクリア

Solution of Fick's Second Law of Diffusion Equation

161 ビュー (過去 30 日間)
Soumyadeep
Soumyadeep 2023 年 1 月 9 日
編集済み: Torsten 2023 年 1 月 11 日
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
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 ?
Soumyadeep
Soumyadeep 2023 年 1 月 9 日
Exactly. Thats correct.

サインインしてコメントする。

採用された回答

Torsten
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.
  7 件のコメント
Soumyadeep
Soumyadeep 2023 年 1 月 11 日
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
Soumyadeep 2023 年 1 月 11 日
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?

サインインしてコメントする。

その他の回答 (1 件)

J. Alex Lee
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:
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.
  6 件のコメント
Soumyadeep
Soumyadeep 2023 年 1 月 11 日
Yes it is definitely helpful Alex.
Torsten
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 ExchangeMathematics and Optimization についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by