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

Diffusion through interface: Fick's Law, Matlab solution with error function

30 Sep. 2012

This Matlab program uses the error function to calculate the concentration on both sides of a material interface, across which the phenomenon of diffusion is taking place.

The program solves Problem 5.12 in the 6th. or 7th. edition of the textbook "Materials Science for Engineers", Shackelford, J F. The original 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)].
___________________________________________________________________

The concentration profile of material A (%) is calculated and plotted, for a block of pure metal A adjacent to a block of pure metal B, at time t. The independent variable x is measured from the original interface, between A and B. The constant of diffusion is D = 10E-14 m^2/s, and the concentration is plotted over a range of 20 micrometers on either side of the original interface.

The error function is the integral of the normal probability distribution; thus it is the cummulative probability distribution, and it depends on the parameters of the problem. Usually a parameter z is defined as the argument of the error function. The error function is available in Matlab as erf.
_________________________________________________________________________________

% error01.m

% Ricardo Ernesto Avila Montoya        Feb. 5, 2007
% Materials Science for Design.        ricardoavila.org
% Universidad Autonoma de Ciudad Juarez, Mexico

% A Matlab program to calculate and plot the error function:
% a solution of the diffusion problem.

% The plot displays the concentration of solutes at time
% t = 1800 s, on both sides of a material interface.

clc;    % clear command window

clf;    % clear figure
clear;  % clear memory

x = linspace(-20e-6, 20e-6, 40)'; % generate x-axis values

D = 1e-14;                        % m^2/s, coefficient of diffusion
t = 1800;  % seconds, time
c=50 * (1 - erf(x/(2*sqrt(D*t)))); % calculate concentration

% Plot the result in figure 1

figure(1)
plot(x,c, 'r')
axis([x(1) x(size(x,1)) 0 100]); % axis limits
grid
xlabel('distance from interface: x, m')
ylabel('% concentration')
title('Plot of error function, concentration at time 1800 s')
hold on

% Display the variables (data in memory at end of run)  

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

% end of file
_______________________________________________________________________

The result of this program is shown in this figure:

ricardoavila.org