Diffusion in Solids

Fick's laws, diffusion mechanisms, and applications in materials processing

Fick's First Law

Fick's first law relates the diffusive flux to the concentration gradient under steady-state conditions (concentration profile does not change with time):

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

where $J$ is the flux (atoms/m²/s or mol/m²/s),$D$ is the diffusion coefficient (m²/s), and $dC/dx$ is the concentration gradient. The negative sign indicates that diffusion occurs down the concentration gradient (from high to low concentration).

Fick's Second Law

For non-steady-state diffusion where the concentration changes with time:

$$\frac{\partial C}{\partial t} = D\frac{\partial^2 C}{\partial x^2}$$

This is a parabolic partial differential equation (the heat equation). For a semi-infinite solid with constant surface concentration $C_s$ and initial uniform concentration $C_0$, the solution involves the complementary error function:

$$C(x,t) = C_0 + (C_s - C_0)\,\text{erfc}\!\left(\frac{x}{2\sqrt{Dt}}\right)$$

The characteristic diffusion length is $x_d = 2\sqrt{Dt}$, which sets the penetration depth of the diffusion front.

Temperature Dependence of Diffusion

Like reaction rates, the diffusion coefficient follows an Arrhenius-type temperature dependence:

$$D = D_0 \exp\!\left(-\frac{Q}{RT}\right)$$

where $D_0$ is the pre-exponential factor (frequency factor),$Q$ is the activation energy for diffusion,$R = 8.314$ J/(mol·K), and $T$ is absolute temperature. Plotting $\ln D$ vs $1/T$ gives a straight line with slope $-Q/R$.

Diffusion Mechanisms

Vacancy Diffusion

An atom jumps into an adjacent vacant lattice site. Requires both a neighboring vacancy and sufficient thermal energy to overcome the migration barrier. Dominant in substitutional alloys. $Q$ is typically 1–3 eV.

Interstitial Diffusion

Small atoms (C, N, H, O) move through the interstices between larger host atoms. No vacancy needed, so it is much faster than vacancy diffusion. Typical $Q \sim 0.5\text{--}1$ eV. Example: carbon in iron.

Kirkendall Effect

In a binary diffusion couple (e.g., Cu-Zn), the two species diffuse at different rates. The faster-diffusing species creates a net vacancy flux in the opposite direction, causing the interface (marked by inert markers) to shift. This was the first direct proof that diffusion in metals occurs by the vacancy mechanism, not by direct atomic interchange.

Applications

Carburization

Diffusing carbon into the surface of steel at 900–950°C to increase surface hardness while maintaining a tough core. Case depth controlled by $2\sqrt{Dt}$.

Semiconductor Doping

Thermal diffusion of dopants (B, P, As) into silicon at controlled temperatures to create p-n junctions. Being replaced by ion implantation for precise depth control.

Oxidation

Growth of oxide layers on metals. Parabolic growth law $x^2 = 2Dt$ when oxygen diffusion through the growing oxide is rate-limiting (Wagner theory).

Python: Fick's Second Law (Finite Difference)

Solves the diffusion equation numerically using an explicit finite difference method (FTCS scheme) for carburization of steel. Compares the numerical solution with the analytic error function result.

Fick's Second Law: Carburization

Python

Finite difference solution of the diffusion equation with error function comparison for steel carburization

script.py101 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.special import erfc

# ──────────────────────────────────────────────────────
# Fick's Second Law: Finite Difference Solution
# Carburization of steel: C(x,t) with fixed surface conc.
# Also compare with the analytic error function solution
# ──────────────────────────────────────────────────────

# Parameters
D = 2.0e-11    # Diffusion coefficient [m^2/s] (C in Fe at ~950 C)
L = 5.0e-3     # Domain length [m]
C_s = 1.0      # Surface concentration (wt% C)
C_0 = 0.2      # Initial bulk concentration (wt% C)

# Spatial grid
Nx = 200
dx = L / (Nx - 1)
x = np.linspace(0, L, Nx)

# Time parameters
dt = 0.4 * dx**2 / D   # stability: dt < dx^2/(2D)
t_max = 3600.0 * 4      # 4 hours
Nt = int(t_max / dt) + 1

# Snapshots to record
snap_times = [0.5, 1.0, 2.0, 4.0]  # hours
snap_steps = [int(t_h * 3600 / dt) for t_h in snap_times]
snap_profiles = {}

# Initialize concentration
C = np.full(Nx, C_0)
C[0] = C_s  # fixed surface BC

# Finite difference time-stepping (explicit FTCS)
for n in range(Nt):
    t_hours = n * dt / 3600.0
    if n in snap_steps:
        snap_profiles[t_hours] = C.copy()

C_new = C.copy()
    for i in range(1, Nx - 1):
        C_new[i] = C[i] + D * dt / dx**2 * (C[i+1] - 2*C[i] + C[i-1])
    # BCs: C[0] = C_s (Dirichlet), C[-1] = C_0 (far field)
    C_new[0] = C_s
    C_new[-1] = C_0
    C = C_new

# Capture final if not already captured
t_final = Nt * dt / 3600.0

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

# Left: Numerical profiles at different times
colors_list = ['cyan', 'lime', 'yellow', 'red']
for idx, (t_h, prof) in enumerate(sorted(snap_profiles.items())):
    ax1.plot(x * 1000, prof, color=colors_list[idx], linewidth=2,
             label=f't = {t_h:.1f} hr (FD)')

# Analytic error function solution for comparison
for idx, t_h in enumerate(snap_times):
    t_sec = t_h * 3600
    C_analytic = C_s - (C_s - C_0) * erfc(x / (2 * np.sqrt(D * t_sec)))
    # erfc solution: C = C_0 + (C_s - C_0) * erfc(x/(2*sqrt(Dt)))
    C_analytic2 = C_0 + (C_s - C_0) * erfc(x / (2 * np.sqrt(D * t_sec)))
    ax1.plot(x * 1000, C_analytic2, '--', color=colors_list[idx],
             linewidth=1.5, alpha=0.7)

ax1.set_xlabel('Depth x (mm)', fontsize=12)
ax1.set_ylabel('Carbon concentration (wt%)', fontsize=12)
ax1.set_title("Fick's 2nd Law: Carburization", fontsize=14)
ax1.legend(fontsize=9)
ax1.set_xlim(0, 2.5)
ax1.grid(True, alpha=0.3)
ax1.text(1.5, 0.85, 'Solid: FD numerical\nDashed: erf analytic',
         fontsize=9, color='gray')

# Right: Penetration depth vs time
times_hr = np.linspace(0.1, 4.0, 100)
times_sec = times_hr * 3600
pen_depth = 2 * np.sqrt(D * times_sec) * 1000  # in mm (approx x where C ~ C_0 + 0.16*(Cs-C0))

ax2.plot(times_hr, pen_depth, 'cyan', linewidth=2.5)
ax2.fill_between(times_hr, 0, pen_depth, alpha=0.15, color='cyan')
ax2.set_xlabel('Time (hours)', fontsize=12)
ax2.set_ylabel('Penetration depth $2\\sqrt{Dt}$ (mm)', fontsize=12)
ax2.set_title('Diffusion Penetration Depth', fontsize=14)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('output.png', dpi=150, bbox_inches='tight')

print(f"Diffusion coefficient D = {D:.2e} m^2/s")
print(f"Grid: Nx={Nx}, dx={dx*1e6:.1f} um, dt={dt:.3f} s")
print(f"Stability parameter r = D*dt/dx^2 = {D*dt/dx**2:.4f} (< 0.5)")
print(f"Penetration depth at 4 hr: 2*sqrt(Dt) = {2*np.sqrt(D*4*3600)*1000:.3f} mm")
print("Plot saved.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran: Diffusion Coefficient & Penetration Depth

Computes the diffusion coefficient at various temperatures and the penetration depth for carburization of steel, using the Arrhenius relation.

Diffusion Coefficient & Penetration Depth

Fortran

Computes D(T) for carbon in iron and carburization penetration depth vs time

diffusion_calc.f9061 lines

program diffusion_calc
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  real(dp), parameter :: R = 8.314_dp   ! J/(mol*K)

! Carbon in alpha-Fe (BCC) and gamma-Fe (FCC)
  real(dp) :: D0_alpha, Q_alpha   ! BCC ferrite
  real(dp) :: D0_gamma, Q_gamma   ! FCC austenite
  real(dp) :: T, D_val, pen_depth
  real(dp) :: t_hours, t_sec
  integer :: i, j

! Diffusion parameters (C in Fe)
  D0_alpha = 6.2e-7_dp   ! m^2/s
  Q_alpha  = 80000.0_dp   ! J/mol (80 kJ/mol)
  D0_gamma = 2.3e-5_dp   ! m^2/s
  Q_gamma  = 148000.0_dp  ! J/mol (148 kJ/mol)

write(*,'(A)') '======================================================'
  write(*,'(A)') '  Diffusion of Carbon in Iron: D(T) and Penetration'
  write(*,'(A)') '======================================================'
  write(*,'(A)') ''

! Part 1: D vs Temperature
  write(*,'(A)') '--- Diffusion Coefficient vs Temperature ---'
  write(*,'(A)') ''
  write(*,'(A8, A15, A15)') 'T (C)', 'D_alpha', 'D_gamma'
  write(*,'(A8, A15, A15)') '', '(m^2/s)', '(m^2/s)'
  write(*,'(A)') '--------------------------------------'

do i = 300, 1100, 100
    T = real(i, dp) + 273.15_dp
    write(*,'(I8, ES15.3, ES15.3)') i, &
      D0_alpha * exp(-Q_alpha / (R * T)), &
      D0_gamma * exp(-Q_gamma / (R * T))
  end do

write(*,'(A)') ''

! Part 2: Carburization penetration depth
  write(*,'(A)') '--- Carburization Penetration Depth (austenite, 950 C) ---'
  write(*,'(A)') ''
  T = 950.0_dp + 273.15_dp
  D_val = D0_gamma * exp(-Q_gamma / (R * T))
  write(*,'(A,ES12.4,A)') 'D at 950 C = ', D_val, ' m^2/s'
  write(*,'(A)') ''
  write(*,'(A12, A18)') 'Time (hr)', 'Depth 2*sqrt(Dt)'
  write(*,'(A12, A18)') '', '(mm)'
  write(*,'(A)') '------------------------------'

do j = 1, 10
    t_hours = real(j, dp)
    t_sec = t_hours * 3600.0_dp
    pen_depth = 2.0_dp * sqrt(D_val * t_sec) * 1000.0_dp  ! convert to mm
    write(*,'(F12.1, F18.4)') t_hours, pen_depth
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Note: Penetration depth ~ sqrt(t), so doubling time'
  write(*,'(A)') '      increases depth by factor of sqrt(2) = 1.414'
end program diffusion_calc

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Video Lectures

Lecture 35: Diffusion I

Lecture 36: Diffusion II

← Part 8 Overview Course Overview →