% Transient 1D diffusion, uniform layer, one side impermeable, % as in a medical drug patch, for example. % Demonstrate a finite-difference method of numerical solution % by Richard K. Herz, www.ReactorLab.net % SUGGSTIONS: Vary n and lambda values to get good approximation to % analytical solution for this case, which has an analytical solution. % Numerical methods are most useful for problems not having available % analytical solution, e.g., conc at layer face varies unpredictably % or complex variation of D or initial conc within layer. % OTHER THINGS TO DO: compute rate of diffusion out of layer vs. time; % compute fraction diffusing component left in layer vs. time. clear all fprintf('------------------- \n') % separates runs in command window % BEGIN USER INPUTS D = 1.0e-5; % (m2/h), diffusion coefficient cI = 1; % (mol/m3), initial concentration c0 = 0; % new conc at layer "face" x = 0 for t>0 L = 0.01; % (m), length of layer, impermeable on one side, "interior" n = 20; % integer number divisions of L tfinal = 10; % (h), final time % lambda = D*dt/dx^2 is a dimensionless diffusion coefficient % require lambda <= 0.25 for numerical stability lambda = 0.05; % END USER INPUTS % compute time step size dt (h) to get desired lambda dt = lambda*(L/n)^2/D % where dx = (L/n), length step % Matlab array indexes are integers starting at 1 % c is 2D concentration array that saves all position and time data % c(i,j), where i is position index, and j is time index % c(1,1) is at x = 0, t = 0 % c(n+1,1) is at x = L, t = 0 % specify initial conditions c(1:n+1,1) = cI; t(1) = 0; j = 1; % time index value at t = 0 jmax = 1+round(tfinal/dt); % max value of time index j tfinal = (jmax-1) * dt; % recompute tfinal using integer jmax m = n+1; % index for center node to match notes % step in time while t(j) <= tfinal % exterior node c(1,j+1) = c0; % interior nodes for i = 2:n c(i,j+1) = c(i,j) + lambda*(c(i+1,j)-2*c(i,j)+c(i-1,j)); end % node at impermeable boundary at m = n+1 c(m,j+1) = c(m,j) + lambda*(2*c(m-1,j)-2*c(m,j)); t(j+1) = t(j) + dt; % update time array j = j+1; % increment time index end fprintf('number of time steps = %i \n',jmax-1) fprintf('interior conc at final time = %6.4f \n',c(m,jmax)) x = linspace(0,L,n+1); % position in layer for plot plot(x,c(1:n+1,jmax)); hold on plot(x,c(1:n+1,round( 0.01 * jmax))); plot(x,c(1:n+1,round( 0.1 * jmax))); plot(x,c(1:n+1,round( 0.5 * jmax))); tt = sprintf('Transient 1D diffusion, conc profiles at time = %3.1f, %3.1f, %3.1f, %3.1f', ... t(round(0.01 * jmax)),t(round(0.1 * jmax)),t(round(0.5 * jmax)),t(jmax)); title(tt); ylabel('concentration (mol/m3)') xlabel('position in layer (m)') hold off figure(2) plot(t,c(m,1:jmax)) axis([0 10 0 1]) title('Transient 1D diffusion, conc at impermeable boundary vs. time') ylabel('interior conc (mol/m3)') xlabel('time (h)')