Material diffusion in time: Fick's law, Matlab program solution

30 Sep. 2012

This Matlab program simulates the diffusion process through a material interface.
The interface separates two pure metals, and the changes of the concentration are calculated and plotted (in %) during a period of time.

The program solves Problem 5.12 in the 6th. or 7th. edition of the textbook
"Materials Science for Engineers", Shackelford, J F. The original stated problem is:
__________________________________________________________________
A diffusion couple is formed when two different materials are allowed ot interdiffuse at an
elevated temperature. For a block of pure metal A adjacent to a block of pure metal B, the concentration profile of A (in at%) after interdiffusion is given by

cx = 50 * (1 - erf(x / (2 * sqrt(D * t)))

where x is measured from the original interface. For a diffusion couple with D = 1E-14 m^2/s, plot the concentration profile of metal A over a range of 20E-6 m, on either side of the original interface (x=0), after a time of 1 hour. [Note that erf(-z) = -erf(z)].
__________________________________________________________________

This simulation results as the solution to the diffusion problem, from Fick's laws of diffusion.

The frames of the "movie" show how concentration changes with time, on both sides of
the material interface.

__________________________________________________________________

% error02.m

% Ricardo E. Avila        Feb. 5, 2007

% Materials for Design.   ricardoavila.org  

% Universidad Autonoma de Ciudad Juarez, Mexico

% This file plots the change of concentration

% with time, on both sides of an interface.

% x = 0 is the position of the interface.

clc;   % clear command window

clf;   % clear figure window

clear; % clear memory

x = linspace(-20e-6, 20e-6, 40)'; % x-axis column array

D = 1e-14;    % m^2/s,   material diffusion coefficient

t_end = 1800; % seconds, end simulation time

% calculate concentration at time t

x1 = x(1);              % left  end of plot

x_end = x(size(x, 1));  % right end of plot

axis([x1 x_end 0 100]); % axis limits

for t = 1 : t_end; 

   c = single(50 * (1 - erf(x/(2*sqrt(D*t))))); % error function

   plot(x, c, 'r')

   axis([x1 x_end 0 100]); % axis limits

   xlabel('x, m')

   ylabel('% concentration')

   title('Plot error function, concentration change with time, s')

   grid

   % create movie frame

   M(t) = getframe;

end

% display variables in memory at end of run: 

whos

disp(' ***** Succesful execution of program ***** ');

disp(' ****************************************** ');

% end of file
_________________________________________________________

One of the "frames" obtained during the simulation is shown here:

ricardoavila.org