I am solving an 1-dimensional diffusion equation with Neumann boundary condition at outlet and constant concentration, C, at the inlet. In the end, I want to observe how the concentration diffuses over time along the x direction.

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}$$

The following is the code written in MATLAB.

function main()
C0 = [500 480 460 440 420 400 380 360 340];
tspan = [0 10];
[t C]  = ode23s(@(t,s) fun(t,s), tspan , C0);

function dC = fun(t,C)
D  = 3000;
retain_rows = 2:8;
tri = -(D/(1/5^2))*full(gallery('tridiag',9,-1,2,-1)); % 5 is del x
tri = tri(retain_rows,:);
dC(1,1) = 0; % constant concentration at inlet
dC(retain_rows,1) = tri*C;
dC(9,1) = -D*((C(9)-C(8))/2); Neumann boundary condition at outlet

I have the solution obtained from ode solver. I'd like to know how to check for mole balance from the solution that is obtained. Also, could someone explain how the above code can be changed to apply constant flux at the inlet?

Any suggestions?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.