testing.py

import scipy.optimize as opt
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import os


class Sand():

    def __init__(self):
       
        self.temperature   = 0.0        # [C]
        self.density       = 0.0        # [kg/m3]
        self.enthalpy      = 0.0        # [J/kg]
        self.specific_heat = 0.0        # [J/kg-K]
        self.conductivity  = 0.0        # [W/m-K]
        self.molar_mass    = 0.0600843  # [kg/mol]
        self.bulk_conductivity = 0.27   # [W/m-K]
        self.packing_fraction  = 0.61   # [-]

    def __repr__(self):
        return (
            f"{'T':.<5}{self.temperature:.>15.3f} [C]\n"
            f"{'rho':.<5}{self.density:.>15.3f} [kg/m3]\n"
            f"{'h':.<5}{self.enthalpy:.>15.3f} [J/kg]\n"
            f"{'c':.<5}{self.specific_heat:.>15.3f} [J/kg-K]\n"
            f"{'M':.<5}{self.molar_mass:.>15.3f} [kg/kmol]\n"
            f"{'gamma':.<5}{self.packing_fraction:.>15.3f} [-]"
        )

    def update(self, temperature=None, enthalpy=None, bounds=(25, 1700)): 
        if isinstance(temperature, (float, int, np.ndarray)) and enthalpy is None: 
            self.temperature = temperature
            self._density()
            self._enthalpy()
            self._specific_heat()
            self._conductivity()
        elif isinstance(enthalpy, (float, int, np.ndarray)) and temperature is None: 
            self.enthalpy = enthalpy
            self._temperature(T_min=bounds[0], T_max=bounds[1])
            self._density()
            self._specific_heat()
            self._conductivity()
        else: raise ValueError('Must specify either temperature or enthalpy, but not both.')
        return self

    def _density(self):
        '''
        Calculates the packed bed density of SiO2. 

        Accepts T [C]
        Returns p [kg/m3]
        '''
        if isinstance(self.temperature, np.ndarray):
            density = np.empty_like(self.temperature, dtype=float)
            density[:] = np.nan  # Default value
           
            # Apply conditions using vectorized masking
            mask1 = self.temperature < 573
            mask2 = (self.temperature >= 573) & (self.temperature < 870)
            mask3 = (self.temperature >= 870) & (self.temperature < 1470)
            mask4 = (self.temperature >= 1470) & (self.temperature < 1705)

            density[mask1] = self.packing_fraction * 2648
            density[mask2] = self.packing_fraction * 2530
            density[mask3] = self.packing_fraction * 2250
            density[mask4] = self.packing_fraction * 2200
           
            self.density = density

        else: # if not vectorized
            if self.temperature < 573: 
                self.density = self.packing_fraction * 2648
            elif self.temperature < 870: 
                self.density = self.packing_fraction * 2530
            elif self.temperature < 1470: 
                self.density = self.packing_fraction * 2250
            elif self.temperature < 1705: 
                self.density = self.packing_fraction * 2200

        # if self.temperature < 573: 
        #     self.density = self.packing_fraction * 2648
        # elif self.temperature < 870: 
        #     self.density = self.packing_fraction * 2530
        # elif self.temperature < 1470: 
        #     self.density = self.packing_fraction * 2250
        # elif self.temperature < 1705: 
        #     self.density = self.packing_fraction * 2200

    def _specific_heat(self):
        """
        Calculates specific heat capacity (cp) for a given temperature T (in Kelvin).
        Uses different sets of coefficients depending on the temperature range.
       
        Accepts T [C]
        Returns c [J/kg-K]
        """
        # Temperature units correction
        self._temperature_K = self.temperature + 273.15
       
        # coefficients for different temperature ranges
        coeffs_lower = {'A': -6.076591, 'B': 251.6755, 'C': -324.7964, 'D': 168.5604, 'E': 0.002548}
        coeffs_upper = {'A':  58.75340, 'B': 10.27925, 'C': -0.131384, 'D': 0.025210, 'E': 0.025601}

        # Check if input is an array or scalar
        if isinstance(self._temperature_K, np.ndarray):
            # Initialize an empty array for specific heat
            self.specific_heat = np.empty_like(self._temperature_K, dtype=float)
           
            # Create masks for temperature ranges
            mask_lower = (298 <= self._temperature_K) & (self._temperature_K < 847)
            mask_upper = (847 <= self._temperature_K) & (self._temperature_K <= 1996)
            mask_invalid = ~ (mask_lower | mask_upper)
           
            if np.any(mask_invalid):
                raise ValueError("Temperature out of valid range (298 - 1996 K)")
           
            # Calculate specific heat for the lower range
            t_lower = self._temperature_K[mask_lower] / 1000
            self.specific_heat[mask_lower] = (1 / self.molar_mass) * (
                coeffs_lower['A'] +
                coeffs_lower['B'] * t_lower +
                coeffs_lower['C'] * t_lower**2 +
                coeffs_lower['D'] * t_lower**3 +
                coeffs_lower['E'] / t_lower**2
            )
           
            # Calculate specific heat for the upper range
            t_upper = self._temperature_K[mask_upper] / 1000
            self.specific_heat[mask_upper] = (1 / self.molar_mass) * (
                coeffs_upper['A'] +
                coeffs_upper['B'] * t_upper +
                coeffs_upper['C'] * t_upper**2 +
                coeffs_upper['D'] * t_upper**3 +
                coeffs_upper['E'] / t_upper**2
            )
        else: # if not vectorized
            if 298 <= self._temperature_K < 847:
                coeffs = coeffs_lower
            elif 847 <= self._temperature_K <= 1996:
                coeffs = coeffs_upper
            else: raise ValueError("Temperature out of valid range (298 - 1996 K)")
           
            t = self._temperature_K / 1000
            self.specific_heat = (1 / self.molar_mass) * (
                coeffs['A'] +
                coeffs['B'] * t +
                coeffs['C'] * t**2 +
                coeffs['D'] * t**3 +
                coeffs['E'] / t**2
            )

        # # Temperature units correction
        # self._temperature_K = self.temperature + 273.15
       
        # # Define coefficients for different temperature ranges
        # coeffs_lower = {'A': -6.076591, 'B': 251.6755, 'C': -324.7964, 'D': 168.5604, 'E': 0.002548}
        # coeffs_upper = {'A':  58.75340, 'B': 10.27925, 'C': -0.131384, 'D': 0.025210, 'E': 0.025601}

        # # Select coefficients based on temperature range
        # if 298 <= self._temperature_K < 847:
        #     coeffs = coeffs_lower
        # elif 847 <= self._temperature_K <= 1996:
        #     coeffs = coeffs_upper
        # else: raise ValueError("Temperature out of valid range (298 - 1996 K)")
       
        # # Compute scaled temperature for specific heat calculation
        # t = self._temperature_K / 1000
       
        # # Calculate specific heat capacity
        # self.specific_heat = (1 / self.molar_mass) * (coeffs['A'] + coeffs['B'] * t + coeffs['C'] * t**2 +
        #             coeffs['D'] * t**3 + coeffs['E'] / t**2)

    def _enthalpy(self):
        """
        Calculates enthalpy relative to 298.15 K for a given temperature T (in Kelvin).
        Uses different sets of coefficients depending on the temperature range.
       
        Accepts T [C]
        Returns h [J/kg]
        """
        self.enthalpy = self._get_enthalpy(self.temperature)

    def _get_enthalpy(self, input) -> float:
        """
        Calculates enthalpy relative to 298.15 K for a given temperature T (in Kelvin).
        Uses different sets of coefficients depending on the temperature range.
       
        Accepts T [C]
        Returns h [J/kg]
        """
        # Temperature units correction
        self._temperature_K = input + 273.15
       
        # Define coefficients for different temperature ranges
        coeffs_lower = {'A': -6.076591, 'B': 251.6755, 'C': -324.7964, 'D': 168.5604, 'E': 0.002548, 'F': -917.6893, 'H': -910.8568}
        coeffs_upper = {'A': 58.75340, 'B': 10.27925, 'C': -0.131384, 'D': 0.025210, 'E': 0.025601, 'F': -929.3292, 'H': -910.8568}

        # Check if input is an array or scalar
        if isinstance(self._temperature_K, np.ndarray):
            # Initialize an empty array for enthalpy
            enthalpy = np.empty_like(self._temperature_K, dtype=float)
           
            # Create masks for temperature ranges
            mask_lower = (298 <= self._temperature_K) & (self._temperature_K < 847)
            mask_upper = (847 <= self._temperature_K) & (self._temperature_K <= 1996)
            mask_invalid = ~ (mask_lower | mask_upper)
           
            # Raise an error if any temperatures are out of the valid range
            if np.any(mask_invalid):
                raise ValueError("Temperature out of valid range (298 - 1996 K)")
           
            # Calculate enthalpy for the lower range
            t_lower = self._temperature_K[mask_lower] / 1000
            enthalpy[mask_lower] = ((1 / self.molar_mass) * (
                coeffs_lower['A'] * t_lower +
                coeffs_lower['B'] * t_lower**2 / 2 +
                coeffs_lower['C'] * t_lower**3 / 3 +
                coeffs_lower['D'] * t_lower**4 / 4 +
                coeffs_lower['E'] / t_lower +
                coeffs_lower['F'] - 
                coeffs_lower['H']
            )) * 1000
           
            # Calculate enthalpy for the upper range
            t_upper = self._temperature_K[mask_upper] / 1000
            enthalpy[mask_upper] = ((1 / self.molar_mass) * (
                coeffs_upper['A'] * t_upper +
                coeffs_upper['B'] * t_upper**2 / 2 +
                coeffs_upper['C'] * t_upper**3 / 3 +
                coeffs_upper['D'] * t_upper**4 / 4 +
                coeffs_upper['E'] / t_upper +
                coeffs_upper['F'] - 
                coeffs_upper['H']
            )) * 1000
           
            self.enthalpy = enthalpy
        else: # if not vectorized
            if 298 <= self._temperature_K < 847:
                coeffs = coeffs_lower
            elif 847 <= self._temperature_K <= 1996:
                coeffs = coeffs_upper
            else: raise ValueError(f"Temperature ({self._temperature_K:.2f} K) is out of valid range (298 - 1996 K)")
           
            t = self._temperature_K / 1000
            enthalpy = ((1 / self.molar_mass) * (
                coeffs['A'] * t +
                coeffs['B'] * t**2 / 2 +
                coeffs['C'] * t**3 / 3 +
                coeffs['D'] * t**4 / 4 +
                coeffs['E'] / t +
                coeffs['F'] - 
                coeffs['H']
            )) * 1000

        return enthalpy

        # # Temperature units correction
        # self._temperature_K = input + 273.15
       
        # # Define coefficients for different temperature ranges
        # coeffs_lower = {'A': -6.076591, 'B': 251.6755, 'C': -324.7964, 'D': 168.5604, 'E': 0.002548, 'F': -917.6893, 'H': -910.8568}
        # coeffs_upper = {'A': 58.75340, 'B': 10.27925, 'C': -0.131384, 'D': 0.025210, 'E': 0.025601, 'F': -929.3292, 'H': -910.8568}
       
        # # Select coefficients based on temperature range
        # if 298 <= self._temperature_K < 847:
        #     coeffs = coeffs_lower
        # elif 847 <= self._temperature_K <= 1996:
        #     coeffs = coeffs_upper
        # else: raise ValueError("Temperature out of valid range (298 - 1996 K)")
       
        # # Compute scaled temperature for specific heat calculation
        # t = self._temperature_K / 1000
       
        # # Calculate enthalpy difference
        # enthalpy = ((1 / self.molar_mass) * (coeffs['A'] * t + coeffs['B'] * t**2 / 2 + coeffs['C'] * t**3 / 3 +
        #                 coeffs['D'] * t**4 / 4 + coeffs['E'] / t + coeffs['F'] - coeffs['H'])) * 1000
        # return enthalpy

    def _temperature(self, T_min=25, T_max=1700):
        """
        Solves for temperature given a target enthalpy value using Brent's method.
       
        Accepts h [J/kg]
        Returns T [C]
        """
       
        # Check if enthalpy is an array
        if isinstance(self.enthalpy, np.ndarray):
            # Initialize an array to store temperatures
            temperatures = np.empty_like(self.enthalpy, dtype=float)
           
            # Loop through each target enthalpy value
            for i, h_target in enumerate(self.enthalpy):
                if self._get_enthalpy(T_min) > h_target or self._get_enthalpy(T_max) < h_target:
                    raise ValueError("Target enthalpy is outside the valid range of temperatures.")
                temperatures[i] = opt.brentq(lambda T: self._get_enthalpy(T) - h_target, T_min, T_max)            
            self.temperature = temperatures
           
        else: # if not vectorized
            if self._get_enthalpy(T_min) > self.enthalpy or self._get_enthalpy(T_max) < self.enthalpy:
                raise ValueError("Target enthalpy is outside the valid range of temperatures.")
           
            # Solve for temperature using Brent's method
            self.temperature = opt.brentq(lambda T: self._get_enthalpy(T) - self.enthalpy, T_min, T_max)
        return self.temperature

        # # Ensure enthalpy function is well-behaved in the given range
        # if self._get_enthalpy(T_min) > self.enthalpy or self._get_enthalpy(T_max) < self.enthalpy:
        #     raise ValueError("Target enthalpy is outside the valid range of temperatures.")
       
        # # Solve for temperature using brentq root-finding
        # self.temperature = opt.brentq(lambda T: self._get_enthalpy(T) - self.enthalpy, T_min, T_max)

    def _conductivity(self): 
        """
        Solves for SiO2 conductivity, given a temperature.
        """

        self._temperature_K = self.temperature + 273.15
        # getting coefficients based on temperature
        A = np.where(self._temperature_K < 597, 0.00144,
                    np.where((self._temperature_K >= 597) & (self._temperature_K < 800), 0.00209,
                            np.where((self._temperature_K >= 800) & (self._temperature_K < 1002), 0.00398, 0.0)))

        B = np.where(self._temperature_K < 597, 0.96928, 
                    np.where((self._temperature_K >= 597) & (self._temperature_K < 800), 0.49472,
                            np.where((self._temperature_K >= 800) & (self._temperature_K < 1002), 0.49472, 2.87)))

        self.conductivity = A * self._temperature_K + B

def shomate(T, A, B, C, D, E):
    t = T / 1000
    return A + B * t + C * t**2 + D * t**3 + E/t**2

def intshomate(T, A, B, C, D, E):
    t = T / 1000
    return A * t + B * t**2 / 2 + C * t**3 / 3 + D * t**4 / 4 - E / t

def sandfit(display=True):

    tspan = np.linspace(600, 1200 + 273.15)
   
    sand  = Sand()
    sand.update(
        temperature = tspan - 273.15
    )
    sifit = sand.specific_heat / 1000

    if display:
        xs_bauxite = [500.0, 700.0, 900.0, 1100., 1300.]
        ys_bauxite = np.array([56.97, 59.21, 60.34, 60.96, 61.34]) / 42.9809
        plt.plot(tspan, sifit, label='sand')
        plt.scatter(xs_bauxite, ys_bauxite, label='AlO2', zorder=3)

def datafit(T, display=True, dataset=1):
    file = f'HSP4070_{dataset}.csv'
    path = os.path.join(os.getcwd(), 'Figures and Data', 'data', file)

    df = pd.read_csv(path)
    nm = df.columns
    df = df[nm]

    temps = df[nm[0]] + 273.15
    scaps = df[nm[1]]

    popt, pcov = curve_fit(shomate, temps, scaps)
    A, B, C, D, E = popt

    print(f"Specific Heat (set {dataset})")
    print("A + B*T + C*T**2 + D*T**3 + E/T**2")
    print(f"{A:>10.5f}")
    print(f"{B:>10.5f}")
    print(f"{C:>10.5f}")
    print(f"{D:>10.5f}")
    print(f"{E:>10.5f}")
    print()

    tspan = T
    scfit = shomate(tspan, A, B, C, D, E)

    if display:
        # plt.plot(temps, scaps, label=f'data set {dataset}')
        plt.plot(tspan, scfit, label=f'HSP 40/70 data')

    return A, B, C, D, E

def tempfit():

    sandfit(False)
    A, B, C, D, E = datafit(False, 1)
    A, B, C, D, E = datafit(False, 2)

    Tmin = 300 + 273.15
    Tmax = 1100 + 273.15
    temps = np.linspace(Tmin, Tmax, 50)
    enths = intshomate(temps, A, B, C, D, E)

    popt, pcov = curve_fit(shomate, enths * 1000, temps)
    A, B, C, D, E = popt

    tsfit = shomate(enths * 1000, A, B, C, D, E)

    print(f"Temperature (set {2})")
    print("A + B*T + C*T**2 + D*T**3 + E/T**2")
    print(f"{A:>10.5f}")
    print(f"{B:>10.5f}")
    print(f"{C:>10.5f}")
    print(f"{D:>10.5f}")
    print(f"{E:>10.5f}")
    print()

    plt.scatter(enths, temps, label='HSP 40/70', zorder=3, s=10)
    plt.plot(enths, tsfit, label='fit', color='orange')
    plt.legend()
    plt.ylabel('Temperature [K]')
    plt.xlabel('Enthalpy [kJ/kg]')
    plt.grid(True)
    plt.margins(x=0)
    plt.show()



def cpAl2O3(T): 
    M = 101.9613    # [g/mol]
    A = 102.4290
    B = 38.74980
    C = -15.9109
    D = 2.628181
    E = -3.00755
    cp_mol = shomate(T, A, B, C, D, E) # [J/mol-K]
    return cp_mol / M

def cpTiO2(T): 
    return 0.683 # [J/g-K]

def cpFeO(T): 
    M = 71.844
    A = 45.75120
    B = 18.78553
    C = -5.95220
    D = 0.852779
    E = -0.08127
    cp_mol = shomate(T, A, B, C, D, E)
    return cp_mol / M


def cpHSP4070(T, a=0.486, b=0.2054, c=0, d=0.3086): 

    Al_oxide = a * cpAl2O3(T)
    Ti_oxide = b * cpTiO2(T)
    Fe_oxide = c * cpFeO(T)
    Si_oxide = d

    sand = Sand()
    sand.update(
        temperature = T - 273.15
    )

    cp = Al_oxide + Ti_oxide + Fe_oxide + (Si_oxide * sand.specific_heat / 1000)
    return cp


def opt_chem_formula():
   
    def routine(params): 

        a, b, c, d = params
        T = np.linspace(600, 1400, 100)
        cp1 = cpHSP4070(T, a, b, c, d)
        A, B, C, D, E = datafit(T, display=False, dataset=2)
        cp2 = shomate(T, A, B, C, D, E)
        return np.sum((cp1 - cp2)**2)

    x0 = [0.70, 0.05, 0.10, 0.15]
    lincon = opt.LinearConstraint([1, 1, 1, 1], [1], [1])
    bounds = [(0, 1), (0, 1), (0, 1), (0, 1)]

    result = opt.minimize(routine, x0, bounds=bounds, method='SLSQP', constraints=[lincon])
    print(result.x)

if __name__ == '__main__':

    Ts = np.linspace(600, 1475, 100)
    cp = cpHSP4070(Ts)
    
    plt.plot(Ts, cp, label='HSP 40/70 estimate')
    plt.xlabel('Temperature [K]')
    plt.ylabel('Specific Heat [J/g-K]')
    plt.margins(x=0)
    plt.grid()
    
    datafit(Ts, display=True, dataset=2)
    
    plt.legend()
    plt.show()

    opt_chem_formula()

#EOF