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 pint`

The problem I'm working on involes the diffusion of nitrogen gas (N_{2}) 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:

$$ J = -D \frac{dC}{dx} $$

We need a concentration gradient $dC$ in units of mass per unit volume like kg/m^{3} or g/cm^{3} 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/m^{3}.

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/m^{3} that works. So let's pick 1 kilogram of iron, and use the density of iron as 7.874 g/cm^{3}.

We set a variable `p`

to equal the density of iron in g/cm^{3} and use the `.ito()`

method to convert the density to units of kg/m^{3}. 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 m^{3}.

```
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/m^{3}, 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/m^{3}. $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}$