Hydrodynamic and Magnetohydrodynamic Simulations of Wire Turbulence

TL;DR

This paper presents simulations of magnetized and non-magnetized supersonic flows through wire meshes, demonstrating the generation of marginally supersonic turbulence with anisotropic structures, relevant for laboratory MHD turbulence studies.

Contribution

It provides the first detailed simulation analysis of wire mesh driven MHD turbulence, highlighting the transition from hypersonic shocks to marginally supersonic turbulence.

Findings

Downstream turbulence is marginally supersonic to transonic.

Turbulence exhibits high anisotropy.

Initial hypersonic shock leads to turbulence with Mach number around 1.

Abstract

We report on simulations of laboratory experiments in which magnetized supersonic flows are driven through a wire mesh. The goal of the study was to investigate the ability of such a configuration to generate supersonic, MHD turbulence. We first report on the morphological structures that develop in both magnetized and non-magnetized cases. We then analyze the flow using a variety of statistical measures, including power spectra and probability distribution functions of the density. Using these results we estimate the sonic mach number in the flows downstream of the wire mesh. We find the initially hypersonic (M=20) planar shock through the wire mesh does lead to downstream turbulent conditions. However, in both magnetized and non-magnetized cases, the resultant turbulence was marginally supersonic to transonic (M~1), and highly anisotropic in structure.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\mathbf{v}+P\mathbf{I}-\mathbf{B}\mathbf{B}\right)=0$$

$$\frac{\partial E}{\partial t}+\nabla \cdot [(E+P)\mathbf{v}-(\mathbf{v}\cdot \mathbf{B})\mathbf{B}]=0$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \cdot (\mathbf{v}\mathbf{B}-\mathbf{v}\cdot \mathbf{B})\mathbf{B})=0$$

$$p(s)=\frac{1}{\sqrt{2\pi {\sigma }_s^2}}\exp \left(1-\frac{(s-s_0{)}^2}{2{\sigma }_s^2}\right).$$

$${\sigma }_s^2=\ln \left(1+b^2{M}^2\frac{\beta }{\beta +1}\right).$$

Full text

††pubyear: 2019††pagerange: Hydrodynamic and Magnetohydrodynamic Simulations of Wire Turbulence–References

1 Introduction

Supersonic, magnetohydrodynamic (MHD) turbulence occurs in many astrophysical settings, from star formation in the interstellar medium (Mac Low & Klessen, 2004; Federrath & Klessen, 2012; Kritsuk et al., 2017; Offner & Liu, 2018), to supernova engines (Couch & Ott, 2015; Fryer et al., 2017; Radice et al., 2018) and remnants (Balsara et al., 2001; Inoue et al., 2009; Roy et al., 2009), to the solar wind (Alexandrova et al., 2008; Bruno & Carbone, 2013; Usmanov et al., 2014). As such, theoretical and simulation-based studies of MHD turbulence have been a robust endeavor within the astrophysical community, articulating important properties of MHD turbulence with regard to turbulent power spectra, decay rates and observational characteristics (Elmegreen & Scalo, 2004; Padoan et al., 2004; Kritsuk et al., 2011; Federrath & Banerjee, 2015; Kritsuk et al., 2017). A problem for simulation-based studies of MHD turbulence, however, has been the limited range of Reynolds numbers (both hydrodynamic and magnetic) achievable with even modern numerical codes. Typically these values are many orders of magnitude smaller than what would be expected for real astrophysical flows (c.f. Elmegreen and Scalo 2004).

Over the last two decades, high energy density laboratory astrophysics (HEDLA) studies have opened new paths to the study of astrophysical phenomena (Remington et al., 2006). The stability of supernova blast waves (Drake et al., 2002), magnetized jets (Lebedev et al., 2002; Ciardi et al., 2007; Bellan et al., 2009; Suzuki-Vidal et al., 2015) and shock-clump interactions (Hartigan et al., 2016) have all been usefully explored in laboratory settings using high energy density platforms, such as high intensity lasers and pulse power machines. Laser-based experiments have recently shown the potential for HEDLA studies to explore issues of MHD turbulence. For example, Meinecke et al. (2015) used colliding plasma jets to study the development of Kolmogorov-like turbulence, and showed that the magnetic field in the flows was amplified by turbulent motions. Earlier work by Meinecke et al. (2014) showed that MHD turbulence could be achieved after a laser driven plasma flow generated from a carbon rod had been passed through a grid. Experiments to generate MHD turbulence could also be performed using magnetized plasma flows generated with pulsed power drivers (Lebedev et al., 2014; Bott-Suzuki et al., 2015; Burdiak et al., 2017; Lebedev et al., ress)

In this paper, we report on simulations that also used a grid to generate turbulence from an initially laminar flow. Using the adaptive mesh refinement (AMR), MHD code AstroBEAR, we tracked the evolution of the flow to explore the conditions under which turbulence could be generated. The work and setup in this paper is intended to provide guidance for future studies of laboratory supersonic turbulence. The paper is organized as follows. The numerical methods and simulations are described in Section 2, Results are presented in Section 3, and a discussion of the findings are given in Section 4.

2 Methods and Simulation Parameters

To explore the possible generation of supersonic turbulence in optically thin plasma in a laboratory setting, we conducted a set of high resolution, 3D simulations of a hypersonic wind colliding with a wire mesh. Our simulations were carried out using the AstroBEAR111See: https://astrobear.pas.rochester.edu/trac/ for a listing of current capabilities. code (Cunningham et al., 2009; Carroll-Nellenback et al., 2013), a state-of-the-art, multiphysics AMR platform for solving the equations of hydrodynamics/MHD in the Eulerian frame. Two simulations were performed of a hypersonic wind being passed through a wire mesh – one in which the wind was hydrodynamic, and the other in which the wind carried a magnetic field. For the MHD case, the simulations solved the ideal MHD equations,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the mass density, $\boldsymbol{v}$ is the velocity, $P$ is the thermal pressure, $\mathbf{B}$ is the magnetic field, and $E$ is the total energy, given by $E=\frac{p}{\rho(\gamma-1)}+\frac{1}{2}\rho v^{2}+\frac{1}{2}B^{2}$ (note, the hydro simulation solved an identical set of equations, but with $\mathbf{B}=0$). To close the system of equations, an ideal gas equation of state was used with an adiabatic index very close to one ($\gamma=1.001$). Setting $\gamma\approx 1$ effectively treats the gas as isothermal. Gases typically used in laboratory experiments of wire turbulence cool efficiently through radiative loses. Thus, treating the gas isothermally allowed us to approximate radiative cooling without explicitly including radiative processes in the calculations.

In each of the simulations, a hypersonic wind was injected into a computational domain filled with a single gas at two different densities. The different densities were used to represent a dense wire mesh embedded within a sparse ambient gas (gas that made up the wires was $10^{4}$ times denser than the surrounding ambient medium). Since the wires simply represented over dense regions in the flow (i.e. did not carry an electrical current), wire material could be ablated over the course of the simulation through interaction with the incident hypersonic flow.

The wires and ambient gas were initialized to be in pressure equilibrium at the start of the simulation in order to minimize expansion of the wires into the surrounding gas. Correspondingly, the wires were a factor of $10^{4}$ colder than the ambient gas. This setup allowed the gas to remain at rest until passage of the wind. The wire mesh was composed of a 2D lattice of cylindrical wires alternating in radius between $r=0.0625$ and $r=.03125$, located at $x=1.25$ (note, the entire computational volume was $10\times 1\times 1$ in $x$, $y$, and $z$, see Fig. 1). The different wire widths were used to impose multiple wavelength perturbations to the flow. The interwire spacing of the mesh was $\Delta=.40625$. A smoothing function of the form $f=\max[0,1-(x^{2}+y^{2})^{4}]$ was applied to the wires to avoid sharp discontinuities in the fluid variables between the wires and ambient medium (c.f. lower-right inset of Fig. 1).

At $t=0$, a hypersonic post-shock wind (corresponding to a planar shock of sonic Mach number $M_{s}=20$, with respect to the pre-shock ambient gas) was continuously injected into the lower $x$ boundary of the computational domain ($x=0$). Once the leading edge of this wind reached the upper $x$ boundary ($x=10$), the given simulation was terminated. Note, a magnetic field was injected into the domain along with the wind in the MHD case (i.e. only the wind contained a magnetic field, pre-shock gas did not). The injected field was oriented perpendicular to the direction of the wind ($\vec{B}=B\hat{y})$. The magnetic pressure in the wind ($P_{mag}\equiv\frac{B^{2}}{\sqrt{4\pi\rho_{wind}}}$) was ten times less than its thermal pressure (i.e. $\beta\equiv\frac{P_{gas}}{P_{mag}}=10$). Given these parameters, the wind was super-Alfvénic, with an Alfvén Mach number $M_{A}\equiv(0.5\gamma\beta M_{s}^{2})^{1/2}\approx 45$.

Both simulations were performed using upwinded Godunov-type integration methods that were second-order accurate in space and time. This included the piece-wise parabolic method (PPM) of Colella & Woodward (1984) for the spatial reconstruction, combined with the Harten-Lax-van Leer Contact (HLLC) Riemann Solver (Toro et al., 1994) for the hydro flux update and the HLLD Riemann solver for the MHD flux update (Miyoshi & Kusano, 2005). The magnetic field was evolved according to a corner transport upwind (CTU) scheme (Cunningham et al., 2009). Time stepping was performed using a two-step Runge-Kutta method (Shu & Osher, 1988).

The computational domain had a base grid composed of $1600x160^{2}$ computing zones. On top of this base grid, the simulations were initialized with $1$ level of refinement, which was centered on the wire mesh. Over the course of the simulation, gradients in the fluid variables triggered refinement in other regions of the computational domain. With the chosen refinement criteria, turbulent substructure was adequately captured at the highest resolution of the simulation, which was $\Delta x_{min}=3.13\times 10^{-3}$. Periodic boundary conditions were imposed on the $x-y$ and $x-z$ planes to enforce slab symmetry of the setup. Inflow/outflow boundaries were used on the remaining faces of the domain, corresponding to the lower/upper $y-z$ planes, respectively.

3 Results

In the following section, we present results of velocity and magnetic field dispersion along the length of the shock tube (Section 3.2), power spectra of the velocity and magnetic fields (Section 3.3), and probability density function fitting of the gas density (Section 3.4). We begin by orienting the reader to the general morphological features of supersonic flows past a wire mesh revealed in our study.

3.1 Flow Morphology

The volume rendering in Figure 2 illustrates typical flow structures found in the simulations. First, intersecting bow shocks are formed along the face of the wire mesh as the injected hypersonic wind collides with wire material or passes through the mesh openings. The flow is strongly sheared as it passes through the bow shocks that wrap around the wires. The shocked flow shows excitation of Kelvin-Helmholtz (KH) instabilities along the shear layer. In addition to KH modes, the thin, isothermal shocks are unstable to Richtmeyer-Meshkov (RM) and Nonlinear Thin Shell (NTS; Vishniac 1994) instabilities. Together, these unstable modes seed perturbations that grow and propagate downstream. Thus, passage of the shock through the wire mesh transforms an initially laminar flow into one rich in multiscale, inhomogeneous structure.

Figure 3 shows density projections perpendicular to the bulk flow direction for the hydro and MHD runs at the final simulation time. Beginning with the top panel of the figure for the hydro case, the flow looks virtually the same irrespective of projection angle. Additionally, the dark ‘knots’ occurring at $x\approx 1.5$ in both the $x-y$ and $x-z$ planes represent the wires (c.f. upper inset of Fig. 1). Just behind the wires at $x\approx 2$, dense post-shock wind material is visible at the nodes of intersecting shocks (i.e. Mach stems), which occur as the gas passes through the ‘cells’ of the wire mesh and collides downstream. In a sense, each cell creates an expanding jet-like flow that interacts with its neighbor. This produces interaction regions that are filamentary in structure and oriented mainly along $x$. Finally, further downstream, the projected density decreases as the flow rebounds from the intersecting shocks and expands. Eventually multiscale, turbulent substructure develops along the remaining length of the shock tube.

The MHD run (bottom panel Fig. 3) looks nearly identical to the hydro case, with two important differences. First, the MHD case exhibits some variation in the flow between the $x-y$ and $x-z$ projected planes. Second, the initial density inhomogeneities formed just beyond the wire mesh are denser in the MHD run compared to the hydro run. The break in symmetry between the $x-y$ and $x-z$ planes is due to the magnetic field being initially orientated along $y$ ($\vec{B_{0}}=B_{0}\hat{y}$). As material streams through the wires, the magnetic field is dragged around the wires leading to extended/elongated loops of field oriented in the $x$ direction. Where the extended loops of field from neighboring mesh cells interact, the gas is compressed to higher densities. Thus, the increased density of the nonlinear structures formed just downstream of the wire mesh, as well as the differences between the projection maps along $y$ and $z$, are due to the distortion of the field as it is dragged around and through the wires.

3.2 Time-averaged Velocity and Magnetic Field Dispersion

We now turn to more quantitative measures of the flow, beginning with velocity and magnetic field dispersion along the length of the shock tube. Standard deviations were calculated for each Cartesian component of the velocity and magnetic field in 10 sequential sub-regions of the domain. These were calculated over the last ten time states of the given simulation (where $\Delta t=0.005~{}t_{final}$), and averaged. The results are plotted in Figure 4.

Beginning with the velocity field (left hand panel, Fig. 4), the flow begins at $\tilde{x}=1$ oriented solely along $x$ ($\vec{v_{0}}=v_{0}\hat{x}$), hence $\sigma_{v_{y}}=\sigma_{v_{z}}=0$ at $\tilde{x}=1$ (note, $v_{x}>>1$, so $\sigma_{v_{x}}\neq 0$ at this position is due to numerical noise). At $\tilde{x}=2$, the wires are contained within the analysis sub-region. This accounts for the increase in the standard deviation of all velocity components as the flow travels through the wires and is deflected. However, the component perturbed * most strongly* upon passage through the wire mesh is $v_{x}$, evidenced by the sharp increase in $\sigma_{v_{x}}$ at this position. This result is consistent with the prominent elongated structures visible in the column density maps along $x$ (see Section 3.1). In the other coordinate directions, the velocity standard deviations are roughly isotropic.

Beyond the wires ($\tilde{x}>3$), the flow relaxes into more coherent velocity structures, visible by the steep decline in $\sigma_{v_{x}}$ and the flattening out of each of the curves. Yet, the velocity continues to be largely anisotropic, remaining preferentially perturbed along the bulk flow direction ($\sigma_{v_{x}}>\sigma_{v_{y}}\sim\sigma_{v_{z}}$). Further downstream still ($\tilde{x}>4$), $\sigma_{v_{x}}$ is lower in the MHD case than in the hydro, which is consistent with the magnetic field resisting motions along this direction. Note, however, that only slight differences exist between the MHD and hydro runs for $\sigma_{v_{y}}$ and $\sigma_{v_{z}}$, which, once again, is consistent with the flow being only weakly perturbed in these directions as it passes through the wire mesh. Thus, the MHD velocity field is slightly more isotropic downstream, given the restriction of fluid motions along $x$ by the magnetic field.

Turning now to the dispersion in the magnetic field (right hand panel, Fig. 4), note that the overall shape and trend seen in the curves are similar to those of the velocity standard deviations. This is not surprising given the field is perfectly coupled to the gas in ideal MHD (except on the grid scale, where numerical dissipation of the field can occur). Thus, magnetic field perturbations are tied to velocity motions and scatter is mirrored in each of the distributions. As described above for velocity, before the flow passes through the wire mesh at $\tilde{x}=1$, the magnetic field is in its initial, unperturbed state ($\vec{B_{0}}=B_{0}\hat{y}$). Once the flow drags the field through the wire mesh at $\tilde{x}=2$, a sharp increase is again visible in the standard deviation of the $x$ component of the field ($\sigma_{B_{x}}>>\sigma_{B_{y}}\sim\sigma_{B_{z}}$). The smaller, concomitant increases in $\sigma_{B_{y}}$ and $\sigma_{B_{z}}$ at this point are due to the flow being only marginally perturbed in $v_{y}$ and $v_{z}$, as described previously. Finally, far away from the wires ($\tilde{x}>2$), a decline in the dispersion of all components is again visible, with the standard deviation in $B_{x}$ remaining the largest downstream.

3.3 Power Spectra

Power spectra were constructed for the Cartesian ($\hat{x}$, $\hat{y}$, $\hat{z}$) and Helmholtz-decomposed solenoidal and compressive components of the velocity and the magnetic fields (satisfying $\nabla\cdot\vec{a}=0$ and $\nabla\times\vec{a}=0$, respectively) using AstroBEAR’s discrete Fourier transform module, which performs fast Fourier transforms (FFT) of the data on an AMR mesh (see Carroll-Nellenback et al. 2014 for details). Spectra were calculated for the last time state in each of the runs, in the last region of the computational domain ($\tilde{x}=10$). Given the location of this region furthest away from the wire mesh, nonlinear perturbations seeded in the flow by the KH, RM, and NTS instabilities have had the most time to homogenize and decay into a turbulent cascade. Before discussing specific details about the power spectra for the different runs, we would like to make a few broad remarks about the trends seen in each of the panels of Figure 5. First, note that the driving scale for energy injection is visible at $k\approx 2$, corresponding to a physical wavelength of $\lambda\approx.5$. This length scale is consistent with the inter-wire spacing of the mesh, $\Delta\approx 0.4$. Second, at the high $k$ end of the spectrum, numerical dissipation on the grid scale limits the inertial range of any turbulent energy cascade that could be captured by the simulation. This is evidenced by the decline in power for each of the curves beyond $k>20$. Thus, we take the inertial range captured by the simulations to lie between $2<k<20$. Lastly, each spectrum exhibits a scaling law over this limited range in $k$ ($P\propto k^{\alpha}$). To quantify the slope of the power-law scaling in each of the cases, best-fit lines are overlaid in the figure.

Turning now to the power spectra of velocity for the hydro case (leftmost panel of Figure 5), the slope of the best-fit line to the $v_{x}^{2}$ spectrum is $\alpha=-1.9$, which is very near the value expected for a turbulent cascade under supersonic, compressible, isothermal conditions (i.e. $\alpha=-2$, see for example, Kritsuk et al. 2007; Schmidt et al. 2009). By comparison, the slopes of the best fit lines in the other dimensions are $\alpha=-1.0$ and $\alpha=-.97$ for $v_{y}^{2}$ and $v_{z}^{2}$, respectively. That the slopes vary for the different components of $v$ indicate that specific kinetic energy was deposited to the flow anisotropically. Moreover, the results support that the flow was most strongly perturbed along the direction of the wind propagation, since the overall power of the $v_{x}$ mode was greater than $v_{y}$ and $v_{z}$ on all scales.

Considering next the total hydro velocity power spectrum, it may be unsurprising that the derived slope lies between $1<\alpha<2$, given the slopes of the individual components. However, it is interesting how close this value ($\alpha=-1.6779$) lies to the classical Kolmogorov (1941) spectrum for incompressible turbulence ($\alpha=-5/3$). Power spectrum of the Helmholtz-decomposed velocity field for the hydro case indicates that more energy was injected into solenoidal modes than compressive (i.e. $v^{2}_{sol}>v^{2}_{div}~{}\forall~{}k$). The ratio of solenoidal to compressive energy can be used to determine the effective sonic Mach number of a turbulent flow (Molina et al., 2012), and is investigated below in Section 3.4.

Turning now to the MHD case, power spectra of the velocity (middle panel, Fig. 5) show similar trends and scaling to those in the hydro simulation, with the largest difference being the slope of the compressive velocity mode, $v^{2}_{div}$, which has decreased to $\alpha=-1.67$. Consequently, even more energy ($E_{k}=\int P_{k}dk$) has been deposited into solenoidal modes than compressive in the MHD case. Power spectra of the magnetic field components (right panel, Fig. 5) mirror the MHD velocity spectra as expected, due to the ideal MHD nature of the simulations.

3.4 Mach Number Estimation

The distribution of gas density in simulations of isothermal, supersonic turbulence is well known to follow a log-normal probability density function (PDF) (Blaisdell et al., 1993; Vazquez-Semadeni, 1994; Padoan et al., 1997; Passot & Vázquez-Semadeni, 1998; Scalo et al., 1998; Kritsuk et al., 2007; Federrath & Klessen, 2012). In the isothermal limit, it is often convenient to define the logarithmic density, $s=\ln{(\rho/\rho_{0})}$ (where $\rho_{0}$ is the mean gas density). This change of variables allows the log-normal PDF to be expressed as a simple Gaussian distribution in $s$:

[TABLE]

Here, $\sigma_{s}$ denotes the standard deviation of $s$, and $s_{0}$ is the mean logarithmic density in the flow, as usual for the normal distribution.

Early studies of hydrodynamic turbulence showed that increasing the RMS sonic Mach number ($M$) of the flow led to a proportional increase in the standard deviation of the gas density, i.e. $\sigma_{\rho}=bM$ (e.g. Padoan et al. 1997; Passot and Vázquez-Semadeni 1998). The proportionality constant $b$ in this relation was later shown to depend on the type of turbulent modes present in the flow, ranging from $b=1/3$ for purely solenoidal to $b=1$ for compressive (Federrath et al., 2008, 2010). For this reason, $b$ is today commonly referred to as the ‘turbulent forcing parameter’.

Assuming the logarithmic gas density follows a Gaussian PDF (Eqn. 5), it can be shown for hydrodynamic turbulence that $\sigma_{s}$, $b$, and $M$ are related through $\sigma_{s}^{2}=\ln({1+b^{2}M^{2}})$ (Padoan et al., 1997). For MHD turbulence, this formula must be corrected to account for the magnetic field. Molina et al. (2012) derived the following semi-analytical formula for $\sigma_{s}$ (again assuming a Gaussian PDF in $s$), as a function of $b$, RMS sonic Mach number $M$, and $\beta$

[TABLE]

Equation 6 has the necessary property that it reduces to the hydrodynamic formula in the limit ($\beta\rightarrow\infty$), and agrees well with 3D simulations of super-Alfvénic MHD turbulence (Molina et al., 2012).

Using Equation 6, one can estimate $M$ of a turbulent flow, given suitable approximations of $\sigma_{s}$, $b$, and $\beta$. To approximate $\sigma_{s}$ for the different runs, we fit Gaussian distributions to volume-weighted density PDFs of the data (Fig. 6). The PDFs were constructed for each of the sub-regions along the shock tube, and averaged over the last ten time states of the simulation. These are shown clockwise in the figure for $\tilde{x}=4$, $6$, $8$, and $10$ (blue points corresponding to hydro data and red to MHD). Best fit Gaussians are overlaid in each of the panels (solid black and green curves for the hydro and MHD fits, respectively). As can be seen in the figure, Gaussian distributions closely approximate the PDFs, and the dispersion of the data is minimal over the last ten time states. Thus, we use the standard deviations of these best fit Gaussians (given in the legend) in Equation 6 to compute the estimated $M$ as described below.

Before discussing our best estimate of $M$ in each of the cases, we draw the readers attention to the evolution of the PDFs shown in Figure 6. Note that the shapes of the PDFs change in traversing the computational domain from $\tilde{x}=4$ to $\tilde{x}=10$; both the hydro and MHD PDFs become more tightly peaked for larger $\tilde{x}$. Federrath et al. (2008, 2010) showed that the width of density PDFs in simulations of fully developed turbulence correlates with the relative strength of compressive vs. solenoidal modes in the flow, with solenoidal flows producing narrower distributions than compressive (see for example fig. 2 and 6 in those papers). Thus, these results might suggest that the flow switches from being more strongly compressive in regions closer to the wires, to more solenoidal in regions further downstream. This behavior would be consistent with Federrath et al. (2010), which shows that even turbulence driven by purely compressive forcing can decay into some fraction of solenoidal modes (about $1/2$ for 3D flows). Power spectra of the flows also supported that the velocity field was predominantly solenoidal at large $\tilde{x}$ (Section 5).

We now turn to Figure 7, which shows the estimated RMS sonic Mach number $M$ along the shock tube in both the hydro and MHD cases. For each case, we used a spatially averaged value for $\beta$ in each of the sub-domains of the shock tube ($\tilde{x}=1-10$), averaged over the last ten time states. For both the hydro (blue regions of the plot) and MHD data (black regions of the plot), a range in potential turbulent forcing parameter $b$ is considered ($1/3\leq b\leq 1$) in calculating $M$ for a given $\tilde{x}$. This is represented in the figure as the shadowed area above and below the solid curves of the same color (the intermediate value of $b$ lies along the solid curve). Thus, for a given $\tilde{x}$, the highest possible $M$ corresponds to $b=1/3$, and the lowest $b=1$. The range of possible $M$ also depends on the best fit $\sigma_{s}$, which varies along $\tilde{x}$ as described previously. As can be seen in the figure, at best ($b=1/3$), $M$ is marginally supersonic for most of the domain beyond the wires ($\tilde{x}>5$) in both the hydro and MHD cases, and transonic in the longest evolving regions of the flow, i.e. $M\approx 1$ at $\tilde{x}=10$ (note that shock compression of the field immediately beyond the wires led to a steep decrease in $\beta$ in that region, thus explaining the high mach numbers immediately adjacent to the wires, i.e. $\tilde{x}=3$).

4 Discussion

We have presented a set of high-resolution, 3D AMR simulations to test whether supersonic MHD turbulence could be generated by passing a hypersonic ($M_{s}=20$), super-Alfvénic ($M_{A}=45$) flow through a wire mesh. Our findings support that marginally supersonic turbulence can be achieved by this experimental setup, irrespective of upstream hydrodynamic or MHD conditions. Further, our results are strongly consistent with the generation of anisotropic turbulence in the downstream (post-shock) flow. We presented a number of analyses in support of these conclusions, including a study of the morphological features of the flows (Section 3.1), time-averaged dispersion of the velocity and magnetic field along the shock tube (Section 3.2), power spectra of the velocity and magnetic fields (Section 3.3), and a sonic Mach number estimate based on the post-shock gas density PDF for each of the runs (Section 3.4).

The prominent morphological features of hypersonic flow past a wire mesh include unstable bow shocks surrounding the wires and multiscale post-shock turbulent substructure oriented mainly along the bulk flow direction ($\hat{x}$). These structures arise as follows. In both the hydro and MHD cases, the incident hypersonic wind is shocked and deflected around the wires. This leads to the excitation of various fluid instabilities in the bow shock layer (such as the KH, RM, and NTS instabilities), which causes ripples and distortions along the bow shock front. This effect is clearly visible in the volume rendering of Figure 2, but is also apparent when viewed in projection (Fig. 3). As each bow shock wraps around its associated wire, the interaction of bow shocks from neighboring wires then leads to the formation of mach stems in the region immediately downstream of the wire mesh. Together, these shock intersections and their associated instabilities seed density perturbations which grow into inhomogeneous clumpy and filamentary substructure further downstream. Throughout the the simulation, the density inhomogeneities remain oriented predominantly along the bulk flow direction in both the hydro and MHD cases.

The standard deviations of the velocity field ($\sigma_{vx}$, $\sigma_{vy}$, $\sigma_{vz}$) reflect the development of the morphological structures seen in the simulation. As shown in Figure 4, the component of the flow perturbed the most upon passage through the wire mesh was $v_{x}$, evidenced by a sharp increase in $\sigma_{v_{x}}$ in the region of the shock tube containing the wires. Further, this component remained the most variable along the remainder of the shock tube, indicative of perturbations being predominantly amplified in the $x$ direction. These trends were largely the same in both the hydro and MHD cases, with only minor differences arising between the runs. The dispersion of the magnetic field for the MHD case was also presented in Figure 4 (right hand panel), with the trends mirroring those found in velocity, as to be expected due to flux-freezing in ideal MHD.

Power spectra were presented of the velocity and magnetic field at the end of the simulation showing a power-law scaling relation in wave number (i.e. $P\propto k^{-\alpha}$), characteristic of turbulent flows. While the scaling power $\alpha$ varied between $1\lesssim\alpha\lesssim 2$ for the different components, they remained comparable between the hydro and MHD cases. In both runs, however, the component closest to a supersonic Burgers-type cascade ($\alpha=2$) was $v^{2}_{x}$. This supports that turbulence was most strongly seeded along the initial flow direction.

PDFs of the gas density indicated a relatively stronger solenoidal component of the flow compared to compressive in the regions furthest away from the wire mesh (Fig. 6). This finding was somewhat qualitative; Federrath et al. (2008, 2010) showed that the width of the gas density PDF for purely solenoidally driven turbulence is significantly narrower than purely compressive turbulence of the same $v_{rms}$ and $M$. In traversing the shock tube, we found the widths (at half maximum) for both the hydro and MHD gas density PDFs to decrease by roughly $50\%$, suggesting that the flow started out more compressive in regions nearest the wire grate, but became more solenoidally-dominated downstream. This could be the result of energy transfer between the modes by nonlinear interactions in the flow. For instance, Federrath et al. (2010) showed that energy transfer between compressive and solenoidal modes readily occurs in high resolution 3D turbulence simulations. As shown in that paper, even turbulence driven by fully compressive forcing decays into a mixture of compressive and solenoidal velocity modes, with the ratio being $\approx 1:1$.

Given that our gas density followed a nearly log-normal distribution, in agreement with many previous studies of supersonic, isothermal turbulence, e.g. Blaisdell et al. (1993); Passot & Vázquez-Semadeni (1998); Kritsuk et al. (2007), we were able to estimate the RMS sonic Mach number of the flows according to Equation 6. We found $M$ to lie between $1\lesssim M\lesssim 3$ for most of the downstream shock tube (in both the hydro and MHD cases). In the region furthest away from the wire mesh, the range in probable $M$ was found to relax to transonic values ($M\approx 1$ for $\tilde{x}=10$).

We would like to finish this discussion with a few remarks on our choice of wire initialization. By treating the wires as an overdense fluid rather than a fixed, solid boundary, some ablation of wire material naturally resulted as the shock interacted with the mesh. In the present simulations, roughly $30\%$ of the downstream gas (located at $\tilde{x}=10$) was ablated wire material. In the current setup, we expect the dominant effect of this wire ablation and entrainment would have been the seeding of turbulence, given the relative temperatures of the two phases of the gas (recall the wire gas was $10^{4}$ times colder than the surrounding post-shock gas which constituted the majority of the shock tube). The difference in sound speeds of these two materials alone would have translated into the wire gas being insufficient in driving significant thermal dynamics in the ambient medium.

While we did not investigate the effect of wire ablation and entrainment on turbulence generation in the present set of simulations, it would be very interesting to explore this in future work. This is especially true since at least some degree of ablation could occur experimentally. While ablated mass in experiments such as Burdiak et al. (2017) would probably have been small on the time-scale of the experiment, the situation could be quite different in the case where the wires are “preconditioned," e.g. by a separate current pulse, and converted into a lower density vapor state. To our knowledge, no physical experiments have been done on this as yet, but we would expect considerable mass entrainment in those cases. One way of varying the degree of wire ablation numerically would be to vary the density contrast between the wires and the ambient medium. Experimentally, this could be done using different materials for the plasma flow and the wire obstacle (solid or with reduced density), and then use spectroscopic diagnostics to observe the presence of the obstacle material in the post-wire plasma flow.

5 Conclusions

Our results suggest that experiments using a shock-tube like configuration with a wire mesh can, in fact, produce turbulent magnetized flows. However, because the turbulence requires passage through an obstacle (the wire mesh), the downstream flows remain close to sonic. Thus, any experiments using this configuration must be tuned to astrophysical environments dominated by such transonic flows, such as those in galaxy clusters. The strongly supersonic turbulence found in star formation regions likely can not be explored using the configuration presented in this paper.

Acknowledgments

This work used the computational and visualization resources in the Center for Integrated Research Computing (CIRC) at the University of Rochester and the computational resources of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, provided through allocation TG-AST120060 from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. Financial support for this project was provided by the Department of Energy grant DE-SC0001063, the National Science Foundation grants AST-1515648 and AST-1813298, and the Space Telescope Science Institute grant HST-AR-12832.01-A. This work was also supported by the US Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration (Contract No. 89233218CNA000001).