Nonlinear Dynamics of Non-uniform Current-Vortex Sheets in Magnetohydrodynamic Flows

A theoretical model is proposed to describe fully nonlinear dynamics of interfaces in two-dimensional MHD flows based on an idea of non-uniform current-vortex sheet. Application of vortex sheet model to MHD flows has a crucial difficulty because of non-conservative nature of magnetic tension. However, it is shown that when a magnetic field is initially parallel to an interface, the concept of vortex sheet can be extended to MHD flows (current-vortex sheet). Two-dimensional MHD flows are then described only by a one-dimensional Lagrange parameter on the sheet. It is also shown that bulk magnetic field and velocity can be calculated from their values on the sheet. The model is tested by MHD Richtmyer–Meshkov instability with sinusoidal vortex sheet strength. Two-dimensional ideal MHD simulations show that the nonlinear dynamics of a shocked interface with density stratification agrees fairly well with that for its corresponding potential flow. Numerical solutions of the model reproduce properly the results of the ideal MHD simulations, such as the roll-up of spike, exponential growth of magnetic field, and its saturation and oscillation. Nonlinear evolution of the interface is found to be determined by the Alfvén and Atwood numbers. Some of their dependence on the sheet dynamics and magnetic field amplification are discussed. It is shown by the model that the magnetic field amplification occurs locally associated with the nonlinear dynamics of the current-vortex sheet. We expect that our model can be applicable to a wide variety of MHD shear flows.

This work was supported by a Grant-in-Aid for Scientific Research (B) (Grant No. 26287147) and (C) (Grant No. 23540453) from the Japan Society for the Promotion of Science, a Grant-in-Aid for Research Promotion, Ehime University, and joint research project of ILE, Osaka University. The authors would like to thank Professor K. Hiraide and Professor S. Yanagi for their mathematical advice. We are deeply grateful to Professor A. Kageyama, Professor K. Kusano, Professor J. G. Wouchuk, and Professor Z. Yoshida for their valuable comments and discussions on plasma physics. We are also particularly indebted to Professor Y. Kaneda for his advice and suggestions.

Authors and Affiliations

1. Laboratory of Applied Mathematics, Graduate School of Engineering, Osaka City University, Sugimoto-cho, Sumiyoshi, Osaka, 558-8585, Japan

   C. Matsuoka

2. Institute of Laser Engineering, Osaka University, Suita, Osaka, 565-0871, Japan

   K. Nishihara & T. Sano

Corresponding author

Correspondence to C. Matsuoka.

Communicated by Andrew Szeri.

Appendix 1: Current Density and Surface Current Density

Here, we mention the relation between the current density \({\varvec{j}}\) and the surface current density \({\varvec{j}}_\mathrm{s}\) (or the current sheet strength \(j_\mathrm{s}\)). From the integral form of the Ampère’s law (2.5), the following relation holds for the (total) current I:

where \({\varvec{j}}= j {\varvec{\hat{e}}}_z\), \(\mathrm{d}{\varvec{A}} = {\varvec{\hat{e}}}_z \mathrm{d}A\) is the area vector associated with a small rectangular area \(\mathrm{d}A = \mathrm{d}\xi \mathrm{d}s\) (\(\mathrm{d}\xi \) is a small element in the direction of the unit normal \({\varvec{n}}\)) encircling the interface and \(d\mathrm{s}\) is the line element of the interface. The current sheet strength \(j_\mathrm{s} \equiv \int j \mathrm{d}\xi \) is given by

From this equation, we obtain the relation

Appendix 2: Velocity Field in Bulk

Now we decompose the velocity field \({\varvec{u}}\) as

where \({\varvec{u}}_\mathrm{s}\) denotes a solenoidal vector field and \(\Xi \) is a irrotational potential field that satisfies the Laplace equation \(\triangle \Xi = 0\) in the bulk and the boundary conditions:

Since the normal component of the fluid velocity is continuous at the interface; \(U_n = q_n\), we have

The potential field that satisfies the conditions (9.1) and (9.3) is \(\Xi = 0\) only. Therefore,

holds in the velocity field. This solenoidal field corresponds to the vortex-induced velocity \({\varvec{q}}\) or \(q^*\) in (2.13). Then, the Birkhoff–Rott equation (2.13) can describe the bulk velocity field \({\varvec{u}}=(u,v)\) in its complex form as

We present this bulk velocity field at \(t=4.0\) for the Atwood number \(A= -0.2\) and the Alfvén (Mach) number \(R_A^2=10^3\) in Fig. 8 (right). Unlike the magnetic field, we see that the velocity field has the normal component at the interface.

Appendix 3: Numerical Methods for the Computation of the Sheet Model

In this appendix, we present numerical methods in order to solve the governing equations in the sheet model (2.12), (2.29), and (6.2) (Matsuoka and Nishihara 2006). Discretized equations to Eq. (6.2) are given as

in which \(X_j \equiv X(\theta _j)\), \(Y_j \equiv Y(\theta _j)\), and \(\gamma _j \equiv \gamma (\theta _j)\) are expanded into discrete Fourier series

with the derivatives

where \(X_{\theta , m} = \left( {\partial X}/{\partial \theta }\right) _m\) and so on. Note that the derivatives (10.3) do not involve errors which necessarily arise in derivative representations by usual difference approximations.

As pointed out by Kerr Kerr (1988), point vortices, i.e., the grid points on the interface, tend to make a cluster around a bubble when the Atwood number is high, which is caused by the fact that the velocity difference between a bubble and spike becomes large for higher Atwood numbers. As a result of that, grid points around a spike decrease, and the calculation fails due to numerical instabilities. In order to avoid that, we use a grid redistribution method developed by Baker et al. Baker and Nachbin (1998) so that grid points are arranged equidistantly. This method is analogous to the node spreading presented by Kerr so as to obtain equally located grid points (Kerr 1988); however, the accuracy of the grid redistribution method is higher than his method, and we can take more grid points. The grid redistribution method is as follows. Now we have a representation for the interface \(\left( X(\theta , t), Y(\theta , t)\right) \) with equally spaced Lagrangian markers \(\theta \) at a time t. Then, we seek a new mapping from \([0, 2\pi ]\) onto itself, \(\theta \rightarrow p\), such that

in which \(s_\theta = \sqrt{X_\theta ^2 + Y_\theta ^2}\) and L is the whole length of the interface at the time t:

Since we want to evenly spaced grid points \(p = mh\) \((m=0, \dots N)\), where \(h=2\pi /N\), N, the number of grid points, we seek the following sequence:

Note that \({\bar{\theta }_m}\) \((m=0, \dots N)\) with \({\bar{\theta }_0} = \theta _0\) and \({\bar{\theta }_N} = \theta _N\) in equation (10.6) mapped to mh in p is not equally divided \(\theta _m = 2\pi m/N\) but a new position in \(\theta \) which are not necessarily evenly spaced. In order to find these new parameterizations successively, Newton’s method is used. Integrals in equations (10.4) and (10.5) are evaluated by the Fourier series of the integrand. Once new marker \({\bar{\theta }_m}\) is given, evenly spaced new position \(\left( X({\bar{\theta }_m},t),Y({\bar{\theta }_m},t)\right) \) and the strength \(\gamma ({\bar{\theta }_m},t)\) are determined by cubic splines using \(\theta _m\), \({\bar{\theta }_m}\), \(\left( X(\theta _m,t),Y(\theta _m,t)\right) \), and \(\gamma (\theta _m,t)\). Thus, the redistribution of grid points at a time t is completed. With these new dependent variables, new velocities \(X_t({\bar{\theta }_m},t)\), \(Y_t({\bar{\theta }_m},t)\), and \(\gamma _t({\bar{\theta }_m},t)\) are evaluated at time t; then, we can regard the discrete variable \(\theta _j\) in the Fourier series \(X_j - {\bar{\theta }_j}\), \(Y_j\), \(\gamma _j\) as the ones in the mapped space p, where the points are distributed with equal interval h. This redistribution is performed every time step.

Appendix 4: Evaluation of Mesh Points for a Multivalued Interface

In order to numerically calculate the bulk magnetic field for an interface that is multivalued as shown in Fig. 9, we need to distinguish that a grid point \({\varvec{x}}_i =(x_i, y_i)\) (\(i \in {\mathbb N}\)) in the plane is on which side (above or below) of the interface. In this appendix, we present a method of evaluation for that. Suppose that y(i, j) (\(i,j \in {\mathbb N}\)) denotes the y coordinate designated by a mesh point (i, j) (i: fixed) and Y(i, k) (\(k \in {\mathbb N}\), i: fixed) is the y coordinate that the line \(x_i =\) constant and the interface intersect, where the integer k satisfies \(k \le K\), K the number of intersection points (refer to Fig. 18).

When the point \({\varvec{x}}_i\) is located below of the interface (\({\varvec{x}}_i \in D_1\)), the integer k that satisfies \(y(i,j) > Y(i,k)\) becomes even (e.g., the point \(y(i,5) \in D_1\) satisfies \(y(i,5) > Y(i,4)\) (\(k=4\)) in Fig. 18). On the other hand, when the point \({\varvec{x}}_i\) is located above of the interface (\({\varvec{x}}_i \in D_2\)), the integer k that satisfies \(y(i,j) > Y(i,k)\) becomes odd (e.g., the point \(y(i,6) \in D_2\) satisfies \(y(i,6) > Y(i,5)\) (\(k=5\)) in Fig. 18).

Schematic figure of the interface, grid points y(i, j), and the crossing points Y(i, k) (\(k=1,2 \ldots K\)) for a fixed \(x_i\), where the cross and the black circle denote y(i, j) and Y(i, k), respectively. In this example, we set \(K=5\)

Full size image

In order to find the intersection points Y(i, k) (\(k = 1,2 \cdots K\)) numerically, we increase the number of grid point of the interface N (\(N=512\) here) to mN (\(m \in {\mathbb N}\)) in advance by the interpolation (we adopt the spline interpolation here) and detect Y(i, k) for a fixed \(x_i\), where we select \(m = 20 \sim 60\). If the condition \(|y(i,j)-Y(i,k)| < \epsilon \) (\(0 < \epsilon \ll 1\)) is satisfied between a point on the interface y(i, j) and an intersection point Y(i, k), we regard the point Y(i, k) as being located on the interface, where we select, e.g., as \(\epsilon = 10^{-4}\).

Appendix 5: Linear Growth Rate of Richtmyer–Meshkov Instability

We consider the same configuration as Richtmyer (Richtmyer 1960) introduced as shown in Fig. 19, where an incident shock propagates through a light fluid (\(i=2\)) in the \(-y\)-direction with shock speed \(-V_{Si}\) in the laboratory system of reference and fluid velocity behind the incident shock is \(-U_2\). When it hits the contact surface between the light and heavy fluids, the reflected shock starts to move in the y-direction and the transmitted one does in the \(-y\)-direction with the shock speeds of \(V_{SR}\) and \(-V_{ST}\), respectively. The contact surface moves to the \(-y\)-direction with a speed of \(-V_C\) and two fluids behind the reflected and transmitted shocks also move with the same velocity \(-U = -V_C\). The mass density of the heavy fluid and that of light fluid after the shock interaction are denoted as \(\rho _1\) and \(\rho _2\), respectively.

Shock flow diagram for a reflected shock case in a laboratory frame. A incident shock propagates from the right to the left with shock speed of \(-V_{Si}\), where \(-U_2\) is the fluid velocity behind the incident shock, \(V_{SR}\) and \(-V_{ST}\) are the reflected and transmitted shock speeds, respectively, and \(-V_C\) is the speed of the contact surface after the shock interaction. Here, \(\rho _1\) and \(\rho _2\) are the mass densities behind the transmitted and reflected shocks, respectively

Full size image

We assume the initial ripple of the interface to be of the form \(\psi _0\cos kx\). Hence, just after the shock refraction at \(t = 0+\), the reflected and transmitted shock fronts are deformed, and those of the initial ripples resemble that at the contact surface, but with different amplitudes. A simple kinematic calculation shows that the amplitudes of initial shock ripples are given by (Richtmyer 1960)

for the reflected and transmitted fronts, respectively. As is well known, the tangential velocity must be continuous across any deformed shock fronts. For the geometry shown in Figs. 1 and 19, the initial tangential velocities can be equal to

These two tangential velocities provide the initial circulation that is distributed along the contact surface ripple proportional to \(\sin kx\), which is actually the initial cause of the growth of the contact surface ripple (Fraley 1986; Samtaney and Zabusky 1993, 1994; Velikovich 1996; Wouchuk and Nishihara 1997; Wouchuk and Cavada 2004; Herrmann et al. 2008). As the shocks separate away from the interface, their corrugation amplitudes oscillate and decrease with time. As the shock ripples evolve in time, pressure perturbations are generated behind the shocks and propagate through the compressed fluids with the local sound speed. As a result, the interface ripple and the shock ripples interact with each other, and they oscillate in time and asymptotically tend to zero. For an ideal gas EOS, the pressure perturbation at the interface decreases asymptotically in time as \(t^{-3/2}\) (Zaidel 1960; Fraley 1986; Wouchuk and Cavada 2004) independent with the shock intensity. It should also be noted that the rippled shock leaves vorticity behind in the bulk, but they are the second order and also decays with time (Wouchuk96). Therefore, no baroclinic effects are required for RMI.

The asymptotic interface growth rate can be obtained within a linear theory and a weak shock limit by integrating tangential component of equation of motion

where we assumed that \(\delta p_i = \delta {\tilde{p}}_i(y,t)\cos kx\) and \(\delta {v}_{xi} = \delta {\tilde{v}}_{xi}(y,t) \sin kx\). This equation is valid at both sides of the contact surface at any time \(t > 0+\). From the pressure continuity at \(y = 0\) and the integration in the time interval \(0+< t < \infty \), we obtain a relation for the asymptotic tangential velocity \(\delta {\tilde{v}}_{xi}^\infty \) in fluid \(i = 1\) or 2 as

Since the normal component of fluid velocity is continuous at the interface and the velocity becomes irrotational asymptotically in a weak shock limit, the linear growth rate \(v_{\mathrm{lin}}\) can be obtained from Eq. (12.4) as

It should be noted that the growth rate becomes smaller than the above growth rate for a strong shock because the bulk vorticity left by the ripple shocks suppresses the growth (Wouchuk and Nishihara 1997). The pressure perturbation at the interface decays with time as \(t^{-3/2}\), independent of the shock Mach number. As a consequence, the plasma velocity perturbations become incompressible and density perturbations approach zero for a very large time, as assumed in the present current-vortex sheet model. The potential flows assumed in the bulk; however, generate the velocity shear at the interface.

Reprints and permissions