A&A 451, 303-317 (2006)
DOI: 10.1051/0004-6361:20054499
M. C. M. Cheung1 - F. Moreno-Insertis2,3 - M. Schüssler1
1 - Max Planck Institute for Solar System Research, 37191
Katlenburg-Lindau, Germany
2 - Institúto de Astrofísica de Canarias, 38200 La Laguna (Tenerife), Spain
3 - Dept. of Astrophysics, Faculty of Physics, University of La
Laguna, 38200 La Laguna (Tenerife), Spain
Received 9 November 2005 / Accepted 13 January 2006
Abstract
Aims. We study the buoyant rise of magnetic flux tubes in a stratified layer over a range of Reynolds numbers (
)
by means of numerical simulations. Special emphasis is placed on studying the fragmentation of the rising tube, its trailing wake and the formation of a vortex street in the high-Reynolds number regime. Furthermore, we evaluate the relevance of the thin flux tube approximation with regard to describing the evolution of magnetic flux tubes in the simulations.
Methods. We used the FLASH code, which has an adaptive mesh refinement (AMR) algorithm, thus allowing the simulations to be carried out at high Reynolds numbers.
Results. The evolution of the magnetic flux tube and its wake depends on the Reynolds number. At Re up to a few hundred, the wake consists of two counter-rotating vortex rolls. At higher Re, the vortex rolls break up and the shedding of flux into the wake occurs in a more intermittent fashion. The amount of flux retained by the central portion of the tube increases with the field line twist (in agreement with previous literature) and with Re. The time evolution of the twist is compatible with a homologous expansion of the tube. The motion of the central portion of the tube in the simulations is very well described by the thin flux tube model whenever the effects of flux loss or vortex forces can be neglected. If the flux tube has an initial net vorticity, it undergoes asymmetric vortex shedding. In this case, the lift force accelerates the tube in such a way that an oscillatory horizontal motion is super-imposed on the vertical rise of the tube, which leaves behind a vortex street. This last result is in accordance with previous simulations reported in the literature, which were carried out at lower Reynolds number.
Key words: magnetohydrodynamics (MHD) - Sun: magnetic fields - Sun: interior
A branch in this undertaking is devoted to the basic magnetohydrodynamics of buoyant flux tubes rising in stratified and (otherwise) unmagnetized media studied by way of 2D or 3D numerical experiments (see review by Fan 2004). One example is an initially horizontal magnetic flux tube embedded in a stratified layer. The tube is endowed with a density deficit with respect to the surroundings, so that it rises and, in doing so, it expands, displaces the surrounding medium and develops a trailing wake. A number of results, obtained mostly for 2D (more precisely, 2.5D) configurations, concern the conversion of the rising magnetic tube into pairs of vortex tubes. These experiments focus on the evolution in a vertical plane normal to the axis of the tube and use the simplifying assumption of independence of all quantities (scalars or vectors) with respect to the coordinate along that axis. When the magnetic field in the horizontal tube is purely longitudinal (i.e., has no components in the plane normal to the tube axis), an initially cylindrical tube, after rising a height equivalent to a few times the tube diameter, splits into two roughly mirror-symmetric vortex rolls. The rolls have vorticity of opposite signs pointing in the direction of the tube axis and separate horizontally from each other. This behavior had been noted in an early paper by Schüssler (1979) and was analyzed by Longcope et al. (1996) who showed that the motion of the resulting vortex tubes could be explained as a result of the combined action of the buoyancy and lift forces on them. On the other hand, Moreno-Insertis & Emonet (1996) demonstrated that a sufficient amount of field line twist around the tube axis can prevent the splitting of the tube: the transverse component of the field (i.e., the component in the normal plane to the tube axis) imparts rigidity to the tube interior against the stresses in their periphery. They could quantify this effect both by deriving approximate criteria and by considering the results of numerical experiments of rising tubes with different levels of twist. Moving magnetic tubes have a trailing wake with two vortex rolls in a wide range of Reynolds numbers, including those relevant for these calculations. Moreno-Insertis & Emonet (1996) showed that the rising tube loses magnetic flux to the wake by an amount that is higher the smaller the pitch angle of the field lines around the tube axis. For an untwisted tube, the major part of the initial magnetic flux is just lost to the wake, a phenomenon that was described by the earlier authors as the conversion of the tube into two vortex rolls. Further to that, the stresses exerted on the periphery of the tube can be integrated to obtain global values for the drag and lift forces on the tube and for the added inertia of the tube caused by the co-acceleration of the surrounding plasma. Those authors showed that the motion of the tube in their experiments could be fitted reasonably using simple estimates for the drag force and enhanced inertia of the tube. Later on, Emonet & Moreno-Insertis (1998) studied the dynamics within the vortex tube and in its surroundings and, in particular, the central role played by the magnetic and viscous boundary layers that serve as interface between them. The boundary layer is the site of generation of the vorticity that accumulates in the wake; by studying the vorticity equation it was shown that the torque of the Lorentz force is a major source of vorticity in the boundary layer. Emonet & Moreno-Insertis (1998) also showed how the behavior of the rising tube concerning its deformation and its wake is intermediate between that of an air bubble rising in water and that of a rigid cylinder in a flow, with those two limits being reached for smaller (air bubbles) and larger (rigid cylinders) levels of twist. Hughes et al. (1998) studied in detail the effect of using different distributions of transverse field within the tube and concluded that it is the strength of the transverse field and not its distribution which is the key factor determining the coherence of the rising tube.
In a wide range of Reynolds numbers, the wake trailing rigid
cylinders is unstable to vortex shedding. Shedding of eddies from
the wake was observed also in rising magnetic flux tubes by
Emonet & Moreno-Insertis
(1998, see
their Fig. 13 and Sect. 6), Hughes & Falle (1998) and, for a
pair of rising tubes in interaction, by Fan et al. (1998a). The
latter authors noted that each time a vortex roll was emitted, the
head of the tube together with the remaining eddy of the wake got
a net non-zero integrated vorticity, say
,
of equal
strength (but opposite sign) to the total circulation around the
eddy shed. The tube therefore becomes subject to a lift force
,
with
the forward
speed of the tube. Fan et al. (1998a) explained in that way the
sideways deviation of the tube trajectory away from a vertical
line observed in their simulation after shedding of the eddies
from the wake. In fact, as vortices were shed alternatively from
the right and left side of the wake with opposite signs of the
vorticity, the result is an oscillatory motion superimposed on the
vertical rising trajectory for the pairs of magnetic tubes.
A further step was undertaken by
Emonet et al. (2001). They studied the
trajectory of a (single) flux tube initially endowed with non-zero
vorticity. The vortex-shedding instability developed also in this
case and led to the formation of a von Kármán vortex street
behind the tube. In the simplest case, the upgoing trajectory of
their tube was oscillatory, like for the pairs of tubes of
Fan et al. (1998a). They attempted to understand the trajectory
in a variety of cases by studying the equation of motion of a
vortex filament with a time-varying value of the integrated
vorticity and subject to buoyancy force, drag and lift. The
solutions depend on the value of a dimensionless parameter,
,
the ratio of the timescales for the drag and lift force,
respectively, to modify significantly the velocity of the tube. In
the drag-dominated asymptotic regime (
), the
trajectory is a straight line which subtends a small angle
to the vertical and the tube moves with the
terminal speed given by the equilibrium between the drag and
buoyancy forces. In the drag-free (or lift- dominated) regime
(i.e.,
), the trajectory is a cycloid with horizontal
drift speed given by the ratio between buoyancy acceleration and
total circulation around the tube. The intermediate cases contain
rising trajectories with epicycles superimposed on them similar to
features observed in the 2D simulations of
Hughes & Falle (1998).
From all the foregoing we deduce that the behavior of moving
magnetic flux tubes is importantly influenced by the interaction
with the surroundings via the stresses (gas pressure plus viscous
and Maxwell stresses) acting at the interface between them. On
the other hand, the Reynolds number in these experiments, even if
still far removed from the actual astrophysical values, must be
high, at least a few hundred, for the wake formation and vortex
shedding phenomena described above to take place. This, in turn,
makes the boundary layer thin (its width varies typically like Re-1/2). So, high resolution is necessary, at least
to resolve the important small-scale features that develop in the
calculation. An adaptive mesh refinement (AMR) algorithm, as used
by Hughes & Falle (1998), is well suited to this end.
Calculations of that sort with adequate resolution are computationally expensive, and, as just seen, necessarily so, even though they are two-dimensional. To model the rise of magnetic tubes through the whole convection zone, typically with a view to understanding the origin of active regions and some of their global features, one has to resort to drastic simplifications. The approach commonly applied makes use of the slender (or thin) flux tube approximation (Roberts & Webb 1978; Spruit 1981), which models a flux tube as a quasi-1D string of mass elements under the assumption that its diameter is small compared to any other relevant length scale of the problem. The comparatively simple 1D equations derived through this approximation allow to study the evolution of a flux tube in a spherical shell with a stratification taken from the mixing-length models of the solar convection zone. Using this kind of approach, it has been possible to explain the emergence latitudes of active regions, their tilt angles as well as the asymmetry between leading and following polarities (Moreno-Insertis et al. 1994; Caligari et al. 1995; Fan et al. 1993,1994; D'Silva & Choudhuri 1991). For details and references, the reader is referred to the reviews by Moreno-Insertis (1997) and Fan (2004). The agreement with those observations provides support for the use of the slender tube approximation in this context. However, important aspects of the flux tube dynamics, like those discussed in the previous paragraphs, are disregarded by the approximation. Furthermore, the condition of small radius compared to the local pressure scale height is not fulfilled as soon as the tubes get close to the photosphere in their rise, typically when they reach levels of about a few Mm below the surface or move above them (as in the 2.5D simulations of Magara 2001). So, it is important to check how the thin tube approximation performs, compared to a more exact solution, when using it close to (or even somewhat beyond) its limits of applicability.
The objective of the present paper is to carry out 2.5D numerical experiments of the rise of a buoyant horizontal magnetic flux tube in a stratified medium using a state-of-the-art AMR code (the FLASH code) and profiting from the unprecedented computational power allowed by today's massively parallel computers. The article is organized as follows. In Sect. 2, we present details of the simulation setup, including the system of equations solved, the numerical method used, initial conditions of the simulation. Additionally, in Sect. 2.4, we present the method used to track the flux tube. In Sect. 3, the results from the 2.5D simulations are presented. The following aspects of the simulations are discussed: dependence on Reynolds number (Sect. 3.1), the dependence of the flux retention on twist (Sect. 3.2) and the evolution of the twist in the tube (Sect. 3.3). In Sect. 4, we use a thin flux tube model to explain the evolution of the flux tube as it rises through the stratified layer. In Sect. 5, we explore the transition from the thin flux tube regime to the thick flux tube regime. Finally, in Sect. 6 we discuss possible implications for understanding real magnetic flux tubes in the solar interior.
The medium is taken to be a
compressible, electrically conducting ideal gas with the equation
of state:
![]() |
(1) |
![]() |
= | ![]() |
(6) |
e | = | ![]() |
(7) |
A hydrostatic, adiabatically stratified polytropic layer of ideal
monatomic gas was chosen as the initial background stratification.
This initial polytropic layer is described by the following
temperature, density and pressure profiles:
We choose to use ,
p0 and Hp0 as units for the
density, pressure and length respectively. RT is used as the
temperature variable. The units for the velocity and time are
(the isothermal sound speed at the top of the
layer) and
t0 = Hp0/c0 respectively. The unit for specific
energies is
and the unit for the magnetic field
is
.
Expressed in these units, the
initial polytropic profiles (8) to (10) become
In the simulations presented in Sect. 3,
the polytropic layer is enclosed in the region
and
.
Therefore the density and pressure contrasts between
the bottom and top of the polytropic layer are 58 and 871respectively. The number of pressure scale heights spanned over
the height of the layer is
.
This is
comparable to the number of pressure scale heights spanned between
the bottom of the solar convection zone (at depth of 200 Mm) and
a depth of about 20 Mm. Thin flux tube simulations of flux tubes
carrying magnetic flux comparable to active regions
(
Mx) are considered to yield reliable
results up to a depth of about 10 Mm. Above that depth, the
radii of these flux tubes become comparable to the local Hp.
The number of pressure scale heights between a depth of 200 Mm
to 10 Mm is
.
The initial flux tube is taken to be axisymmetric. In cylindrical
coordinates, the longitudinal and azimuthal components of its
magnetic field have the form:
For each simulation, a flux tube was inserted near the bottom
of the stratified layer at time t=0. We have carried out a number of simulations with
different values of R0 and .
As mentioned in Sect. 2.1, the MHD Eqs. (2)-(5) do not take into account the effects of thermal, viscous and Ohmic diffusion. In practice, however, such non-ideal effects are always present in simulations as a result of numerical diffusion. Of course, if diffusive effects are actually important for the problem of interest, diffusive terms can be added to the MHD equations to capture the relevant physics. In astrophysical problems, one often encounters situations with very large dynamic and magnetic Reynolds numbers. This is also the case for the solar interior. As such, we have chosen not to impose explicit diffusive terms in the MHD equations. Diffusive effects in the following simulations are purely numerical. The amount of numerical diffusion present in the simulation is dependent on the numerical resolution used. The higher the numerical resolution, the smaller is the amount of numerical diffusion and the larger are the effective Reynolds numbers. Since the initial state of the background atmosphere in our problem is uniform (except for small pressure perturbations) in the horizontal direction and smoothly varying in the vertical direction, only relatively large grid spacing is required to resolve regions far away from the flux tube. On the other hand, high spatial resolution is needed to resolve the small-scale features at the interface between the tube and its surroundings (e.g., in regions where the flux tube fragmented).
The numerical resolution we can use is limited by the
computational resources available. To get the highest numerical
resolution where we need it, we made use of the Adaptive Mesh
Refinement (AMR) feature in FLASH. The Cartesian domain is
comprised of adjacent square blocks, each consisting of grid cells. At each time-step, all the blocks are checked to
determine if the block should be refined. If the normalized
second-order spatial derivative of the absolute field strength,
|B|, exceeds some fixed threshold in any grid cell, the
resolution of the corresponding block is doubled by interpolation
and the original block is split into four sub-blocks, increasing
the "refinement level'' of the original block by one.
The reverse process (coarsening) occurs when the normalized
second-order spatial derivative of |B| is smaller than some
threshold for all four sub-blocks. Then the resolution of each
sub-block is halved and the sub-blocks are merged.
For further details on the FLASH code, the reader is referred to the FLASH user manual (ASCI/ Alliance Center for Thermonuclear Flashes 2003).
In the simulations discussed here, the
initial flux tube does not remain a single, monolithic structure
as it rises to the top. From the results of the literature
(Fan et al. 1998a; Emonet & Moreno-Insertis 1998; Hughes & Falle 1998), we expect the tube to fragment
and lose flux by means of vortex shedding. However, for a
sufficiently high level of field line twist, a central portion of
the tube retains its identity throughout the simulation. We refer
to this central flux filament as the "main tube''. To track the
main tube, we make use of the flux function:
![]() |
(16) |
If we arbitrarily chose a value of
to define the main
tube, we have no guarantee that at a later time, the structure we
track remains coherent. In order to define a coherent main tube,
we reverse the aforementioned procedure. At the end of a
simulation (i.e. when the main flux structure has reached the top
of the domain), t=t1, we calculate
.
We test
different contour levels
.
If the contours
corresponding to a particular value of
consists of more
than one closed loop, we dismiss them. For the remaining values of
(each of which has only one corresponding closed curve),
we pick the one that encloses the maximum amount of flux (
). This is defined as the main tube for this
particular simulation. To back-track the main tube at an earlier
time, we simply calculate
for that time. Then the contour
which encloses a flux equal to
represents the main
tube at that time.
The centre of the main tube is located at the extremum of (maximum or minimum depending on whether the tube has right or
left-handed twist). Although the main tube is, in general, not
circular, we can define an effective radius,
We have carried out a number of simulations in order to study different aspects of the problem of the buoyant rise of magnetic flux tubes. Various aspects of the results are discussed in the following sections.
From the same initial setup, we have carried out simulations with
different levels of grid refinement in order to study how the
numerical resolution influences the outcome of the simulation. The
initial condition is as follows: a flux tube was inserted near the
bottom of the polytropic atmosphere at
(x0,y0) = (0.0,2.5) at
t=0. The flux tube has R0=0.5 (corresponding to
of the
local pressure scale height),
and
.
The material inside the tube has the same entropy as the external
atmosphere, so that it is buoyant. We carried out four runs from
this initial condition. Run A1 has the lowest effective
resolution. If the simulation domain was fully refined, the domain
would be spanned by
grid cells in the x and y directions. Runs A2, A3 and A4 have 2, 4 and 8 times the effective
resolution of A1 respectively. Table 1 gives
the effective resolution and effective Reynolds numbers (Re) for
each of the runs. The latter is defined as
Table 1: Simulation runs carried out to study the dependence of the simulation result on numerical resolution (and hence Reynolds number).
Figure 1 shows the distribution of the
longitudinal field (Bz) over the entire domain at t=200 for
all four runs. Figure 2 shows the zcomponent of the vorticity ()
at the same time. To
emphasize the difference in resolution between the runs, the axes
are labelled in terms of grid-points. To calculate Re, we
examined the profiles of Bz and
near the vertical
x=0 for each of the runs at t=200. From the Bz profile we
can find the tube diameter D. From the corresponding
profile, we identify a thin boundary layer near the apex of the
tube. The left and right halves (about x=0) of this boundary
layer have opposite sign. The thickness of this boundary layer -
which is the site of vorticity generation - corresponds to
.
Inspection of the vertical profile of the magnetic field
along x=0 gives the distance over which the magnetic field goes
to zero above the tube apex. This gives the thickness of the
magnetic boundary layer. For the simulations in this study, we
found that the magnetic and viscous boundary layers have similar
thickness, about 6 grid cells, indicating that the magnetic
Reynolds number
.
This is not a coincidence, since
the viscous and magnetic diffusion stem from diffusion inherent in
the numerical scheme.
![]() |
Figure 1: The structure of the wake below the rising flux tube depends on the Reynolds number of the flow. The four panels show the distribution of the longitudinal magnetic field at Reynolds numbers ranging from 25 to 2600. |
Open with DEXTER |
![]() |
Figure 2:
Same as Fig. 1 but for the zcomponent of the vorticity. At low Reynolds numbers - see cases
with ![]() ![]() ![]() ![]() |
Open with DEXTER |
The fraction of magnetic flux retained by the head of the flux
tube (i.e. the main tube) also depends on the Reynolds number.
Figure 3 shows the flux retained in the
main tube at t=280 for the four different runs (diamonds).
Clearly, with increasing Re, the percentage of flux retained by
the main tube increases. At
,
the main tube retains
of the original flux of the initial tube. Although we
cannot conclude from these simulations that the flux retained
converges to some value in the limit
,
Fig. 3 does seems to suggest that the
curve levels off for increasing Re.
![]() |
Figure 3:
Magnetic flux retained in the main tube as a function of
the effective Reynolds number. The diamonds plot the values from
simulations A1 to A4. If the amount of flux lost scaled as
![]() |
Open with DEXTER |
Emonet & Moreno-Insertis (1998)
demonstrated that in a rising, twisted flux tube, vorticity is
generated in the magnetic boundary layer between the tube and the
surrounding flow. The material in this boundary is then advected
towards the wake, leading to a loss of magnetic flux from the
tube. We can estimate the flux loss per unit time as
![]() |
(19) |
![]() |
(20) |
In ideal MHD, the ratio of mass and longitudinal magnetic flux
enclosed in the main tube, ,
should remain constant. In
numerical simulations, however, the ratio always increases with
time because of some mass diffusion across field lines. The size
of this change tells us how well the simulation approximates the
ideal MHD case. Figure 4 shows the
percentage change of this quantity for the main tube between t=0and t=280, as a function of Re. At
,
the ratio
increased by
.
This increase in the ratio diminishes for
higher resolution. At
,
the change is only on the order
of
.
Consequently, in order to compare the results from
numerical simulations with predictions in the approximation of
thin flux tubes (which assumes ideal MHD), we should take the
results from runs with the highest values of Re.
![]() |
Figure 4:
The change of ![]() |
Open with DEXTER |
For magnetic flux tubes with non-zero twist, we were able to track
a main tube. Figure 5 shows the amount
of flux retained in the main tube at t=280. It is a
monotonically increasing function of ,
a result
consistent with the previous work
of Moreno-Insertis & Emonet (1996). Thus, given
that a main tube can be tracked, the amount of flux it retains
increases with Re.
![]() |
Figure 5:
The dependence of the fraction of flux retained in the
main tube as a function of the twist parameter ![]() |
Open with DEXTER |
As Parker (1979,1974) pointed out, the radial expansion (compression) of a twisted flux tube leads to an increase (decrease) of the pitch angle of the field lines. This is a consequence of magnetic flux conservation (Fan et al. 1998b).
The radial profiles of the longitudinal and transverse field we
have chosen for the initial flux tube yield a pitch angle of the
field lines that depends on radial distance from the axis. A more
appropriate measure for the amount of twist in the tube is ,
which is dimensionless and constant over the initial
tube. If the flux tube undergoes a homologous expansion (or
compression),
evolves according to
In Fig. 6, the relationship given by this equation
is shown as a solid line. Overplotted (as diamonds) are values of
the mean twist
of the main tube in run
A4 against its effective radius
.
Following (15), we define the mean twist as
![]() |
(22) |
![]() |
Figure 6:
The relationship between the tube's effective radius and
its twist ![]() ![]() |
Open with DEXTER |
One of the main aims of this paper is to evaluate the relevance of the thin flux tube approximation with regards to describing the behavior of flux tubes in 2.5D simulations. In the following, we consider how a rising magnetic flux tube behaves in the context of this approximation (Spruit 1981; Roberts & Webb 1978).
The basic assumption of the thin flux tube approximation is that the radius of the flux tube is much smaller than any other characteristic length scale in the system (e.g., local pressure scale height and the radius of curvature of the tube axis). We assume that the quantities are uniform over the tube cross-section, so that their values at the tube axis are representative of their off-axis values. This assumption corresponds to retaining only the zeroth-order term in the axis-centered Taylor-expansion of the quantities in the tube. Higher-order treatments can also be derived (Ferriz-Mas et al. 1989; Roberts & Webb 1978). In the following, we develop a model based on the zeroth-order approximation, which is already sufficient for modelling how the physical properties near the tube centre evolve. To model the evolution of the twist in the main tube, it is necessary to extend to a second-order approach.
Instantaneous pressure balance (
)
between the tube and its surroundings is assumed.
For this assumption to hold, we require that the sound-crossing
time over the tube diameter be much smaller than the time required
for the tube to transverse a distance comparable to its diameter.
This means the adiabatic sound speed
is much larger than the
speed of the tube. Taking the terminal
velocity (Parker 1975; Moreno-Insertis & Emonet 1996):
![]() |
(23) |
For a uniform horizontal flux tube, conservation of longitudinal
magnetic flux leads to
![]() |
Figure 7: Comparison between the simulation (run A4) and the thin flux tube model. Diamonds indicate values of the physical quantities at the tube centre in the simulation and the solid lines show the predictions from the thin tube model (Eqs. (29)-(33)). |
Open with DEXTER |
Of the four simulation runs A1 to A4, we chose to compare the thin
flux tube model with results from Run A4 because the effect of
magnetic diffusion is smallest for this case.
Figure 7 shows the dependence of |B|, T,
and the tube radius as functions of
.
The values
of |B|, T and
in the simulation were taken at the
centre of the main tube, and are plotted as diamonds. The
effective radius
of the main tube is defined by
Eq. (17). The solid lines show these
quantities as calculated with Eqs. (29) to (33). Since
is the pressure contrast,
corresponds heights above the original position of the
tube. For this simulation, we have tracked the main tube until
,
corresponding to over 5 pressure
scale heights. The thin tube predictions agree well with the
simulation results over this wide range of heights, even at the
lowest values of
where the radius becomes comparable to
or larger than the local pressure scale height. The thin flux tube
predictions of temperature, density and |B| deviate from their
actual values in the main tube centre by less than
.
The
values of
calculated with the thin tube model deviate from
the actual values by
at most and the effective radius of the
main tube differs from the theoretical value by less than
.
The comparison we have made here shows that Eqs. (29) to (33) accurately describe the height dependence of the physical properties in the tube centre. In order for our thin tube model to be a dynamical model, we must also solve the equation of motion for a thin flux tube. This will then allow us to model the motion of the main tube as well as the time-dependence of its physical properties in the thin flux tube framework. This comparison is carried out in the following section.
As explained in the introduction, the motion of a twisted magnetic
flux tube in an unmagnetized environment shares a number of
features with the motion of a rigid cylinder in a flow
(Emonet et al. 2001; Emonet & Moreno-Insertis 1998). Under a number of
simplifying assumptions, the equation of motion of the magnetic
tube can be written in a simple way; among them we count: (1) zero
circulation and small Mach number of the surrounding flow; (2) not
too low Reynolds numbers (
); (3) small
length-scales and long timescales of change for the flux tube
compared to the intrinsic length and time scales, respectively, of
the flow. Under those assumptions, the integrated effect of the
fluid stresses on the periphery of the tube can be simply
described by a drag force given by the classical aerodynamic
formula and an enhancement of the inertia of the tube because of
the co-acceleration of the external medium in the vicinity of the
tube (Batchelor 1967). For a cylinder or flux tube
driven by its own buoyancy in rectilinear motion this would yield:
Figure 8 shows the rise velocity of the
main tube (upper panel) and its height (lower panel), both as
functions of time. Values from the simulation are indicated as
diamonds. The solid lines indicates the theoretical profiles found
by integrating Eq. (34). The mean density
deficit of the main tube,
,
and its initial radius, R=0.64, were used as initial
conditions for the path integration. The values of the drag
coefficient and enhanced inertia factor used are CD=2.0 and
I=2.0 respectively. To take into account the effect of tube
expansion on the buoyancy and drag forces,
Eqs. (31) and (33) were used to
update the tube radius and density at each time step of the path
integration.
![]() |
Figure 8: Rise velocity ( upper panel) of the main tube and its height ( lower panel), both as functions of time. Diamonds indicate values from the MHD simulation (run A4). The solid line shows the velocity profile calculated with the thin flux tube model, with CD=2.0 and I=2.0. |
Open with DEXTER |
The velocity and height profiles from the thin tube approximation
are in general agreement with the motion of the main tube between
t=0 and t=200. Between t=0 and t=20, the main tube
approximately undergoes a free-fall acceleration in accordance
with its own buoyancy. The motion of the main tube during this
time interval is well matched by the solid lines. The time taken
for the main tube to rise a height difference of
(corresponding to 4.1 pressure scale heights) is
.
The corresponding rise time predicted by the thin flux tube
mode is
,
which is within
of the actual
value.
One feature of the motion of the main tube which is not predicted by the thin tube calculations is the temporary deceleration of the tube between t=20 and t=30. As reported by Moreno-Insertis & Emonet (1996), this is a result of the differential acceleration experienced by the different parts of the flux tube. Since, initially, the core of the tube has the largest density deficit, it undergoes a stronger acceleration than the parts of the tube above it. This differential acceleration leads to a compression at the apex of the tube, which enhances the transverse magnetic field there. The enhanced magnetic tension near the apex of the tube eventually decelerates the core of the flux tube (our main tube) and induces an internal oscillation.
This effect, which is not described by Eq. (34) for the motion of a thin flux tube, is a source of discrepancy between the simulation results and the thin flux tube calculations. Another discrepancy is the deceleration of the main tube after t=200, which is not predicted by the thin flux tube result. This deceleration is due to the closed top boundary condition used. Equation (34) does not take this into account.
Near the top of the simulation domain, the main tube is so large that it can no longer be considered a thin flux tube. At t=253, the main tube has an effective radius equal to Hp. In Sect. 5, we examine in more detail how the limit of the thin flux tube model is reached as the flux tube increases in size.
In the simulations discussed thus far, the background atmosphere
was initially plane-parallel and the flux tube initially
axisymmetric and stationary. This confines the flux tube to a
purely vertical trajectory. To study the asymmetric rise of flux
tubes, we carried out an additional simulation run. Run B has
essentially the same initial condition as Run A4 (
,
,
R0=0.5), with the exception that the flux tube
is initially rotating solidly about its axis with angular velocity
.
The vorticity integrated over the initial flux
tube (
)
is
.
As already mentioned in the introduction, a flux tube with a net
vorticity
travelling with forward velocity
with respect to the external medium experiences a lift
force equal to
.
The lift force
causes a sideways acceleration of the flux tube so that its motion
deviates from the vertical. In simulation runs A1 to A4, we have
seen that a flux tube rising purely vertically sheds equal but
opposite amounts of vorticity to the left and right halves of its
wake. When the tube motion is no longer purely vertical, the two
sides of the tube shed unequal amounts of vorticity. Each time a
vortex roll is emitted from the flux tube, the remaining tube and
wake structure is left with a net
vorticity.
Emonet et al. (2001) found
that the quasi-periodic shedding of vorticity of alternating sign by a rising
flux tube leaves this tube and wake structure with a net
circulation that reverses its sign periodically in time. Thus
the horizontal component of the lift acceleration also
alternates periodically. This results in an oscillatory,
horizontal motion of the flux tube super-imposed on the general
vertical rise of the tube, so that it traces out a zigzag path. By adding the lift acceleration
to the equation of motion of the thin flux tube and assuming a
sinusoidal time-varying vorticity for the tube and wake, they
could model the zigzag motion of the tube. We found that
gives a good
agreement between the integrated path from the thin flux tube
model (plotted as a dashed line in
Fig. 9) and the actual path of the main
tube. The sequence of circles indicate the positions and effective
radii of the main tube as it rises and expands.
![]() |
Figure 9: The trajectory of the main tube in run B. The circles indicate the position of the main tube at different times during the simulation. The effective radius of the main tube at different instances is given by the size of the circles in the plot. The dashed line shows the trajectory from a thin flux tube calculation, taking into account the aerodynamic lift force. |
Open with DEXTER |
![]() |
Figure 10: The distribution of the vorticity at t=280 for run B. In this simulation, the flux tube had an initial net vorticity. Aerodynamic lift causes the flux tube to rise in a zigzag fashion, leaving behind a vortex street in its wake. |
Open with DEXTER |
Figure 10 shows the vorticity distribution at t=280. The three vortex rolls of alternating sign constitute a pattern reminiscent of a von Kármán vortex street. The first vortex roll shed by the flux tube is centered at (x,y) = (-3,12)and has negative sign. This means that, as this roll was being shed, the remaining tube and wake system was gaining a net positive vorticity. The lift force then acts to steer the tube and wake system towards the right. During this time, vortex rolls with positive vorticity are preferentially shed off the tube. This continues until the tube and wake system has a net negative vorticity. At this point, the lift force pushes the tube back towards the left.
In this section, we address the question: How relevant are the
predictions of the thin flux tube model in describing the
quasi-static structure of a rising flux tube as its radius
approaches and exceeds the local pressure scale height? To explore
the transition from the thin flux tube regime to the "thick'' flux
tube regime, we performed a simulation (Run C) of a rising tube,
starting with a flux tube with an initial flux 100 times larger
than in the previous cases. Here, the domain spans
and
.
The flux tube is initially centered at
(x0,y0)=(0,20) with
B0 = 12.3, R0 = 5 (corresponding to 1/4 of the local pressure scale height) and
.
Figure 11 shows the flux tube at t=0 and t=220. The five concentric green contours in the left panel are the planar projection of different field lines. Each of these contours defines a flux roll. By tracking these contours, we follow the evolution of these flux rolls. The longitudinal flux and initial radius of each of these rolls is given in Table 2.
![]() |
Figure 11: Evolution of a large flux tube (Run C). Shown in color is the longitudinal field. Each of the five circular green contours represent a magnetic field line projected onto the plane. Every contour encloses a certain amount of flux, which defines a flux roll. Table 2 gives the flux and initial radius of each flux roll. |
Open with DEXTER |
In order to evaluate the validity of the thin flux tube
approximation, we compare its predictions with the simulation
results for all five flux rolls. For the thin flux tube
predictions, Eqs. (30)-(33) are used
to calculate how the physical quantities evolve. The initial
values (B0, ,
etc.) are taken as averages inside the
rolls. The values from the thin flux tube calculations are then
compared with the average values measured in the flux rolls at
later times in the simulation.
Figure 12 shows the relative discrepancies
between the model predictions and the average values from the
simulation as a function of the effective radius of largest flux
roll (flux roll 5). The size of this flux roll is representative
of the "true'' size of the rising flux tube. As expected, the
discrepancy grows as the effective radius of the flux tube becomes
comparable to the local pressure scale height. In this regime, the
flux tube interior is sufficiently stratified that average values
of its physical properties do not match with the corresponding
thin flux tube values, so the thin flux tube approximation is no
longer appropriate for describing the state of the flux tube. The
discrepancies in temperature and density do not grow as much as
the discrepancies in field strength and
for the following
reason: for an ideal thin flux tube in pressure balance with its
surroundings, the relative difference of T and
between
the tube and its surroundings is always of order
.
Thus relative differences in T and
decrease with
increasing
(the case of an expanding tube). Consequently,
discrepancies in T and
between the thin tube predictions
and the simulation results do not tell whether the thin flux tube
approximation is good at describing the average properties of the
flux tube. In contrast, discrepancies in
and |B| show
clearly the transition between the thin flux tube regime to the
"fat'' flux tube regime, when the radius of the tube is comparable
to the pressure scale height.
Table 2: The initial radius and flux of the flux rolls as shown in Fig. 11a.
We have also examined how the average twist,
,
of each flux roll evolves as they expand.
Figure 13 shows
as a function of the effective radii of the flux rolls.
Again, there is good agreement between the data points and the
solid curve, which plots the relation given by Eq. (21).
![]() |
Figure 12: Comparison between thin flux tube calculations with average quantities inside the flux rolls (Run C). The different symbols correspond to quantities in the different flux rolls (see Table 2). As the effective radius of the flux tube approaches the local pressure scale height Hp, the discrepancy between the simulation results and the thin flux tube predictions grows. |
Open with DEXTER |
![]() |
Figure 13: Variation of the average twist of the flux rolls as a function of their effective radii. The different symbols show the average twist of the five flux rolls (Run C). The solid line shows the relation given by Eq. (21). The good match between the simulation results and relation (21) indicates that the flux rolls expand homologously. |
Open with DEXTER |
We have carried out idealized 2.5D MHD simulations of buoyant
magnetic flux tubes rising in a stratified layer over a range of
Reynolds numbers (
). Our simulations
confirm previous results in the literature. Additionally, we have
analyzed the dependence of the results on the Reynolds number. We
found that the detailed structure of the wake, as well as the
amount of flux retained in the main tube, varies with the Reynolds
number. At sufficiently high Reynolds number (
),
the vortex pair in the wake breaks into secondary rolls. The
amount of flux retained in the tube also increases with Re.
We have studied how the twist in a flux tube varies as it expands.
In particular, the dimensionless twist
,
averaged over the main tube, scales approximately linearly with
the tube radius. If flux tubes originating from the bottom of the
solar convection zone have any amount of initial twist, this twist
will be amplified upon the rise of the tube. The twist will be
maximum at the apex of the rising loop, where the cross-sectional
radius is largest.
We derived thin flux tube equations (Eqs. (30)-(33)) to model the evolution of the properties of a horizontal flux tube rising adiabatically through the atmosphere. Using these equations to model the expansion of the tube, and using Eq. (34), the motion of the main tube can be reproduced. For a tube that undergoes asymmetric vortex shedding, the lift force can be included into the equation of motion to explain the zigzag motion of the tube. The vortex shedding associated with this zigzag motion leaves behind a vorticity distribution resembling a von Kármán vortex street.
Furthermore, we studied the transition from the regime of thin flux tubes to the regime of fat flux tubes (Run C). We found that the discrepancy between the thin tube calculations and the average quantities of flux rolls increases as the flux tube expands. When its radius is comparable to one pressure scale height, the discrepancy between the thin tube calculations and the simulation results can be at least of order O(1), meaning the thin flux tube approximation is no longer valid. This result is in accordance with previous expectations. In the solar convection zone, rising toroidal flux tubes approach this limit at a depth of about 10 Mm, so it is no longer appropriate to continue thin flux tube simulations above those depths. On the other hand, our results point in the direction that below such depths, the thin flux tube approximation is useful for studying the evolution of flux tubes.
Acknowledgements
We would like to thank the anonymous referee for helpful comments which improved the presentation of the results. MCMC acknowledges financial support from the International Max Planck Research School at the Max Planck Institute for Solar System Research.