Simulation of High-Current Intersecting Plasma Beams by MHD and Monte Carlo Methods

This work was carried out in collaboration between all authors. Authors LZ and SE simulated the model by MC and MHD methods, respectively, performed the statistical analysis, and wrote the first draft of the manuscript and managed literature searches. Authors HAG and DB managed the analyses of the model and simulations. Author EB provided MCNPX software and gave some advice. All authors read and approved the final manuscript. ABSTRACT This paper deals with the 3D time-dependent intersecting plasma beams model using Magnetohydrodynamics and Monte Carlo methods under the conditions of high pressures (from 0.01 MPa up to 0.1 MPa) and high current (100 kA). After the detailed presentation of model, two methods have been fully analyzed in terms of intersecting plasma beams properties in the focal region. Here, we have compared the results of MHD time-dependent numerical simulation with MC stochastic and statistical particles simulation. Through success of these comparisons, we have demonstrated that MHD and MC methods provide practical tools to capture essential physics of intersecting plasma beams.


INTRODUCTION
A plasma is a collection of electrons and ions in the fourth fundamental state [1]. Plasmas can be regarded as a truly multi-physics phenomenon, which generally exhibit strongly coupled interaction between electron and ion transport, electromagnetic fields, heat transfer, and fluid mechanics [2]. Due to the complexity of plasma physics, experiences are often insufficient to predict plasma behavior to the desired level of accuracy and detailed experimental characterization of plasmas is typically difficult to be achieved due to the extreme operating conditions. For these two reasons, computational models that can provide insight into the plasma dynamics in a complex nonlinear system and quantify effects of design changes are valuable to understand factors. As a result, computational modeling has been increasingly employed in the design of a variety of plasma-based technologies.
Most of MHD models have focused on atmospheric or sub-atmospheric pressure plasma discharges and under high-current conditions [3][4]. Besides, Lebouvier et al. [5] successfully reported MHD modelling under lowcurrent high-voltage conditions at atmospheric pressure. The low current condition implies numerical instabilities which make the model difficult to converge. In addition, the problems associated with very high-pressure conditions have made the implementation of this model even more challenging. The original of MC method for computations can back to hundreds of years ago, e.g., the famous Buffon's needle experiment. With the emergence of the modern computer in 1940s and pioneer works by John von Neumann, Stanislaw Ulam and Nicholas Metropolis, MC method brings amounts of practical applications. The famous and classical paper is said [6] mainly contributed by the hero of plasma physics, M. N. Rosenbluth. Different computational modeling techniques are used for phenomena at different spatiotemporal scales: at the smallest scales one can use MC method based on kinetic theory, while at large scales MHD method based on fluid theory can be applied. In this paper, author uses two methods to model and simulate in-sphere focal region of intersecting plasma beams. The experiment is used to demonstrate HOPE Innovations Inc's alternative approach to generate fusion energy based on the concept of high-current plasma beams passing through a common intersection point called focal region. In order to accomplish the experiment, the 3D focal region of intersecting plasma beams model is required. By analyzing and comparing simulation results, important parameters will be presented in various forms. The presented results will be compared with the corresponding experimental data to verify the model in the next steps.

Assumptions
The 3D MHD model studied is time dependent and based on following main assumptions [7]: • The plasma is treated as a single conducting fluid using lumped macroscopic variables and their corresponding hydrodynamic conservation equations. • The single-fluid equations describe the very-low-frequency and large-scale fluidlike behavior of plasma. • The spatial and temporal scales of the variations of the fluids and fields are substantially longer than the corresponding scales of the heaviest component of the plasma, i.e. ions.
The aim of MHD method is to get general behavior of intersecting plasm beam under high pressures and high current in the focal region.
The single quantities that we are interested in are plasma current (I), current density (J), magnetic field (B), pressure (P), pressure gradient (∇P).

MHD flowchart
The 1D, 2D and 3D profiles of the abovementioned parameters have been successfully modelled. The 2D profiles are constructed in the polar coordinates (r, θ). The 3D profiles are generated in the cylindrical coordinate system (r, θ, z). The flowchart of MHD simulation is shown in Fig. 1. The equations for mass density ρ , velocity v, and charge density q are obtained by summing corresponding multi-fluid equations [8] over all species:

a. Conservation of Mass
We start from the equation of mass for particles of species α and sum over all species α .
Equation (1) is the continuity equation for particles of species α . It tells us that the particle number density, mass density and charge density remain unchanged in the absence of any interaction processes which can create or annihilate particles. The equation of total mass conservation is It is convenient to introduce the total time derivative, Eulerian time derivative or Lagrangian time derivative with respect to the flow of the plasma as a whole. The equation of total mass conservation is

b. Conservation of Charge
The equation for conservation of charge is

c. Conservation of Momentum
The momentum conservation equation is On the right side of Equation (5), they mean pressure tensor, gravitational force, viscous stress, and Lorentz force, sequentially.
From Equation (1) to Equation (5), n is particle number density, v is velocity, α is species for particles, t is time, ρ is mass density, J is current density, q is charge density, g is gravitational constant, B is magnetic field, and P is pressure, respectively.

Fig. 1. Flowchart of MHD simulation
The primary sources of nuclear data are evaluations from the Evaluated Nuclear Data File (ENDF) system, and the Evaluated Nuclear Data Library (ENDL), etc., and evaluations from the Applied Nuclear Science (T-2) Group at Los Alamos. MCNP5's generalized user-input source capability allows user to specify a various source conditions without having to make a code modification. Independent or dependent probability distributions can specified for the source variables of energy, time, and position, and for other parameters such as starting cell(s). User can use MCNP5 to make various tallies related to particle current, particle flux, and energy deposition. MCNP5 tallies are normalized and printed in the output accompanied by the R, which is the estimated relative error.

B. Kinetic theory
Kinetic theory averages out microscopic information to obtain statistical and kinetic equations. When we wish to deal with a plasma with many particles but it has too many to calculate individual orbits, then we can take a statistical approach and define a distribution of identical particles [10] such that the number of particles in phase space volume 3 3 d xd v r r at time t is So the total number of particles N (where N can be very large) at time t is Various average quantities of the plasma can be calculated by integrating the distribution function over velocity. For example, the zeroth order moment is the particle number density n at position x and time t. One introduces the plasma aspects via the force term, returning to the governing Lorentz equation for a charged particle. Therefore, the Vlasov equation is From Equation (6) to Equation (8), f is distribution function, x is particle position in the sphere, v is velocity, t is time, N is total number of particles, q ' is electric charge, m is mass, E is electrical field, and B is magnetic field, respectively.

MCNP5 flowchart
In Fig. 2, there are six steps to finish the simulation of high-current plasma beams by MC method. The detailed steps are to develop specifications, geometry, material, variancereduction techniques, neutron source and output visualization.

Comparison between MHD and MC Methods
MC method is different from MHD method. Firstly, they are applied on different conditions. MC method is used to provide a microscopic description of plasma phenomena, since a plasma consists of a very large number of interacting particles. By contrast, MHD method consists in treating whole plasma as a single conducting fluid. Secondly, they are based on

Purpose
The main purpose of Stern experiment is to test formation of a high density pinched plasma beam, which will then be used to guide the construction of an apparatus in which four such beams are arranged into a balanced tetrahedral structure so that all four beams will pass through a central focal region [13][14]. According to HOPE Innovation Inc's fusion theory, the plasma at focal region will attain a stable condition that, at sufficiently high current, shall meet or exceed Lawson Criterion (density, temperature, and confinement time) which is necessary for fusion to happen.

Description of procedures
To initiate a plasma, a voltage is applied to both ends of a single carbon tube passing through the chamber, as shown in Fig. 3. A notch or a length of thinner tube is made at the center point of the carbon tube. The tube is expected to vaporize starting from the center thinned-down location, and thus a plasma arc is generated and maintained by the supplied voltage and current. When the current passing through the plasma media is high enough, a pinched plasma beam forms under the influence of the Lorentz force. The remaining halves of the tube continue to act as electrodes to feed current through the plasma beam. The high temperature plasma as well as the current continue to vaporize the ends of the carbon graphite tube and widen the gap until the tube is consumed or until its temperature can no longer increase due to the cooling effect of the terminal blocks on which the tube was mounted. The plasma beam may is extinguished when the gap between the electrodes exceeds the length.

Simulations
Objective: To obtain 1D, 2D, and 3D modelling of the beams and in-sphere focal region profiles for the plasma beams. In Table 2, the input specifications of plasma beams are shown in details.

Simulation based on MHD method
Each beam interacts with the others in an intersecting point called the focal region. The interaction between the beams in the common focal region occurs between the vector components of plasma parameters in spherical coordinates. In order to study the interaction between the beams, the self-induced force density (or the Lorentz force F = J × B) of the plasma beams has to be calculated from the vector components of J and B. The profiles of J and B are first constructed. Afterwards, the 1D & 2D profiles of Lorentz force are calculated from the J and B profiles which then are used to construct a reference plasma beam in the cylindrical coordinate system. An in-sphere is generated from the reference beam and converted from cylindrical to spherical coordinates. The in-spheres of the three remaining beams are obtained by rotating the reference in-sphere in 3D space. The common focal region is then constructed by superimposing the in-spheres. Fig. 4 shows the steps of constructing the reference in-sphere and converting it to the spherical coordinate system (r, θ, φ).

Simulation based on MC method
The Visual Editor [16] is developed to assist the user in easily displaying geometries and in the creation of MCNP5 input file. The Visual Editor has many powerful features to help the user create and display MCNP5 geometries. These features include the ability to display multiple cross-sectional views of the geometry, optional displays of the geometry in 3D using either wire mesh or ray tracing, plotting of the source, and optional displays of particle tracks during the random walks. In Fig. 5, for the (a), the green and blue sphere mean two deuterium neutron sources models and the transparent sphere means the focal region. For the (b) and (c), they are colored according to different rules. The (b) is colored by the type of materials and the (c) is colored by the number of cell.
Tally 2 represents the average surface flux tally, F2 and tally 4 represents the average cell flux tally, F4. The formulas are 2 ( , , ) where A, E p , r, t, Φ and V mean surface area (cm 2 ), particle energy (MeV), radial distance, time, particle flux and volume (cm 3 ), respectively.
In Table 3, the "value at NPS (number of particle histories) " column shows the TFC bin values of the current history, while the "value at NPS+1" column shows the results after the largest previous history has been added to the tally. The last column shows the relative change of the TFC bin values from the NPS value to the NPS+1 value. The relative error increased by 0.1435% while the figure of merit decreased by 0.2864%. One negative effect is that the VOV increased by 0.1593%, however, it is beneath the required value of 0.1.
In Table 4, MCNP5 passed 100% of the ten TFC bin statistical checks. The RE is less than 10% and the VOV is below the required 0.1 maximum and is decreasing as 1/NPS. The probability density function (PDF) slope is greater than 3and it is 10. Both indicate that the problem is sampled 7 In Table 4, MCNP5 passed 100% of the ten TFC bin statistical checks. The RE is less than 10% and the VOV is below the required 0.1 maximum and is decreasing as 1/NPS. The probability PDF) slope is greater than 3and it is 10. Both indicate that the problem is sampled adequately. These ten statistical checks do not ensure a totally reliable result and they can provide a more rigorous check of the tally reliability. For tally 2 and tally 4, the specific values of four TFC bin statistical checks are shown in Fig. 6.
Constructing and converting the reference in-sphere from cylindrical to spherical coordinate system

MHD Method
The radial profiles of plasma parameters in a single plasma beam are shown in Fig. 7. The profiles are calculated using Equations (3), (4), and (6) which are identical for all four plasma beams. These profiles are used to construct the polar profiles required for 3D modelling of plasma beams. Fig. 8 shows the polar profiles of plasma parameters of single plasma beam. The color bar in each plot represents the variation of the plasma parameter over the beam radius.
The numerical values of plasma parameters (as calculated by MATLAB codes) are listed in An example of 3D model for one of the plasma parameters is shown in Fig. 9. The model represents the Lorentz force [17] (or the force density) acting on the plasma beam. The length of model is about 10 mm which corresponds to the diameter of the focal region. This model is used for constructing a reference in-sphere which can be rotated to obtain the in-spheres of other plasma beams. Fig. 10 shows the inspheres of Lorentz force for the plasma beams after rotating the reference in-sphere [18]. The force in-spheres for the plasma beams are calculated by integrating the force density over the cylinder volume and rotating the corresponding reference in-sphere.

MC Method
In MCNP5 output tables, the results provide a description of the TFC bin checks that test the tally for its reliability. This problem dramatically illustrates the importance of the VOV and the PDF slop checks in determining the reliability of the results. The materials are assumed to be uniform throughout all spheres. Computer running is terminated when 100000 particle histories were done. The problem summary table provides an accounting of particle track, weight, and energy creation and loss. For this problem, there is total 36,726 collisions for 10,000 source histories. In Fig. 11, the weight per escaping source particle is 1.0035, meaning that the flux on the surface 1 of radius 0.5cm is approximately The unnormalized probability density for tally 2 is a log-log plot of the PDF that is shown in Fig. 12, along with the central mean (denoted by the red dash line). For this tally 2, the slope of the PDF must be greater than or equal to three in order to achieve a reliable confidence interval. The tally 2 successes this criterion indicating that the Central Limit theorem is satisfied.

Comparison between Simulations of MHD and MC Methods
In Fig. 13, two focal region modelling of intersecting plasma beams are shown. For the part (a), the lines with arrows represent the 2D profiles of Lorentz force, which come from a reference plasma beam in the cylindrical coordinate system [19]. For the part (b), two blue circles which are neutron sources represent the in-sphere part of the intersecting plasma beams in the focal region (the red circle). However, they have some similar features, for example, they regard the focal point as a sphere in the cylindrical coordinate system. The simulation results are evaluated and summarized in Table 6, illustrating the performances, running times and limitations of the two modeling methods. Note that there exist other modeling approaches in literatures as well, such as hybrid [20] MC methods for fluid and plasma dynamics, but they are not included in this table.

CONCLUSION
For the Stern experiment, the model simulates an assumed experiment based in the focal region for fusion reactor. As a first approach in terms of MHD, the initial work focuses on developing analytical models for the intersecting plasma beams. Numerical simulations are successfully conducted using MATLAB codes to obtain various plasma parameters in the plasma beams. The second approach, MC method, is used to simulate the neutron emission from the focal region of intersecting plasma beams. The neutron flux averaged over the surface 1 is calculated using F2 tally function based on the virtual sphere (cell 1), which includes one neutron point source (initial energy 2.45MeV). The data from MCNP5 output provide a description of the TFC bins checks to determine the reliability of the tally 2. The unified viewpoint by combining two different modelling methods in same conditions provides us a totally new picture.
In this new view, we can see that difficult problems naturally become the simple ones. For example, MC method is new for kinetic problems and MHD method can be seen an advanced tool for fluid problems.