1#!/usr/bin/env python3
2"""
3MHD Energy Principle
4
5Evaluates the potential energy δW for a given plasma perturbation ξ.
6
7The energy principle states that an equilibrium is stable if δW > 0
8for all perturbations ξ.
9
10δW = δW_field + δW_compression + δW_pressure
11
12where:
13- δW_field: field line bending energy
14- δW_compression: magnetic compression energy
15- δW_pressure: pressure-driven destabilization
16
17Author: Claude
18Date: 2026-02-14
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23from scipy.integrate import simpson
24
25
26class EnergyPrinciple:
27 """
28 Class for evaluating the MHD energy principle.
29
30 Attributes
31 ----------
32 r : ndarray
33 Radial grid
34 """
35
36 def __init__(self, r_max=0.01, n_points=200):
37 """
38 Initialize with radial grid.
39
40 Parameters
41 ----------
42 r_max : float
43 Maximum radius (m)
44 n_points : int
45 Number of radial points
46 """
47 self.r = np.linspace(0, r_max, n_points)
48 self.r_max = r_max
49 self.n_points = n_points
50
51 def equilibrium_fields(self, B_z, I_total):
52 """
53 Compute equilibrium magnetic fields for a Z-pinch.
54
55 Parameters
56 ----------
57 B_z : float
58 Axial magnetic field (T)
59 I_total : float
60 Total current (A)
61
62 Returns
63 -------
64 B_theta : ndarray
65 Azimuthal field
66 p : ndarray
67 Pressure profile
68 """
69 mu_0 = 4 * np.pi * 1e-7
70
71 # Current density (parabolic profile)
72 J_z = 2 * I_total / (np.pi * self.r_max**2) * (1 - (self.r/self.r_max)**2)
73 J_z = np.where(self.r <= self.r_max, J_z, 0.0)
74
75 # Azimuthal field from Ampere's law
76 B_theta = np.zeros_like(self.r)
77 for i, ri in enumerate(self.r):
78 if ri > 1e-10:
79 I_enclosed = np.pi * ri**2 * 2 * I_total / (np.pi * self.r_max**2) * (1 - 0.5*(ri/self.r_max)**2)
80 B_theta[i] = mu_0 * I_enclosed / (2 * np.pi * ri)
81
82 # Pressure from force balance
83 p = np.zeros_like(self.r)
84 p_edge = 100.0 # Pa
85 p[-1] = p_edge
86
87 for i in range(len(self.r)-2, -1, -1):
88 dr = self.r[i+1] - self.r[i]
89 p[i] = p[i+1] - J_z[i] * B_theta[i] * dr
90
91 return B_theta, p, J_z
92
93 def perturbation_radial(self, m, amplitude=1e-4):
94 """
95 Generate radial displacement perturbation.
96
97 ξ_r(r) = A * r^m * (1 - r/r_max)
98
99 Parameters
100 ----------
101 m : int
102 Mode number
103 amplitude : float
104 Perturbation amplitude (m)
105
106 Returns
107 -------
108 xi_r : ndarray
109 Radial displacement
110 """
111 xi_r = amplitude * (self.r/self.r_max)**m * (1 - self.r/self.r_max)
112 return xi_r
113
114 def field_bending_energy(self, xi_r, B_z, B_theta, m):
115 """
116 Compute field line bending energy.
117
118 δW_bend = (1/2μ₀) ∫ |δB_⊥|² dV
119
120 Parameters
121 ----------
122 xi_r : ndarray
123 Radial displacement
124 B_z : float
125 Axial field
126 B_theta : ndarray
127 Azimuthal field
128 m : int
129 Mode number
130
131 Returns
132 -------
133 W_bend : float
134 Field bending energy (J/m)
135 """
136 mu_0 = 4 * np.pi * 1e-7
137
138 # Perturbed field components (simplified)
139 # δB ~ B·∇ξ
140 dxi_dr = np.gradient(xi_r, self.r)
141
142 # Perpendicular field perturbation
143 delta_B_perp_sq = B_z**2 * (m * xi_r / (self.r + 1e-10))**2 + B_theta**2 * dxi_dr**2
144
145 # Integrate over volume (per unit length in z)
146 integrand = delta_B_perp_sq * self.r
147 W_bend = (np.pi / mu_0) * simpson(integrand, x=self.r)
148
149 return W_bend
150
151 def compression_energy(self, xi_r, B_z, B_theta):
152 """
153 Compute magnetic compression energy.
154
155 δW_comp = (1/2μ₀) ∫ |B|² |∇·ξ|² dV
156
157 Parameters
158 ----------
159 xi_r : ndarray
160 Radial displacement
161 B_z : float
162 Axial field
163 B_theta : ndarray
164 Azimuthal field
165
166 Returns
167 -------
168 W_comp : float
169 Compression energy (J/m)
170 """
171 mu_0 = 4 * np.pi * 1e-7
172
173 # Divergence of displacement (cylindrical): ∇·ξ = (1/r) d(r*ξ_r)/dr
174 div_xi = np.gradient(self.r * xi_r, self.r) / (self.r + 1e-10)
175
176 # Total field squared
177 B_squared = B_z**2 + B_theta**2
178
179 # Integrate
180 integrand = B_squared * div_xi**2 * self.r
181 W_comp = (np.pi / mu_0) * simpson(integrand, x=self.r)
182
183 return W_comp
184
185 def pressure_energy(self, xi_r, p):
186 """
187 Compute pressure-driven energy (destabilizing).
188
189 δW_p = -∫ ξ·∇p (∇·ξ) dV
190
191 Parameters
192 ----------
193 xi_r : ndarray
194 Radial displacement
195 p : ndarray
196 Pressure profile
197
198 Returns
199 -------
200 W_p : float
201 Pressure energy (J/m)
202 """
203 # Pressure gradient
204 dp_dr = np.gradient(p, self.r)
205
206 # Divergence
207 div_xi = np.gradient(self.r * xi_r, self.r) / (self.r + 1e-10)
208
209 # This term is typically destabilizing (negative)
210 integrand = -xi_r * dp_dr * div_xi * self.r
211 W_p = 2 * np.pi * simpson(integrand, x=self.r)
212
213 return W_p
214
215 def total_energy(self, xi_r, B_z, B_theta, p, m):
216 """
217 Compute total potential energy δW.
218
219 Parameters
220 ----------
221 xi_r : ndarray
222 Radial displacement
223 B_z : float
224 Axial field
225 B_theta : ndarray
226 Azimuthal field
227 p : ndarray
228 Pressure
229 m : int
230 Mode number
231
232 Returns
233 -------
234 dict
235 Dictionary with energy components
236 """
237 W_bend = self.field_bending_energy(xi_r, B_z, B_theta, m)
238 W_comp = self.compression_energy(xi_r, B_z, B_theta)
239 W_p = self.pressure_energy(xi_r, p)
240
241 W_total = W_bend + W_comp + W_p
242
243 return {
244 'W_bend': W_bend,
245 'W_comp': W_comp,
246 'W_p': W_p,
247 'W_total': W_total,
248 'stable': W_total > 0
249 }
250
251
252def scan_mode_numbers(ep, B_z, B_theta, p, m_max=5, n_max=5):
253 """
254 Scan over mode numbers (m, n) and compute δW.
255
256 Parameters
257 ----------
258 ep : EnergyPrinciple
259 Energy principle instance
260 B_z : float
261 Axial field
262 B_theta : ndarray
263 Azimuthal field
264 p : ndarray
265 Pressure
266 m_max, n_max : int
267 Maximum mode numbers
268
269 Returns
270 -------
271 W_map : ndarray
272 2D array of δW values
273 """
274 m_values = np.arange(0, m_max + 1)
275 n_values = np.arange(1, n_max + 1)
276
277 W_map = np.zeros((len(m_values), len(n_values)))
278
279 for i, m in enumerate(m_values):
280 for j, n in enumerate(n_values):
281 if m == 0:
282 m_eff = 1 # Avoid m=0 singularity
283 else:
284 m_eff = m
285
286 xi_r = ep.perturbation_radial(m_eff, amplitude=1e-4)
287 result = ep.total_energy(xi_r, B_z, B_theta, p, m_eff)
288 W_map[i, j] = result['W_total']
289
290 return W_map, m_values, n_values
291
292
293def plot_energy_decomposition(ep, B_z, B_theta, p, m=2):
294 """
295 Plot energy decomposition for a given mode.
296
297 Parameters
298 ----------
299 ep : EnergyPrinciple
300 Energy principle instance
301 B_z : float
302 Axial field
303 B_theta : ndarray
304 Azimuthal field
305 p : ndarray
306 Pressure
307 m : int
308 Mode number
309 """
310 xi_r = ep.perturbation_radial(m, amplitude=1e-4)
311 result = ep.total_energy(xi_r, B_z, B_theta, p, m)
312
313 fig, axes = plt.subplots(1, 2, figsize=(12, 5))
314
315 # Plot perturbation
316 axes[0].plot(ep.r * 100, xi_r * 1e6, 'b-', linewidth=2)
317 axes[0].set_xlabel('Radius (cm)', fontsize=12)
318 axes[0].set_ylabel('ξ_r (μm)', fontsize=12)
319 axes[0].set_title(f'Radial Displacement (m={m})', fontsize=13, fontweight='bold')
320 axes[0].grid(True, alpha=0.3)
321
322 # Plot energy components
323 components = ['W_bend', 'W_comp', 'W_p', 'W_total']
324 values = [result[key] for key in components]
325 colors = ['blue', 'green', 'red', 'purple']
326 labels = ['Field Bending', 'Compression', 'Pressure', 'Total']
327
328 bars = axes[1].bar(labels, values, color=colors, alpha=0.7, edgecolor='black')
329 axes[1].axhline(y=0, color='black', linestyle='-', linewidth=1)
330 axes[1].set_ylabel('Energy δW (J/m)', fontsize=12)
331 axes[1].set_title('Energy Decomposition', fontsize=13, fontweight='bold')
332 axes[1].grid(True, alpha=0.3, axis='y')
333
334 # Add value labels on bars
335 for bar, val in zip(bars, values):
336 height = bar.get_height()
337 axes[1].text(bar.get_x() + bar.get_width()/2., height,
338 f'{val:.2e}',
339 ha='center', va='bottom' if height > 0 else 'top',
340 fontsize=9)
341
342 # Stability text
343 stability_text = "STABLE" if result['stable'] else "UNSTABLE"
344 stability_color = "green" if result['stable'] else "red"
345 axes[1].text(0.5, 0.95, stability_text,
346 transform=axes[1].transAxes,
347 fontsize=14, fontweight='bold',
348 color=stability_color,
349 ha='center', va='top',
350 bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
351
352 plt.suptitle(f'MHD Energy Principle Analysis (m={m})',
353 fontsize=14, fontweight='bold')
354 plt.tight_layout()
355 plt.savefig('energy_principle.png', dpi=150, bbox_inches='tight')
356 plt.show()
357
358
359def plot_stability_map(W_map, m_values, n_values):
360 """
361 Plot 2D stability map.
362
363 Parameters
364 ----------
365 W_map : ndarray
366 Energy values
367 m_values, n_values : ndarray
368 Mode numbers
369 """
370 fig, ax = plt.subplots(figsize=(10, 8))
371
372 # Plot heatmap
373 im = ax.imshow(W_map, cmap='RdBu_r', aspect='auto',
374 extent=[n_values[0]-0.5, n_values[-1]+0.5,
375 m_values[-1]+0.5, m_values[0]-0.5],
376 vmin=-np.max(np.abs(W_map)), vmax=np.max(np.abs(W_map)))
377
378 # Contour at δW = 0
379 ax.contour(n_values, m_values, W_map, levels=[0],
380 colors='black', linewidths=3)
381
382 # Labels
383 ax.set_xlabel('Toroidal mode number n', fontsize=13)
384 ax.set_ylabel('Poloidal mode number m', fontsize=13)
385 ax.set_title('MHD Stability Map (δW > 0: stable, δW < 0: unstable)',
386 fontsize=14, fontweight='bold')
387
388 # Colorbar
389 cbar = plt.colorbar(im, ax=ax)
390 cbar.set_label('δW (J/m)', fontsize=12)
391
392 # Add text annotations for stable/unstable regions
393 for i, m in enumerate(m_values):
394 for j, n in enumerate(n_values):
395 if W_map[i, j] > 0:
396 marker = '✓'
397 color = 'green'
398 else:
399 marker = '✗'
400 color = 'red'
401 ax.text(n, m, marker, ha='center', va='center',
402 color=color, fontsize=16, fontweight='bold')
403
404 plt.tight_layout()
405 plt.savefig('stability_map.png', dpi=150, bbox_inches='tight')
406 plt.show()
407
408
409def main():
410 """Main execution function."""
411 # Initialize
412 r_max = 0.01 # 1 cm
413 ep = EnergyPrinciple(r_max=r_max, n_points=200)
414
415 # Equilibrium parameters
416 B_z = 0.5 # T
417 I_total = 50e3 # 50 kA
418
419 print("MHD Energy Principle Analysis")
420 print("=" * 60)
421 print(f"Configuration: Z-pinch")
422 print(f" Radius: {r_max*100:.1f} cm")
423 print(f" Axial field: {B_z:.2f} T")
424 print(f" Current: {I_total/1e3:.1f} kA")
425 print()
426
427 # Compute equilibrium
428 B_theta, p, J_z = ep.equilibrium_fields(B_z, I_total)
429
430 print("Equilibrium computed:")
431 print(f" Peak B_θ: {np.max(B_theta)*1e3:.2f} mT")
432 print(f" Peak pressure: {np.max(p)/1e3:.2f} kPa")
433 print()
434
435 # Analyze single mode
436 m = 2
437 print(f"Analyzing m={m} mode...")
438 xi_r = ep.perturbation_radial(m, amplitude=1e-4)
439 result = ep.total_energy(xi_r, B_z, B_theta, p, m)
440
441 print(f" Field bending energy: {result['W_bend']:.3e} J/m")
442 print(f" Compression energy: {result['W_comp']:.3e} J/m")
443 print(f" Pressure energy: {result['W_p']:.3e} J/m")
444 print(f" Total δW: {result['W_total']:.3e} J/m")
445 print(f" Stability: {'STABLE' if result['stable'] else 'UNSTABLE'}")
446 print()
447
448 # Plot energy decomposition
449 plot_energy_decomposition(ep, B_z, B_theta, p, m=2)
450 print("Energy decomposition plot saved as 'energy_principle.png'")
451
452 # Scan mode numbers
453 print("Scanning mode numbers (m, n)...")
454 W_map, m_values, n_values = scan_mode_numbers(ep, B_z, B_theta, p,
455 m_max=5, n_max=5)
456
457 # Plot stability map
458 plot_stability_map(W_map, m_values, n_values)
459 print("Stability map saved as 'stability_map.png'")
460
461
462if __name__ == '__main__':
463 main()