I was working through a diffusion problem and thought that Python and a package for dealing with units and unit conversions called pint would be usefull.
I'm using the Anaconda distribution of Python, which comes with the Anaconda Prompt already installed. For help installing Anaconda, see a previous blog post: Installing Anaconda on Windows 10.
To use the pint package, I needed to install pint using the Anaconda Prompt:
> pip install pintThe problem I'm working on involes the diffusion of nitrogen gas (N2) into a thin plate of iron.
Given:¶
When α-iron is put in a nitrogen atmosphere, the concentration of nitrogen in the α-iron, $C_{N}$ (in units of wt%) is a function of the nitrogen pressure $P_{N_2}$ and temperature $T$ accoding to the relationship:
$$C_{N} = 4.9 \times 10^{-6} \sqrt{P_{N_2}} exp{\frac{Q_n}{RT}} $$Where:
$Q_n = 37,600 \frac{J}{mol}$
$R=8.31 \frac{J}{mol-K}$
$T$ is the temperature in Kelvin.
At 300 °C the nitrogen gas pressure on one side of an iron plate is 0.10 MPa. On the other side of the iron plate, the nitrogen gas pressure is 5.0 MPa. The iron plate is a 1.5 mm thick. Assume the pre-exponential term $D_0$ and the activation energy of diffusion of nitrogen in carbon, $Q_d$ are equal to the values below:
$D_0 = 5 \times 10^{-7} \frac{m^2}{s}$
$Q_d = 77,000 \frac{J}{mol} $
Find:¶
Calculate the diffusion flux, J through the plate using Fick's First Law of Diffusion:
$$ J = -D \frac{dC}{dx} $$Solution:¶
We have a couple different quantities and a couple of different units to handle to solve this problem. We'll start out importing pint and creating a UnitRegistry object. We'll also need the exp (e raised to a power) and sqrt (square root) functions from the math module, part of the Python standard library.
import pint
from math import exp, sqrt
u = pint.UnitRegistry()
Let's start with the temperature, T = 300 °C.
Temperature units in °F and °C are relative units with an off-set scale. °C and °F are not multiplicative units. Non-multiplicatve units are handled by pint a little differently compared to regular multiplicative units.
To create a variable including a unit of degrees C, we instantiate a Quantity object and pass in the temperature in °C along with the unit (u.degC). We can convert the temperature to Kelvin (K) using the .ito method.
Since we want to do some mulipication, division and exponentiation with our temperature, we need to convert the temperature to a multiplicative unit. Pint has two versions of the temperature unit in Kelvin (K). There is the non-multiplicative type degK and the multiplicative type kelvin.
We convert the temperature variable T to the multiplicative type kelvin by pulling out the magnitude (the number part without the degK unit) from the T variable and multiplying it by the kevlin unit from pint.
Q_ = u.Quantity
T = Q_(300, u.degC)
print('T = {}'.format(T))
T.ito('degK')
print('T = {}'.format(T))
T = T.magnitude * u.kelvin
print(T)
Next we'll create variables for $Q_n = 37,600 \frac{J}{mol}$ and the universal gas contant $R=8.31 \frac{J}{mol-K}$
Qn = 37600 * u.J/u.mol
R = 8.31 * u.J/(u.mol*u.kelvin)
Our first nitrogen pressure is 0.10 MPa and our second nitrogen pressure is 5.0 MPa, we'll make variables for both:
PN1 = 0.10
PN2 = 5.0
Now we can calculate the two nitrogren concentrations in wt% using the equation:
$$C_{N} = 4.9 \times 10^{-6} \sqrt{P_{N_2}} exp{\frac{Q_n}{RT}}$$where $P_{N_2}$ = 0.10 for one side of the iron plate and $P_{N_2}$ = 5.0 for the other side of the iron plate
CN1 = (4.9e-3)*sqrt(PN1)*exp(-Qn/(R*T))
print(CN1)
CN2 = (4.9e-3)*sqrt(PN2)*exp(-Qn/(R*T))
print(CN2)
These values CN1 and CN2 are in units of wt% N in an iron-nitrogen "alloy" where almost all of the alloy is iron with only a small amount of nitrogen. To use Fick's First Law of Diffusion:
We need a concentration gradient $dC$ in units of mass per unit volume like kg/m3 or g/cm3 not in units of wt %. Therefore we need to convert the two concentrations of nitrogen in iron, CN1 and CN2 from units of wt% to units of kg/m3.
To make the conversion between wt% and mass per unit volume we have to pick a sample mass of iron. This mass of iron will contain a mass of nitrogen (based on wt%). We can divide this mass of nitrogen by the volume of iron that corresponds to the mass of iron we picked. As long as we divide the mass of nitrogen by the volume of iron that contains that mass of nitrogen, we will end up with a unit conversion from wt% to kg/m3 that works. So let's pick 1 kilogram of iron, and use the density of iron as 7.874 g/cm3.
We set a variable p to equal the density of iron in g/cm3 and use the .ito() method to convert the density to units of kg/m3. Then we divide the mass of iron that we picked (1 kg) and convert it to volume of iron using the density p. This will give use the volume of 1kg of iron in units of m3.
p=7.874*u.g/u.cm**3
p.ito(u.kg/u.m**3)
mFe = 1*u.kg
vFe = mFe/p
Now we'll determine how many kg of nitrogen there are in 1 kg of iron given our concentrations CN1 and CN2 in wt%. Note that we have to multiply CN1 and CN2 by 0.01 because CN1 and CN2 are in units of %.
When we divide the mass of nitrogen by the volume of iron, we get a concentration of nitrogen in iron in units of kg/m3, which is the concentration units we need to use the Fick's First Law of Diffusion.
mN1 = mFe*CN1*0.01
CN1 = mN1/vFe
print(CN1)
mN2 = mFe*CN2*0.01
CN2 = mN2/vFe
print(CN2)
Back to Fick's Fist Law of Diffusion:
$$ J = -D \frac{dC}{dx} $$The difference in concentration $dC$, is just the difference between the two concentrations CN1 and CN2 now that they are both in units of kg/m3. $dx$, the change in distance is the thickness of the plate, 1.5 mm. We'll convert the change in distance, $dx$ to units of meters using the ito() method.
dC = CN2-CN1
dx = 1.5 *u.mm
dx.ito(u.m)
Next we need to find the diffusion coefficient $D$. To do this, we need the pre-exponential term $D_0$ and the activating envery of diffusion $Q_d$.
From the beginning of the problem:
$D_0 = 5 \times 10^{-7} \frac{m^2}{s}$
$Q_d = 77,000 \frac{J}{mol} $
Let's assign these to variables with the applicable units.
D0 = 5e-7 * u.m**2/u.s
Qd = 77000 * u.J/u.mol
To calculate diffusion constant $D$, we use the equation which relates diffusion coefficient, $D$ to temperature, $T$ according to:
$$ D = D_0e^{\frac{-Q_d}{RT}} $$D = D0 * exp(-Qd/(R*T))
print(D)
Now that we have $D$, $dC$ and $dx$, we can finally calculate diffusion flux, $J$ through the plate using Fick's First Law of Diffusion:
$$ J = -D \frac{dC}{dx} $$J = -D*(dC/dx)
J
Final Answer:¶
So the final answer rounded to 3 sig figs is:
$J = -8.77 \times 10^{-15} \frac{kg}{m^{2}s}$