1#!/usr/bin/env python3
2"""
3Ponomarenko Dynamo
4
5This module implements the Ponomarenko dynamo, which demonstrates that a simple
6helical flow in a cylinder can generate a self-sustaining magnetic field.
7
8The flow configuration is:
9 v_z = V (for r < a) - axial flow
10 v_θ = Ωr (for r < a) - azimuthal rotation
11 v = 0 (for r > a) - no flow outside
12
13This creates a helical motion that can amplify magnetic fields when the
14magnetic Reynolds number Rm = Va/η exceeds a critical value.
15
16The problem is formulated as an eigenvalue problem for the growth rate γ.
17
18Author: MHD Course Examples
19License: MIT
20"""
21
22import numpy as np
23import matplotlib.pyplot as plt
24from scipy.special import jv, yv, iv, kv # Bessel functions
25from scipy.optimize import fsolve, brentq
26
27
28class PonomarenkoDynamo:
29 """
30 Ponomarenko dynamo eigenvalue solver.
31
32 Attributes:
33 a (float): Cylinder radius
34 V (float): Axial velocity
35 Omega (float): Angular velocity
36 eta (float): Magnetic diffusivity
37 k (float): Axial wavenumber
38 m (int): Azimuthal mode number
39 """
40
41 def __init__(self, a=1.0, V=1.0, Omega=1.0, eta=0.1, k=1.0, m=1):
42 """
43 Initialize Ponomarenko dynamo parameters.
44
45 Parameters:
46 a (float): Cylinder radius
47 V (float): Axial velocity
48 Omega (float): Angular velocity
49 eta (float): Magnetic diffusivity
50 k (float): Axial wavenumber
51 m (int): Azimuthal mode number
52 """
53 self.a = a
54 self.V = V
55 self.Omega = Omega
56 self.eta = eta
57 self.k = k
58 self.m = m
59
60 # Magnetic Reynolds numbers
61 self.Rm_axial = V * a / eta
62 self.Rm_rot = Omega * a**2 / eta
63 self.Rm = np.sqrt(self.Rm_axial**2 + self.Rm_rot**2)
64
65 def dispersion_relation(self, gamma_guess, Rm):
66 """
67 Dispersion relation for Ponomarenko dynamo.
68
69 This is a simplified version. The full problem requires matching
70 boundary conditions at r=a.
71
72 Parameters:
73 gamma_guess (float): Growth rate guess
74 Rm (float): Magnetic Reynolds number
75
76 Returns:
77 complex: Residual of dispersion relation
78 """
79 # Simplified dispersion relation (approximate)
80 # Full version requires numerical solution of Bessel function equation
81 k = self.k
82 m = self.m
83
84 # Effective wavenumber including growth
85 q_squared = k**2 + 1j * gamma_guess / self.eta
86
87 # For simplicity, use approximate relation
88 # Real Ponomarenko requires matching at r=a
89 lhs = q_squared * self.a**2
90 rhs = Rm * k * self.a # Coupling term
91
92 return abs(lhs - rhs**2) - 1.0
93
94 def find_critical_Rm(self, k_value):
95 """
96 Find critical magnetic Reynolds number for given wavenumber.
97
98 Parameters:
99 k_value (float): Axial wavenumber
100
101 Returns:
102 float: Critical Rm
103 """
104 self.k = k_value
105
106 # Critical Rm is where growth rate becomes positive
107 # For Ponomarenko, Rm_c ≈ 17.7 for optimal k and m=1
108 def growth_rate_sign(Rm):
109 # Approximate growth rate
110 gamma_approx = self.eta * (Rm * self.k * self.a / self.a**2 - self.k**2)
111 return gamma_approx
112
113 # Search for Rm where gamma = 0
114 Rm_range = np.linspace(10, 30, 100)
115 growth_rates = [growth_rate_sign(Rm) for Rm in Rm_range]
116
117 # Find zero crossing
118 for i in range(len(growth_rates) - 1):
119 if growth_rates[i] < 0 and growth_rates[i+1] > 0:
120 return Rm_range[i]
121
122 return 17.7 # Theoretical value for m=1
123
124 def compute_growth_rate(self, Rm_value):
125 """
126 Compute growth rate for given magnetic Reynolds number.
127
128 Parameters:
129 Rm_value (float): Magnetic Reynolds number
130
131 Returns:
132 float: Growth rate γ
133 """
134 # Approximate growth rate formula
135 # Derived from perturbation analysis
136 k = self.k
137 a = self.a
138 eta = self.eta
139
140 # Growth rate scales as Rm above critical
141 Rm_c = 17.7
142
143 if Rm_value < Rm_c:
144 gamma = -eta * k**2 # Decay
145 else:
146 # Approximate linear growth above threshold
147 gamma = eta * (k * a) * (Rm_value / Rm_c - 1.0)
148
149 return gamma
150
151 def compute_eigenfunction(self, r_array, Rm_value):
152 """
153 Compute magnetic field eigenfunction B(r).
154
155 Parameters:
156 r_array (ndarray): Radial positions
157 Rm_value (float): Magnetic Reynolds number
158
159 Returns:
160 tuple: (B_r, B_theta, B_z) components
161 """
162 k = self.k
163 m = self.m
164 a = self.a
165
166 gamma = self.compute_growth_rate(Rm_value)
167 q = np.sqrt(k**2 + 1j * gamma / self.eta)
168
169 B_r = np.zeros_like(r_array, dtype=complex)
170 B_theta = np.zeros_like(r_array, dtype=complex)
171 B_z = np.zeros_like(r_array, dtype=complex)
172
173 for i, r in enumerate(r_array):
174 if r < a:
175 # Inside cylinder: Bessel J function
176 arg = q * r
177 B_r[i] = jv(m, arg)
178 B_theta[i] = 1j * m * jv(m, arg) / (q * r) if r > 0 else 0
179 B_z[i] = k * jv(m, arg) / q
180 else:
181 # Outside cylinder: Modified Bessel K function (decay)
182 arg = k * (r - a)
183 B_r[i] = jv(m, q * a) * np.exp(-arg)
184 B_theta[i] = 0
185 B_z[i] = k * jv(m, q * a) * np.exp(-arg) / q
186
187 # Normalize
188 B_max = np.max(np.abs(B_z))
189 if B_max > 0:
190 B_r /= B_max
191 B_theta /= B_max
192 B_z /= B_max
193
194 return np.real(B_r), np.real(B_theta), np.real(B_z)
195
196 def plot_growth_rate_vs_Rm(self):
197 """
198 Plot growth rate as a function of magnetic Reynolds number.
199 """
200 Rm_array = np.linspace(10, 30, 100)
201 gamma_array = [self.compute_growth_rate(Rm) for Rm in Rm_array]
202
203 fig, ax = plt.subplots(figsize=(10, 6))
204
205 ax.plot(Rm_array, gamma_array, 'b-', linewidth=2, label='Growth rate γ')
206 ax.axhline(y=0, color='k', linestyle='--', alpha=0.5)
207 ax.axvline(x=17.7, color='r', linestyle='--', alpha=0.5,
208 label=r'$Rm_c \approx 17.7$')
209
210 ax.set_xlabel('Magnetic Reynolds Number Rm', fontsize=12)
211 ax.set_ylabel('Growth Rate γ', fontsize=12)
212 ax.set_title('Ponomarenko Dynamo: Growth Rate vs Rm\n(m=1 mode)',
213 fontsize=14)
214 ax.grid(True, alpha=0.3)
215 ax.legend(fontsize=11)
216
217 # Mark regions
218 ax.fill_between(Rm_array, -1, 1, where=(Rm_array < 17.7),
219 alpha=0.2, color='red', label='Stable (decay)')
220 ax.fill_between(Rm_array, -1, 1, where=(Rm_array >= 17.7),
221 alpha=0.2, color='green', label='Unstable (growth)')
222
223 ax.set_ylim([-0.5, 0.5])
224
225 plt.tight_layout()
226 return fig
227
228 def plot_eigenfunction(self, Rm_value):
229 """
230 Plot magnetic field eigenfunction B(r).
231
232 Parameters:
233 Rm_value (float): Magnetic Reynolds number
234 """
235 r_array = np.linspace(0, 2 * self.a, 200)
236 B_r, B_theta, B_z = self.compute_eigenfunction(r_array, Rm_value)
237
238 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
239
240 # Field components
241 ax1.plot(r_array / self.a, B_r, 'b-', label=r'$B_r$', linewidth=2)
242 ax1.plot(r_array / self.a, B_theta, 'g-', label=r'$B_\theta$', linewidth=2)
243 ax1.plot(r_array / self.a, B_z, 'r-', label=r'$B_z$', linewidth=2)
244 ax1.axvline(x=1.0, color='k', linestyle='--', alpha=0.5,
245 label='Cylinder boundary')
246 ax1.set_xlabel('Radius r/a', fontsize=12)
247 ax1.set_ylabel('Magnetic Field (normalized)', fontsize=12)
248 ax1.set_title(f'Ponomarenko Dynamo Eigenfunction (Rm = {Rm_value:.1f})',
249 fontsize=14)
250 ax1.grid(True, alpha=0.3)
251 ax1.legend(fontsize=11)
252
253 # Total field magnitude
254 B_total = np.sqrt(B_r**2 + B_theta**2 + B_z**2)
255 ax2.plot(r_array / self.a, B_total, 'k-', linewidth=2, label='|B|')
256 ax2.axvline(x=1.0, color='k', linestyle='--', alpha=0.5,
257 label='Cylinder boundary')
258 ax2.fill_between(r_array / self.a, 0, B_total,
259 where=(r_array <= self.a), alpha=0.3, color='blue',
260 label='Inside (helical flow)')
261 ax2.set_xlabel('Radius r/a', fontsize=12)
262 ax2.set_ylabel('|B|', fontsize=12)
263 ax2.set_title('Total Magnetic Field Magnitude', fontsize=14)
264 ax2.grid(True, alpha=0.3)
265 ax2.legend(fontsize=11)
266
267 plt.tight_layout()
268 return fig
269
270
271def main():
272 """
273 Main function demonstrating the Ponomarenko dynamo.
274 """
275 print("=" * 60)
276 print("Ponomarenko Dynamo Simulation")
277 print("=" * 60)
278
279 # Create dynamo instance
280 dynamo = PonomarenkoDynamo(
281 a=1.0,
282 V=1.0,
283 Omega=1.0,
284 eta=0.05,
285 k=2.0,
286 m=1
287 )
288
289 print(f"\nParameters:")
290 print(f" Cylinder radius a = {dynamo.a}")
291 print(f" Axial velocity V = {dynamo.V}")
292 print(f" Angular velocity Ω = {dynamo.Omega}")
293 print(f" Magnetic diffusivity η = {dynamo.eta}")
294 print(f" Axial wavenumber k = {dynamo.k}")
295 print(f" Azimuthal mode m = {dynamo.m}")
296 print(f"\nMagnetic Reynolds numbers:")
297 print(f" Rm_axial = Va/η = {dynamo.Rm_axial:.2f}")
298 print(f" Rm_rot = Ωa²/η = {dynamo.Rm_rot:.2f}")
299 print(f" Total Rm = {dynamo.Rm:.2f}")
300 print(f"\nCritical Rm ≈ 17.7 for m=1 mode")
301
302 # Plot growth rate vs Rm
303 fig1 = dynamo.plot_growth_rate_vs_Rm()
304
305 # Plot eigenfunction for Rm above critical
306 fig2 = dynamo.plot_eigenfunction(Rm_value=20.0)
307
308 plt.savefig('/tmp/ponomarenko_dynamo.png', dpi=150, bbox_inches='tight')
309 print("\nPlots saved to /tmp/ponomarenko_dynamo.png")
310
311 plt.show()
312
313
314if __name__ == "__main__":
315 main()