1#!/usr/bin/env python3
2"""
3Z-Pinch Equilibrium Analysis
4
5This script solves the Z-pinch equilibrium equation:
6 p'(r) = -J_z(r) * B_θ(r)
7
8for a given current profile and verifies the Bennett relation:
9 I² = 8π N k T / μ₀
10
11The Z-pinch is a cylindrical plasma configuration where an axial current
12creates an azimuthal magnetic field that confines the plasma.
13
14Author: Claude
15Date: 2026-02-14
16"""
17
18import numpy as np
19import matplotlib.pyplot as plt
20from scipy.integrate import odeint, simpson
21from scipy.constants import mu_0, k as k_B, e
22
23# Physical constants
24MU_0 = mu_0 # Permeability of free space
25
26
27def current_density_profile(r, r_max, I_total, profile_type='uniform'):
28 """
29 Define current density profile J_z(r).
30
31 Parameters
32 ----------
33 r : array_like
34 Radial positions
35 r_max : float
36 Maximum radius of the pinch
37 I_total : float
38 Total current (A)
39 profile_type : str
40 Type of profile: 'uniform', 'parabolic', or 'skin'
41
42 Returns
43 -------
44 J_z : ndarray
45 Current density at each radius
46 """
47 if profile_type == 'uniform':
48 # Uniform current density
49 J_0 = I_total / (np.pi * r_max**2)
50 J_z = np.where(r <= r_max, J_0, 0.0)
51
52 elif profile_type == 'parabolic':
53 # Parabolic profile: J ~ (1 - (r/r_max)²)
54 J_0 = 2 * I_total / (np.pi * r_max**2)
55 J_z = np.where(r <= r_max, J_0 * (1 - (r/r_max)**2), 0.0)
56
57 elif profile_type == 'skin':
58 # Skin current: J ~ exp(-(r/δ)²)
59 delta = r_max / 3 # Skin depth
60 J_0 = I_total / (np.pi * delta**2 * (1 - np.exp(-r_max**2/delta**2)))
61 J_z = J_0 * np.exp(-(r/delta)**2)
62
63 else:
64 raise ValueError(f"Unknown profile type: {profile_type}")
65
66 return J_z
67
68
69def azimuthal_field(r, r_max, I_total, profile_type='uniform'):
70 """
71 Compute azimuthal magnetic field B_θ(r) from Ampere's law.
72
73 B_θ(r) = μ₀ I(r) / (2π r)
74 where I(r) is the current enclosed within radius r.
75
76 Parameters
77 ----------
78 r : array_like
79 Radial positions
80 r_max : float
81 Maximum radius
82 I_total : float
83 Total current
84 profile_type : str
85 Current profile type
86
87 Returns
88 -------
89 B_theta : ndarray
90 Azimuthal magnetic field
91 """
92 r = np.asarray(r)
93 B_theta = np.zeros_like(r)
94
95 for i, ri in enumerate(r):
96 if ri < 1e-10:
97 B_theta[i] = 0.0
98 else:
99 # Integrate current density to get enclosed current
100 r_int = np.linspace(0, ri, 100)
101 J_int = current_density_profile(r_int, r_max, I_total, profile_type)
102 I_enclosed = 2 * np.pi * simpson(r_int * J_int, x=r_int)
103 B_theta[i] = MU_0 * I_enclosed / (2 * np.pi * ri)
104
105 return B_theta
106
107
108def solve_pressure_balance(r, r_max, I_total, p_edge, profile_type='uniform'):
109 """
110 Solve pressure balance equation: dp/dr = -J_z * B_θ
111
112 Parameters
113 ----------
114 r : array_like
115 Radial grid
116 r_max : float
117 Pinch radius
118 I_total : float
119 Total current
120 p_edge : float
121 Pressure at edge
122 profile_type : str
123 Current profile type
124
125 Returns
126 -------
127 pressure : ndarray
128 Pressure profile
129 """
130 J_z = current_density_profile(r, r_max, I_total, profile_type)
131 B_theta = azimuthal_field(r, r_max, I_total, profile_type)
132
133 # Integrate from edge inward
134 pressure = np.zeros_like(r)
135 pressure[-1] = p_edge
136
137 for i in range(len(r)-2, -1, -1):
138 dr = r[i+1] - r[i]
139 # dp/dr = -J_z * B_θ, integrate backward
140 pressure[i] = pressure[i+1] - J_z[i] * B_theta[i] * dr
141
142 return pressure
143
144
145def safety_factor(r, r_max, I_total, profile_type='uniform'):
146 """
147 Compute safety factor q(r) = r * B_z / (R * B_θ)
148
149 For a Z-pinch with small B_z (from external field),
150 we assume B_z = constant for this calculation.
151
152 Parameters
153 ----------
154 r : array_like
155 Radial positions
156 r_max : float
157 Pinch radius
158 I_total : float
159 Total current
160 profile_type : str
161 Current profile type
162
163 Returns
164 -------
165 q : ndarray
166 Safety factor
167 """
168 B_z = 0.1 # Assume small axial field (Tesla)
169 B_theta = azimuthal_field(r, r_max, I_total, profile_type)
170 R_major = 1.0 # Assume major radius = 1 m for calculation
171
172 # Avoid division by zero
173 q = np.where(np.abs(B_theta) > 1e-10,
174 r * B_z / (R_major * B_theta),
175 np.inf)
176 return q
177
178
179def bennett_relation(I_total, n_avg, T_avg):
180 """
181 Verify Bennett relation: I² = 8π N k T / μ₀
182
183 where N is the line density (particles per unit length).
184
185 Parameters
186 ----------
187 I_total : float
188 Total current (A)
189 n_avg : float
190 Average particle density (m⁻³)
191 T_avg : float
192 Average temperature (eV)
193
194 Returns
195 -------
196 dict
197 Bennett relation verification results
198 """
199 # Convert temperature to Joules
200 T_J = T_avg * e
201
202 # For a pinch of radius a, N ~ n * π * a²
203 # We'll compute the expected current from Bennett relation
204 a = 0.01 # Assume 1 cm radius
205 N_line = n_avg * np.pi * a**2
206
207 I_bennett = np.sqrt(8 * np.pi * N_line * k_B * T_J / MU_0)
208
209 return {
210 'I_actual': I_total,
211 'I_bennett': I_bennett,
212 'ratio': I_total / I_bennett,
213 'N_line': N_line
214 }
215
216
217def plot_equilibrium(r, pressure, B_theta, J_z, q, profile_type):
218 """
219 Plot Z-pinch equilibrium profiles.
220
221 Parameters
222 ----------
223 r : ndarray
224 Radial grid
225 pressure : ndarray
226 Pressure profile
227 B_theta : ndarray
228 Azimuthal field
229 J_z : ndarray
230 Current density
231 q : ndarray
232 Safety factor
233 profile_type : str
234 Current profile type
235 """
236 fig, axes = plt.subplots(2, 2, figsize=(12, 10))
237
238 # Pressure
239 axes[0, 0].plot(r * 100, pressure / 1e3, 'b-', linewidth=2)
240 axes[0, 0].set_xlabel('Radius (cm)', fontsize=12)
241 axes[0, 0].set_ylabel('Pressure (kPa)', fontsize=12)
242 axes[0, 0].set_title('Pressure Profile', fontsize=13, fontweight='bold')
243 axes[0, 0].grid(True, alpha=0.3)
244
245 # Azimuthal field
246 axes[0, 1].plot(r * 100, B_theta * 1e3, 'r-', linewidth=2)
247 axes[0, 1].set_xlabel('Radius (cm)', fontsize=12)
248 axes[0, 1].set_ylabel('B_θ (mT)', fontsize=12)
249 axes[0, 1].set_title('Azimuthal Magnetic Field', fontsize=13, fontweight='bold')
250 axes[0, 1].grid(True, alpha=0.3)
251
252 # Current density
253 axes[1, 0].plot(r * 100, J_z / 1e6, 'g-', linewidth=2)
254 axes[1, 0].set_xlabel('Radius (cm)', fontsize=12)
255 axes[1, 0].set_ylabel('J_z (MA/m²)', fontsize=12)
256 axes[1, 0].set_title('Current Density', fontsize=13, fontweight='bold')
257 axes[1, 0].grid(True, alpha=0.3)
258
259 # Safety factor (plot only finite values)
260 q_plot = np.where(np.isfinite(q) & (q < 100), q, np.nan)
261 axes[1, 1].plot(r * 100, q_plot, 'm-', linewidth=2)
262 axes[1, 1].set_xlabel('Radius (cm)', fontsize=12)
263 axes[1, 1].set_ylabel('Safety factor q', fontsize=12)
264 axes[1, 1].set_title('Safety Factor', fontsize=13, fontweight='bold')
265 axes[1, 1].grid(True, alpha=0.3)
266
267 plt.suptitle(f'Z-Pinch Equilibrium ({profile_type.capitalize()} Current Profile)',
268 fontsize=14, fontweight='bold')
269 plt.tight_layout()
270 plt.savefig(f'z_pinch_{profile_type}.png', dpi=150, bbox_inches='tight')
271 plt.show()
272
273
274def main():
275 """Main execution function."""
276 # Parameters
277 r_max = 0.01 # Pinch radius: 1 cm
278 I_total = 100e3 # Total current: 100 kA
279 p_edge = 100.0 # Edge pressure: 100 Pa
280 profile_type = 'parabolic' # 'uniform', 'parabolic', or 'skin'
281
282 # Plasma parameters for Bennett relation
283 n_avg = 1e20 # Average density: 10²⁰ m⁻³
284 T_avg = 100.0 # Average temperature: 100 eV
285
286 # Create radial grid
287 r = np.linspace(0, r_max, 200)
288
289 # Compute profiles
290 print(f"Computing Z-pinch equilibrium with {profile_type} current profile...")
291 J_z = current_density_profile(r, r_max, I_total, profile_type)
292 B_theta = azimuthal_field(r, r_max, I_total, profile_type)
293 pressure = solve_pressure_balance(r, r_max, I_total, p_edge, profile_type)
294 q = safety_factor(r, r_max, I_total, profile_type)
295
296 # Print key results
297 print(f"\nKey Parameters:")
298 print(f" Pinch radius: {r_max*100:.2f} cm")
299 print(f" Total current: {I_total/1e3:.1f} kA")
300 print(f" Peak current density: {np.max(J_z)/1e6:.2f} MA/m²")
301 print(f" Peak B_θ: {np.max(B_theta)*1e3:.2f} mT")
302 print(f" Peak pressure: {np.max(pressure)/1e3:.2f} kPa")
303 print(f" Central pressure: {pressure[0]/1e3:.2f} kPa")
304
305 # Bennett relation
306 bennett = bennett_relation(I_total, n_avg, T_avg)
307 print(f"\nBennett Relation Verification:")
308 print(f" Actual current: {bennett['I_actual']/1e3:.1f} kA")
309 print(f" Bennett current: {bennett['I_bennett']/1e3:.1f} kA")
310 print(f" Ratio I_actual/I_bennett: {bennett['ratio']:.3f}")
311 print(f" Line density N: {bennett['N_line']:.2e} m⁻¹")
312
313 if np.abs(bennett['ratio'] - 1.0) < 0.2:
314 print(" ✓ Bennett relation approximately satisfied")
315 else:
316 print(" ✗ Bennett relation not satisfied (adjust n or T)")
317
318 # Plot results
319 plot_equilibrium(r, pressure, B_theta, J_z, q, profile_type)
320 print(f"\nPlot saved as 'z_pinch_{profile_type}.png'")
321
322
323if __name__ == '__main__':
324 main()