Acoustic Transmission at Low Gas Pressures

acoustic_transmission_in_vac.py

import numpy as np
import matplotlib.pyplot as plt

def plot_acoustic_transfer_function_n2():
    """
    Calculates and plots the Acoustic Transfer Function (P_recv / P_source) 
    for Nitrogen gas across Continuum, Transition, and Free-Molecular regimes.
    
    Pedagogy & Physics:
    -------------------
    1. GENERATION (Impedance Mismatch):
       The first barrier is getting energy from the wall (Steel) into the gas (N2).
       Transmission Coeff (Pressure Amplitude) T_12 = 2*Z2 / (Z1 + Z2).
       Since Z_gas << Z_wall, T_12 ~ 2 * Z_gas / Z_wall.
       Because Z_gas ~ rho ~ Pressure, this coupling scales linearly with P.
       Ref: Kinsler, Frey, Coppens, & Sanders, "Fundamentals of Acoustics", Ch 6.

    2. PROPAGATION (Attenuation):
       a. Continuum Regime (Kn << 1): 
          Stokes-Kirchhoff attenuation alpha ~ f^2 / P.
          Includes bulk viscosity correction for N2 rotational relaxation.
       b. Transition/Free-Molecular Regime (Kn > 1):
          The gas ceases to be a fluid. Attenuation per wavelength saturates 
          as the wave becomes a ballistic molecular beam.
          Model uses Greenspan's scaling parameter (f/P) to bridge regimes.
       Ref: Greenspan, M. (1956). J. Acoust. Soc. Am. 28, 644.
       Ref: Meyer, E., & Sessler, G. (1957). Z. Phys. 149, 15.

    Target Audience: LIGO / Precision Measurement Physicists.
    """

    # --- Constants for Nitrogen (N2) at 293 K ---
    T = 293.15          
    M = 28.0134e-3      # kg/mol
    R_gas = 8.31446     
    R_specific = R_gas / M
    gamma = 1.4         
    eta = 1.78e-5       # Pa s
    Pr = 0.72           
    
    # Adiabatic Sound Speed (Continuum limit)
    c0 = np.sqrt(gamma * R_specific * T) 
    
    # Impedance of the Source Wall (Stainless Steel)
    # Z_wall = rho_steel * c_steel ~ 8000 * 5700 ~ 45 MRayls
    Z_WALL = 45e6 

    # --- Derived Viscosity Parameters ---
    # N2 has rotational degrees of freedom (diatomic). 
    # Bulk viscosity eta_B accounts for energy lost to internal modes.
    # For N2, eta_B / eta ~ 0.8.
    eta_bulk = 0.8 * eta
    
    # Effective dissipation parameter (Stokes + Kirchhoff + Bulk)
    # This combines shear viscosity, thermal conductivity, and bulk viscosity losses.
    dissipation_eff = (4/3)*eta + eta_bulk + (gamma - 1)*(eta/Pr)

    # --- Kinetic Theory Limits ---
    # As Kn -> infinity (UHV), attenuation per wavelength saturates.
    # Experimental limit for diatomic gases is approx 2.5 Nepers/lambda.
    LIMIT_ALPHA_LAMBDA = 2.5 

    def get_transfer_function(f, P, dist=1.0):
        """
        Returns P_received / P_wall_displacement_stress
        Combines Coupling Loss (Impedance) * Propagation Loss (Attenuation).
        """
        
        # --- 1. Impedance Matching (The "Interface" Loss) ---
        rho = P / (R_specific * T)
        Z_gas = rho * c0
        
        # Pressure Transmission Coefficient (Normal Incidence)
        # T_p = P_trans / P_inc = 2 * Z_gas / (Z_gas + Z_wall)
        # Approximation valid since Z_gas (400 Rayls @ 1atm) << Z_wall (45 MRayl)
        coupling_efficiency = 2 * Z_gas / Z_WALL
        
        # --- 2. Attenuation (The "Path" Loss) ---
        omega = 2 * np.pi * f
        
        # A. Classical Stokes-Kirchhoff (Valid for f/P < 10^4 Hz/Pa)
        # alpha_cl = (omega^2 / (2 rho c^3)) * dissipation_eff
        alpha_cl_np = (omega**2 / (2 * rho * c0**3)) * dissipation_eff
        
        # B. Greenspan/Meyer Bridge to Free-Molecular
        # Normalize to attenuation per wavelength: alpha * lambda
        alpha_lambda_cl = alpha_cl_np * (c0 / f)
        
        # Apply saturation function: y = x / (1 + x/limit)
        # This smoothly transitions from x (linear) to limit (constant)
        alpha_lambda_unified = alpha_lambda_cl / (1 + alpha_lambda_cl / LIMIT_ALPHA_LAMBDA)
        
        # Convert back to spatial attenuation [Nepers/m]
        alpha_unified_np = alpha_lambda_unified * (f / c0)
        
        propagation_loss = np.exp(-alpha_unified_np * dist)
        
        return coupling_efficiency * propagation_loss

    # --- Plotting Routine ---
    fig, ax = plt.subplots(figsize=(10, 7))
    
    # Pressures of interest for UHV/Vacuum systems
    pressures = [101325, 10, 1, 1e-5] 
    labels = [
        "1 Bar (Laboratory Air)", 
        "10 Pa (Roughing Line)", 
        "1 Pa (Turbomolecular Transition)", 
        "1e-5 Pa (LIGO Beam Tube / UHV)"
    ]
    
    # Colors for visibility against dark/light backgrounds
    colors = ['#2C3E50', '#00FFFF', '#FF00FF', '#FF4500'] # Dark Blue, Cyan, Magenta, OrangeRed

    # Frequency range: 10 Hz to 1 MHz
    f = np.logspace(1, 6, 500) 

    for P, label, color in zip(pressures, labels, colors):
        tf = get_transfer_function(f, P)
        ax.loglog(f, tf, label=label, color=color, linewidth=2.5)

    # --- Plot Styling ---
    ax.set_title("Acoustic Transfer Function in $N_2$ (Path $L=1.0$ m)\nCombining Impedance Mismatch & Rarefied Gas Attenuation", 
                 fontsize=14, fontweight='bold', pad=15)
    
    ax.set_xlabel(r'Frequency $f$ [Hz]', fontsize=12)
    ax.set_ylabel(r'Pressure Transfer Function $|P_{rx} / P_{wall}|$', fontsize=12)
    
    # Grid: Major (solid) and Minor (dotted) for reading log scales
    ax.grid(True, which="major", linestyle='-', alpha=0.5, color='gray')
    ax.grid(True, which="minor", linestyle=':', alpha=0.3, color='gray')

    # Legend
    ax.legend(loc='upper right', fontsize=10, framealpha=0.9, edgecolor='black')
    
    # Y-Limits: Show from unity down to the "thermal noise floor" of reality
    ax.set_ylim(1e-18, 1e-1) 
    
    # --- "Pedagogical" Code Commentary (Mental Noise Budget) ---
    # 1 Bar:  Coupling is ~1e-5 (-100dB). You shout, the wall barely moves. Standard acoustics.
    # 10 Pa:  Coupling drops by 10^4. Viscosity kills the ultrasonic band. 
    #         If you are hearing pings here, it's structure-borne, not gas-borne.
    # 1 Pa:   The "Cutoff" (Kn ~ 1) is audible. High freq is just heat.
    # UHV:    Transfer function is 1e-16. That's -320 dB. 
    #         The mean free path is kilometers. 
    #         If you see signal here, check your grounding or blame the backscattered light.
    
    plt.tight_layout() # tight!
    plt.show()

if __name__ == "__main__":
    plot_acoustic_transfer_function_n2()