Structure, Functionality, and Bioceramics - Design Principles for High Performance
Upon completing this chapter, you will be able to explain:
Structural ceramics are ceramic materials with excellent mechanical properties (high strength, high hardness, heat resistance) that are used as structural components in harsh environments. They can be used in high-temperature or corrosive environments where metallic materials cannot, with important applications including:
Structural ceramics are indispensable in the aerospace, automotive, and medical fields. Advanced ceramic materials account for roughly 60% of the global ceramics market (over $230B as of 2023). The reasons include:
High-strength ceramics are represented by the following three major materials:
flowchart LR
A[Al₂O₃
Alumina] --> B[High Hardness
Hv 2000]
C[ZrO₂
Zirconia] --> D[High Toughness
10-15 MPa√m]
E[Si₃N₄
Silicon Nitride] --> F[High-Temperature Strength
Usable to 1400°C]
style A fill:#e3f2fd
style C fill:#fff3e0
style E fill:#e8f5e9
style B fill:#f3e5f5
style D fill:#fce4ec
style F fill:#fff9c4
While ceramics have high strength and high hardness, their greatest drawback is brittleness (low toughness). Microscopic defects (pores, cracks) act as stress concentration points and trigger sudden fracture (Griffith theory). Fracture toughness is less than 1/10 that of metals, which is why toughening technology remains an important research topic.
This is the most effective toughening mechanism for zirconia (ZrO₂):
Toughening mechanism:
Implementation: Y₂O₃ (3-8 mol%) or MgO (9-15 mol%) is added to stabilize the tetragonal phase at room temperature (PSZ: Partially Stabilized Zirconia)
This technique combines high-strength fibers with a ceramic matrix:
Toughening mechanism:
Applications: SiC/SiC composites (aircraft engine components), C/C composites (brake disks)
The piezoelectric effect is a phenomenon in which electrical polarization is generated when mechanical stress is applied (the direct piezoelectric effect), and conversely, mechanical strain is generated when an electric field is applied (the converse piezoelectric effect).
Characteristics of PZT (Lead Zirconate Titanate):
Because PZT contains more than 60 wt% lead (Pb), its use is restricted under the European RoHS directive. Lead-free alternatives such as BaTiO₃-based, (K,Na)NbO₃-based, and BiFeO₃-based materials are being researched, but none yet match PZT's performance (d₃₃ = 100-300 pC/N). Piezoelectric devices are currently exempt for applications such as medical equipment, but development of alternative materials remains necessary in the long term.
The piezoelectric effect appears only in materials with a non-centrosymmetric crystal structure:
Dielectric ceramics are used as capacitor materials with a high dielectric constant (εᵣ) that store electrical energy.
Origin of the high dielectric constant:
Modern MLCCs have been miniaturized and improved to an extreme degree:
BaTiO₃-based MLCCs are a key material enabling the miniaturization and performance improvement of electronic devices.
Ferrites are oxide-based magnetic materials with low loss at high frequency, and are widely used in transformers, inductors, and electromagnetic wave absorbers.
Characteristics of spinel ferrites:
Characteristics of hexagonal ferrites (hard ferrites):
Ferrite magnetism arises from the antiparallel alignment of the magnetic moments of the ions at the A-sites (tetrahedral positions) and B-sites (octahedral positions) of the spinel structure (AB₂O₄) — that is, ferrimagnetism. In Mn-Zn ferrite, the magnetic moments of Mn²⁺ and Fe³⁺ partially cancel, so the net magnetization is small, but this is precisely what makes the high permeability possible.
Bioceramics are ceramic materials that do not trigger a rejection response when in contact with living tissue (biocompatibility) and that can bond directly with bone tissue (osteoconductivity).
Characteristics of hydroxyapatite (HAp):
Unlike HAp, β-TCP (tricalcium phosphate) is gradually resorbed within the body:
Because of this bioresorbability, no permanent foreign material remains in the body, enabling ideal bone regeneration in which the implant is fully replaced by the patient's own bone tissue.
# ===================================
# Example 1: Arrhenius Equation Simulation
# ===================================
import numpy as np
import matplotlib.pyplot as plt
# Physical constant
R = 8.314 # J/(mol·K)
# Diffusion parameters for the BaTiO3 system (literature values)
D0 = 5e-4 # m²/s (pre-exponential factor)
Ea = 300e3 # J/mol (activation energy, 300 kJ/mol)
def diffusion_coefficient(T, D0, Ea):
"""Calculate the diffusion coefficient using the Arrhenius equation
Args:
T (float or array): Temperature [K]
D0 (float): Pre-exponential factor [m²/s]
Ea (float): Activation energy [J/mol]
Returns:
float or array: Diffusion coefficient [m²/s]
"""
return D0 * np.exp(-Ea / (R * T))
# Temperature range 800-1400°C
T_celsius = np.linspace(800, 1400, 100)
T_kelvin = T_celsius + 273.15
# Calculate the diffusion coefficient
D = diffusion_coefficient(T_kelvin, D0, Ea)
# Plot
plt.figure(figsize=(10, 6))
# Logarithmic plot (Arrhenius plot)
plt.subplot(1, 2, 1)
plt.semilogy(T_celsius, D, 'b-', linewidth=2)
plt.xlabel('Temperature (°C)', fontsize=12)
plt.ylabel('Diffusion Coefficient (m²/s)', fontsize=12)
plt.title('Arrhenius Plot', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
# 1/T vs ln(D) plot (linear relationship)
plt.subplot(1, 2, 2)
plt.plot(1000/T_kelvin, np.log(D), 'r-', linewidth=2)
plt.xlabel('1000/T (K⁻¹)', fontsize=12)
plt.ylabel('ln(D)', fontsize=12)
plt.title('Linearized Arrhenius Plot', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('arrhenius_plot.png', dpi=300, bbox_inches='tight')
plt.show()
# Display the diffusion coefficient at key temperatures
key_temps = [1000, 1100, 1200, 1300]
print("Comparison of temperature dependence:")
print("-" * 50)
for T_c in key_temps:
T_k = T_c + 273.15
D_val = diffusion_coefficient(T_k, D0, Ea)
print(f"{T_c:4d}°C: D = {D_val:.2e} m²/s")
# Example output:
# Comparison of temperature dependence:
# --------------------------------------------------
# 1000°C: D = 1.89e-12 m²/s
# 1100°C: D = 9.45e-12 m²/s
# 1200°C: D = 4.01e-11 m²/s
# 1300°C: D = 1.48e-10 m²/s
# ===================================
# Example 2: Conversion Calculation with the Jander Equation
# ===================================
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
def jander_equation(alpha, k, t):
"""Jander equation
Args:
alpha (float): Conversion (0-1)
k (float): Rate constant [s⁻¹]
t (float): Time [s]
Returns:
float: Left-hand side of the Jander equation minus k*t
"""
return (1 - (1 - alpha)**(1/3))**2 - k * t
def calculate_conversion(k, t):
"""Calculate the conversion at time t
Args:
k (float): Rate constant
t (float): Time
Returns:
float: Conversion (0-1)
"""
# Solve the Jander equation for alpha numerically
alpha0 = 0.5 # Initial guess
alpha = fsolve(lambda a: jander_equation(a, k, t), alpha0)[0]
return np.clip(alpha, 0, 1) # Restrict to the 0-1 range
# Parameter settings
D = 1e-11 # m²/s (diffusion coefficient at 1200°C)
C0 = 10000 # mol/m³
r0_values = [1e-6, 5e-6, 10e-6] # Particle radius [m]: 1μm, 5μm, 10μm
# Time array (0-50 hours)
t_hours = np.linspace(0, 50, 500)
t_seconds = t_hours * 3600
# Plot
plt.figure(figsize=(12, 5))
# Effect of particle size
plt.subplot(1, 2, 1)
for r0 in r0_values:
k = D * C0 / r0**2
alpha = [calculate_conversion(k, t) for t in t_seconds]
plt.plot(t_hours, alpha, linewidth=2,
label=f'r₀ = {r0*1e6:.1f} μm')
plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Conversion (α)', fontsize=12)
plt.title('Effect of Particle Size', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 1])
# Effect of temperature (fixed particle size)
plt.subplot(1, 2, 2)
r0_fixed = 5e-6 # Fixed at 5μm
temperatures = [1100, 1200, 1300] # °C
for T_c in temperatures:
T_k = T_c + 273.15
D_T = diffusion_coefficient(T_k, D0, Ea)
k = D_T * C0 / r0_fixed**2
alpha = [calculate_conversion(k, t) for t in t_seconds]
plt.plot(t_hours, alpha, linewidth=2,
label=f'{T_c}°C')
plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Conversion (α)', fontsize=12)
plt.title('Effect of Temperature (r₀ = 5 μm)', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 1])
plt.tight_layout()
plt.savefig('jander_simulation.png', dpi=300, bbox_inches='tight')
plt.show()
# Calculate the time required to reach 50% conversion
print("\nTime required to reach 50% conversion:")
print("-" * 50)
for r0 in r0_values:
k = D * C0 / r0**2
t_50 = fsolve(lambda t: jander_equation(0.5, k, t), 10000)[0]
print(f"r₀ = {r0*1e6:.1f} μm: t₅₀ = {t_50/3600:.1f} hours")
# Example output:
# Time required to reach 50% conversion:
# --------------------------------------------------
# r₀ = 1.0 μm: t₅₀ = 1.9 hours
# r₀ = 5.0 μm: t₅₀ = 47.3 hours
# r₀ = 10.0 μm: t₅₀ = 189.2 hours
# ===================================
# Example 3: Calculating Activation Energy with the Kissinger Method
# ===================================
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import linregress
# Kissinger method: Ea is obtained from the slope of ln(β/Tp²) vs 1/Tp
# β: heating rate [K/min]
# Tp: peak temperature [K]
# slope = -Ea/R
# Experimental data (DSC peak temperatures at different heating rates)
heating_rates = np.array([5, 10, 15, 20]) # K/min
peak_temps_celsius = np.array([1085, 1105, 1120, 1132]) # °C
peak_temps_kelvin = peak_temps_celsius + 273.15
def kissinger_analysis(beta, Tp):
"""Calculate the activation energy using the Kissinger method
Args:
beta (array): Heating rate [K/min]
Tp (array): Peak temperature [K]
Returns:
tuple: (Ea [kJ/mol], A [min⁻¹], R²)
"""
# Left-hand side of the Kissinger equation
y = np.log(beta / Tp**2)
# 1/Tp
x = 1000 / Tp # Scaled by 1000/T for readability
# Linear regression
slope, intercept, r_value, p_value, std_err = linregress(x, y)
# Calculate the activation energy
R = 8.314 # J/(mol·K)
Ea = -slope * R * 1000 # J/mol → kJ/mol
# Pre-exponential factor
A = np.exp(intercept)
return Ea, A, r_value**2
# Calculate the activation energy
Ea, A, R2 = kissinger_analysis(heating_rates, peak_temps_kelvin)
print("Results of the Kissinger analysis:")
print("=" * 50)
print(f"Activation energy Ea = {Ea:.1f} kJ/mol")
print(f"Pre-exponential factor A = {A:.2e} min⁻¹")
print(f"Coefficient of determination R² = {R2:.4f}")
print("=" * 50)
# Plot
plt.figure(figsize=(10, 6))
# Kissinger plot
y_data = np.log(heating_rates / peak_temps_kelvin**2)
x_data = 1000 / peak_temps_kelvin
plt.plot(x_data, y_data, 'ro', markersize=10, label='Experimental data')
# Fitted line
x_fit = np.linspace(x_data.min()*0.95, x_data.max()*1.05, 100)
slope = -Ea * 1000 / (R * 1000)
intercept = np.log(A)
y_fit = slope * x_fit + intercept
plt.plot(x_fit, y_fit, 'b-', linewidth=2, label=f'Fit: Ea = {Ea:.1f} kJ/mol')
plt.xlabel('1000/Tp (K⁻¹)', fontsize=12)
plt.ylabel('ln(β/Tp²)', fontsize=12)
plt.title('Kissinger Plot for Activation Energy', fontsize=14, fontweight='bold')
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
# Show the results in a text box
textstr = f'Ea = {Ea:.1f} kJ/mol\nR² = {R2:.4f}'
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
plt.text(0.05, 0.95, textstr, transform=plt.gca().transAxes, fontsize=11,
verticalalignment='top', bbox=props)
plt.tight_layout()
plt.savefig('kissinger_plot.png', dpi=300, bbox_inches='tight')
plt.show()
# Example output:
# Results of the Kissinger analysis:
# ==================================================
# Activation energy Ea = 287.3 kJ/mol
# Pre-exponential factor A = 2.45e+12 min⁻¹
# Coefficient of determination R² = 0.9956
# ==================================================
The temperature profile used in a solid-state reaction is the single most important control parameter determining whether the reaction succeeds. The following three elements must be designed carefully:
flowchart TD
A[Temperature Profile Design] --> B[Heating Rate]
A --> C[Holding Time]
A --> D[Cooling Rate]
B --> B1[Too fast: thermal stress → cracking]
B --> B2[Too slow: unwanted phase transformations]
C --> C1[Too short: incomplete reaction]
C --> C2[Too long: excessive grain growth]
D --> D1[Too fast: thermal stress → cracking]
D --> D2[Too slow: unfavorable phases]
style A fill:#f093fb
style B fill:#e3f2fd
style C fill:#e8f5e9
style D fill:#fff3e0
Typical recommended value: 2-10°C/min
Factors to consider:
In BaTiO₃ synthesis, BaCO₃ → BaO + CO₂ decomposition occurs at 800-900°C. At heating rates above 20°C/min, CO₂ can be released so rapidly that the sample ruptures. A heating rate of 5°C/min or below is recommended.
How to determine it: Estimate from the Jander equation, then optimize experimentally
The required holding time can be estimated with the following equation:
Typical holding times:
Typical recommended value: 1-5°C/min (slower than the heating rate)
Why it matters:
# ===================================
# Example 4: Temperature Profile Optimization
# ===================================
import numpy as np
import matplotlib.pyplot as plt
def temperature_profile(t, T_target, heating_rate, hold_time, cooling_rate):
"""Generate a temperature profile
Args:
t (array): Time array [min]
T_target (float): Holding temperature [°C]
heating_rate (float): Heating rate [°C/min]
hold_time (float): Holding time [min]
cooling_rate (float): Cooling rate [°C/min]
Returns:
array: Temperature profile [°C]
"""
T_room = 25 # Room temperature
T = np.zeros_like(t)
# Heating duration
t_heat = (T_target - T_room) / heating_rate
# Time cooling begins
t_cool_start = t_heat + hold_time
for i, time in enumerate(t):
if time <= t_heat:
# Heating phase
T[i] = T_room + heating_rate * time
elif time <= t_cool_start:
# Holding phase
T[i] = T_target
else:
# Cooling phase
T[i] = T_target - cooling_rate * (time - t_cool_start)
T[i] = max(T[i], T_room) # Never goes below room temperature
return T
def simulate_reaction_progress(T, t, Ea, D0, r0):
"""Calculate reaction progress based on a temperature profile
Args:
T (array): Temperature profile [°C]
t (array): Time array [min]
Ea (float): Activation energy [J/mol]
D0 (float): Pre-exponential factor [m²/s]
r0 (float): Particle radius [m]
Returns:
array: Conversion
"""
R = 8.314
C0 = 10000
alpha = np.zeros_like(t)
for i in range(1, len(t)):
T_k = T[i] + 273.15
D = D0 * np.exp(-Ea / (R * T_k))
k = D * C0 / r0**2
dt = (t[i] - t[i-1]) * 60 # min → s
# Simple integration (reaction progress over a small time step)
if alpha[i-1] < 0.99:
dalpha = k * dt / (2 * (1 - (1-alpha[i-1])**(1/3)))
alpha[i] = min(alpha[i-1] + dalpha, 1.0)
else:
alpha[i] = alpha[i-1]
return alpha
# Parameter settings
T_target = 1200 # °C
hold_time = 240 # min (4 hours)
Ea = 300e3 # J/mol
D0 = 5e-4 # m²/s
r0 = 5e-6 # m
# Comparison across heating rates
heating_rates = [2, 5, 10, 20] # °C/min
cooling_rate = 3 # °C/min
# Time array
t_max = 800 # min
t = np.linspace(0, t_max, 2000)
# Plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10))
# Temperature profiles
for hr in heating_rates:
T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
ax1.plot(t/60, T_profile, linewidth=2, label=f'{hr}°C/min')
ax1.set_xlabel('Time (hours)', fontsize=12)
ax1.set_ylabel('Temperature (°C)', fontsize=12)
ax1.set_title('Temperature Profiles', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, t_max/60])
# Reaction progress
for hr in heating_rates:
T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
alpha = simulate_reaction_progress(T_profile, t, Ea, D0, r0)
ax2.plot(t/60, alpha, linewidth=2, label=f'{hr}°C/min')
ax2.axhline(y=0.95, color='red', linestyle='--', linewidth=1, label='Target (95%)')
ax2.set_xlabel('Time (hours)', fontsize=12)
ax2.set_ylabel('Conversion', fontsize=12)
ax2.set_title('Reaction Progress', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim([0, t_max/60])
ax2.set_ylim([0, 1])
plt.tight_layout()
plt.savefig('temperature_profile_optimization.png', dpi=300, bbox_inches='tight')
plt.show()
# Calculate the time to reach 95% conversion for each heating rate
print("\nComparison of time to reach 95% conversion:")
print("=" * 60)
for hr in heating_rates:
T_profile = temperature_profile(t, T_target, hr, hold_time, cooling_rate)
alpha = simulate_reaction_progress(T_profile, t, Ea, D0, r0)
# Time at which 95% is reached
idx_95 = np.where(alpha >= 0.95)[0]
if len(idx_95) > 0:
t_95 = t[idx_95[0]] / 60
print(f"Heating rate {hr:2d}°C/min: t₉₅ = {t_95:.1f} hours")
else:
print(f"Heating rate {hr:2d}°C/min: reaction incomplete")
# Example output:
# Comparison of time to reach 95% conversion:
# ============================================================
# Heating rate 2°C/min: t₉₅ = 7.8 hours
# Heating rate 5°C/min: t₉₅ = 7.2 hours
# Heating rate 10°C/min: t₉₅ = 6.9 hours
# Heating rate 20°C/min: t₉₅ = 6.7 hours
pycalphad is a Python library for phase diagram calculations based on the CALPHAD (CALculation of PHAse Diagrams) method. It calculates equilibrium phases from a thermodynamic database and is useful for designing reaction pathways.
# ===================================
# Example 5: Phase Diagram Calculation with pycalphad
# ===================================
# Note: requires pycalphad to be installed
# pip install pycalphad
from pycalphad import Database, equilibrium, variables as v
import matplotlib.pyplot as plt
import numpy as np
# Load a TDB database (a simplified example is used here)
# In practice an appropriate TDB file is required
# Example: the BaO-TiO2 system
# Simplified TDB string (the real one is far more complex)
tdb_string = """
$ BaO-TiO2 system (simplified)
ELEMENT BA BCC_A2 137.327 !
ELEMENT TI HCP_A3 47.867 !
ELEMENT O GAS 15.999 !
FUNCTION GBCCBA 298.15 +GHSERBA; 6000 N !
FUNCTION GHCPTI 298.15 +GHSERTI; 6000 N !
FUNCTION GGASO 298.15 +GHSERO; 6000 N !
PHASE LIQUID:L % 1 1.0 !
PHASE BAO_CUBIC % 2 1 1 !
PHASE TIO2_RUTILE % 2 1 2 !
PHASE BATIO3 % 3 1 1 3 !
"""
# Note: an actual calculation requires a formal TDB file
# Here we stop at a conceptual explanation
print("Concept of phase diagram calculation with pycalphad:")
print("=" * 60)
print("1. Load a TDB database (thermodynamic data)")
print("2. Set the temperature and composition ranges")
print("3. Run the equilibrium calculation")
print("4. Visualize the stable phases")
print()
print("Real-world applications:")
print("- BaO-TiO2 system: formation temperature and composition range of BaTiO3")
print("- Si-N system: stability region of Si3N4")
print("- Phase relationships in multicomponent ceramics")
# Conceptual plot (an illustration based on real data)
fig, ax = plt.subplots(figsize=(10, 7))
# Temperature range
T = np.linspace(800, 1600, 100)
# Stability regions of each phase (conceptual diagram)
# BaO + TiO2 → BaTiO3 reaction
BaO_region = np.ones_like(T) * 0.3
TiO2_region = np.ones_like(T) * 0.7
BaTiO3_region = np.where((T > 1100) & (T < 1400), 0.5, np.nan)
ax.fill_between(T, 0, BaO_region, alpha=0.3, color='blue', label='BaO + TiO2')
ax.fill_between(T, BaO_region, TiO2_region, alpha=0.3, color='green',
label='BaTiO3 stable')
ax.fill_between(T, TiO2_region, 1, alpha=0.3, color='red', label='Liquid')
ax.axhline(y=0.5, color='black', linestyle='--', linewidth=2,
label='BaTiO3 composition')
ax.axvline(x=1100, color='gray', linestyle=':', linewidth=1, alpha=0.5)
ax.axvline(x=1400, color='gray', linestyle=':', linewidth=1, alpha=0.5)
ax.set_xlabel('Temperature (°C)', fontsize=12)
ax.set_ylabel('Composition (BaO mole fraction)', fontsize=12)
ax.set_title('Conceptual Phase Diagram: BaO-TiO2', fontsize=14, fontweight='bold')
ax.legend(fontsize=10, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_xlim([800, 1600])
ax.set_ylim([0, 1])
# Text annotation
ax.text(1250, 0.5, 'BaTiO₃\nformation\nregion',
fontsize=11, ha='center', va='center',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.7))
plt.tight_layout()
plt.savefig('phase_diagram_concept.png', dpi=300, bbox_inches='tight')
plt.show()
# Example of real usage (commented out)
"""
# Example of real pycalphad usage
db = Database('BaO-TiO2.tdb') # Load the TDB file
# Equilibrium calculation
eq = equilibrium(db, ['BA', 'TI', 'O'], ['LIQUID', 'BATIO3'],
{v.X('BA'): (0, 1, 0.01),
v.T: (1000, 1600, 50),
v.P: 101325})
# Plot the results
eq.plot()
"""
Design of Experiments (DOE) is a statistical method for finding optimal conditions with the minimum number of experiments in a system where multiple parameters interact.
Key parameters to optimize in a solid-state reaction:
# ===================================
# Example 6: Optimizing Conditions with DOE
# ===================================
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize
# A hypothetical conversion model (a function of temperature and time)
def reaction_yield(T, t, noise=0):
"""Calculate conversion from temperature and time (hypothetical model)
Args:
T (float): Temperature [°C]
t (float): Time [hours]
noise (float): Noise level
Returns:
float: Conversion [%]
"""
# Optimum: T=1200°C, t=6 hours
T_opt = 1200
t_opt = 6
# Quadratic (Gaussian) model
yield_val = 100 * np.exp(-((T-T_opt)/150)**2 - ((t-t_opt)/3)**2)
# Add noise
if noise > 0:
yield_val += np.random.normal(0, noise)
return np.clip(yield_val, 0, 100)
# Experimental design points (central composite design)
T_levels = [1000, 1100, 1200, 1300, 1400] # °C
t_levels = [2, 4, 6, 8, 10] # hours
# Arrange experimental points on a grid
T_grid, t_grid = np.meshgrid(T_levels, t_levels)
yield_grid = np.zeros_like(T_grid, dtype=float)
# Measure the conversion at each experimental point (simulated)
for i in range(len(t_levels)):
for j in range(len(T_levels)):
yield_grid[i, j] = reaction_yield(T_grid[i, j], t_grid[i, j], noise=2)
# Display the results
print("Optimizing reaction conditions with Design of Experiments")
print("=" * 70)
print(f"{'Temperature (°C)':<20} {'Time (hours)':<20} {'Yield (%)':<20}")
print("-" * 70)
for i in range(len(t_levels)):
for j in range(len(T_levels)):
print(f"{T_grid[i, j]:<20} {t_grid[i, j]:<20} {yield_grid[i, j]:<20.1f}")
# Find the conditions giving the maximum conversion
max_idx = np.unravel_index(np.argmax(yield_grid), yield_grid.shape)
T_best = T_grid[max_idx]
t_best = t_grid[max_idx]
yield_best = yield_grid[max_idx]
print("-" * 70)
print(f"Optimal conditions: T = {T_best}°C, t = {t_best} hours")
print(f"Maximum conversion: {yield_best:.1f}%")
# 3D plot
fig = plt.figure(figsize=(14, 6))
# 3D surface plot
ax1 = fig.add_subplot(121, projection='3d')
T_fine = np.linspace(1000, 1400, 50)
t_fine = np.linspace(2, 10, 50)
T_mesh, t_mesh = np.meshgrid(T_fine, t_fine)
yield_mesh = np.zeros_like(T_mesh)
for i in range(len(t_fine)):
for j in range(len(T_fine)):
yield_mesh[i, j] = reaction_yield(T_mesh[i, j], t_mesh[i, j])
surf = ax1.plot_surface(T_mesh, t_mesh, yield_mesh, cmap='viridis',
alpha=0.8, edgecolor='none')
ax1.scatter(T_grid, t_grid, yield_grid, color='red', s=50,
label='Experimental points')
ax1.set_xlabel('Temperature (°C)', fontsize=10)
ax1.set_ylabel('Time (hours)', fontsize=10)
ax1.set_zlabel('Yield (%)', fontsize=10)
ax1.set_title('Response Surface', fontsize=12, fontweight='bold')
ax1.view_init(elev=25, azim=45)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=5)
# Contour plot
ax2 = fig.add_subplot(122)
contour = ax2.contourf(T_mesh, t_mesh, yield_mesh, levels=20, cmap='viridis')
ax2.contour(T_mesh, t_mesh, yield_mesh, levels=10, colors='black',
alpha=0.3, linewidths=0.5)
ax2.scatter(T_grid, t_grid, c=yield_grid, s=100, edgecolors='red',
linewidths=2, cmap='viridis')
ax2.scatter(T_best, t_best, color='red', s=300, marker='*',
edgecolors='white', linewidths=2, label='Optimum')
ax2.set_xlabel('Temperature (°C)', fontsize=11)
ax2.set_ylabel('Time (hours)', fontsize=11)
ax2.set_title('Contour Map', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
fig.colorbar(contour, ax=ax2, label='Yield (%)')
plt.tight_layout()
plt.savefig('doe_optimization.png', dpi=300, bbox_inches='tight')
plt.show()
In an actual solid-state reaction, DOE is applied following these steps:
Results from a research group that used DOE to optimize LiCoO₂ synthesis conditions:
# ===================================
# Example 7: Fitting Reaction Rate Curves
# ===================================
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Experimental data (time vs. conversion)
# Example: BaTiO3 synthesis @ 1200°C
time_exp = np.array([0, 1, 2, 3, 4, 6, 8, 10, 12, 15, 20]) # hours
conversion_exp = np.array([0, 0.15, 0.28, 0.38, 0.47, 0.60,
0.70, 0.78, 0.84, 0.90, 0.95])
# Jander equation model
def jander_model(t, k):
"""Calculate conversion using the Jander equation
Args:
t (array): Time [hours]
k (float): Rate constant
Returns:
array: Conversion
"""
# Solve [1 - (1-α)^(1/3)]² = kt for α
kt = k * t
alpha = 1 - (1 - np.sqrt(kt))**3
alpha = np.clip(alpha, 0, 1) # Restrict to the 0-1 range
return alpha
# Ginstling-Brounshtein equation (an alternative diffusion model)
def gb_model(t, k):
"""Ginstling-Brounshtein equation
Args:
t (array): Time
k (float): Rate constant
Returns:
array: Conversion
"""
# 1 - 2α/3 - (1-α)^(2/3) = kt
# This must be solved numerically; an approximate form is used here
kt = k * t
alpha = 1 - (1 - kt/2)**(3/2)
alpha = np.clip(alpha, 0, 1)
return alpha
# Power law (empirical formula)
def power_law_model(t, k, n):
"""Power-law model
Args:
t (array): Time
k (float): Rate constant
n (float): Exponent
Returns:
array: Conversion
"""
alpha = k * t**n
alpha = np.clip(alpha, 0, 1)
return alpha
# Fit each model
# Jander equation
popt_jander, _ = curve_fit(jander_model, time_exp, conversion_exp, p0=[0.01])
k_jander = popt_jander[0]
# Ginstling-Brounshtein equation
popt_gb, _ = curve_fit(gb_model, time_exp, conversion_exp, p0=[0.01])
k_gb = popt_gb[0]
# Power law
popt_power, _ = curve_fit(power_law_model, time_exp, conversion_exp, p0=[0.1, 0.5])
k_power, n_power = popt_power
# Generate the predicted curves
t_fit = np.linspace(0, 20, 200)
alpha_jander = jander_model(t_fit, k_jander)
alpha_gb = gb_model(t_fit, k_gb)
alpha_power = power_law_model(t_fit, k_power, n_power)
# Calculate the residuals
residuals_jander = conversion_exp - jander_model(time_exp, k_jander)
residuals_gb = conversion_exp - gb_model(time_exp, k_gb)
residuals_power = conversion_exp - power_law_model(time_exp, k_power, n_power)
# Calculate R²
def r_squared(y_true, y_pred):
ss_res = np.sum((y_true - y_pred)**2)
ss_tot = np.sum((y_true - np.mean(y_true))**2)
return 1 - (ss_res / ss_tot)
r2_jander = r_squared(conversion_exp, jander_model(time_exp, k_jander))
r2_gb = r_squared(conversion_exp, gb_model(time_exp, k_gb))
r2_power = r_squared(conversion_exp, power_law_model(time_exp, k_power, n_power))
# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
# Fitting results
ax1.plot(time_exp, conversion_exp, 'ko', markersize=8, label='Experimental data')
ax1.plot(t_fit, alpha_jander, 'b-', linewidth=2,
label=f'Jander (R²={r2_jander:.4f})')
ax1.plot(t_fit, alpha_gb, 'r-', linewidth=2,
label=f'Ginstling-Brounshtein (R²={r2_gb:.4f})')
ax1.plot(t_fit, alpha_power, 'g-', linewidth=2,
label=f'Power law (R²={r2_power:.4f})')
ax1.set_xlabel('Time (hours)', fontsize=12)
ax1.set_ylabel('Conversion', fontsize=12)
ax1.set_title('Kinetic Model Fitting', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, 20])
ax1.set_ylim([0, 1])
# Residual plot
ax2.plot(time_exp, residuals_jander, 'bo-', label='Jander')
ax2.plot(time_exp, residuals_gb, 'ro-', label='Ginstling-Brounshtein')
ax2.plot(time_exp, residuals_power, 'go-', label='Power law')
ax2.axhline(y=0, color='black', linestyle='--', linewidth=1)
ax2.set_xlabel('Time (hours)', fontsize=12)
ax2.set_ylabel('Residuals', fontsize=12)
ax2.set_title('Residual Plot', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kinetic_fitting.png', dpi=300, bbox_inches='tight')
plt.show()
# Results summary
print("\nReaction kinetics model fitting results:")
print("=" * 70)
print(f"{'Model':<25} {'Parameter':<30} {'R²':<10}")
print("-" * 70)
print(f"{'Jander':<25} {'k = ' + f'{k_jander:.4f} h⁻¹':<30} {r2_jander:.4f}")
print(f"{'Ginstling-Brounshtein':<25} {'k = ' + f'{k_gb:.4f} h⁻¹':<30} {r2_gb:.4f}")
print(f"{'Power law':<25} {'k = ' + f'{k_power:.4f}, n = {n_power:.4f}':<30} {r2_power:.4f}")
print("=" * 70)
print(f"\nBest model: {'Jander' if r2_jander == max(r2_jander, r2_gb, r2_power) else 'GB' if r2_gb == max(r2_jander, r2_gb, r2_power) else 'Power law'}")
# Example output:
# Reaction kinetics model fitting results:
# ======================================================================
# Model Parameter R²
# ----------------------------------------------------------------------
# Jander k = 0.0289 h⁻¹ 0.9953
# Ginstling-Brounshtein k = 0.0412 h⁻¹ 0.9867
# Power law k = 0.2156, n = 0.5234 0.9982
# ======================================================================
#
# Best model: Power law
In solid-state reactions, holding at high temperature for a long time causes undesirable grain growth. Strategies to suppress it include:
Mechanochemical methods (high-energy ball milling) can also drive solid-state reactions to proceed near room temperature:
# ===================================
# Example 8: Simulating Grain Growth
# ===================================
import numpy as np
import matplotlib.pyplot as plt
def grain_growth(t, T, D0, Ea, G0, n):
"""Time evolution of grain growth
Burke-Turnbull equation: G^n - G0^n = k*t
Args:
t (array): Time [hours]
T (float): Temperature [K]
D0 (float): Pre-exponential factor
Ea (float): Activation energy [J/mol]
G0 (float): Initial grain size [μm]
n (float): Grain-growth exponent (usually 2-4)
Returns:
array: Grain size [μm]
"""
R = 8.314
k = D0 * np.exp(-Ea / (R * T))
G = (G0**n + k * t * 3600)**(1/n) # hours → seconds
return G
# Parameter settings
D0_grain = 1e8 # μm^n/s
Ea_grain = 400e3 # J/mol
G0 = 0.5 # μm
n = 3
# Effect of temperature
temps_celsius = [1100, 1200, 1300]
t_range = np.linspace(0, 12, 100) # 0-12 hours
plt.figure(figsize=(12, 5))
# Temperature dependence
plt.subplot(1, 2, 1)
for T_c in temps_celsius:
T_k = T_c + 273.15
G = grain_growth(t_range, T_k, D0_grain, Ea_grain, G0, n)
plt.plot(t_range, G, linewidth=2, label=f'{T_c}°C')
plt.axhline(y=1.0, color='red', linestyle='--', linewidth=1,
label='Target grain size')
plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Grain Size (μm)', fontsize=12)
plt.title('Grain Growth at Different Temperatures', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 5])
# Effect of two-step sintering
plt.subplot(1, 2, 2)
# Conventional sintering: 1300°C, 6 hours
t_conv = np.linspace(0, 6, 100)
T_conv = 1300 + 273.15
G_conv = grain_growth(t_conv, T_conv, D0_grain, Ea_grain, G0, n)
# Two-step: 1300°C 1h → 1200°C 5h
t1 = np.linspace(0, 1, 20)
G1 = grain_growth(t1, 1300+273.15, D0_grain, Ea_grain, G0, n)
G_intermediate = G1[-1]
t2 = np.linspace(0, 5, 80)
G2 = grain_growth(t2, 1200+273.15, D0_grain, Ea_grain, G_intermediate, n)
t_two_step = np.concatenate([t1, t2 + 1])
G_two_step = np.concatenate([G1, G2])
plt.plot(t_conv, G_conv, 'r-', linewidth=2, label='Conventional (1300°C)')
plt.plot(t_two_step, G_two_step, 'b-', linewidth=2, label='Two-step (1300°C→1200°C)')
plt.axvline(x=1, color='gray', linestyle=':', linewidth=1, alpha=0.5)
plt.xlabel('Time (hours)', fontsize=12)
plt.ylabel('Grain Size (μm)', fontsize=12)
plt.title('Two-Step Sintering Strategy', fontsize=14, fontweight='bold')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.ylim([0, 5])
plt.tight_layout()
plt.savefig('grain_growth_control.png', dpi=300, bbox_inches='tight')
plt.show()
# Compare the final grain size
G_final_conv = grain_growth(6, 1300+273.15, D0_grain, Ea_grain, G0, n)
G_final_two_step = G_two_step[-1]
print("\nComparison of grain growth:")
print("=" * 50)
print(f"Conventional (1300°C, 6h): {G_final_conv:.2f} μm")
print(f"Two-step (1300°C 1h + 1200°C 5h): {G_final_two_step:.2f} μm")
print(f"Grain size suppression: {(1 - G_final_two_step/G_final_conv)*100:.1f}%")
# Example output:
# Comparison of grain growth:
# ==================================================
# Conventional (1300°C, 6h): 4.23 μm
# Two-step (1300°C 1h + 1200°C 5h): 2.87 μm
# Grain size suppression: 32.2%
Upon completing this chapter, you should be able to explain:
In the synthesis reaction BaCO₃ + TiO₂ → BaTiO₃ + CO₂ for BaTiO₃, which step is the slowest (rate-limiting)?
a) Release of CO₂
b) Nucleation of BaTiO₃
c) Diffusion of Ba²⁺ ions through the product layer
d) Chemical reaction at the interface
Correct answer: c) Diffusion of Ba²⁺ ions through the product layer
Explanation:
In a solid-state reaction, the product layer physically separates the reactants, so the process of ions diffusing through that layer is the slowest step.
Key point: Because the diffusion coefficient increases exponentially with temperature, choosing the right reaction temperature is critically important.
In the diffusion coefficient D(T) = D₀ exp(-Eₐ/RT), how does the sensitivity of the diffusion coefficient to temperature change as Eₐ (the activation energy) increases?
a) It increases (temperature dependence becomes stronger)
b) It decreases (temperature dependence becomes weaker)
c) It stays the same
d) There is no relationship
Correct answer: a) It increases (temperature dependence becomes stronger)
Explanation:
Because the activation energy Eₐ sits in the exponent of exp(-Eₐ/RT), a larger Eₐ means D changes more steeply with temperature.
Numerical example:
This is why temperature control is especially critical in systems with large activation energies.
According to the Jander equation k = D·C₀/r₀², if the particle radius r₀ is halved, by what factor does the rate constant k increase?
a) 2x
b) 4x
c) 1/2x
d) 1/4x
Correct answer: b) 4x
Calculation:
k ∝ 1/r₀²
When r₀ → r₀/2, k → k/(r₀/2)² = k/(r₀²/4) = 4k
Practical significance:
This is exactly why "grinding/refining" is so important in solid-state reactions.
In BaTiO₃ synthesis, the heating rate was changed from 20°C/min to 5°C/min. What is the most appropriate reason for this change?
a) To speed up the reaction rate
b) To prevent the sample from rupturing due to rapid release of CO₂
c) To save on electricity costs
d) To lower the crystallinity
Correct answer: b) To prevent the sample from rupturing due to rapid release of CO₂
Detailed explanation:
In the reaction BaCO₃ + TiO₂ → BaTiO₃ + CO₂, barium carbonate decomposes and releases CO₂ at 800-900°C.
Practical advice: In syntheses involving a decomposition reaction, slow the heating rate specifically through the relevant temperature range to control the rate of gas release (e.g., pass through 750-950°C at 2°C/min).
The following data were obtained from DSC measurements. Use the Kissinger method to determine the activation energy.
Heating rate β (K/min): 5, 10, 15
Peak temperature Tp (K): 1273, 1293, 1308
Kissinger equation: slope of ln(β/Tp²) vs. 1/Tp = -Eₐ/R
Answer:
Step 1: Organize the data
| β (K/min) | Tp (K) | ln(β/Tp²) | 1000/Tp (K⁻¹) |
|---|---|---|---|
| 5 | 1273 | -11.558 | 0.7855 |
| 10 | 1293 | -11.171 | 0.7734 |
| 15 | 1308 | -10.932 | 0.7645 |
Step 2: Linear regression
Plot y = ln(β/Tp²) vs. x = 1000/Tp
slope = Δy/Δx = (-10.932 - (-11.558)) / (0.7645 - 0.7855) = 0.626 / (-0.021) ≈ -29.8
Step 3: Calculate Eₐ
slope = -Eₐ / (R × 1000) (divided by 1000 because 1000/Tp was used)
Eₐ = -slope × R × 1000
Eₐ = 29.8 × 8.314 × 1000 = 247,757 J/mol ≈ 248 kJ/mol
Answer: Eₐ ≈ 248 kJ/mol
Physical interpretation:
This value falls within the typical range of activation energies (250-350 kJ/mol) for solid-state reactions in BaTiO₃ systems, and is thought to correspond to the solid-state diffusion of Ba²⁺ ions.
Using Design of Experiments, you examine two factors: temperature (1100, 1200, 1300°C) and time (4, 6, 8 hours). How many total experiments are needed? Also list two advantages of this approach compared to the conventional method of changing one factor at a time.
Answer:
Number of experiments:
3 levels × 3 levels = 9 runs (full factorial design)
Advantages of DOE (compared to the conventional method):
Additional benefits:
Design a temperature profile for synthesizing Li₁.₂Ni₀.₂Mn₀.₆O₂ (a lithium-rich cathode material) under the following conditions:
Explain the temperature profile (heating rate, holding temperature/time, cooling rate) and the reasoning behind the design.
Recommended temperature profile:
Phase 1: Preheating (Li₂CO₃ decomposition)
Phase 2: Intermediate heating (precursor formation)
Phase 3: Main firing (target-phase synthesis)
Phase 4: Cooling
Key design considerations:
Total time required: Approximately 30 hours (12h heating + 18h holding)
Alternative methods to consider:
From the following data, infer the reaction mechanism and calculate the activation energy.
Experimental data:
| Temperature (°C) | Time to reach 50% conversion, t₅₀ (hours) |
|---|---|
| 1000 | 18.5 |
| 1100 | 6.2 |
| 1200 | 2.5 |
| 1300 | 1.2 |
Assuming the Jander equation: [1-(1-0.5)^(1/3)]² = k·t₅₀
Answer:
Step 1: Calculate the rate constant k
For the Jander equation at α=0.5:
[1-(1-0.5)^(1/3)]² = [1-0.794]² = 0.206² = 0.0424
Therefore k = 0.0424 / t₅₀
| T (°C) | T (K) | t₅₀ (h) | k (h⁻¹) | ln(k) | 1000/T (K⁻¹) |
|---|---|---|---|---|---|
| 1000 | 1273 | 18.5 | 0.00229 | -6.080 | 0.7855 |
| 1100 | 1373 | 6.2 | 0.00684 | -4.985 | 0.7284 |
| 1200 | 1473 | 2.5 | 0.01696 | -4.077 | 0.6788 |
| 1300 | 1573 | 1.2 | 0.03533 | -3.343 | 0.6357 |
Step 2: Arrhenius plot
Plot ln(k) vs. 1/T (linear regression)
Linear fit: ln(k) = A - Eₐ/(R·T)
slope = -Eₐ/R
Linear regression calculation:
slope = Δ(ln k) / Δ(1000/T)
= (-3.343 - (-6.080)) / (0.6357 - 0.7855)
= 2.737 / (-0.1498)
= -18.27
Step 3: Calculate the activation energy
slope = -Eₐ / (R × 1000)
Eₐ = -slope × R × 1000
Eₐ = 18.27 × 8.314 × 1000
Eₐ = 151,899 J/mol ≈ 152 kJ/mol
Step 4: Discuss the reaction mechanism
Step 5: Suggested verification methods
Final conclusion:
Activation energy Eₐ = 152 kJ/mol
Inferred mechanism: interface-reaction control, or diffusion control in a fine-particle system
Additional experiments are recommended.
In this chapter, we studied the fundamentals of advanced ceramic materials (structural, functional, and bioceramics). In the next chapter, Chapter 2, we turn to advanced polymer materials (high-performance engineering plastics, functional polymers, and biodegradable polymers).