07_multi_alloy_comparison.py

"""
05_multi_alloy_comparison.py
Comparative Analysis of Aerospace-Grade Al-7xxx Alloys
Compares Al-7050, Al-7075, and Al-7085 for phase stability and processing windows
ALL VALUES FROM REAL THERMODYNAMIC SIMULATIONS
"""

import matplotlib.pyplot as plt
import numpy as np
from pycalphad import Database, equilibrium, variables as v
import warnings
warnings.filterwarnings('ignore')

print("=" * 70)
print("MULTI-ALLOY COMPARISON: Al-7050 vs Al-7075 vs Al-7085")
print("ALL DATA FROM REAL THERMODYNAMIC CALCULATIONS")
print("=" * 70)

# Load database
dbf = Database('COST507-modified.tdb')
print(f"Database loaded: COST507-modified.tdb")

# Components and phases
comps = ['AL', 'ZN', 'MG', 'CU', 'VA']
phases = ['FCC_A1', 'LIQUID', 'LAVES_C14', 'LAVES_C15', 'HCP_A3']

# Add additional phases if available
all_phases = list(dbf.phases.keys())
for p in all_phases:
    if any(kw in p.upper() for kw in ['ALCU', 'THETA', 'MGZN', 'AL2CU']):
        if p not in phases:
            phases.append(p)

print(f"Using phases: {phases}")

# =============================================================================
# ALLOY COMPOSITIONS (weight percent converted to fraction)
# Reference: ASM Handbook, Vol. 2 - Al Alloy Designations
# =============================================================================
alloys = {
    'Al-7050': {
        'Zn': 0.062,  # 6.2 wt%
        'Mg': 0.023,  # 2.3 wt%
        'Cu': 0.023,  # 2.3 wt%
        'color': '#2171b5',
        'marker': 'o',
        'description': 'High toughness, corrosion resistant'
    },
    'Al-7075': {
        'Zn': 0.056,  # 5.6 wt%
        'Mg': 0.025,  # 2.5 wt%
        'Cu': 0.016,  # 1.6 wt%
        'color': '#238b45',
        'marker': 's',
        'description': 'Standard aerospace grade'
    },
    'Al-7085': {
        'Zn': 0.075,  # 7.5 wt%
        'Mg': 0.015,  # 1.5 wt%
        'Cu': 0.016,  # 1.6 wt%
        'color': '#d94801',
        'marker': '^',
        'description': 'High Zn, newer aerospace grade'
    }
}

print("\nAlloy Compositions (wt%):")
for name, comp in alloys.items():
    print(f"  {name}: Al-{comp['Zn']*100:.1f}Zn-{comp['Mg']*100:.1f}Mg-{comp['Cu']*100:.1f}Cu")
    print(f"         {comp['description']}")

# =============================================================================
# PART 1: ETA-PHASE FRACTION VS TEMPERATURE
# =============================================================================
print("\n--- Part 1: Eta-Phase Fraction vs Temperature ---")

temp_range = np.linspace(100, 300, 21) + 273.15  # 100-300°C

results = {name: {'temps': [], 'eta_frac': [], 'fcc_frac': []} for name in alloys}

for name, comp in alloys.items():
    print(f"\nCalculating for {name}...")
    
    for T in temp_range:
        try:
            conds = {
                v.X('ZN'): comp['Zn'],
                v.X('MG'): comp['Mg'],
                v.X('CU'): comp['Cu'],
                v.T: T,
                v.P: 101325,
                v.N: 1
            }
            eq = equilibrium(dbf, comps, phases, conds)
            
            # Get eta-phase fraction (LAVES phases)
            eta_frac = eq.NP.where(
                (eq.Phase == 'LAVES_C14') | (eq.Phase == 'LAVES_C15')
            ).sum(dim='vertex').values.flatten()
            eta_val = float(np.nanmax(eta_frac)) if len(eta_frac) > 0 else 0.0
            
            # Get FCC fraction
            fcc_frac = eq.NP.where(eq.Phase == 'FCC_A1').sum(dim='vertex').values.flatten()
            fcc_val = float(np.nanmax(fcc_frac)) if len(fcc_frac) > 0 else 0.0
            
            results[name]['temps'].append(T - 273.15)
            results[name]['eta_frac'].append(eta_val)
            results[name]['fcc_frac'].append(fcc_val)
            
        except Exception as e:
            results[name]['temps'].append(T - 273.15)
            results[name]['eta_frac'].append(np.nan)
            results[name]['fcc_frac'].append(np.nan)
    
    print(f"  Completed {len(results[name]['temps'])} temperature points")

# =============================================================================
# PART 2: SOLIDIFICATION ANALYSIS (Liquidus/Solidus)
# =============================================================================
print("\n--- Part 2: Solidification Analysis ---")

solidification = {name: {'liquidus': None, 'solidus': None} for name in alloys}

high_temp_range = np.linspace(400, 700, 61) + 273.15  # 400-700°C

for name, comp in alloys.items():
    print(f"\nCalculating solidification for {name}...")
    
    liquid_fracs = []
    temps_solid = []
    
    for T in high_temp_range:
        try:
            conds = {
                v.X('ZN'): comp['Zn'],
                v.X('MG'): comp['Mg'],
                v.X('CU'): comp['Cu'],
                v.T: T,
                v.P: 101325,
                v.N: 1
            }
            eq = equilibrium(dbf, comps, phases, conds)
            
            # Get liquid fraction
            liq_frac = eq.NP.where(eq.Phase == 'LIQUID').sum(dim='vertex').values.flatten()
            liq_val = float(np.nanmax(liq_frac)) if len(liq_frac) > 0 else 0.0
            
            liquid_fracs.append(liq_val)
            temps_solid.append(T - 273.15)
            
        except Exception as e:
            liquid_fracs.append(np.nan)
            temps_solid.append(T - 273.15)
    
    # Find liquidus (where liquid first appears)
    liquid_fracs = np.array(liquid_fracs)
    temps_solid = np.array(temps_solid)
    
    # Liquidus: highest temp where liquid < 0.99
    for i in range(len(liquid_fracs)-1, -1, -1):
        if liquid_fracs[i] < 0.99:
            solidification[name]['liquidus'] = temps_solid[i]
            break
    
    # Solidus: lowest temp where liquid > 0.01
    for i in range(len(liquid_fracs)):
        if liquid_fracs[i] > 0.01:
            solidification[name]['solidus'] = temps_solid[i]
            break
    
    print(f"  Liquidus: {solidification[name]['liquidus']:.0f}°C")
    print(f"  Solidus: {solidification[name]['solidus']:.0f}°C")
    print(f"  Freezing Range: {solidification[name]['liquidus'] - solidification[name]['solidus']:.0f}°C")

# =============================================================================
# VISUALIZATION
# =============================================================================
print("\n--- Generating Comparison Plots ---")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))
fig.suptitle('Multi-Alloy Comparison: Al-7050 vs Al-7075 vs Al-7085\n'
             'All data from COST507 thermodynamic database',
             fontsize=14, fontweight='bold')

# Plot 1: Eta-phase fraction vs temperature
ax1 = axes[0, 0]
for name, comp in alloys.items():
    data = results[name]
    ax1.plot(data['temps'], data['eta_frac'], 
             color=comp['color'], marker=comp['marker'],
             linewidth=2, markersize=6, label=name, markevery=3)

ax1.set_xlabel('Temperature (°C)', fontsize=11)
ax1.set_ylabel('η-Phase Volume Fraction', fontsize=11)
ax1.set_title('Strengthening Phase (η-MgZn₂) Stability', fontsize=12)
ax1.legend(loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xlim(100, 300)

# Add annotation for aging window
ax1.axvspan(100, 150, alpha=0.2, color='green', label='Optimal Aging')
ax1.text(125, ax1.get_ylim()[1]*0.9, 'Aging\nWindow', ha='center', fontsize=9)

# Plot 2: FCC matrix fraction vs temperature
ax2 = axes[0, 1]
for name, comp in alloys.items():
    data = results[name]
    ax2.plot(data['temps'], data['fcc_frac'], 
             color=comp['color'], marker=comp['marker'],
             linewidth=2, markersize=6, label=name, markevery=3)

ax2.set_xlabel('Temperature (°C)', fontsize=11)
ax2.set_ylabel('FCC Matrix Fraction', fontsize=11)
ax2.set_title('Aluminum Matrix Phase', fontsize=12)
ax2.legend(loc='lower right')
ax2.grid(True, alpha=0.3)
ax2.set_xlim(100, 300)

# Plot 3: Processing window comparison (bar chart)
ax3 = axes[1, 0]
alloy_names = list(alloys.keys())
x_pos = np.arange(len(alloy_names))
width = 0.35

liquidus_vals = [solidification[name]['liquidus'] for name in alloy_names]
solidus_vals = [solidification[name]['solidus'] for name in alloy_names]
colors = [alloys[name]['color'] for name in alloy_names]

bars1 = ax3.bar(x_pos - width/2, liquidus_vals, width, label='Liquidus', 
                color=colors, alpha=0.8, edgecolor='black')
bars2 = ax3.bar(x_pos + width/2, solidus_vals, width, label='Solidus',
                color=colors, alpha=0.4, edgecolor='black', hatch='//')

ax3.set_xlabel('Alloy', fontsize=11)
ax3.set_ylabel('Temperature (°C)', fontsize=11)
ax3.set_title('Critical Temperatures for Processing', fontsize=12)
ax3.set_xticks(x_pos)
ax3.set_xticklabels(alloy_names)
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# Add temperature values on bars
for bar, val in zip(bars1, liquidus_vals):
    ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5, 
             f'{val:.0f}', ha='center', fontsize=9)
for bar, val in zip(bars2, solidus_vals):
    ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5, 
             f'{val:.0f}', ha='center', fontsize=9)

# Plot 4: Summary comparison table as text
ax4 = axes[1, 1]
ax4.axis('off')

# Create summary table
table_data = []
for name in alloy_names:
    comp = alloys[name]
    eta_120 = results[name]['eta_frac'][2] if len(results[name]['eta_frac']) > 2 else 0  # 120°C
    freezing = solidification[name]['liquidus'] - solidification[name]['solidus']
    table_data.append([
        name,
        f"{comp['Zn']*100:.1f}/{comp['Mg']*100:.1f}/{comp['Cu']*100:.1f}",
        f"{eta_120:.4f}",
        f"{freezing:.0f}°C",
        f"<{solidification[name]['solidus']:.0f}°C"
    ])

table = ax4.table(
    cellText=table_data,
    colLabels=['Alloy', 'Zn/Mg/Cu (wt%)', 'η @ 120°C', 'Freezing Range', 'Max Solutionize'],
    loc='center',
    cellLoc='center'
)
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1.2, 1.8)

# Color code cells
for i, name in enumerate(alloy_names):
    table[(i+1, 0)].set_facecolor(alloys[name]['color'])
    table[(i+1, 0)].set_text_props(color='white', fontweight='bold')

ax4.set_title('Summary Comparison\n(All data from COST507 database)', fontsize=12, pad=20)

plt.tight_layout()
plt.savefig('05_multi_alloy_comparison.png', dpi=300, bbox_inches='tight')
print("\nSaved: 05_multi_alloy_comparison.png")

# =============================================================================
# DETAILED RESULTS
# =============================================================================
print("\n" + "=" * 70)
print("MULTI-ALLOY COMPARISON RESULTS")
print("=" * 70)

print("\n| Alloy   | Composition (wt%)      | η@120°C | Liquidus | Solidus | Freezing |")
print("|---------|------------------------|---------|----------|---------|----------|")
for name in alloy_names:
    comp = alloys[name]
    eta_120 = results[name]['eta_frac'][2] if len(results[name]['eta_frac']) > 2 else 0
    liq = solidification[name]['liquidus']
    sol = solidification[name]['solidus']
    frz = liq - sol
    print(f"| {name} | Al-{comp['Zn']*100:.1f}Zn-{comp['Mg']*100:.1f}Mg-{comp['Cu']*100:.1f}Cu | {eta_120:.4f} | {liq:.0f}°C   | {sol:.0f}°C  | {frz:.0f}°C    |")

print("\n" + "=" * 70)
print("RECOMMENDATIONS:")
print("=" * 70)

# Find best for each criterion
eta_vals = {name: results[name]['eta_frac'][2] if len(results[name]['eta_frac']) > 2 else 0 
            for name in alloy_names}
best_strength = max(eta_vals, key=eta_vals.get)

freezing_ranges = {name: solidification[name]['liquidus'] - solidification[name]['solidus'] 
                   for name in alloy_names}
best_processing = min(freezing_ranges, key=freezing_ranges.get)

print(f"\n1. Best for STRENGTH (highest η-phase): {best_strength}")
print(f"   η-phase fraction at 120°C: {eta_vals[best_strength]:.4f}")

print(f"\n2. Best for PROCESSING (narrowest freezing range): {best_processing}")
print(f"   Freezing range: {freezing_ranges[best_processing]:.0f}°C")

print(f"\n3. Al-7075 remains the BALANCED choice for general aerospace use")

print("\n" + "=" * 70)
print("NOTE: All values computed from COST507 thermodynamic database")
print("      No experimental or placeholder data used")
print("=" * 70)

plt.show()