Chapter 5: Fundamentals of Additive Manufacturing

flowchart LR A[1986
SLAInvention
Chuck Hull] -->B[1988
SLSEmergence
Carl Deckard] B -->C[1992
FDMPatent
Stratasys] C -->D[2005
RepRap
Open Source] D -->E[2012
MetalAMAdoption
EBM/SLM] E -->F[2023
industryaddfast
Large-scale・Highfast] style A fill:#e3f2fd style B fill:#fff3e0 style C fill:#e8f5e9 style D fill:#f3e5f5 style E fill:#fce4ec style F fill:#fff9c4

flowchart TD AM[Additive Manufacturing
7 Processes] -->MEX[Material Extrusion
Material Extrusion] AM -->VPP[Vat Photopolymerization
Vat Photopolymerization] AM -->PBF[Powder Bed Fusion
Powder Bed Fusion] AM -->MJ[Material Jetting
Material Jetting] AM -->BJ[Binder Jetting
Binder Jetting] AM -->SL[Sheet Lamination
Sheet Lamination] AM -->DED[Directed Energy Deposition
Directed Energy Deposition] MEX -->MEX_EX[FDM/FFF
Lowcost・Adoptiontype] VPP -->VPP_EX[SLA/DLP
Highprecision・Highsurfacequality] PBF -->PBF_EX[SLS/SLM/EBM
Highstrength・metalcompatible] style AM fill:#f093fb style MEX fill:#e3f2fd style VPP fill:#fff3e0 style PBF fill:#e8f5e9 style MJ fill:#f3e5f5 style BJ fill:#fce4ec style SL fill:#fff9c4 style DED fill:#fce4ec

1.3 STL File Format and Data Processing

1.3.1 Structure of STL Files

STL (STereoLithography) isthe most widely used 3D model file format in AM, developed by 3D Systems in 1987.STL files represent object surfaces asa collection of triangle meshes.

Basic Structure of STL Files

ASCII STL format example:

solid cube facet normal 0 0 1 outer loop vertex 0 0 10 vertex 10 0 10 vertex 10 10 10 endloop endfacet facet normal 0 0 1 outer loop vertex 0 0 10 vertex 10 10 10 vertex 0 10 10 endloop endfacet ... endsolid cube

Two types of STL format:

1. ASCII STL: Human-readable text format. Large file size (10-20 times binary for same model). Useful for debugging and verification.
2. Binary STL: Binary format, small file size, fast processing. Standard in industrial applications. Structure: 80-byte header + 4 bytes (triangle count) + 50 bytes per triangle (normal 12B + vertices 36B + attribute 2B).

1.3.2 Key Concepts of STL Files

1. Normal Vector

Each triangle face has anormal vector (outward direction)defined to distinguish object "inside" and "outside". Normal direction is determined byright-hand rule:

Vertex order rule:Vertices v1, v2, v3 are arranged counter-clockwise (CCW), and when viewed from outside, the counter-clockwise order makes the normal point outward.

2. Manifold Conditions

For STL mesh to be 3D printable, it must bemanifold:

• Edge sharing: Every edge is shared by exactly two triangles
• Vertex sharing: Every vertex belongs to a continuous triangle fan
• Closed surface: Forms a completely closed surface without holes or openings
• No self-intersection: Triangles do not intersect or penetrate each other

1.3.3 Quality Metrics for STL Files

STL mesh quality is evaluated by the following metrics:

1. Triangle Count: Typically 10,000-500,000. Avoid too few (coarse model) or too many (large file size, processing delay).
2. Edge Length Uniformity: Mixed extremely large and small triangles degrade build quality. Ideally 0.1-1.0 mm range.
3. Aspect Ratio: Elongated triangles (high aspect ratio) cause numerical errors. Ideally aspect ratio< 10.
4. Normal Consistency: All normals unified outward. Mixed inverted normals cause inside/outside determination errors.

1.3.4 STL Processing with Python Libraries

Major libraries for handling STL files in Python:

1. numpy-stl: Fast STL read/write, volume/surface area calculation, normal vector manipulation. Simple and lightweight.
2. trimesh: Comprehensive 3D mesh processing library. Mesh repair, Boolean operations, ray casting, collision detection. Feature-rich but many dependencies.
3. PyMesh: Advanced mesh processing (remeshing, subdivision, feature extraction). Installation somewhat complex.

Basic usage of numpy-stl:

from stl import mesh import numpy as np # Load STL file your_mesh = mesh.Mesh.from_file('model.stl') # Basic geometric information volume, cog, inertia = your_mesh.get_mass_properties() print(f"Volume: {volume:.2f} mm³") print(f"Center of Gravity: {cog}") print(f"Surface Area: {your_mesh.areas.sum():.2f} mm²") # Number of triangles print(f"Number of Triangles: {len(your_mesh.vectors)}")

1.4 Slicing and Toolpath Generation

The process of converting STL files into commands (G-code) that 3D printers can understand is calledslicing. This section covers the basic principles of slicing, toolpath strategies, and G-code fundamentals.

1.4.1 Basic Principles of Slicing

Slicing is the process of horizontally cutting a 3D model at constant height (layer height) and extracting the contour of each layer:

flowchart TD A[3D Model
STL File] -->B[Slice in layers
along Z-axis] B -->C[Extract contour of each layer
Contour Detection] C -->D[Generate shell
Perimeter Path] D -->E[Generate infill
Infill Path] E -->F[Add support
Support Structure] F -->G[Optimal toolpath
Retraction/Travel] G -->H[G-codeoutput] style A fill:#e3f2fd style H fill:#e8f5e9

Layer Height Selection

Layer height is the most important parameter determining the trade-off between build quality and build time:

1.4.2 Shell and Infill Strategies

Shell (Perimeter) Generation

Shell (Perimeter)is the path forming the outer perimeter of each layer:

• Perimeter Count: Typically 2-4. Affects external quality and strength.
  • 1: Very weak, high transparency, decorative only
  • 2 lines: Standard (good balance)
  • 3-4: High strength, improved surface quality, improved airtightness
• Shell order: Inside-out is common. Outside-in used when surface quality is priority.

Infill Pattern

Infillforms internal structure and controls strength and material usage:

1.4.3 Generation of Support Structures

Parts with overhang angle exceeding 45 degrees requiresupport structures:

Support Types

• Linear Support (straight support): Vertical pillar-like support. Simple and easy to remove, but uses more material。
• Tree Support (tree support): Tree-like branching support. Material usage reduced by 30-50%, easy to remove. Standard in Cura and PrusaSlicer。
• Interface Layers: Thin interface layer on support top. Easy to remove, improves surface quality. Typically 2-4 layers.

Important Support Parameters

1.4.4 Fundamentals of G-code

G-code is the standard numerical control language for controlling 3D printers and CNC machines. Each line represents one command:

Major G-code Commands

G-codeExample（Build startpart minutes）

; === Start G-code === M140 S60 ; Start heating bed to 60°C (non-blocking) M104 S210 ; Start heating nozzle to 210°C (non-blocking) G28 ; Home all axes G29 ; Auto-leveling (bed mesh measurement) M190 S60 ; Wait for bed temperature to reach target M109 S210 ; Wait for nozzle temperature to reach target G92 E0 ; Reset extrusion to zero G1 Z2.0 F3000 ; Raise Z-axis 2mm (safety) G1 X10 Y10 F5000 ; Move to prime position G1 Z0.3 F3000 ; Lower Z-axis to 0.3mm (first layer height) G1 X100 E10 F1500 ; Draw prime line (clear nozzle clog) G92 E0 ; Reset extrusion to zero again ; === Build start ===

1.4.5 Major Slicing Software

1.4.6 Toolpath Optimization Strategies

Efficient toolpath improves build time, quality, and material usage:

• Retraction: Retracting filament during movement to prevent stringing。
  • Distance: 1-6mm (Bowden tube 4-6mm, direct 1-2mm)
  • Speed: 25-45 mm/s
  • Excessive retraction causes nozzle clogging
• Z-hop (Z-axis hopping): Raising nozzle during movement to avoid collision with built object。0.2-0.5mm raise. Slightly increases build time but improves surface quality。
• Combing: Restricting movement paths above infill to reduce surface movement marks。Effective when appearance is important。
• Seam Position: Strategy for aligning layer start/end points.
  • Random: Random placement (inconspicuous)
  • Aligned: Aligned placement (easy to remove seam in post-processing)
  • Sharpest Corner: Place at sharpest corner (less noticeable)

Example 1: Loading STL Files and Obtaining Basic Information

# =================================== # Example 1: Load STL file and obtain basic information # =================================== import numpy as np from stl import mesh # Load STL file your_mesh = mesh.Mesh.from_file('model.stl') # Get basic geometric information volume, cog, inertia = your_mesh.get_mass_properties() print("=== STL File Basic Information ===") print(f"Volume: {volume:.2f} mm³") print(f"Surface Area: {your_mesh.areas.sum():.2f} mm²") print(f"Center of Gravity: [{cog[0]:.2f}, {cog[1]:.2f}, {cog[2]:.2f}] mm") print(f"Number of Triangles: {len(your_mesh.vectors)}") # Calculate bounding box min_coords = your_mesh.vectors.min(axis=(0, 1)) max_coords = your_mesh.vectors.max(axis=(0, 1)) dimensions = max_coords - min_coords print(f"\n=== Bounding Box ===") print(f"X: {min_coords[0]:.2f} to {max_coords[0]:.2f} mm (Width: {dimensions[0]:.2f} mm)") print(f"Y: {min_coords[1]:.2f} to {max_coords[1]:.2f} mm (Depth: {dimensions[1]:.2f} mm)") print(f"Z: {min_coords[2]:.2f} to {max_coords[2]:.2f} mm (Height: {dimensions[2]:.2f} mm)") # Build time simple estimation (Layer Height 0.2mm, speed 50mm/s assumed) layer_height = 0.2 # mm print_speed = 50 # mm/s num_layers = int(dimensions[2] / layer_height) # Simple calculation: estimation based on surface area estimated_path_length = your_mesh.areas.sum() / layer_height # mm estimated_time_seconds = estimated_path_length / print_speed estimated_time_minutes = estimated_time_seconds / 60 print(f"\n=== Build Estimation ===") print(f"Number of layers（0.2mm/layer）: {num_layers} layer") print(f"Estimated build time: {estimated_time_minutes:.1f} minutes ({estimated_time_minutes/60:.2f} hours)") # Output example: # === STL File Basic Information === # Volume: 12450.75 mm³ # Surface Area: 5832.42 mm² # Center of Gravity: [25.34, 18.92, 15.67] mm # Number of Triangles: 2456 # # === Bounding Box === # X: 0.00 to 50.00 mm (Width: 50.00 mm) # Y: 0.00 to 40.00 mm (Depth: 40.00 mm) # Z: 0.00 to 30.00 mm (Height: 30.00 mm) # # === Build Estimation === # Number of layers（0.2mm/layer）: 150 layer # Estimated build time: 97.2 minutes (1.62 hours)

Example 2: Mesh Normal Vector Verification

# =================================== # Example 2: Mesh normal vector verification # =================================== import numpy as np from stl import mesh def check_normals(mesh_data): """Check STL mesh normal vector consistency Args: mesh_data: numpy-stl Mesh object Returns: tuple: (flipped_count, total_count, percentage) """ # Confirm normal direction according to right-hand rule flipped_count = 0 total_count = len(mesh_data.vectors) for i, facet in enumerate(mesh_data.vectors): v0, v1, v2 = facet # Calculate edge vectors edge1 = v1 - v0 edge2 = v2 - v0 # Calculate normal by cross product (right-hand system) calculated_normal = np.cross(edge1, edge2) # Normalize norm = np.linalg.norm(calculated_normal) if norm >1e-10: # Confirm not zero vector calculated_normal = calculated_normal / norm else: continue # Skip degenerate triangles # Compare with stored normal in file stored_normal = mesh_data.normals[i] stored_norm = np.linalg.norm(stored_normal) if stored_norm >1e-10: stored_normal = stored_normal / stored_norm # Check direction match with dot product dot_product = np.dot(calculated_normal, stored_normal) # If dot product is negative, opposite orientation if dot_product< 0: flipped_count += 1 percentage = (flipped_count / total_count) * 100 if total_count >0 else 0 return flipped_count, total_count, percentage # Load STL file your_mesh = mesh.Mesh.from_file('model.stl') # Execute normal check flipped, total, percent = check_normals(your_mesh) print("=== Normal Vector Verification Results ===") print(f"Total triangles: {total}") print(f"Inverted normals: {flipped}") print(f"Inversion rate: {percent:.2f}%") if flipped == 0: print("\n✅ All normals are correctly oriented") print(" This mesh is 3D printable") elif percent< 5: print("\n⚠️ Some normals are inverted (minor)") print(" High possibility slicer can auto-correct") else: print("\n❌ Many normals are inverted (critical)") print(" Repair with mesh repair tools (Meshmixer, netfabb) recommended") # Output example: # === Normal Vector Verification Results === # Total triangles: 2456 # Inverted normals: 0 # Inversion rate: 0.00% # # ✅ All normals are correctly oriented # This mesh is 3D printable

Example 3: Manifold Check

# =================================== # Example 3: Manifold property (Watertight) check # =================================== import trimesh # Load STL file (trimesh attempts automatic repair) mesh = trimesh.load('model.stl') print("=== Mesh Quality Diagnostics ===") # Basic information print(f"Vertex count: {len(mesh.vertices)}") print(f"Face count: {len(mesh.faces)}") print(f"Volume: {mesh.volume:.2f} mm³") # Manifold property check print(f"\n=== 3D Printability Check ===") print(f"Is watertight: {mesh.is_watertight}") print(f"Is winding consistent: {mesh.is_winding_consistent}") print(f"Is valid: {mesh.is_valid}") # Detailed problem diagnosis if not mesh.is_watertight: # Detect number of holes try: edges = mesh.edges_unique edges_sorted = mesh.edges_sorted duplicate_edges = len(edges_sorted) - len(edges) print(f"\n⚠️ Problem detected:") print(f" - Mesh has holes") print(f" - Duplicate edge count: {duplicate_edges}") except: print(f"\n⚠️ Mesh structure problem exists") # Attempt repair if not mesh.is_watertight or not mesh.is_winding_consistent: print(f"\n🔧 Executing automatic repair...") # Correct normals trimesh.repair.fix_normals(mesh) print(" ✓ Fixed normal vectors") # Fill holes trimesh.repair.fill_holes(mesh) print(" ✓ Filled holes") # Remove degenerate triangles mesh.remove_degenerate_faces() print(" ✓ Removed degenerate faces") # Merge duplicate vertices mesh.merge_vertices() print(" ✓ Merged duplicate vertices") # Confirm state after repair print(f"\n=== Post-Repair Status ===") print(f"Is watertight: {mesh.is_watertight}") print(f"Is winding consistent: {mesh.is_winding_consistent}") # Save repaired mesh if mesh.is_watertight: mesh.export('model_repaired.stl') print(f"\n✅ Repair complete! Saved as model_repaired.stl") else: print(f"\n❌ Automatic repair failed. Specialized tools like Meshmixer recommended") else: print(f"\n✅ This mesh is 3D printable") # Output example: # === Mesh Quality Diagnostics === # Vertex count: 1534 # Face count: 2456 # Volume: 12450.75 mm³ # # === 3D Printability Check === # Is watertight: True # Is winding consistent: True # Is valid: True # # ✅ This mesh is 3D printable

Example 4: Basic Slicing Algorithm

# =================================== # Example 4: Basic slicing algorithm # =================================== import numpy as np from stl import mesh def slice_mesh_at_height(mesh_data, z_height): """Generate temperature profile Args: t (array): Time array [min] T_target (float): Hold temperature [°C] heating_rate (float): Heating rate [°C/min] hold_time (float): Hold 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 time t_heat = (T_target - T_room) / heating_rate # Cooling start time 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) # Not below room temperature return T def simulate_reaction_progress(T, t, Ea, D0, r0): """Calculate reaction progress based on temperature profile Args: T (array): Temperature profile [°C] t (array): Time array [min] Ea (float): Activation energy [J/mol] D0 (float): Frequency factor [m²/s] r0 (float): Particle radius [m] Returns: array: Reaction rate """ 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 (infinitesimal time reaction progress) 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 setting T_target = 1200 # °C hold_time = 240 # min (4 hours) Ea = 300e3 # J/mol D0 = 5e-4 # m²/s r0 = 5e-6 # m # Comparison at different 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 profile 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 time to reach 95% reaction for each heating rate print("\nComparison of time to reach 95% reaction:") 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 to reach 95% 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: Incomplete reaction") # Output example: # Comparison of time to reach 95% reaction: # ============================================================ # 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

Exercises

1.5.1 What is pycalphad

pycalphadis a Python library for phase diagram calculations based on the CALPHAD (CALculation of PHAse Diagrams) method. It calculates equilibrium phases from thermodynamic databases and is useful for reaction pathway design.

1.5.2 Binary Phase Diagram Calculation Example

# =================================== # Example 5: Phase diagram calculation with pycalphad # =================================== # Note: pycalphad installation required # pip install pycalphad from pycalphad import Database, equilibrium, variables as v import matplotlib.pyplot as plt import numpy as np # Load TDB database (simplified example here) # Actually requires proper TDB file # Example: BaO-TiO2 system # Simplified TDB string (more complex in practice) 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: Actual calculations require official TDB file # Limited to conceptual explanation here print("Concept of phase diagram calculation with pycalphad:") print("=" * 60) print("1. Load TDB database (thermodynamic data)") print("2. Set temperature and composition range") print("3. Execute equilibrium calculation") print("4. Visualize stable phases") print() print("Practical application examples:") print("- BaO-TiO2 system: Temperature and composition range for BaTiO3 formation") print("- Si-Nsystem: Si3N4 stability region") print("- Phase relationships in multicomponent ceramics") # Conceptual plot (image based on actual database) fig, ax = plt.subplots(figsize=(10, 7)) # Temperature range T = np.linspace(800, 1600, 100) # Each phase stability region (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() # Practical usage example (commented out) """ # Practical example using pycalphad db = Database('BaO-TiO2.tdb') # Load 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 result eq.plot() """

1.6 Condition Optimization with Design of Experiments (DOE)

1.6.1 What is DOE

Design of Experiments (DOE) is a statistical technique to find optimal conditions with minimum number of experiments in systems with multiple interacting parameters。

Main parameters to optimize in solid-state reactions:

• reactionTemperature（T）
• Hold time (t)
• Particle size (r)
• Raw material ratio (molar ratio)
• Atmosphere (air, nitrogen, vacuum, etc.)

1.6.2 Response Surface Methodology

# =================================== # Example 6: Condition optimization with DOE # =================================== import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy.optimize import minimize # Virtual reaction yield model (Temperature and time function) def reaction_yield(T, t, noise=0): """Calculate reaction yield from temperature and time (virtual model) Args: T (float): Temperature [°C] t (float): Time [hours] noise (float): Noise level Returns: float: Reaction rate [%] """ # Optimal value: T=1200°C, t=6 hours T_opt = 1200 t_opt = 6 # 2nd order model (Gaussian type) 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) # Experiment point placement (central composite design) T_levels = [1000, 1100, 1200, 1300, 1400] # °C t_levels = [2, 4, 6, 8, 10] # hours # Grid experiment point placement T_grid, t_grid = np.meshgrid(T_levels, t_levels) yield_grid = np.zeros_like(T_grid, dtype=float) # Measure reaction yield at each experiment point (simulation) 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 results print("Reaction Condition Optimization 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 condition with highest reaction yield 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 condition: T = {T_best}°C, t = {t_best} hours") print(f"Highest reaction yield: {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()

1.6.3 Practical Approach to Experimental Design

In actual solid-state reactions, apply DOE with the following procedure:

1. Screening experiment (two-level factorial design): Identify parameters with large influence
2. Response surface method (central composite design): Search for optimal conditions
3. Confirmation experiment: Conduct experiment at predicted optimal condition to verify model

1.7 Fitting of Reaction Rate Curves

1.7.1 Determination of Rate Constants from Experimental Data

# =================================== # Example 7: Reaction kinetics curve fitting # =================================== 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 reaction rate using Jander equation Args: t (array): Time [hours] k (float): Rate constant Returns: array: Reaction rate """ # [1 - (1-α)^(1/3)]² = kt solve for α kt = k * t alpha = 1 - (1 - np.sqrt(kt))**3 alpha = np.clip(alpha, 0, 1) # Limit to 0-1 range return alpha # Ginstling-Brounshtein equation (alternative diffusion model) def gb_model(t, k): """Ginstling-Brounshtein equation Args: t (array): hours k (float): Rate constant Returns: array: Reaction rate """ # 1 - 2α/3 - (1-α)^(2/3) = kt # Requires numerical solution, using approximate form here kt = k * t alpha = 1 - (1 - kt/2)**(3/2) alpha = np.clip(alpha, 0, 1) return alpha # Power law (empirical equation) def power_law_model(t, k, n): """Power law model Args: t (array): hours k (float): Rate constant n (float): Exponent Returns: array: Reaction rate """ 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 prediction 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 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() # Result summary print("\nKinetic 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'}") # Output example: # Kinetic 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

1.8 Advanced Topics: Microstructure Control

1.8.1 Suppression of Grain Growth

In solid-state reactions, high temperature and long hold times can cause undesirable grain growth. Strategies to suppress this:

• Two-step sintering: Short hold at high temperature followed by long hold at low temperature
• Use of additives: Add small amounts of grain growth inhibitors (e.g., MgO, Al₂O₃)
• Spark Plasma Sintering (SPS): Rapid heating and short sintering time

1.8.2 Mechanochemical Activation of Reactions

Using mechanochemical methods (high-energy ball milling), solid-state reactions can also proceed near room temperature:

# =================================== # Example 8: Grain growth simulation # =================================== import numpy as np import matplotlib.pyplot as plt def grain_growth(t, T, D0, Ea, G0, n): """Grain size time evolution Burke-Turnbull equation: G^n - G0^n = k*t Args: t (array): Time [hours] T (float): Temperature [K] D0 (float): Frequency factor Ea (float): Activation energy [J/mol] G0 (float): Initial grain size [μm] n (float): particleachievelongExponent（ usually2-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 setting D0_grain = 1e8 # μm^n/s Ea_grain = 400e3 # J/mol G0 = 0.5 # μm n = 3 # Temperature influence temps_celsius = [1100, 1200, 1300] t_range = np.linspace(0, 12, 100) # 0-12 hours plt.figure(figsize=(12, 5)) # Temperature dependency 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]) # Two-step sintering effect 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() # Final grain size comparison G_final_conv = grain_growth(6, 1300+273.15, D0_grain, Ea_grain, G0, n) G_final_two_step = G_two_step[-1] print("\nGrain growth comparison:") 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 effect: {(1 - G_final_two_step/G_final_conv)*100:.1f}%") # Output example: # Grain growth comparison: # ================================================== # Conventional (1300°C, 6h): 4.23 μm # Two-step (1300°C 1h + 1200°C 5h): 2.87 μm # Grain size suppression effect: 32.2%

Learning Objectivesconfirm

Upon completing this chapter, you will be able to explain:

Basic Understanding

• ✅ Can explain the three rate-limiting stages of solid-state reactions (nucleation, interface reaction, diffusion)
• ✅ Understands the physical meaning of the Arrhenius equation and temperature dependency
• ✅ Can explain the difference between Jander and Ginstling-Brounshtein equations
• ✅ Understands the three essential elements of temperature profiles (heating rate, hold time, cooling rate))

Practical Skills

• ✅ Can simulate temperature dependency of diffusion coefficient with Python
• ✅ Can predict reaction progress using Jander equation
• ✅ Can calculate activation energy from DSC/TG data using Kissinger method
• ✅ Can optimize reaction conditions using DOE (Design of Experiments)
• ✅ Understands basics of phase diagram calculation using pycalphad

Applied Expertise

• ✅ Can design synthesis process for new ceramic materials
• ✅ Can infer reaction mechanism from experimental data and select appropriate kinetic equation
• ✅ Can propose optimization strategy for industrial process conditions
• ✅ Can propose grain growth control strategies (two-step sintering, etc.)

Exercises

Easy (Fundamentals)

Medium (Application)

import trimesh mesh = trimesh.load('model.stl') # Add code here: Check manifold property, # if problem exists, automatically repair, then # save repaired mesh as 'model_fixed.stl'

import trimesh mesh = trimesh.load('model.stl') # Manifold property check print(f"Is watertight: {mesh.is_watertight}") print(f"Is winding consistent: {mesh.is_winding_consistent}") # Repair if problem exists if not mesh.is_watertight or not mesh.is_winding_consistent: print("Executing mesh repair...") # Correct normals trimesh.repair.fix_normals(mesh) # Fill holes trimesh.repair.fill_holes(mesh) # Remove degenerate triangles mesh.remove_degenerate_faces() # Merge duplicate vertices mesh.merge_vertices() # Confirm repair result print(f"Watertight after repair: {mesh.is_watertight}") # Save repaired mesh if mesh.is_watertight: mesh.export('model_fixed.stl') print("Repair complete: Saved as model_fixed.stl") else: print("⚠️ Automatic repair failed. Please use Meshmixer or other tools") else: print("✓ Mesh is 3D printable")

DOE Advantages (Compared to Conventional Methods):

1. Can Detect Interactions
   • Conventional method: Evaluate temperature influence and time influence separately
   • DOE: Can determine interactions like "shorter time works at high temperature"
   • Example: 1300°C4hours10 minutes、1100°C8hoursnecessary、etc.
2. experimenttimesnumberreduction
   • Conventional method (OFAT: One Factor At a Time):
     • Temperature study: 3 times (fixed time)
     • Time study: 3 times (fixed temperature)
     • confirmexperiment: multipletimes
     • Total: 10+ times
   • DOE: Complete in 9 times (all conditions covered + interaction analysis)
   • Furthermore, using central composite design can reduce to 7 times

Additional Advantages:

• Obtain statistically significant conclusions (can evaluate variance)
• Construct response surface, can predict untested conditions
• Can detect cases where optimal condition is outside experimental range

Hard (Advanced)

nextstep

In Chapter 5, we learned the fundamentals of additive manufacturing (AM), including the seven process categories according to ISO/ASTM 52900, STL file format structure, slicing, and G-code basics. In Chapter 2, we will learn about the detailed build process, material properties, and process parameter optimization for Material Extrusion (FDM/FFF)。

References

1. Gibson, I., Rosen, D., & Stucker, B. (2015).Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping, and Direct Digital Manufacturing(2nd ed.). Springer. pp. 1-35, 89-145, 287-334. - Comprehensive AM technology textbook with detailed explanation of seven process categories and STL data processing
2. ISO/ASTM 52900:2021.Additive manufacturing — General principles — Fundamentals and vocabulary. International Organization for Standardization. - AM terminology and process category international standard specification, widely referenced in manufacturing industry
3. Kruth, J.P., Leu, M.C., & Nakagawa, T. (1998). "Progress in Additive Manufacturing and Rapid Prototyping."CIRP Annals - Manufacturing Technology, 47(2), 525-540. - Theoretical basis of selective laser sintering and binding mechanisms
4. Hull, C.W. (1986).Apparatus for production of three-dimensional objects by stereolithography. US Patent 4,575,330. - World's first AM technology (SLA) patent, important document that became the origin of AM industry
5. Wohlers, T. (2023).Wohlers Report 2023: 3D Printing and Additive Manufacturing Global State of the Industry. Wohlers Associates, Inc. pp. 15-89, 156-234. - Latest statistical report on AM market trends and industrial applications, updated annually as industry standard reference
6. 3D Systems, Inc. (1988).StereoLithography Interface Specification. - STL file format official specification document, defining ASCII/Binary STL structure
7. numpy-stl Documentation. (2024).Python library for working with STL files.https://numpy-stl.readthedocs.io/- Python library for STL file loading and volume calculation
8. trimesh Documentation. (2024).Python library for loading and using triangular meshes.https://trimsh.org/- Comprehensive library for mesh repair, Boolean operations, and quality evaluation

usefortoolandlibrary

• NumPy(v1.24+): numbervaluecalculationlibrary -https://numpy.org/
• numpy-stl(v3.0+): STLfileprocessinglibrary -https://numpy-stl.readthedocs.io/
• trimesh(v4.0+): 3D mesh processing library (repair, verification, Boolean operations) -https://trimsh.org/
• Matplotlib(v3.7+): Data visualization library -https://matplotlib.org/
• SciPy(v1.10+): Scientific computing library (optimization, interpolation) -https://scipy.org/