1#!/usr/bin/env python3
2"""
3Solar Flux Tube Buoyancy Model
4
5This module simulates the rise of a magnetic flux tube through the solar
6convection zone due to magnetic buoyancy.
7
8Physical principles:
9- Magnetic pressure reduces internal gas pressure and density
10- Buoyant force: F_b = (ρ_ext - ρ_int)g × Volume
11- Drag force: F_d = -C_D × ρ_ext × v² × Area
12- Flux tube expands as it rises (pressure balance)
13
14The model tracks:
15- Trajectory of flux tube from base of convection zone to surface
16- Rise velocity evolution
17- Magnetic field strength variation with height
18- Comparison with observed solar active region emergence
19
20Author: MHD Course Examples
21License: MIT
22"""
23
24import numpy as np
25import matplotlib.pyplot as plt
26from scipy.integrate import odeint
27
28
29class SolarFluxTube:
30 """
31 Thin flux tube model in stratified solar atmosphere.
32
33 Attributes:
34 R_sun (float): Solar radius (m)
35 g_sun (float): Solar surface gravity (m/s²)
36 depth_base (float): Depth of convection zone base (m)
37 T_base (float): Temperature at base (K)
38 rho_base (float): Density at base (kg/m³)
39 B0 (float): Initial magnetic field strength (T)
40 tube_radius (float): Initial flux tube radius (m)
41 """
42
43 def __init__(self, R_sun=6.96e8, g_sun=274.0,
44 depth_base=2.0e8, T_base=2e6, rho_base=200.0,
45 B0=1e5, tube_radius=1e7):
46 """
47 Initialize solar flux tube model.
48
49 Parameters:
50 R_sun (float): Solar radius (m)
51 g_sun (float): Surface gravity (m/s²)
52 depth_base (float): Convection zone depth (m)
53 T_base (float): Temperature at base (K)
54 rho_base (float): Density at base (kg/m³)
55 B0 (float): Initial field strength (Gauss → Tesla)
56 tube_radius (float): Initial tube radius (m)
57 """
58 self.R_sun = R_sun
59 self.g_sun = g_sun
60 self.depth_base = depth_base
61 self.T_base = T_base
62 self.rho_base = rho_base
63 self.B0 = B0 * 1e-4 # Convert Gauss to Tesla
64 self.tube_radius = tube_radius
65
66 # Gas constant
67 self.mu = 0.6 # Mean molecular weight
68 self.m_p = 1.67e-27 # Proton mass (kg)
69 self.k_B = 1.38e-23 # Boltzmann constant (J/K)
70 self.R_gas = self.k_B / (self.mu * self.m_p)
71
72 # Stratification scale height
73 self.H = self.R_gas * self.T_base / self.g_sun
74
75 # Drag coefficient
76 self.C_D = 1.0
77
78 # History
79 self.time_history = []
80 self.height_history = []
81 self.velocity_history = []
82 self.B_history = []
83 self.radius_history = []
84
85 def external_density(self, z):
86 """
87 External atmospheric density as function of height.
88
89 Assumes exponential stratification: ρ(z) = ρ₀ exp(-z/H)
90
91 Parameters:
92 z (float): Height above base (m)
93
94 Returns:
95 float: Density (kg/m³)
96 """
97 return self.rho_base * np.exp(-z / self.H)
98
99 def external_pressure(self, z):
100 """
101 External atmospheric pressure.
102
103 Parameters:
104 z (float): Height above base (m)
105
106 Returns:
107 float: Pressure (Pa)
108 """
109 rho = self.external_density(z)
110 T = self.T_base * np.exp(-z / (2 * self.H)) # Simplified temp profile
111 return rho * self.R_gas * T
112
113 def internal_density(self, z, B):
114 """
115 Internal flux tube density from pressure balance.
116
117 Total pressure balance: p_int + B²/(2μ₀) = p_ext
118
119 Parameters:
120 z (float): Height (m)
121 B (float): Magnetic field (T)
122
123 Returns:
124 float: Internal density (kg/m³)
125 """
126 mu_0 = 4 * np.pi * 1e-7 # Permeability
127 p_ext = self.external_pressure(z)
128 p_mag = B**2 / (2 * mu_0)
129
130 # Internal gas pressure
131 p_int = p_ext - p_mag
132
133 if p_int <= 0:
134 # Flux tube has expanded to very low density
135 return 1e-3
136
137 # Assuming same temperature inside and outside
138 T = self.T_base * np.exp(-z / (2 * self.H))
139 rho_int = p_int / (self.R_gas * T)
140
141 return max(rho_int, 1e-3)
142
143 def buoyancy_force(self, z, B):
144 """
145 Compute buoyancy force per unit mass.
146
147 Parameters:
148 z (float): Height (m)
149 B (float): Magnetic field (T)
150
151 Returns:
152 float: Buoyancy acceleration (m/s²)
153 """
154 rho_ext = self.external_density(z)
155 rho_int = self.internal_density(z, B)
156
157 # Buoyancy acceleration
158 g_eff = self.g_sun * (1 - rho_int / rho_ext)
159
160 return g_eff
161
162 def drag_force(self, z, v):
163 """
164 Compute drag force per unit mass.
165
166 Parameters:
167 z (float): Height (m)
168 v (float): Velocity (m/s)
169
170 Returns:
171 float: Drag acceleration (m/s²)
172 """
173 rho_ext = self.external_density(z)
174
175 # Drag per unit mass (assumes cylindrical tube)
176 # F_d/m_tube = -C_D × (ρ_ext/ρ_int) × v²
177 # Simplified: proportional to velocity
178 drag_coeff = self.C_D * rho_ext / self.rho_base
179 a_drag = -drag_coeff * v
180
181 return a_drag
182
183 def magnetic_field_evolution(self, z, B):
184 """
185 Magnetic field evolves due to tube expansion.
186
187 Flux conservation: B × A = constant
188 For thin tube: B ∝ ρ_int (from pressure balance)
189
190 Parameters:
191 z (float): Height (m)
192 B (float): Current field (T)
193
194 Returns:
195 float: Updated field (T)
196 """
197 rho_int = self.internal_density(z, B)
198 B_new = self.B0 * np.sqrt(rho_int / self.rho_base)
199
200 return B_new
201
202 def rhs(self, state, t):
203 """
204 Right-hand side for flux tube motion.
205
206 state = [z, v]
207 dz/dt = v
208 dv/dt = g_buoyancy + a_drag
209
210 Parameters:
211 state (list): [height, velocity]
212 t (float): Time
213
214 Returns:
215 list: Time derivatives
216 """
217 z, v = state
218
219 # Current magnetic field
220 B = self.magnetic_field_evolution(z, self.B0)
221
222 # Forces
223 a_buoy = self.buoyancy_force(z, B)
224 a_drag = self.drag_force(z, v)
225
226 dz_dt = v
227 dv_dt = a_buoy + a_drag
228
229 return [dz_dt, dv_dt]
230
231 def run_simulation(self, t_end=1e6):
232 """
233 Simulate flux tube rise.
234
235 Parameters:
236 t_end (float): End time (s)
237 """
238 print("Simulating solar flux tube rise...")
239 print(f"Initial depth: {self.depth_base/1e8:.1f} × 10^8 m")
240 print(f"Initial field: {self.B0*1e4:.0f} Gauss")
241 print(f"Scale height: {self.H/1e8:.2f} × 10^8 m")
242
243 # Initial state: at base, zero velocity
244 state0 = [0.0, 0.0]
245
246 # Time array
247 t_array = np.linspace(0, t_end, 1000)
248
249 # Integrate
250 solution = odeint(self.rhs, state0, t_array)
251
252 self.time_history = t_array
253 self.height_history = solution[:, 0]
254 self.velocity_history = solution[:, 1]
255
256 # Compute B and radius history
257 self.B_history = []
258 self.radius_history = []
259
260 for z in self.height_history:
261 B = self.magnetic_field_evolution(z, self.B0)
262 self.B_history.append(B)
263
264 # Radius from flux conservation
265 rho_int = self.internal_density(z, B)
266 r = self.tube_radius * np.sqrt(self.rho_base / rho_int)
267 self.radius_history.append(r)
268
269 self.B_history = np.array(self.B_history)
270 self.radius_history = np.array(self.radius_history)
271
272 # Find emergence time
273 emergence_idx = np.argmax(self.height_history >= self.depth_base)
274 if emergence_idx > 0:
275 t_emerge = self.time_history[emergence_idx]
276 print(f"\nFlux tube emerges at t = {t_emerge/86400:.1f} days")
277 print(f"Emergence velocity: {self.velocity_history[emergence_idx]:.1f} m/s")
278 print(f"Surface field: {self.B_history[emergence_idx]*1e4:.0f} Gauss")
279
280 def plot_trajectory(self):
281 """
282 Plot flux tube trajectory and properties.
283 """
284 fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))
285
286 # Convert to convenient units
287 t_days = self.time_history / 86400
288 z_Mm = self.height_history / 1e6
289 B_gauss = self.B_history * 1e4
290 r_km = self.radius_history / 1e3
291
292 # Height vs time
293 ax1.plot(t_days, z_Mm, 'b-', linewidth=2)
294 ax1.axhline(y=self.depth_base/1e6, color='r', linestyle='--',
295 label='Surface', linewidth=2)
296 ax1.set_xlabel('Time (days)', fontsize=11)
297 ax1.set_ylabel('Height (Mm)', fontsize=11)
298 ax1.set_title('Flux Tube Rise Trajectory', fontsize=13)
299 ax1.grid(True, alpha=0.3)
300 ax1.legend(fontsize=10)
301
302 # Velocity vs height
303 ax2.plot(z_Mm, self.velocity_history, 'g-', linewidth=2)
304 ax2.set_xlabel('Height (Mm)', fontsize=11)
305 ax2.set_ylabel('Velocity (m/s)', fontsize=11)
306 ax2.set_title('Rise Velocity vs Height', fontsize=13)
307 ax2.grid(True, alpha=0.3)
308
309 # Magnetic field vs height
310 ax3.plot(z_Mm, B_gauss, 'r-', linewidth=2)
311 ax3.set_xlabel('Height (Mm)', fontsize=11)
312 ax3.set_ylabel('Magnetic Field (Gauss)', fontsize=11)
313 ax3.set_title('Magnetic Field Evolution', fontsize=13)
314 ax3.grid(True, alpha=0.3)
315
316 # Tube radius vs height
317 ax4.plot(z_Mm, r_km, 'm-', linewidth=2)
318 ax4.set_xlabel('Height (Mm)', fontsize=11)
319 ax4.set_ylabel('Tube Radius (km)', fontsize=11)
320 ax4.set_title('Flux Tube Expansion', fontsize=13)
321 ax4.grid(True, alpha=0.3)
322
323 plt.tight_layout()
324 return fig
325
326
327def main():
328 """
329 Main function demonstrating solar flux tube rise.
330 """
331 print("=" * 60)
332 print("Solar Flux Tube Buoyancy Simulation")
333 print("=" * 60)
334
335 # Create flux tube model
336 tube = SolarFluxTube(
337 R_sun=6.96e8,
338 g_sun=274.0,
339 depth_base=2.0e8, # 200 Mm convection zone
340 T_base=2e6, # 2 MK at base
341 rho_base=200.0, # kg/m³
342 B0=1e5, # 100,000 Gauss
343 tube_radius=1e7 # 10,000 km
344 )
345
346 # Run simulation (about 30 days)
347 tube.run_simulation(t_end=3e6)
348
349 # Plot results
350 tube.plot_trajectory()
351
352 plt.savefig('/tmp/solar_flux_tube.png', dpi=150, bbox_inches='tight')
353 print("\nPlot saved to /tmp/solar_flux_tube.png")
354
355 plt.show()
356
357
358if __name__ == "__main__":
359 main()