Chapman McDaniel and Venmathy Rajarathinam

Diffusivity: Gases, Liquids, and Solids

Matlab Code

DabCalc

function Dab = dabcalc(index,p,T)

 

%  Comprehensive program for calculating diffusivities.  Estimates the

%  phases of components involved then calculates Dab using the Wilke-Chang

%  equation for liquids and the Chapman-Enskog equation for gases.

%

%  function D = dabcalc(index,p,T)

%

%  Argument List:

%    index [=] index of 2 compounds in cnms whose Dab is to be found

%       Note:

%        - dabcalc can also be used to determine solid diffusivities

%          in very select cases.  enter -'s'- as the index to use this

%          function. 

%        - for liquids, the solvent should be entered first

%    p     [=] pressure in kPa

%    T     [=] temperature in the units of Tdeg

%   

%              

%  Returns:

%    D [=] mass diffusivity in units of cm2/s

%

%  Ex:  >> start301 (specify new session, E&M balances, ceng301 data base,

%         K as temp, 2 compounds:  CO2 and helium)

%       >> dabcalc([1 2],101.33,276)

%          The system is vapor/vapor.  If this is not correct, use Dabforce

%          Diffusivity of carbon dioxide and air          

%          ans =

%              0.1316     

%  If the phase prediction portion of the program causes errors, the

%  program Dabforce can be used to force calculation in a certain phase.

 

global TbpK

global TmpK

global mw

global cnms

global LJones

global lrho

global lvp

global Tdeg

global checklj

 

Tk = at(T);

for j=1:2;

if index == 's'

    k = menu('Choose the desired solid:','SiO2','Pyrex');

    if k == 1

            j = menu('Choose the gas involved:','He','H2');

            if j == 1

                disp('Diffusivity of Helium and SiO2')

                if Tk < 500

                    Dab = 6.3E-6*exp(-5600/1.987/Tk);

                    break

                else

                    Dab = 3.4E-6*exp(-5600/1.987/Tk);

                    break

                end

            elseif j == 2

                disp('Diffusivity of Hydrogen and SiO2')

                if Tk > 500

                    Dab = 14.9E-6*exp(-10100/1.987/Tk);

                    break

                else

                    Dab = 5.5E-6*exp(-10100/1.987/Tk);

                    break

                end

            end

    else

        disp('For Pyrex, data is only available for Helium Diffusivity')

        disp('Diffusivity of Helium and Pyrex')

            if Tk < 500

                Dab = 1.3E-6*exp(-8700/1.987/Tk);

                break

            else

                Dab = 5.5E-6*exp(-8700/1.987/Tk);

                break

            end

        end

    end

 

 

 

 

len = length(TbpK);

phases = char(2);

xs = zeros(1,len);

 

% Determine the phase of each component

% the melting point of a component is assumed

% independent of pressure

for i=1:2;

    if i == 1

        k = 2;

    else

        k = 1;

    end

    for j=1:len;

        if j == index(i);

            xs(1,j) = .9999;

        elseif j == index(k)

            xs(1,j) = .0001;

        else

            xs(1,j) = 0;

        end

    end

   

    Tb = bubpt(p,xs);

    Tb = at(Tb);

    if Tk < TmpK(index(i))

        phases(i) = 's';

    elseif Tk < Tb;

        phases(i) = 'l';

    else

        phases(i) = 'v';

    end

end

 

Dabtype = [phases(1),phases(2)];

if Dabtype ~= 'll'

    if Dabtype ~= 'vv'

        Disp('The system appears to be two-phase at these conditions')

        p = input('Should the system be liquid or vapor?(l/v)','s')

        if p == 'v'

            Dabtype = 'vv'

        elseif p == 'l'

            Dabtype = 'll'

        else

            break

        end

    end

end

 

 

% The Wilke-Chang equation for diffusivities in liquids at infinite

% dilution. 

if Dabtype == 'll'

    disp('The system is liquid/liquid.  If this is not correct, use Dabforce.')

    % assigns psi values based on compound names in cnms

    if length(cnms(index(1),:)) == 10

        if cnms(index(1),:) == 'water     '

            psi = 2.6;

        elseif cnms(index(1),:) == 'n-propanol'

            psi = 1.2;

        else

            psi = 1;

        end

    elseif length(cnms(index(1),:)) == 8

        if cnms(index(1),:) == 'methanol'

            psi = 1.9;

        else

            psi = 1;

        end

    elseif length(cnms(index(1),:)) == 7

        if cnms(index(1),:) == 'ethanol'

            psi = 1.5;

        else

            psi = 1;

        end

    else

        psi = 1.0;

    end

    % Checks to see if the density at the boiling point is available

    % if not, prompts the user to enter the molar volume

    check = isnan(lrho(index(2),2));

    if check == 1

        disp('Start301 does not have the required density data');

        Va = input('Please input the molar volume of the solute(cm^3/gmol):');

    else

        Va = 1/(lrho(index(2),2)/1000/mw(index(2)))*100^3;

    end

    % checks to see if liqmucalc will work.  if not, uses mucalc.

    mu = liqmucalc(T,index(1));

    check = isnan(mu);

    if check == 1

        exp1=['Lvp data does not exist for ',cnms(index(1),:),'.'];

        disp(exp1);

        disp('The Chapman-Enskog formula has been used to estimate the viscosity.');

        mu = mucalc(T,index(1));

    end

    % finally, calculates diffusivity

    exp1 = ['Diffusivity of ',cnms(index(1),:),' and ',cnms(index(2),:)];

    disp(exp1)   

    Dab = 7.4E-11*sqrt(psi*mw(index(1)))*T/mu/(Va^.6);

    break

end

 

% Chapman-Enskog formula for the diffusivity of light gases

if Dabtype == 'vv'

    disp('The system is vapor/vapor.  If this is not correct, use Dabforce')

    sigmaAB = .5*(LJones(index(1),1)+LJones(index(2),1));

    epsilon = sqrt(LJones(index(1),2)*LJones(index(2),2));

    kTe = Tk/epsilon;

   

    [mu,k,omega] = omegacalc(kTe);

    exp1 = ['Diffusivity of ',cnms(index(1),:),' and ',cnms(index(2),:)];

    disp(exp1)

    Dab = 1.8583E-7*101.325*sqrt(Tk^3*(1/mw(index(1))+1/mw(index(2))))/(p*sigmaAB^2*omega)*100^2;

    break

end

 

end