Chapter 1: Boltzmann Equation and H-Theorem

import numpy as np import matplotlib.pyplot as plt from scipy.stats import maxwell class MaxwellDistribution: """Analysis and visualization of Maxwell velocity distribution""" def __init__(self, temperature, mass): """ Parameters: temperature: Temperature [K] mass: Particle mass [kg] """ self.T = temperature self.m = mass self.kB = 1.380649e-23 # Boltzmann constant # Characteristic velocities self.v_mean = np.sqrt(8 * self.kB * self.T / (np.pi * self.m)) self.v_rms = np.sqrt(3 * self.kB * self.T / self.m) self.v_prob = np.sqrt(2 * self.kB * self.T / self.m) def distribution(self, v): """Maxwell velocity distribution function""" a = np.sqrt(self.kB * self.T / self.m) return maxwell.pdf(v, scale=a) def plot_distribution(self): """Plot distribution function""" v = np.linspace(0, 4*self.v_rms, 1000) f_v = self.distribution(v) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # Velocity distribution ax1.plot(v, f_v, 'b-', linewidth=2, label='Maxwell distribution') ax1.axvline(self.v_prob, color='r', linestyle='--', label=f'Most probable velocity = {self.v_prob:.1f} m/s') ax1.axvline(self.v_mean, color='g', linestyle='--', label=f'Mean velocity = {self.v_mean:.1f} m/s') ax1.axvline(self.v_rms, color='purple', linestyle='--', label=f'RMS velocity = {self.v_rms:.1f} m/s') ax1.set_xlabel('Velocity v [m/s]', fontsize=12) ax1.set_ylabel('Probability density f(v)', fontsize=12) ax1.set_title(f'Maxwell Velocity Distribution (T={self.T}K)', fontsize=14, fontweight='bold') ax1.legend() ax1.grid(True, alpha=0.3) # Temperature dependence temperatures = [100, 300, 500, 1000] for T in temperatures: md = MaxwellDistribution(T, self.m) v = np.linspace(0, 3000, 1000) ax2.plot(v, md.distribution(v), label=f'T={T}K') ax2.set_xlabel('Velocity v [m/s]', fontsize=12) ax2.set_ylabel('Probability density f(v)', fontsize=12) ax2.set_title('Temperature Dependence', fontsize=14, fontweight='bold') ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() return fig # Execution example: Nitrogen molecule (room temperature) m_N2 = 28.014 * 1.66054e-27 # kg maxwell_N2 = MaxwellDistribution(temperature=300, mass=m_N2) fig = maxwell_N2.plot_distribution() plt.show() print(f"Most probable velocity: {maxwell_N2.v_prob:.2f} m/s") print(f"Mean velocity: {maxwell_N2.v_mean:.2f} m/s") print(f"RMS velocity: {maxwell_N2.v_rms:.2f} m/s")

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class FreeBoltzmannSolver: """Boltzmann equation solver for collision-free flow""" def __init__(self, Nx=100, Nv=100, Lx=10.0, v_max=3.0): """ Parameters: Nx: Number of position grid points Nv: Number of velocity grid points Lx: Spatial domain size v_max: Maximum value of velocity domain """ self.Nx = Nx self.Nv = Nv self.Lx = Lx self.v_max = v_max # Grid generation self.x = np.linspace(0, Lx, Nx) self.v = np.linspace(-v_max, v_max, Nv) self.dx = self.x[1] - self.x[0] self.dv = self.v[1] - self.v[0] # Distribution function (2D array: x, v) self.f = np.zeros((Nx, Nv)) def initialize_gaussian(self, x0, v0, sigma_x, sigma_v): """Gaussian initial distribution""" X, V = np.meshgrid(self.x, self.v, indexing='ij') self.f = np.exp(-((X - x0)**2 / (2*sigma_x**2) + (V - v0)**2 / (2*sigma_v**2))) self.f /= np.sum(self.f) * self.dx * self.dv # Normalization def step(self, dt): """Time evolution (one step) - Upwind method""" f_new = self.f.copy() for i in range(self.Nx): for j in range(self.Nv): v_ij = self.v[j] if v_ij > 0: # Right-moving advection if i > 0: f_new[i, j] = self.f[i, j] - v_ij * dt / self.dx * \ (self.f[i, j] - self.f[i-1, j]) elif v_ij < 0: # Left-moving advection if i < self.Nx - 1: f_new[i, j] = self.f[i, j] - v_ij * dt / self.dx * \ (self.f[i+1, j] - self.f[i, j]) self.f = f_new def get_density(self): """Number density n(x)""" return np.sum(self.f, axis=1) * self.dv def get_mean_velocity(self): """Mean velocity u(x)""" n = self.get_density() u = np.sum(self.f * self.v[np.newaxis, :], axis=1) * self.dv return u / (n + 1e-10) def plot_state(self, t): """Plot current state""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Distribution function f(x,v) ax = axes[0, 0] im = ax.contourf(self.x, self.v, self.f.T, levels=20, cmap='viridis') ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Velocity v', fontsize=12) ax.set_title(f'Distribution function f(x,v) at t={t:.2f}', fontsize=14, fontweight='bold') plt.colorbar(im, ax=ax) # Number density n(x) ax = axes[0, 1] n = self.get_density() ax.plot(self.x, n, 'b-', linewidth=2) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Number density n(x)', fontsize=12) ax.set_title('Number Density Distribution', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Mean velocity u(x) ax = axes[1, 0] u = self.get_mean_velocity() ax.plot(self.x, u, 'r-', linewidth=2) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Mean velocity u(x)', fontsize=12) ax.set_title('Velocity Field', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Velocity distribution (at specific positions) ax = axes[1, 1] for x_pos in [self.Nx//4, self.Nx//2, 3*self.Nx//4]: ax.plot(self.v, self.f[x_pos, :], label=f'x={self.x[x_pos]:.2f}') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Distribution f(v)', fontsize=12) ax.set_title('Velocity Distribution (at each position)', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Simulation execution solver = FreeBoltzmannSolver(Nx=100, Nv=100) solver.initialize_gaussian(x0=2.0, v0=1.0, sigma_x=0.5, sigma_v=0.3) dt = 0.01 for step in range(100): solver.step(dt) fig = solver.plot_state(t=1.0) plt.show()

import numpy as np import matplotlib.pyplot as plt class HTheoremDemo: """Demonstration of H-theorem""" def __init__(self, N_particles=1000, N_bins=50): """ Parameters: N_particles: Number of particles N_bins: Number of velocity histogram bins """ self.N = N_particles self.N_bins = N_bins # Velocity range self.v_min, self.v_max = -5.0, 5.0 self.bins = np.linspace(self.v_min, self.v_max, N_bins+1) self.bin_centers = 0.5 * (self.bins[:-1] + self.bins[1:]) def initialize_two_maxwellians(self, T1=1.0, T2=4.0): """Mixture of two Maxwell distributions at different temperatures""" N1 = self.N // 2 N2 = self.N - N1 v1 = np.random.normal(0, np.sqrt(T1), N1) v2 = np.random.normal(0, np.sqrt(T2), N2) self.velocities = np.concatenate([v1, v2]) def compute_H(self, velocities): """Compute H-function""" # Approximate distribution function with histogram hist, _ = np.histogram(velocities, bins=self.bins, density=True) # H = ∫ f ln(f) dv ≈ Σ f_i ln(f_i) Δv H = 0.0 dv = self.bins[1] - self.bins[0] for f_i in hist: if f_i > 1e-10: # Avoid division by zero H += f_i * np.log(f_i) * dv return H def maxwell_boltzmann_collision(self, dt, cross_section=0.1): """Simplified collision process""" N_collisions = int(cross_section * self.N * dt) for _ in range(N_collisions): # Randomly select two particles i, j = np.random.choice(self.N, size=2, replace=False) v_i = self.velocities[i] v_j = self.velocities[j] # Center-of-mass and relative velocities v_cm = 0.5 * (v_i + v_j) v_rel = v_i - v_j # Scatter at random angle theta = np.random.uniform(0, 2*np.pi) v_rel_new = v_rel * np.array([np.cos(theta), np.sin(theta)])[0] # Post-collision velocities self.velocities[i] = v_cm + 0.5 * v_rel_new self.velocities[j] = v_cm - 0.5 * v_rel_new def simulate(self, T_max=10.0, dt=0.1): """Time evolution simulation""" times = np.arange(0, T_max, dt) H_values = [] for t in times: H = self.compute_H(self.velocities) H_values.append(H) self.maxwell_boltzmann_collision(dt) return times, H_values def plot_evolution(self, times, H_values): """Time evolution of H-function and distribution function changes""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Time evolution of H-function ax = axes[0, 0] ax.plot(times, H_values, 'b-', linewidth=2) ax.set_xlabel('Time t', fontsize=12) ax.set_ylabel('H-function', fontsize=12) ax.set_title('H-theorem: dH/dt ≤ 0', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) ax.axhline(y=H_values[-1], color='r', linestyle='--', label='Equilibrium value') ax.legend() # Time derivative of H-function ax = axes[0, 1] dH_dt = np.gradient(H_values, times[1] - times[0]) ax.plot(times, dH_dt, 'r-', linewidth=2) ax.axhline(y=0, color='k', linestyle='--', alpha=0.5) ax.set_xlabel('Time t', fontsize=12) ax.set_ylabel('dH/dt', fontsize=12) ax.set_title('Entropy Production Rate', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Initial distribution ax = axes[1, 0] demo_init = HTheoremDemo(self.N, self.N_bins) demo_init.initialize_two_maxwellians() ax.hist(demo_init.velocities, bins=self.bins, density=True, alpha=0.7, color='blue', label='Initial state') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Probability density f(v)', fontsize=12) ax.set_title('Initial Velocity Distribution (Two-temperature mixture)', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) # Final distribution ax = axes[1, 1] ax.hist(self.velocities, bins=self.bins, density=True, alpha=0.7, color='green', label='Equilibrium state') # Theoretical Maxwell distribution T_eq = np.var(self.velocities) v_theory = np.linspace(self.v_min, self.v_max, 200) f_maxwell = (1/np.sqrt(2*np.pi*T_eq)) * np.exp(-v_theory**2/(2*T_eq)) ax.plot(v_theory, f_maxwell, 'r--', linewidth=2, label='Maxwell distribution') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Probability density f(v)', fontsize=12) ax.set_title('Equilibrium Velocity Distribution', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Simulation execution demo = HTheoremDemo(N_particles=10000, N_bins=50) demo.initialize_two_maxwellians(T1=1.0, T2=4.0) times, H_values = demo.simulate(T_max=20.0, dt=0.1) fig = demo.plot_evolution(times, H_values) plt.show() print(f"Initial H value: {H_values[0]:.4f}") print(f"Final H value: {H_values[-1]:.4f}") print(f"ΔH = {H_values[-1] - H_values[0]:.4f} < 0 (entropy increase)")

import numpy as np import matplotlib.pyplot as plt class BGKSolver: """BGK equation solver""" def __init__(self, Nx=100, Nv=64, Lx=10.0, v_max=5.0, tau=1.0): """ Parameters: Nx: Number of spatial grid points Nv: Number of velocity grid points Lx: Spatial domain size v_max: Maximum value of velocity domain tau: Relaxation time """ self.Nx = Nx self.Nv = Nv self.Lx = Lx self.v_max = v_max self.tau = tau # Grid self.x = np.linspace(0, Lx, Nx) self.v = np.linspace(-v_max, v_max, Nv) self.dx = self.x[1] - self.x[0] self.dv = self.v[1] - self.v[0] # Distribution function self.f = np.zeros((Nx, Nv)) def local_maxwell(self, n, u, T): """Local Maxwell distribution""" f_eq = np.zeros((self.Nx, self.Nv)) for i in range(self.Nx): for j in range(self.Nv): v_rel = self.v[j] - u[i] f_eq[i, j] = n[i] / np.sqrt(2*np.pi*T[i]) * \ np.exp(-v_rel**2 / (2*T[i])) return f_eq def compute_moments(self): """Compute macroscopic quantities""" # Number density n = np.sum(self.f, axis=1) * self.dv # Mean velocity u = np.sum(self.f * self.v[np.newaxis, :], axis=1) * self.dv / (n + 1e-10) # Temperature v_rel_sq = np.zeros(self.Nx) for i in range(self.Nx): v_rel_sq[i] = np.sum(self.f[i, :] * (self.v - u[i])**2) * self.dv T = v_rel_sq / (n + 1e-10) return n, u, T def step_BGK(self, dt): """Time evolution of BGK equation""" # Compute macroscopic quantities n, u, T = self.compute_moments() # Local Maxwell distribution f_eq = self.local_maxwell(n, u, T) # BGK collision term collision_term = -(self.f - f_eq) / self.tau # Advection term (Upwind method) advection_term = np.zeros_like(self.f) for j in range(self.Nv): v_j = self.v[j] if v_j > 0: advection_term[1:, j] = -v_j * (self.f[1:, j] - self.f[:-1, j]) / self.dx else: advection_term[:-1, j] = -v_j * (self.f[1:, j] - self.f[:-1, j]) / self.dx # Time evolution self.f += dt * (advection_term + collision_term) def initialize_shock(self, x_shock=5.0): """Initial conditions for shock wave problem""" n_L, u_L, T_L = 2.0, 1.0, 1.0 # Left side n_R, u_R, T_R = 1.0, 0.0, 0.5 # Right side for i in range(self.Nx): if self.x[i] < x_shock: n, u, T = n_L, u_L, T_L else: n, u, T = n_R, u_R, T_R for j in range(self.Nv): v_rel = self.v[j] - u self.f[i, j] = n / np.sqrt(2*np.pi*T) * np.exp(-v_rel**2/(2*T)) def plot_solution(self, t): """Plot solution""" n, u, T = self.compute_moments() fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Distribution function ax = axes[0, 0] im = ax.contourf(self.x, self.v, self.f.T, levels=20, cmap='plasma') ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Velocity v', fontsize=12) ax.set_title(f'Distribution function f(x,v) at t={t:.2f}', fontsize=14, fontweight='bold') plt.colorbar(im, ax=ax) # Number density ax = axes[0, 1] ax.plot(self.x, n, 'b-', linewidth=2) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Number density n(x)', fontsize=12) ax.set_title('Density Profile', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Flow velocity ax = axes[1, 0] ax.plot(self.x, u, 'r-', linewidth=2) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Flow velocity u(x)', fontsize=12) ax.set_title('Velocity Profile', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Temperature ax = axes[1, 1] ax.plot(self.x, T, 'g-', linewidth=2) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Temperature T(x)', fontsize=12) ax.set_title('Temperature Profile', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Simulation execution solver = BGKSolver(Nx=100, Nv=64, tau=0.5) solver.initialize_shock(x_shock=5.0) dt = 0.01 for step in range(200): solver.step_BGK(dt) fig = solver.plot_solution(t=2.0) plt.show()

import numpy as np import matplotlib.pyplot as plt class DSMCSimulator: """Direct Simulation Monte Carlo method""" def __init__(self, N_particles=10000, Lx=10.0, N_cells=50): """ Parameters: N_particles: Number of particles Lx: Spatial domain size N_cells: Number of cells """ self.N = N_particles self.Lx = Lx self.N_cells = N_cells self.dx = Lx / N_cells # Particle positions and velocities self.x = np.random.uniform(0, Lx, N_particles) self.v = np.zeros(N_particles) def initialize_two_stream(self): """Two-stream initial conditions""" N_half = self.N // 2 self.v[:N_half] = np.random.normal(2.0, 0.5, N_half) self.v[N_half:] = np.random.normal(-2.0, 0.5, self.N - N_half) def advect(self, dt): """Particle advection""" self.x += self.v * dt # Periodic boundary conditions self.x = self.x % self.Lx def collide(self, dt, sigma_collision=0.5): """Particle collision processing""" # Classify particles by cell cell_indices = (self.x / self.dx).astype(int) for cell in range(self.N_cells): # Get particles in this cell particles_in_cell = np.where(cell_indices == cell)[0] N_cell = len(particles_in_cell) if N_cell < 2: continue # Estimate number of collisions N_collisions = int(0.5 * N_cell * sigma_collision * dt / self.dx + 0.5) for _ in range(N_collisions): # Randomly select two particles if N_cell >= 2: i, j = np.random.choice(particles_in_cell, size=2, replace=False) # Relative velocity v_rel = self.v[i] - self.v[j] # Post-collision velocities (isotropic scattering) v_cm = 0.5 * (self.v[i] + self.v[j]) theta = np.random.uniform(0, 2*np.pi) v_rel_new = abs(v_rel) * np.cos(theta) self.v[i] = v_cm + 0.5 * v_rel_new self.v[j] = v_cm - 0.5 * v_rel_new def compute_distribution(self, N_bins=50): """Compute phase space distribution""" x_bins = np.linspace(0, self.Lx, N_bins+1) v_bins = np.linspace(self.v.min()-1, self.v.max()+1, N_bins+1) H, xedges, vedges = np.histogram2d(self.x, self.v, bins=[x_bins, v_bins]) return H.T, xedges, vedges def run_simulation(self, T_max=5.0, dt=0.01): """Run simulation""" times = np.arange(0, T_max, dt) entropy = [] for t in times: self.advect(dt) self.collide(dt) # Compute entropy H, _, _ = self.compute_distribution(N_bins=30) H_norm = H / (H.sum() + 1e-10) S = -np.sum(H_norm * np.log(H_norm + 1e-10)) entropy.append(S) return times, entropy def plot_phase_space(self, t): """Plot phase space""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Phase space scatter plot ax = axes[0, 0] ax.scatter(self.x, self.v, alpha=0.3, s=1) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Velocity v', fontsize=12) ax.set_title(f'Phase space t={t:.2f}', fontsize=14, fontweight='bold') ax.set_xlim(0, self.Lx) # Phase space density ax = axes[0, 1] H, xedges, vedges = self.compute_distribution(N_bins=50) im = ax.contourf(xedges[:-1], vedges[:-1], H, levels=20, cmap='viridis') ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Velocity v', fontsize=12) ax.set_title('Distribution Function Density', fontsize=14, fontweight='bold') plt.colorbar(im, ax=ax) # Spatial density ax = axes[1, 0] x_hist, x_bins = np.histogram(self.x, bins=50, density=True) ax.bar(x_bins[:-1], x_hist, width=x_bins[1]-x_bins[0], alpha=0.7) ax.set_xlabel('Position x', fontsize=12) ax.set_ylabel('Density', fontsize=12) ax.set_title('Spatial Density Distribution', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Velocity distribution ax = axes[1, 1] v_hist, v_bins = np.histogram(self.v, bins=50, density=True) ax.bar(v_bins[:-1], v_hist, width=v_bins[1]-v_bins[0], alpha=0.7) # Maxwell distribution (theoretical) T = np.var(self.v) v_theory = np.linspace(self.v.min(), self.v.max(), 200) f_maxwell = (1/np.sqrt(2*np.pi*T)) * np.exp(-v_theory**2/(2*T)) ax.plot(v_theory, f_maxwell, 'r--', linewidth=2, label='Maxwell distribution') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Probability density', fontsize=12) ax.set_title('Velocity Distribution', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Simulation execution dsmc = DSMCSimulator(N_particles=20000, Lx=10.0, N_cells=50) dsmc.initialize_two_stream() times, entropy = dsmc.run_simulation(T_max=5.0, dt=0.01) # Plot time evolution fig_phase = dsmc.plot_phase_space(t=5.0) plt.show() # Entropy increase fig_entropy, ax = plt.subplots(figsize=(10, 6)) ax.plot(times, entropy, 'b-', linewidth=2) ax.set_xlabel('Time t', fontsize=12) ax.set_ylabel('Entropy S', fontsize=12) ax.set_title('Entropy Increase by DSMC Method', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) plt.show() print(f"Initial entropy: {entropy[0]:.4f}") print(f"Final entropy: {entropy[-1]:.4f}") print(f"Entropy increase: {entropy[-1] - entropy[0]:.4f}")

import numpy as np import matplotlib.pyplot as plt class TemperatureRelaxation: """Simulation of temperature relaxation process""" def __init__(self, N1=5000, N2=5000, T1=2.0, T2=0.5): """ Mixing of two gases at different temperatures Parameters: N1, N2: Number of particles in each gas T1, T2: Initial temperatures """ self.N1 = N1 self.N2 = N2 # Initial velocity distributions self.v1 = np.random.normal(0, np.sqrt(T1), N1) self.v2 = np.random.normal(0, np.sqrt(T2), N2) def compute_temperatures(self): """Compute temperatures of each gas""" T1 = np.var(self.v1) T2 = np.var(self.v2) return T1, T2 def collide_step(self, dt, collision_rate=0.1): """Temperature relaxation by collisions""" N_collisions = int(collision_rate * min(self.N1, self.N2) * dt) for _ in range(N_collisions): # Randomly select one particle from each gas i = np.random.randint(0, self.N1) j = np.random.randint(0, self.N2) # Collision (energy and momentum conservation) v_cm = 0.5 * (self.v1[i] + self.v2[j]) v_rel = self.v1[i] - self.v2[j] # Scattering angle theta = np.random.uniform(0, 2*np.pi) v_rel_new = abs(v_rel) * np.cos(theta) # Post-collision velocities self.v1[i] = v_cm + 0.5 * v_rel_new self.v2[j] = v_cm - 0.5 * v_rel_new def run_simulation(self, T_max=50.0, dt=0.1): """Run simulation""" times = np.arange(0, T_max, dt) T1_history = [] T2_history = [] for t in times: T1, T2 = self.compute_temperatures() T1_history.append(T1) T2_history.append(T2) self.collide_step(dt) return times, np.array(T1_history), np.array(T2_history) def plot_relaxation(self, times, T1_hist, T2_hist): """Plot relaxation process""" fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # Temperature time evolution ax = axes[0, 0] ax.plot(times, T1_hist, 'r-', linewidth=2, label='Gas 1') ax.plot(times, T2_hist, 'b-', linewidth=2, label='Gas 2') # Equilibrium temperature (theoretical) T_eq = (self.N1 * T1_hist[0] + self.N2 * T2_hist[0]) / (self.N1 + self.N2) ax.axhline(y=T_eq, color='k', linestyle='--', label=f'Equilibrium temperature={T_eq:.3f}') ax.set_xlabel('Time t', fontsize=12) ax.set_ylabel('Temperature T', fontsize=12) ax.set_title('Temperature Relaxation Process', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) # Temperature difference time evolution ax = axes[0, 1] dT = np.abs(T1_hist - T2_hist) ax.semilogy(times, dT, 'g-', linewidth=2) # Exponential decay fit from scipy.optimize import curve_fit def exp_decay(t, A, tau): return A * np.exp(-t/tau) try: popt, _ = curve_fit(exp_decay, times[10:], dT[10:], p0=[dT[10], 10.0]) ax.semilogy(times, exp_decay(times, *popt), 'r--', label=f'Fit: τ={popt[1]:.2f}') ax.legend() except: pass ax.set_xlabel('Time t', fontsize=12) ax.set_ylabel('Temperature difference |T1 - T2|', fontsize=12) ax.set_title('Temperature Difference Relaxation (log scale)', fontsize=14, fontweight='bold') ax.grid(True, alpha=0.3) # Initial velocity distribution ax = axes[1, 0] v1_init = np.random.normal(0, np.sqrt(T1_hist[0]), self.N1) v2_init = np.random.normal(0, np.sqrt(T2_hist[0]), self.N2) ax.hist(v1_init, bins=50, alpha=0.5, density=True, color='red', label='Gas 1 (high temp)') ax.hist(v2_init, bins=50, alpha=0.5, density=True, color='blue', label='Gas 2 (low temp)') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Probability density', fontsize=12) ax.set_title('Initial Velocity Distribution', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) # Final velocity distribution ax = axes[1, 1] ax.hist(self.v1, bins=50, alpha=0.5, density=True, color='red', label='Gas 1') ax.hist(self.v2, bins=50, alpha=0.5, density=True, color='blue', label='Gas 2') # Equilibrium Maxwell distribution v_range = np.linspace(-5, 5, 200) f_eq = (1/np.sqrt(2*np.pi*T_eq)) * np.exp(-v_range**2/(2*T_eq)) ax.plot(v_range, f_eq, 'k--', linewidth=2, label='Maxwell distribution') ax.set_xlabel('Velocity v', fontsize=12) ax.set_ylabel('Probability density', fontsize=12) ax.set_title('Equilibrium Velocity Distribution', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() return fig # Simulation execution relaxation = TemperatureRelaxation(N1=10000, N2=10000, T1=3.0, T2=0.5) times, T1_hist, T2_hist = relaxation.run_simulation(T_max=50.0, dt=0.1) fig = relaxation.plot_relaxation(times, T1_hist, T2_hist) plt.show() print(f"Initial temperatures: T1={T1_hist[0]:.3f}, T2={T2_hist[0]:.3f}") print(f"Final temperatures: T1={T1_hist[-1]:.3f}, T2={T2_hist[-1]:.3f}") print(f"Theoretical equilibrium temperature: {(10000*T1_hist[0] + 10000*T2_hist[0])/20000:.3f}")

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation class RarefiedGasFlow: """Comprehensive simulation of rarefied gas flow""" def __init__(self, Nx=100, Nv=64, Lx=20.0, v_max=5.0, tau=0.5): self.Nx = Nx self.Nv = Nv self.Lx = Lx self.v_max = v_max self.tau = tau self.x = np.linspace(0, Lx, Nx) self.v = np.linspace(-v_max, v_max, Nv) self.dx = self.x[1] - self.x[0] self.dv = self.v[1] - self.v[0] self.f = np.zeros((Nx, Nv)) # Obstacle (wall blocking flow) self.obstacle_start = Nx // 3 self.obstacle_end = Nx // 3 + Nx // 10 def initialize_flow(self, n0=1.0, u0=2.0, T0=1.0): """Initial conditions for flow""" for i in range(self.Nx): for j in range(self.Nv): v_rel = self.v[j] - u0 self.f[i, j] = n0 / np.sqrt(2*np.pi*T0) * \ np.exp(-v_rel**2/(2*T0)) def compute_moments(self): """Macroscopic quantities""" n = np.sum(self.f, axis=1) * self.dv u = np.sum(self.f * self.v[np.newaxis, :], axis=1) * self.dv / (n + 1e-10) v_rel_sq = np.zeros(self.Nx) for i in range(self.Nx): v_rel_sq[i] = np.sum(self.f[i, :] * (self.v - u[i])**2) * self.dv T = v_rel_sq / (n + 1e-10) return n, u, T def apply_boundary_conditions(self): """Boundary conditions (reflection at obstacle)""" # Specular reflection at obstacle for i in range(self.obstacle_start, self.obstacle_end): for j in range(self.Nv): if self.v[j] > 0 and i == self.obstacle_start: # Reflect right-moving velocity to left-moving j_reflect = self.Nv - 1 - j self.f[i, j_reflect] = self.f[i, j] self.f[i, j] = 0 def step(self, dt): """Time evolution (BGK + advection)""" n, u, T = self.compute_moments() # Local Maxwell distribution f_eq = np.zeros_like(self.f) for i in range(self.Nx): for j in range(self.Nv): v_rel = self.v[j] - u[i] f_eq[i, j] = n[i] / np.sqrt(2*np.pi*T[i]) * \ np.exp(-v_rel**2/(2*T[i])) # BGK collision collision = -(self.f - f_eq) / self.tau # Advection (Upwind) advection = np.zeros_like(self.f) for j in range(self.Nv): v_j = self.v[j] if v_j > 0: advection[1:, j] = -v_j * (self.f[1:, j] - self.f[:-1, j]) / self.dx advection[0, j] = -v_j * (self.f[0, j] - self.f[-1, j]) / self.dx # Periodic boundary else: advection[:-1, j] = -v_j * (self.f[1:, j] - self.f[:-1, j]) / self.dx advection[-1, j] = -v_j * (self.f[0, j] - self.f[-1, j]) / self.dx # Time evolution self.f += dt * (advection + collision) # Apply boundary conditions self.apply_boundary_conditions() def plot_solution(self, t): """Solution visualization""" n, u, T = self.compute_moments() fig = plt.figure(figsize=(16, 12)) gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3) # Distribution function ax1 = fig.add_subplot(gs[0, :]) im = ax1.contourf(self.x, self.v, self.f.T, levels=30, cmap='viridis') # Display obstacle ax1.axvspan(self.x[self.obstacle_start], self.x[self.obstacle_end], alpha=0.3, color='red', label='Obstacle') ax1.set_xlabel('Position x', fontsize=12) ax1.set_ylabel('Velocity v', fontsize=12) ax1.set_title(f'Phase space distribution t={t:.2f}', fontsize=14, fontweight='bold') ax1.legend(loc='upper right') plt.colorbar(im, ax=ax1, label='f(x,v)') # Number density ax2 = fig.add_subplot(gs[1, 0]) ax2.plot(self.x, n, 'b-', linewidth=2) ax2.axvspan(self.x[self.obstacle_start], self.x[self.obstacle_end], alpha=0.2, color='red') ax2.set_xlabel('Position x', fontsize=12) ax2.set_ylabel('Number density n', fontsize=12) ax2.set_title('Density Field', fontsize=14, fontweight='bold') ax2.grid(True, alpha=0.3) # Flow velocity ax3 = fig.add_subplot(gs[1, 1]) ax3.plot(self.x, u, 'r-', linewidth=2) ax3.axvspan(self.x[self.obstacle_start], self.x[self.obstacle_end], alpha=0.2, color='red') ax3.set_xlabel('Position x', fontsize=12) ax3.set_ylabel('Flow velocity u', fontsize=12) ax3.set_title('Velocity Field', fontsize=14, fontweight='bold') ax3.grid(True, alpha=0.3) # Temperature ax4 = fig.add_subplot(gs[2, 0]) ax4.plot(self.x, T, 'g-', linewidth=2) ax4.axvspan(self.x[self.obstacle_start], self.x[self.obstacle_end], alpha=0.2, color='red') ax4.set_xlabel('Position x', fontsize=12) ax4.set_ylabel('Temperature T', fontsize=12) ax4.set_title('Temperature Field', fontsize=14, fontweight='bold') ax4.grid(True, alpha=0.3) # Pressure ax5 = fig.add_subplot(gs[2, 1]) p = n * T # Ideal gas ax5.plot(self.x, p, 'purple', linewidth=2) ax5.axvspan(self.x[self.obstacle_start], self.x[self.obstacle_end], alpha=0.2, color='red') ax5.set_xlabel('Position x', fontsize=12) ax5.set_ylabel('Pressure p', fontsize=12) ax5.set_title('Pressure Field', fontsize=14, fontweight='bold') ax5.grid(True, alpha=0.3) return fig # Simulation execution flow = RarefiedGasFlow(Nx=150, Nv=80, Lx=20.0, v_max=6.0, tau=0.3) flow.initialize_flow(n0=1.0, u0=3.0, T0=1.0) dt = 0.005 for step in range(1000): flow.step(dt) if step % 100 == 0: print(f"Step {step}/1000") fig = flow.plot_solution(t=5.0) plt.show() print("\n=== Rarefied gas flow simulation complete ===")