512 python

# main.py

import numpy as np
from RHjump import calculate_rh_conditions
from species import calculate_thermodynamic_properties
from mixture import compute_cp, get_mixture_properties  # Importing necessary functions

def initialize_composition():
    """
    Initialize the gas mixture composition (N2/N Mixture).
    """
    # Example composition: 100% N2, 0% N
    return 1.0  # Return y (mole fraction of N2)

def simulate_shock_flow(x_max, post_shock_conditions):
    """
    Simulate the flow downstream of the shock.
    """
    # Placeholder for flow field simulation
    return {
        "x": np.linspace(0, x_max, 100),
        "rho_s": None,
        "T": None,
        "Tv": None,
        "U": None
    }

def main():
    """
    The main program computes the flowfield in a shock tube.
    """
    # Initialize parameters
    P1 = 5.0        # Pre-shock pressure [Pa]
    T1 = 200.0      # Pre-shock temperature [K]
    TV1 = T1        # Pre-shock vibrational temperature [K]
    U1 = 7000.0     # Shock speed [m/s]
    x_max = 0.1     # Standoff distance [m]

    # Initialize composition
    y = initialize_composition()  # Mole fraction of N2

    # Load species data
    species_data = np.loadtxt("data2.txt", skiprows=1)

    print("Shock-Tube Simulation")
    print("========================================")
    print("Mixture: N2/N")

    # Calculate properties for N (species_idx=0) at T=7000 K
    properties_N = calculate_thermodynamic_properties(
        temperature=7000,
        species_idx=0,
        species_data=species_data
    )

    # Access and display the results
    print(f"Q_in: {properties_N['Q_in']}")
    print(f"Q_tr: {properties_N['Q_tr']}")
    print(f"Q_total: {properties_N['Q_total']}")
    print(f"Entropy: {properties_N['S']} J/mol/K")
    print(f"H-H0: {properties_N['H_H0']} J/mol")

    # Compute Cp for each species at T = 7000 K
    Cp_N = compute_cp(T=7000, species_idx=0, species_data=species_data)
    Cp_N2 = compute_cp(T=7000, species_idx=1, species_data=species_data)

    print(f"Cp (N): {Cp_N} J/mol/K")
    print(f"Cp (N2): {Cp_N2} J/mol/K")

    # Compute mixture properties including Cp
    mixture_props = get_mixture_properties(T=7000, y=y, species_data=species_data)
    print(f"Mixture Cp: {mixture_props['Cp']} J/mol/K")

    print("========================================")

    # Compute initial conditions (Post-Shock)
    print("Computing Rankine-Hugoniot jump relations for hot gas.")
    post_shock_conditions = calculate_rh_conditions(P1, T1, U1, y, species_data)
    print("Post-Shock Conditions (Hot Gas):")
    print(f"P2: {post_shock_conditions['P2']} Pa")
    print(f"T2: {post_shock_conditions['T2']} K")
    print(f"rho2: {post_shock_conditions['rho2']} kg/m³")
    print(f"u2: {post_shock_conditions['u2']} m/s")
    print(f"a2: {post_shock_conditions['a2']} m/s")
    print(f"M2: {post_shock_conditions['M2']}")
    print(f"gamma2: {post_shock_conditions['gamma2']}")
    print("========================================")

    # Simulate the flowfield downstream of the shock
    flow_field = simulate_shock_flow(x_max, post_shock_conditions)
    print("Shock-tube simulation complete.")

if __name__ == "__main__":
    main()



# RHjump.py

import numpy as np
from mixture import get_mixture_properties, get_sound_speed, R_UNIVERSAL

def calculate_rh_conditions(P1: float, T1: float, u1: float, y: float, species_data: np.ndarray) -> dict:
    """
    Calculate Rankine-Hugoniot jump conditions.

    Args:
        P1 (float): Initial pressure (Pa)
        T1 (float): Initial temperature (K)
        u1 (float): Initial velocity (m/s)
        y (float): Mole fraction of N2
        species_data (np.ndarray): Array containing species data

    Returns:
        dict: Dictionary containing post-shock conditions
              Keys: 'P2', 'T2', 'rho2', 'u2', 'a2', 'M2', 'gamma2'
    """
    # Step 1: Get initial mixture properties
    props1 = get_mixture_properties(T1, y, species_data)
    gamma1 = props1['gamma']
    M_mix = props1['M']

    # Step 2: Calculate specific gas constant (R_specific)
    # R_UNIVERSAL is in J/mol/K, M_mix should be in kg/mol
    R_specific = R_UNIVERSAL / M_mix  # J/kg/K

    # Step 3: Calculate initial density (rho1)
    rho1 = P1 / (R_specific * T1)  # kg/m³

    # Step 4: Calculate sound speed (a1)
    a1 = get_sound_speed(T1, y, species_data)  # m/s

    # Step 5: Calculate Mach number (M1)
    M1 = u1 / a1

    # Step 6: Apply Rankine-Hugoniot jump conditions
    P2 = P1 * (2 * gamma1 * M1**2 - (gamma1 - 1)) / (gamma1 + 1)  # Post-shock pressure (Pa)
    rho2 = rho1 * ((gamma1 + 1) * M1**2) / ((gamma1 - 1) * M1**2 + 2)  # Post-shock density (kg/m³)
    u2 = u1 * (((gamma1 - 1) * M1**2 + 2) / ((gamma1 + 1) * M1**2))  # Post-shock velocity (m/s)

    # Step 7: Calculate temperature ratio and absolute temperature (T2)
    T2_T1 = (P2 / P1) * (rho1 / rho2)  # Temperature ratio T2/T1
    T2 = T1 * T2_T1  # Post-shock temperature (K)

    # Step 8: Assume frozen flow; sound speed remains the same
    a2 = a1  # Post-shock sound speed (m/s)

    # Step 9: Calculate post-shock Mach number (M2)
    M2 = u2 / a1  # Assuming a2 = a1

    # Step 10: Assume gamma remains the same across the shock
    gamma2 = gamma1

    # Step 11: Compile results into a dictionary
    return {
        'P2': P2,          # Post-shock pressure (Pa)
        'T2': T2,          # Post-shock temperature (K)
        'rho2': rho2,      # Post-shock density (kg/m³)
        'u2': u2,          # Post-shock velocity (m/s)
        'a2': a2,          # Post-shock sound speed (m/s)
        'M2': M2,          # Post-shock Mach number
        'gamma2': gamma2   # Post-shock specific heat ratio
    }



import numpy as np
import math
from dataclasses import dataclass

# Physical constants
R = 0.25  # O2 N2 ratio
PRESSURE = 1  # pressure
KB = 1.380649E-23  # Boltzmann constant
H = 6.62607015E-34  # Planck constant
C = 3.0E+10  # speed of light, cm/s
DT = 1.0E-8  # T step for finite differences
NA = 6.02214076E+23  # Avogadro number
NU = 1.0  # number of moles
P_ATM = 101325.0  # atmospheric pressure
EV = 1.60218E-19  # eV -> J translation coefficient


@dataclass
class Species:
    constants: np.ndarray

def calculate_thermodynamic_properties(temperature: float, species_idx: int, species_data: np.ndarray) -> dict:
    """
    Calculate thermodynamic properties for a given species at specified temperature.

    Args:
        temperature: Temperature in Kelvin
        species_idx: Species index (0 for N, 1 for N2)
        species_data: Array containing species data

    Returns:
        Dictionary containing:
        - Q_in: Internal partition function
        - Q_tr: Translational partition function
        - Q_total: Total partition function
        - S: Entropy (J/mol/K)
        - H_H0: Enthalpy difference from 0K (kJ/mol)
    """
    # Create Species object for the selected species
    species = Species(constants=species_data[species_idx])

    # Calculate Q_in
    if species.constants[1] == 0:
        Q_rot = 1.0
    else:
        Q_rot = (8.0 * math.pi**2 * species.constants[1] * KB * temperature) / (H**2 * species.constants[0])

    if species.constants[2] == 0:
        Q_vib = 1.0
    else:
        Q_vib = 1.0 / (1.0 - math.exp((-H * C * species.constants[2]) / (KB * temperature)))

    if species.constants[3] == 0:
        Q_el = 1.0
    else:
        Q_el = sum(
            species.constants[3+i] * math.exp(-H * C * species.constants[22+i] / (KB * temperature))
            for i in range(19)
        )

    Q_in = Q_rot * Q_vib * Q_el

    # Calculate Q_tr
    Q_tr = (NU * KB * NA * temperature / (PRESSURE * P_ATM)) * \
           (2 * math.pi * species.constants[41] * KB * temperature / (H * H))**(1.5)

    # Calculate Q_total
    Q_total = Q_in * Q_tr

    # Calculate Entropy [J/(mol*K)]
    S = KB * NA * (
        math.log(Q_in) + math.log(Q_tr)
    ) + KB * NA * temperature * (
        (Q_tr_calc(species, temperature + DT) - Q_tr_calc(species, temperature - DT)) /
        (2.0 * DT * Q_tr) +
        (Q_in_calc(species, temperature + DT) - Q_in_calc(species, temperature - DT)) /
        (2.0 * DT * Q_in)
    ) - KB * NA * math.log(NA)

    # Calculate H-H0 [J/mol]
    H_H0 = KB * NA * temperature * (
        (temperature / Q_in) *
        (Q_in_calc(species, temperature + DT) - Q_in_calc(species, temperature - DT)) /
        (2.0 * DT) + 2.5
    )

    return {
        'Q_in': Q_in,
        'Q_tr': Q_tr,
        'Q_total': Q_total,
        'S': S,
        'H_H0': H_H0
    }

# Helper functions for finite difference calculations
def Q_in_calc(species: Species, T: float) -> float:
    """Calculate Q_in at a specific temperature (helper function for derivatives)"""
    if species.constants[1] == 0:
        Q_rot = 1.0
    else:
        Q_rot = (8.0 * math.pi**2 * species.constants[1] * KB * T) / (H**2 * species.constants[0])

    if species.constants[2] == 0:
        Q_vib = 1.0
    else:
        Q_vib = 1.0 / (1.0 - math.exp((-H * C * species.constants[2]) / (KB * T)))

    if species.constants[3] == 0:
        Q_el = 1.0
    else:
        Q_el = sum(
            species.constants[3+i] * math.exp(-H * C * species.constants[22+i] / (KB * T))
            for i in range(19)
        )

    return Q_rot * Q_vib * Q_el

def Q_tr_calc(species: Species, T: float) -> float:
    """Calculate Q_tr at a specific temperature (helper function for derivatives)"""
    return (NU * KB * NA * T / (PRESSURE * P_ATM)) * \
           (2 * math.pi * species.constants[41] * KB * T / (H * H))**(1.5)


# mixture.py

import numpy as np
from species import calculate_thermodynamic_properties

# Constants
R_UNIVERSAL = 8.314462  # J/mol/K
Na = 6.02214076E23 # Avogadro
M_N = 0.014             # kg/mol
M_N2 = 0.028            # kg/mol
EV_TO_J = 1.60218E-19  # eV to J conversion
H_F_N = 4.87933027 * EV_TO_J * Na  # Convert eV to J/mol
H_F_N2 = 0.0            # Formation enthalpy for N2 in J/mol


def compute_cp(T: float, species_idx: int, species_data: np.ndarray, delta_T: float = 1e-8) -> float:
    """
    Compute the specific heat at constant pressure (Cp) for a given species using the provided formula.

    Args:
        T: Temperature in Kelvin
        species_idx: Species index (0 for N, 1 for N2)
        species_data: Array containing species data
        delta_T: Temperature increment for finite differences (default: 1e-3 K)

    Returns:
        Cp in J/mol/K
    """
    # Calculate Q_in at T
    props_T = calculate_thermodynamic_properties(T, species_idx, species_data)
    Q_in_T = props_T['Q_in']

    # Calculate Q_in at T + delta_T
    props_T_plus = calculate_thermodynamic_properties(T + delta_T, species_idx, species_data)
    Q_in_T_plus = props_T_plus['Q_in']

    # Calculate Q_in at T - delta_T
    props_T_minus = calculate_thermodynamic_properties(T - delta_T, species_idx, species_data)
    Q_in_T_minus = props_T_minus['Q_in']

    # Compute ln(Q_in) at T, T + delta_T, T - delta_T
    ln_Q_in_T = np.log(Q_in_T)
    ln_Q_in_T_plus = np.log(Q_in_T_plus)
    ln_Q_in_T_minus = np.log(Q_in_T_minus)

    # First derivative d(lnQ_in)/dT using central difference
    dlnQin_dT = (ln_Q_in_T_plus - ln_Q_in_T_minus) / (2 * delta_T)

    # Second derivative d²(lnQ_in)/dT² using central difference
    d2lnQin_dT2 = (ln_Q_in_T_plus - 2 * ln_Q_in_T + ln_Q_in_T_minus) / (delta_T ** 2)

    # Calculate Cp_int using Capitelli
    Cp_int = R_UNIVERSAL * (2 * T * dlnQin_dT + T**2 * d2lnQin_dT2)

    # Total Cp
    Cp = Cp_int + (5.0 / 2.0) * R_UNIVERSAL

    return Cp

def get_mixture_properties(T: float, y: float, species_data: np.ndarray) -> dict:
    """
    Calculate mixture properties including formation enthalpies and Cp.

    Args:
        T: Temperature in K
        y: Mole fraction of N2 (1-y is mole fraction of N)
        species_data: Array containing thermodynamic data for species

    Returns:
        Dictionary containing gamma, molar mass, sound speed coefficient, mixture enthalpy, and Cp
    """
    # Get thermodynamic properties for both species
    N_props = calculate_thermodynamic_properties(T, 0, species_data)
    N2_props = calculate_thermodynamic_properties(T, 1, species_data)

    # Calculate mixture enthalpy including formation enthalpies
    H_mix = (1 - y) * (N_props['H_H0'] + H_F_N) + y * (N2_props['H_H0'] + H_F_N2)

    # Compute Cp for each species using the new compute_cp function
    Cp_int_N = compute_cp(T, 0, species_data)
    Cp_int_N2 = compute_cp(T, 1, species_data)

    # Calculate mixture Cp using mole fractions
    Cp_mix = (1 - y) * Cp_int_N + y * Cp_int_N2

    # Calculate gamma using Cp and Cv = Cp - R
    Cv_mix = Cp_mix - R_UNIVERSAL
    gamma = Cp_mix / Cv_mix

    # Calculate mixture molar mass
    M_mix = (1 - y) * M_N + y * M_N2

    # Calculate sound speed coefficient
    sound_coeff = np.sqrt(gamma * R_UNIVERSAL / M_mix)

    return {
        'gamma': gamma,
        'M': M_mix,
        'sound_coeff': sound_coeff,
        'H_mix': H_mix,
        'Cp': Cp_mix
    }

def get_sound_speed(T: float, y: float, species_data: np.ndarray) -> float:
    """
    Calculate mixture sound speed.

    Args:
        T: Temperature in K
        y: Mole fraction of N2
        species_data: Array containing thermodynamic data for species

    Returns:
        Sound speed in m/s
    """
    props = get_mixture_properties(T, y, species_data)
    return props['sound_coeff'] * np.sqrt(T)