Three-dimensional decaying magnetic field
belonging to Beltrami flow

Antuono (2020) obtained a periodic three-dimensional analytic solution shown below using the method of Ross Ethier and Steinman (1994). Subscripts 1 and 2 represent two solutions.

	$\displaystyle u_{1,2}$	$\displaystyle=$	$\displaystyle\alpha U\left[\sin(kx+\theta_{1.2})\cos(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right.$		(6)
			$\displaystyle\quad\left.-\cos(kz+\theta_{1.2})\sin(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right]e^{-3\nu k^{2}t},$
	$\displaystyle v_{1,2}$	$\displaystyle=$	$\displaystyle\alpha U\left[\sin(ky+\theta_{1.2})\cos(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right.$		(7)
			$\displaystyle\quad\left.-\cos(kx+\theta_{1.2})\sin(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right]e^{-3\nu k^{2}t},$
	$\displaystyle w_{1,2}$	$\displaystyle=$	$\displaystyle\alpha U\left[\sin(kz+\theta_{1.2})\cos(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right.$		(8)
			$\displaystyle\quad\left.-\cos(ky+\theta_{1.2})\sin(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right]e^{-3\nu k^{2}t},$
	$\displaystyle p_{1,2}$	$\displaystyle=$	$\displaystyle-\rho\frac{|{\bf u}_{1,2}|^{2}}{2},$		(9)

where $\alpha=4\sqrt{2}/(3\sqrt{3})$. $k$ is the wavenumber and $\lambda=2\pi/k$ is the wavelength. The phase $\theta$, $\phi$, and $\psi$ are given as

$$\theta_{1}=\psi_{1}-\frac{5\pi}{6},\quad\phi_{1}=\psi_{1}-\frac{\pi}{6},\quad\psi_{1}=\cos^{-1}\left(\frac{R}{\sqrt{1+R^{2}}}\right),$$

(10)

$$\theta_{2}=\phi_{1},\quad\phi_{2}=\theta_{1},\quad\psi_{2}=\psi_{1},$$

(11)

where $R$ is a parameter, and the value excluding the singularity value $R=\pm 1/\sqrt{3}$ is set. In this study, $R=0$ is set as in the existing study Antuono (2020). The two solutions show similar distributions. The variables in the fundamental equations are non-dimensionalized using the reference values of length $\lambda$, velocity $U$, time $\lambda/U$, and pressure $\rho U^{2}$. The dimensionless exact solution is as follows:

	$\displaystyle u_{1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(kx+\theta_{1.2})\cos(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right.$		(12)
			$\displaystyle\,\left.-\cos(kz+\theta_{1.2})\sin(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re}t},$
	$\displaystyle v_{1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(ky+\theta_{1.2})\cos(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right.$		(13)
			$\displaystyle\,\left.-\cos(kx+\theta_{1.2})\sin(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re}t},$
	$\displaystyle w_{1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(kz+\theta_{1.2})\cos(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right.$		(14)
			$\displaystyle\,\left.-\cos(ky+\theta_{1.2})\sin(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re}t},$
	$\displaystyle p_{1,2}$	$\displaystyle=$	$\displaystyle-\frac{{\bf u}_{1,2}|^{2}}{2}.$		(15)

where $k=2\pi$ is the non-dimensionalized wavenumber. The Reynolds number is defined as $Re=UL/\nu$.

This flow of the three-dimensional Taylor decaying vortex belongs to the Beltrami flow. Therefore, the velocity vector ${\bf u}$ and the vorticity vector ${\bf\omega}$ are parallel, and ${\bf u}\times{\bf\omega}={\bf u}\times(\nabla\times{\bf u})=0$. In addition, in the velocity and vorticity fields of this Beltrami flow, because the relationship of ${\bf\omega}_{1,2}=\pm\sqrt{3}k{\bf u}_{1,2}$ holds, the pressure is associated with the kinetic energy $\rho{\bf u}\cdot{\bf u}/2$ and enstrophy ${\bf\omega}\cdot{\bf\omega}/2$, respectively. Stagnation point is where the pressure is $p=0$. Antuono (2020) clarified that the high-pressure region including the stagnation point has a Y-shaped structure. The stagnation point occurs at the next position.

$$(x_{l},y_{m},z_{n})=-\frac{\psi}{k}(1,1,1)-\frac{\pi}{k}(l,m,n),$$

(16)

where $l$, $m$, and $n$ are integers.

As in the two-dimensional case (Taylor, 1923), the three-dimensional Taylor decaying vortex flow is also periodic, so the volume integral of the velocity is zero. Since the total amount of pressure corresponds to the volume integral of kinetic energy, it is a constant value of $-1/2$ for inviscid fluids. The volume-integrated kinetic energy is explained later.

$$\int_{V}udV=0,\quad\int_{V}vdV=0,\quad\int_{V}wdV=0,\quad\int_{V}pdV=-\frac{1}{2},$$

(17)

where $V=1\times 1\times 1$.

The kinetic energy $K$ is given as

$$K=\frac{1}{2}|{\bf u}|^{2}.$$

(18)

The volume-integrated kinetic energy $\langle K\rangle$ and the average kinetic energy $K_{av}$ are given as

$$\langle K\rangle=\int_{V}KdV=\frac{1}{2}e^{-\frac{6k^{2}}{Re}t},$$

(19)

$$K_{av}=\frac{1}{V}\langle K\rangle=\frac{1}{2}e^{-\frac{6k^{2}}{Re}t},$$

(20)

The time when the average kinetic energy $K_{av}$ is halved is $t=Re\ln(2)/(6k^{2})$, and the time when $K_{av}$ decays to 10 % of the initial kinetic energy is $t=Re\ln(10)/(6k^{2})$.

Assuming that the velocity vector and magnetic flux density vector represent the Beltrami flow, the exact solution of the velocity and pressure of the three-dimensional Taylor decaying vortex is expressed by the equation (15) even under an applied magnetic field. Using the three-dimensional Taylor decaying vortex model proposed by Antuono (2020), the magnetic flux density in an electromagnetic fluid under an applied magnetic field is given as

	$\displaystyle B_{x1,2}$	$\displaystyle=$	$\displaystyle\alpha B\left[\sin(kx+\theta_{1.2})\cos(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right.$		(21)
			$\displaystyle\quad\left.-\cos(kz+\theta_{1.2})\sin(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right]e^{-3\nu_{m}k^{2}t},$
	$\displaystyle B_{y1,2}$	$\displaystyle=$	$\displaystyle\alpha B\left[\sin(ky+\theta_{1.2})\cos(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right.$		(22)
			$\displaystyle\quad\left.-\cos(kx+\theta_{1.2})\sin(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right]e^{-3\nu_{m}k^{2}t},$
	$\displaystyle B_{z1,2}$	$\displaystyle=$	$\displaystyle\alpha B\left[\sin(kz+\theta_{1.2})\cos(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right.$		(23)
			$\displaystyle\quad\left.-\cos(ky+\theta_{1.2})\sin(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right]e^{-3\nu_{m}k^{2}t}.$

The above equations satisfy the basic equations (3) and (4). As reference values dimensionless variables, the length is $\lambda$, velocity is $U$, time is $\lambda/U$, pressure is $\rho U^{2}$, and magnetic density is $B$. The exact solution is non-dimensionalized as follows:

	$\displaystyle B_{x1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(kx+\theta_{1.2})\cos(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right.$		(24)
			$\displaystyle\,\left.-\cos(kz+\theta_{1.2})\sin(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re_{m}}t},$
	$\displaystyle B_{y1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(ky+\theta_{1.2})\cos(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right.$		(25)
			$\displaystyle\,\left.-\cos(kx+\theta_{1.2})\sin(ky+\phi_{1.2})\sin(kz+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re_{m}}t},$
	$\displaystyle B_{z1,2}$	$\displaystyle=$	$\displaystyle\alpha\left[\sin(kz+\theta_{1.2})\cos(kx+\phi_{1.2})\sin(ky+\psi_{1.2})\right.$		(26)
			$\displaystyle\,\left.-\cos(ky+\theta_{1.2})\sin(kz+\phi_{1.2})\sin(kx+\psi_{1.2})\right]e^{-3\frac{k^{2}}{Re_{m}}t},$

where $k=2\pi$ is the non-dimensionalized wavenumber. The magnetic Reynolds number is defined as $Re_{m}=UL/\nu_{m}$.

Since this three-dimensional magnetic flux density flow is the same as the Beltrami flow, the magnetic flux density vector ${\bf B}$ and the current density vector ${\bf J}$ are parallel, and ${\bf B}\times{\bf J}={\bf B}\times(\nabla\times{\bf B})=0$. Therefore, the Lorentz force does not act, and work does not occur. Assuming that the kinematic viscosity and magnetic diffusivity are zero, the kinetic and magnetic energies are conserved in a periodic flow, and the total energy is also conserved. Therefore, it is useful to use this decaying vortex model as a benchmark test problem for verifying the validity of the energy conservation characteristics in a calculation method.

Since the flow of the three-dimensional Taylor decaying vortex is periodic, the volume integral of the magnetic flux density is zero like the velocity.

$$\int_{V}B_{x}dV=0,\quad\int_{V}B_{y}dV=0,\quad\int_{V}B_{z}dV=0,$$

(27)

where $V=1\times 1\times 1$.

The magnetic energy $M$ is given as

$$M=\frac{1}{2Al^{2}}|{\bf B}|^{2}.$$

(28)

This magnetic energy $M$ corresponds to dimensionless magnetic pressure. The volume-integrated magnetic energy $\langle M\rangle$ and the average magnetic energy $M_{av}$ are obtained as follows:

$$\langle M\rangle=\int_{V}MdV=\frac{1}{2Al^{2}}e^{-\frac{6k^{2}}{Re_{m}}t},$$

(29)

$$M_{av}=\frac{1}{V}\langle M\rangle=\frac{1}{2Al^{2}}e^{-\frac{6k^{2}}{Re_{m}}t},$$

(30)

The total energy $E_{t}$ is $E_{t}=K+M$. The volume-integrated total energy $\langle E_{t}\rangle$ and the average total energy $E_{tav}$ are obtained as follows:

$$\langle E_{t}\rangle=\frac{1}{2}\left(e^{-\frac{6k^{2}}{Re}t}+\frac{1}{Al^{2}}e^{-\frac{6k^{2}}{Re_{m}}t}\right),$$

(31)

$$E_{tav}=K_{av}+M_{av}=\frac{1}{2}\left(e^{-\frac{6k^{2}}{Re}t}+\frac{1}{Al^{2}}e^{-\frac{6k^{2}}{Re_{m}}t}\right).$$

(32)

Cross helicity $H_{c}$ is given as

$$H_{c}=\frac{1}{Al}{\bf u}\cdot{\bf B}.$$

(33)

The volume-integrated cross helicity $\langle H_{c}\rangle$ and the average cross helicity $H_{cav}$ are obtained as follows:

$$\langle H_{c}\rangle=\int_{V}H_{c}dV=\frac{1}{Al}e^{-\frac{3k^{2}}{Re}t}e^{-\frac{3k^{2}\pi^{2}}{Re_{m}}t},$$

(34)

$$H_{cav}=\frac{1}{V}\langle H_{c}\rangle=\frac{1}{Al}e^{-\frac{3k^{2}}{Re}t}e^{-\frac{3k^{2}\pi^{2}}{Re_{m}}t}.$$

(35)

For stable calculations, transport quantities such as kinetic energy must be conserved discretely in periodic flows of inviscid fluids. In addition, the analytical characteristics of the governing equation must hold discretely, and the compatibility between the conservative and non-conservative forms of convection term must be satisfied discretely. Similar to the previous research (Ham et al., 2002; Morinishi, 2009), we use the implicit midpoint law to discretize the time derivatives in the governing equations and perform time marching. The second-order central difference scheme is used for the discretization of the space derivatives. The simplified marker and cell method (Amsden and Harlow, 1970) is applied to solve the discretized equations. In this study, the Lorentz force is discretized using compact interpolation (Yanaoka, 2023), which can maintain the compatibility between conservative and non-conservative forms of the Lorentz force. We describe the Lorentz force discretization method in Appendix A.

The calculation area is a cube with one side $\lambda$. A uniform grid of $41^{3}$ is used for the calculation. For error evaluation and grid dependency verification, we also use the grids with dimensions of $11^{3}$, $21^{3}$, and $81^{3}$. Exact solutions ${\bf u}_{1}$ and $p_{1}$ are given as initial conditions, and periodic boundary conditions are set as boundary conditions.

First, we verify the validity and stability of the present calculation method. To confirm whether non-physical kinetic energy is generated, we analyse the inviscid flow at $Re=\infty$. The calculation is performed up to the time $t/(L/U)=10$ under the condition that the Courant number is CFL=0.5.

Regarding the calculation conditions of the decaying vortex, the Reynolds numbers are $Re=10$, 50, $10^{2}$, and $10^{3}$. We calculate until the average kinetic energy is half and compare the calculation result with the exact solution. For each Reynolds number, the time step is given so that the Courant number CFL is constant regardless of the grid. The Courant number for Reynolds numbers other than $Re=10$ is CFL=0.25. At $Re=10$, the time until the average kinetic energy is half is short, so the Courant number was set to CFL=0.05.

Similar to the earlier study (Antuono, 2020), this study performs a long-time calculation to confirm whether the transition to turbulent flow is observed. As in the previous study, the Reynolds number is $Re=10^{3}$. In this calculation, we also use a grid finer than the number of grid points used in the existing research. The Courant number is set to CFL=0.1 and 0.5 to verify the influence of time increments when the number of grid points is 41. In other grids, the Courant number is CFL=0.5. In addition, we investigate the effect of initial velocity disturbance on the flow transition. Similar to the previous study(Morinishi, 2009), we use the vector potential with uniform random numbers to generate an initial disturbance that satisfies the continuity equation. The disturbance levels are values corresponding to 1of the initial average kinetic energy. As an initial condition, we use the value obtained by adding this initial disturbance to the exact solution of velocity.

Next, we consider the analysis of a three-dimensional Taylor decaying vortex under an applied magnetic field. The grids and boundary conditions used are similar to the analytical model described in the previous subsection. The exact solutions ${\bf u}_{1}$, $p_{1}$, and ${\bf B}_{1}$ are given as initial conditions.

Regarding the calculation conditions, the Reynolds numbers are set to $Re=10^{2}$ and $10^{4}$, referring to the existing research on two-dimensional Taylor decaying vortex (Liu and Wang, 2001; Yanaoka, 2023). The magnetic Reynolds numbers are $Re_{m}=1$ and 50, and the Alfvén number is $Al=1$. We calculate under the condition that the Courant number CFL=0.25.

We explain the method of extracting low- or high-pressure regions. If pressure distribution is concentric around a vortex tube, we can identify the vortex tube by displaying the isosurface of the pressure. However, when the vortex tube and shear layer coexist, the pressure changes due to the two structures, so we cannot extract only the vortex tube. Because the radius of a thin vortex tube is small, we can identify the thin vortex tube by visualising a vortex tube with a large curvature. Therefore, by calculating the curvature of an equipressure surface and displaying the isosurface with a large curvature, we can identify a vortex tube with a small radius of curvature. Now we consider the case where pressure is high at the centre of a concentric circle and low at the periphery. The curvature of the pressure isosurface can be defined as follows:

$$\kappa_{p}=-\nabla\cdot\hat{{\bf n}},\quad\hat{{\bf n}}=\frac{{\bf n}}{|{\bf n}|}=\frac{\nabla p}{|\nabla p|},$$

(36)

where $\hat{{\bf n}}$ is a unit normal vector on the isosurface. Therefore, we can visualise a high-pressure region by displaying the high value of $\kappa_{p}>0$. Conversely, we can extract a low-pressure region by showing the low value of $\kappa_{p}<0$. $1/|\kappa_{p}|$ is the radius of curvature of the vortex tube.

Figure 1 shows the flow field of the analytic solution at $t=0$. The velocity vectors, streamlines, and the contour and isosurface of pressure are shown. The second invariant of velocity gradient tensor and the curvature of an equipressure surface are also displayed in an isosurface form. The red isosurface of pressure shows the distribution of dimensionless pressure $p=0.48$, and we can see the pressure field near a stagnation point. The second invariant $Q$ of the velocity gradient tensor is a quantity that expresses the magnitude relationship between the strain rate tensor and the vorticity tensor. Structures with $Q<0$ represent regions of high shear rates, where the high viscous dissipation rate of kinetic energy occurs. The second invariant of the velocity gradient tensor is a negative value of $Q=-52$ and represents a tubular high-shear region where the strain rate tensor increases. Along the direction of the vector $(1,1,1)$, the tubular structure exists a little away from the stagnation point. The high-pressure region with the stagnation point is a Y-shaped structure (Antuono, 2020). To extract the structure of a high-pressure area near a stagnation point, we calculate the curvature of the isosurface of the pressure. The green isosurface represents its curvature, and the magnitude of the curvature is $\kappa_{p}=48$. The high-pressure regions, which indicate the low-velocity areas with stagnation points, are connected in a mesh pattern, and a distorted cube structure appears. We can see that the isosurface of curvature buries the isosurface of pressure, and the cube structure represents the structure of the pressure field. The pressure is low at the centre of this cube structure, and when looking at the streamline, a swirling flow occurs so as to go around the low-pressure region. The tubular high-shear structure passes through the high-pressure region, and the pressure at the centre of the tubular high-shear structure is higher than the central pressure of the decaying vortex. No clear rotational flow occurs around the tubular high-shear structure.

Figure 2 shows the time variation of the total amount of transport quantity in an inviscid analysis. $\langle u\rangle$, $\langle v\rangle$, and $\langle w\rangle$ are the volume integrals of velocities of $x$-, $y$-, and $z$-directions, respectively, and $\langle K\rangle$ are the volume integral of kinetic energy. $\langle K\rangle_{e}$ is the exact solution. Since the flow is periodic, each total amount of the transport quantity is preserved. The total amounts of the velocities change at the level of rounding error, indicating that the conservation characteristics of the present calculation method are excellent. In addition, the kinetic energy is kept constant and consistent with the exact solution, and no non-physical energy is generated. Therefore, since this calculation method has excellent stability and does not generate non-physical kinetic energy, it is considered that this method can simulate a decaying vortex in a viscous fluid without destroying its characteristics.

Figure 3 shows the flow field of $Re=10^{2}$ at time $t/(L/U)=0.3$. The isosurface values were set so that the size of the structure represented by the isosurface was the same as the structure size shown in figure 1. The red isosurface expresses the distribution of the dimensionless pressure $p=0.24$. The pressure near the stagnation point becomes lower than the initial value shown in figure 1. The second invariant of the velocity gradient tensor is $Q=-26$, and we find that the tubular high-shear structure with the initial strength of $Q=-52$ diffuses and decays. Since the vorticity of the Taylor decaying vortex decreases with time, the magnitude of the strain rate tensor inside a high-shear structure generated by the interaction between vortices also decreases. The curvature of the pressure isosurface is $\kappa_{p}=48$, which is the same as the initial value. Over time, the cube structure in which the high-pressure regions are connected in a mesh pattern is maintained. The velocity of the swirling flow inside this cubic structure becomes slow.

Figure 4 shows the time variations of the total amounts of velocities and kinetic energy at $Re=10^{2}$. The total amount of each velocity is the level of rounding error. The kinetic energy decays over time. This calculation results well agree with the exact solution, indicating that this calculation method accurately captures the characteristics of the decaying vortex. The three-dimensional periodic decaying vortex proposed by Antuono (2020) is an effective model for investigating the conservation characteristics of calculation methods and changes in kinetic energy.

For $Re=10^{2}$, the relative error $|\varepsilon_{K}|$ of kinetic energy at the time $t/(L/U)=0.3$ after the average kinetic energy is half is shown in figure 5. The error is defined as $\varepsilon_{K}=(\langle K\rangle-\langle K\rangle_{e})/\langle K\rangle_{e}$. This error also corresponds to the relative error of pressure. As the number of grid points increases, the error decreases with the slope of $-2$, indicating that the calculation is the second-order accuracy.

Figure 6 (a) shows the time variation of the total amount of kinetic energy at each Reynolds number. Regardless of $Re$, the time variation of kinetic energy well agrees with the exact solution. Also, as $Re$ increases, the attenuation of kinetic energy slows down. Ham et al. (2002) performed a calculation in which the time progress was inverted in an inviscid periodic flow field and showed the reversibility of the calculation method. This study set a negative time step at $Re=10^{3}$ and used the result of $t/(L/U)=3$ as an initial value. Then, the time marching was performed in the opposite direction. The result is plotted in the figure. Intriguingly, we can calculate the phenomenon that the decaying vortex returns to its original state even in the viscous fluid. It was possible to calculate until the time near $t/(L/U)=1.5$. However, as time passed, the calculation diverged on the way, and the decaying vortex could not be completely restored to its original state. In the case of Courant number 0.25, it was possible to calculate up to near $t/(L/U)=1$, but the calculation diverged on the way. Although this phenomenon is physically impossible, this time-reversal simulation may be an effective tool as a method for verifying the stability and conservation characteristics of a calculation method. Figure 6 (b) shows the relative error when the Reynolds number varies. The error maintains constant regardless of $Re$. Since we can confirm the calculation accuracy by changing the number of grid points and the Reynolds number, it is considered that the three-dimensional Taylor decaying vortex model is valuable as a benchmark test.

Next, for $Re=10^{3}$, the time variation of the kinetic energy up to the time $t/(L/U)=10$ is shown in figure 7 (a) using three grids. In the analysis using the number of grid points of $41^{3}$, we varied the Courant number and investigated the influence of time increment. At $t/(L/U)=10$, the kinetic energy sufficiently attenuates. In the case of the number of grid points $11^{3}$, there is a slight difference between the calculation result and the analytic solution. On the other hand, the results obtained using finer grids well agree with the exact solution. Even if we perform the calculation over a long time, this calculated value well agrees with the exact solution, and there is no difference with time increments. In the previous study (Antuono, 2020), the difference between the calculation result and the analytic solution increased with time, and the transition to turbulent flow was observed. This study could not confirm such a transition phenomenon at $Re=10^{3}$. Figure 7 (b) shows the calculation result when an initial disturbance is added to the initial condition of the velocity. From around time $t/(L/U)=5$, the difference between the calculation result and the exact solution occurs, and the same trend as the existing research (Antuono, 2020) appears. We can see that increasing the initial disturbance accelerates the transition. For investigating the transition phenomenon of a decaying vortex, the analysis of a decaying vortex with an initial disturbance is an intriguing investigation.

Figure 8 shows the time variations of the total amounts of velocities and kinetic energy at $Re=10^{2}$. In addition, figure 9 shows the time variations of the total amounts of magnetic flux density, magnetic energy, total energy, and cross helicity at $Re=10^{2}$. Since the magnetic flux density is periodic, the total amount is conserved and analytically zero. The total amount of calculated magnetic flux density is the level of rounding error. The kinetic and magnetic energies decay with time, and the magnetic energy decays sharply. Since no Lorentz force occurs, we can confirm from the comparison of figures 4 and 8 that the kinetic energy is consistent with the result with no applied magnetic field. If a non-physical Lorentz force occurs, the attenuation of kinetic energy is overestimated. The total energy and cross helicity agree well with each exact solution. Since the magnetic Reynolds number is low, the magnetic energy decays quickly, and most of the total energy at $t/(L/U)=0.3$ corresponds to the kinetic energy.

For $Re=10^{2}$, the relative error $|\varepsilon_{E}|$ of the total energy at the time $t/(L/U)=0.3$ after the average kinetic energy is half is shown in figure 10. The error decreases with a slope of $-2$, indicating that the calculation method is the second-order accuracy.

Next, we investigate the trend of the decaying vortex when the calculation conditions vary. Figure 11 shows the time variations in the kinetic, magnetic, and total energies and cross helicity at $Re=10^{4}$ and $Re_{m}=50$. For comparison, the results of $Re=10^{2}$ and $Re_{m}=1$ are also included. This calculation result agrees well with the exact solution, and the present analysis can capture the energy decay process. Under $Re=10^{4}$ and $Re_{m}=50$, there is almost no attenuation of the kinetic energy. On the other hand, the magnetic energy, total energy, and cross helicity decay over time.

Figure 12 (a) shows the distribution of the magnetic flux density at time $t/(L/U)=0.3$. The $y$-$z$, $z$-$x$, and $x$-$y$ cross-sections show the distribution of magnetic flux densities in the $x$-, $y$-, and $z$-directions, respectively. At this time, the distinct periodicity of the magnetic flux density remains. In addition, figures 12 (b) and (c) show the magnetic flux lines, the magnetic pressure, the magnitude of the current density vector, the second invariant of the velocity gradient tensor, and the curvature of the isosurface of magnetic pressure. Here, we calculated the curvature of the isosurface of magnetic pressure to visualise the region of low magnetic pressure. The red isosurfaces show the distributions of magnetic pressure $P=0.01$ and $P=0.28$, and the light blue isosurface expresses the current density magnitude $J=7.5$, indicating the occurrence of high current density. The silver isosurface shows the second invariant $Q=-53$ of the velocity gradient tensor and represents a tubular high-shear region. The green isosurface shows the curvature $\kappa=-50$ of the magnetic pressure isosurface, and we can confirm the region of low magnetic pressure. The high magnetic pressure is shown in figure (b), and the low magnetic pressure and the magnitude of high current density are shown in figure (c). In figure (b), the magnetic pressure distribution of a distorted cubic structure appears so that the regions of low magnetic pressure are connected in a mesh pattern and surround the areas of high magnetic pressure. At this time, the attenuation of the velocity field is small, so the clear vortex structure is present. The magnetic flux density decays, but the magnetic flux lines are similar to the streamlines of the velocity field. The distribution of magnetic pressure shown in figure (c) is Y-shaped and the same as the shape of the pressure distribution in figure 3. The dimensionless magnetic pressure corresponds to the dimensionless magnetic energy. Therefore, the magnetic energy becomes high in the region, where the magnetic pressure is high. The high current density occurs in a grid pattern, and the magnetic flux lines swirl so as to surround the high current density region. In the area of high current density, the magnetic pressure, that is, the magnetic energy becomes high.

Figure 13 shows the time variations in the kinetic and total energies at $Re=10^{4}$ and $Re_{m}=50$. Over time, the difference between the result obtained using each grid and the exact solution appears, and the kinetic and total energies decay sharply. An existing study (Antuono, 2020) reported that the difference from the analytic solution suggested a transition to turbulent flow. When the flow transition occurred, the magnetic energy had sufficiently decreased. Therefore, under this condition, the influence of the magnetic field on the flow transition is considered small. As the number of grid points increases, the transition points approach the dimensionless time $t/(L/U)=16$.

The vorticity distribution at time $t/(L/U)=0$ and 16.85 is shown in figure 14. The $y$-$z$, $z$-$x$, and $x$-$y$ cross-sections show the vorticity distributions in the $x$-, $y$-, and $z$-directions, respectively. In the initial state, a large-scale vortex exists, but at $t/(L/U)=16.85$, it can be seen from the figure that it is converted to small-scale vortex structures by a non-linear effect and attenuated.

Figure 15 shows the time variation of the kinetic energy $K_{av}$ and its dissipation rate $\varepsilon_{K}$. The dissipation rate is defined as $\varepsilon_{K}=-dK_{av}/dt$. As the kinetic energy decays sharply, the dissipation rate increases and a maximum value appears. The pressure and vorticity distributions when the dissipation rate is maximum were shown in figures 14 (c) and (d). The vortex structure disappears with time due to viscous dissipation, the induced magnetic field also disappears, and the flow field asymptotically approaches the stationary state.