1#!/usr'bin/env python3
2"""
3Rayleigh-Taylor Instability in MHD
4
5Analyzes the 2D MHD Rayleigh-Taylor instability with magnetic field.
6
7Setup: Heavy fluid (ρ₂) over light fluid (ρ₁) with transverse B field.
8Growth rate: γ² = gk - (k·B)²/(μ₀ρ)
9
10Shows magnetic stabilization of short wavelengths.
11
12Author: Claude
13Date: 2026-02-14
14"""
15
16import numpy as np
17import matplotlib.pyplot as plt
18
19
20def dispersion_relation(k, g, rho_1, rho_2, B, theta):
21 """
22 Compute growth rate for RT instability with magnetic field.
23
24 γ² = A * g * k - (k·B)²/(μ₀ρ_avg)
25
26 where A = (ρ₂ - ρ₁)/(ρ₂ + ρ₁) is the Atwood number.
27
28 Parameters
29 ----------
30 k : float or array_like
31 Wave number (m⁻¹)
32 g : float
33 Gravitational acceleration (m/s²)
34 rho_1, rho_2 : float
35 Fluid densities (kg/m³), rho_2 > rho_1
36 B : float
37 Magnetic field strength (T)
38 theta : float
39 Angle between k and B (radians)
40
41 Returns
42 -------
43 gamma : float or ndarray
44 Growth rate (s⁻¹)
45 """
46 mu_0 = 4 * np.pi * 1e-7
47
48 # Atwood number
49 A = (rho_2 - rho_1) / (rho_2 + rho_1)
50
51 # Average density
52 rho_avg = (rho_1 + rho_2) / 2
53
54 # Gravitational term (destabilizing)
55 term_grav = A * g * np.abs(k)
56
57 # Magnetic stabilization
58 k_parallel = np.abs(k) * np.cos(theta) # Component parallel to B
59 term_mag = (B * k_parallel)**2 / (mu_0 * rho_avg)
60
61 # Growth rate squared
62 gamma_sq = term_grav - term_mag
63
64 # Return growth rate (zero if stable)
65 gamma = np.where(gamma_sq > 0, np.sqrt(gamma_sq), 0.0)
66
67 return gamma
68
69
70def critical_wavelength(g, rho_1, rho_2, B, theta):
71 """
72 Compute critical wavelength for stabilization.
73
74 At λ_c, growth rate goes to zero: γ(k_c) = 0
75
76 Parameters
77 ----------
78 g : float
79 Gravity
80 rho_1, rho_2 : float
81 Densities
82 B : float
83 Magnetic field
84 theta : float
85 Angle
86
87 Returns
88 -------
89 lambda_c : float
90 Critical wavelength (m)
91 """
92 mu_0 = 4 * np.pi * 1e-7
93 A = (rho_2 - rho_1) / (rho_2 + rho_1)
94 rho_avg = (rho_1 + rho_2) / 2
95
96 if np.abs(np.cos(theta)) < 1e-10:
97 # Perpendicular: no stabilization
98 return np.inf
99
100 # k_c from γ²(k_c) = 0
101 k_c = A * g * mu_0 * rho_avg / (B * np.cos(theta))**2
102
103 lambda_c = 2 * np.pi / k_c if k_c > 0 else np.inf
104
105 return lambda_c
106
107
108def plot_growth_vs_wavenumber(g, rho_1, rho_2, B_values, theta=0):
109 """
110 Plot growth rate vs wavenumber for different B fields.
111
112 Parameters
113 ----------
114 g : float
115 Gravity
116 rho_1, rho_2 : float
117 Densities
118 B_values : array_like
119 Magnetic field strengths
120 theta : float
121 Angle (default 0 = parallel)
122 """
123 k = np.logspace(-2, 3, 200) # Wave numbers from 0.01 to 1000 m⁻¹
124
125 fig, ax = plt.subplots(figsize=(10, 7))
126
127 colors = plt.cm.viridis(np.linspace(0, 1, len(B_values)))
128
129 for i, B in enumerate(B_values):
130 gamma = dispersion_relation(k, g, rho_1, rho_2, B, theta)
131
132 label = f'B = {B:.2f} T'
133 if B > 0:
134 lambda_c = critical_wavelength(g, rho_1, rho_2, B, theta)
135 if np.isfinite(lambda_c):
136 label += f' (λ_c={lambda_c*100:.1f} cm)'
137
138 ax.loglog(k, gamma, color=colors[i], linewidth=2.5, label=label)
139
140 # No field case
141 gamma_0 = dispersion_relation(k, g, rho_1, rho_2, 0, theta)
142 ax.loglog(k, gamma_0, 'k--', linewidth=2, label='B = 0 (no field)')
143
144 ax.set_xlabel('Wave number k (m⁻¹)', fontsize=13)
145 ax.set_ylabel('Growth rate γ (s⁻¹)', fontsize=13)
146 ax.set_title('Rayleigh-Taylor Growth Rate: Magnetic Stabilization',
147 fontsize=14, fontweight='bold')
148 ax.legend(fontsize=10)
149 ax.grid(True, alpha=0.3, which='both')
150
151 # Add annotations
152 ax.text(0.05, 0.95,
153 f'ρ₁ = {rho_1:.1e} kg/m³\nρ₂ = {rho_2:.1e} kg/m³\ng = {g:.1f} m/s²\nθ = {np.degrees(theta):.0f}°',
154 transform=ax.transAxes, fontsize=10,
155 verticalalignment='top',
156 bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.7))
157
158 plt.tight_layout()
159 plt.savefig('rt_growth_vs_k.png', dpi=150, bbox_inches='tight')
160 plt.show()
161
162
163def plot_growth_vs_angle(g, rho_1, rho_2, B, k_values):
164 """
165 Plot growth rate vs angle between k and B.
166
167 Parameters
168 ----------
169 g : float
170 Gravity
171 rho_1, rho_2 : float
172 Densities
173 B : float
174 Magnetic field
175 k_values : array_like
176 Wave numbers to plot
177 """
178 theta = np.linspace(0, np.pi, 200) # 0 to 180 degrees
179
180 fig, ax = plt.subplots(figsize=(10, 7))
181
182 colors = plt.cm.plasma(np.linspace(0, 1, len(k_values)))
183
184 for i, k in enumerate(k_values):
185 gamma = np.array([dispersion_relation(k, g, rho_1, rho_2, B, t)
186 for t in theta])
187
188 lambda_val = 2 * np.pi / k
189 ax.plot(np.degrees(theta), gamma, color=colors[i], linewidth=2.5,
190 label=f'λ = {lambda_val*100:.1f} cm')
191
192 ax.set_xlabel('Angle θ between k and B (degrees)', fontsize=13)
193 ax.set_ylabel('Growth rate γ (s⁻¹)', fontsize=13)
194 ax.set_title('RT Growth Rate vs Field Orientation',
195 fontsize=14, fontweight='bold')
196 ax.legend(fontsize=10)
197 ax.grid(True, alpha=0.3)
198 ax.set_xlim([0, 180])
199
200 # Mark key angles
201 ax.axvline(x=0, color='r', linestyle='--', alpha=0.5, label='θ=0° (parallel)')
202 ax.axvline(x=90, color='b', linestyle='--', alpha=0.5, label='θ=90° (perpendicular)')
203
204 plt.tight_layout()
205 plt.savefig('rt_growth_vs_angle.png', dpi=150, bbox_inches='tight')
206 plt.show()
207
208
209def plot_2d_mode_structure(k_x, k_z, g, rho_1, rho_2, B, t_max=0.1):
210 """
211 Visualize 2D mode structure evolution.
212
213 Parameters
214 ----------
215 k_x, k_z : float
216 Wave vector components
217 g : float
218 Gravity
219 rho_1, rho_2 : float
220 Densities
221 B : float
222 Magnetic field (in x-direction)
223 t_max : float
224 Maximum time
225 """
226 # Grid
227 x = np.linspace(0, 2*np.pi/k_x, 100)
228 z = np.linspace(-0.05, 0.05, 100)
229 X, Z = np.meshgrid(x, z)
230
231 # Wave vector and angle
232 k = np.sqrt(k_x**2 + k_z**2)
233 theta = np.arctan2(k_z, k_x)
234
235 # Growth rate
236 gamma = dispersion_relation(k, g, rho_1, rho_2, B, theta)
237
238 # Time snapshots
239 times = np.linspace(0, t_max, 4)
240
241 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
242 axes = axes.flatten()
243
244 for i, t in enumerate(times):
245 ax = axes[i]
246
247 # Displacement: ξ ~ exp(γt) * sin(k·x)
248 amplitude = np.exp(gamma * t)
249 xi = amplitude * np.sin(k_x * X + k_z * Z)
250
251 # Plot
252 im = ax.contourf(X*100, Z*100, xi, levels=20, cmap='RdBu_r')
253 ax.contour(X*100, Z*100, xi, levels=10, colors='black',
254 linewidths=0.5, alpha=0.3)
255
256 ax.set_xlabel('x (cm)', fontsize=11)
257 ax.set_ylabel('z (cm)', fontsize=11)
258 ax.set_title(f't = {t*1e3:.2f} ms, amplitude = {amplitude:.2f}',
259 fontsize=12, fontweight='bold')
260 ax.set_aspect('equal')
261 plt.colorbar(im, ax=ax, label='Displacement ξ')
262
263 plt.suptitle(f'RT Instability Evolution (γ = {gamma:.2f} s⁻¹, B = {B:.2f} T)',
264 fontsize=14, fontweight='bold')
265 plt.tight_layout()
266 plt.savefig('rt_mode_evolution.png', dpi=150, bbox_inches='tight')
267 plt.show()
268
269
270def plot_stability_diagram(g, rho_1, rho_2):
271 """
272 Create stability diagram in (B, λ) space.
273
274 Parameters
275 ----------
276 g : float
277 Gravity
278 rho_1, rho_2 : float
279 Densities
280 """
281 # Grid
282 B_arr = np.linspace(0, 0.5, 100)
283 lambda_arr = np.logspace(-3, 0, 100) # 1 mm to 1 m
284
285 B_grid, lambda_grid = np.meshgrid(B_arr, lambda_arr)
286 k_grid = 2 * np.pi / lambda_grid
287
288 # Growth rate (parallel field)
289 gamma_grid = np.zeros_like(B_grid)
290
291 for i in range(len(B_arr)):
292 for j in range(len(lambda_arr)):
293 gamma_grid[j, i] = dispersion_relation(k_grid[j, i], g, rho_1, rho_2,
294 B_grid[j, i], theta=0)
295
296 # Plot
297 fig, ax = plt.subplots(figsize=(10, 8))
298
299 # Filled contours
300 im = ax.contourf(B_grid, lambda_grid*100, gamma_grid,
301 levels=20, cmap='hot_r')
302
303 # Stability boundary (γ = 0)
304 cs = ax.contour(B_grid, lambda_grid*100, gamma_grid,
305 levels=[0], colors='blue', linewidths=3)
306 ax.clabel(cs, inline=True, fontsize=11, fmt='γ=0 (stable)')
307
308 ax.set_xlabel('Magnetic field B (T)', fontsize=13)
309 ax.set_ylabel('Wavelength λ (cm)', fontsize=13)
310 ax.set_yscale('log')
311 ax.set_title('RT Stability Diagram',
312 fontsize=14, fontweight='bold')
313
314 cbar = plt.colorbar(im, ax=ax)
315 cbar.set_label('Growth rate γ (s⁻¹)', fontsize=12)
316
317 # Add text regions
318 ax.text(0.4, 50, 'UNSTABLE', fontsize=16, fontweight='bold',
319 color='white', ha='center',
320 bbox=dict(boxstyle='round', facecolor='red', alpha=0.6))
321 ax.text(0.4, 0.5, 'STABLE', fontsize=16, fontweight='bold',
322 color='black', ha='center',
323 bbox=dict(boxstyle='round', facecolor='cyan', alpha=0.6))
324
325 plt.tight_layout()
326 plt.savefig('rt_stability_diagram.png', dpi=150, bbox_inches='tight')
327 plt.show()
328
329
330def main():
331 """Main execution function."""
332 # Parameters
333 g = 10.0 # m/s² (gravity or equivalent acceleration)
334 rho_1 = 0.1 # kg/m³ (light fluid, top initially)
335 rho_2 = 1.0 # kg/m³ (heavy fluid, bottom)
336
337 print("Rayleigh-Taylor Instability with Magnetic Field")
338 print("=" * 60)
339 print(f"Configuration:")
340 print(f" Light fluid density ρ₁: {rho_1:.2f} kg/m³")
341 print(f" Heavy fluid density ρ₂: {rho_2:.2f} kg/m³")
342 print(f" Atwood number A: {(rho_2-rho_1)/(rho_2+rho_1):.3f}")
343 print(f" Acceleration g: {g:.1f} m/s²")
344 print()
345
346 # Test different field strengths
347 B_values = [0, 0.05, 0.1, 0.2, 0.3]
348
349 print("Growth rates vs wavenumber:")
350 plot_growth_vs_wavenumber(g, rho_1, rho_2, B_values, theta=0)
351 print(" Saved as 'rt_growth_vs_k.png'")
352
353 # Growth vs angle
354 B_test = 0.2 # T
355 k_test = [10, 50, 100, 200] # m⁻¹
356
357 print("\nGrowth rate vs field orientation:")
358 plot_growth_vs_angle(g, rho_1, rho_2, B_test, k_test)
359 print(" Saved as 'rt_growth_vs_angle.png'")
360
361 # 2D mode evolution
362 print("\n2D mode structure evolution:")
363 k_x, k_z = 100, 0 # Parallel to field
364 plot_2d_mode_structure(k_x, k_z, g, rho_1, rho_2, B=0.1, t_max=0.05)
365 print(" Saved as 'rt_mode_evolution.png'")
366
367 # Stability diagram
368 print("\nStability diagram in (B, λ) space:")
369 plot_stability_diagram(g, rho_1, rho_2)
370 print(" Saved as 'rt_stability_diagram.png'")
371
372 # Critical wavelength examples
373 print("\nCritical wavelengths (parallel field):")
374 for B in [0.05, 0.1, 0.2, 0.3]:
375 lambda_c = critical_wavelength(g, rho_1, rho_2, B, theta=0)
376 if np.isfinite(lambda_c):
377 print(f" B = {B:.2f} T: λ_c = {lambda_c*100:.2f} cm")
378
379
380if __name__ == '__main__':
381 main()