Explicit high-order noncanonical symplectic algorithms for ideal two-fluid systems
An explicit high-order noncanonical symplectic algorithm for ideal two-fluid systems is developed. The fluid is discretized as particles in the Lagrangian description, while the electromagnetic fields and internal energy are treated as discrete differential form fields on a fixed mesh. With the assistance of Whitney interpolating forms Whitney (1957); Desbrun et al. (2008); Xiao et al. (2015a), this scheme preserves the gauge symmetry of the electromagnetic field, and the pressure field is naturally derived from the discrete internal energy. The whole system is solved using the Hamiltonian splitting method discovered by He et al. He et al. (2015a), which was been successfully adopted in constructing symplectic particle-in-cell schemes Xiao et al. (2015a). Because of its structure preserving and explicit nature, this algorithm is especially suitable for large-scale simulations for physics problems that are multi-scale and require long-term fidelity and accuracy. The algorithm is verified via two tests: studies of the dispersion relation of waves in a two-fluid plasma system and the oscillating two-stream instability.
The ideal two-fluid model, a basic non-dissipative model of plasma physics, has been widely used to study fusion and astrophysical plasmas. In this model, the electrons and ions are treated as ideal fluids separately, with coupling to the electromagnetic fields through the charge and current carried by them. Although this system is easily generalized to any number of different charged species, the terminology “two-fluid” will be used here in lieu of “multi-fluid”, as is typically done. Because the ideal two-fluid system has noncanonical Hamiltonian form Morrison (1982, 1998), as was shown in Ref. Spencer and Kaufman (1982), its dynamics preserves geometric structure and there is no dissipation of invariants such as the total energy and momentum in the system. Conventional numerical algorithms for the ideal two-fluid system generally do not preserve geometric structure and thus the truncation error can accumulate coherently over simulation time-steps. This is a serious drawback when solving most electron-ion systems whose behaviors are naturally multi-scale. For example, the ion cyclotron period is thousands of times longer than that of the electron.
Symplectic methods, discovered in the 1980s Lee (1982); Ruth (1983); Feng (1985, 1986); Lee (1987); Veselov (1988); Yoshida (1990); Forest and Ruth (1990); Channell and Scovel (1990); Candy and Rozmus (1991); Tang (1993); Sanz-Serna and Calvo (1994); Shang (1999); Marsden and West (2001); Hairer et al. (2002); Feng and Qin (2010), have proven to be efficient for solving finite-dimensional canonical Hamiltonian systems. Such methods preserve the symplectic geometric structure (2-form) associated with the original canonical Hamiltonian system, and the numerical error of all invariants can be globally bounded by small values throughout simulations Hairer et al. (2006). In plasma physics, accelerator physics, and fluid dynamics, many of the finite-dimensional Hamiltonian systems and most of the infinite-dimensional Hamiltonian systems are noncanonical; for example, this is the case for guiding center dynamics Littlejohn (1979, 1981, 1983), the Euler fluid and magnetohydrodynamics (MHD) equations Morrison and Greene (1980), the Vlasov-Maxwell and Vlasov-Poisson systems Morrison (1980, 1982); Weinstein and Morrison (1981); Marsden and Weinstein (1982); Morrison (1998), and drift and gyrokinetic theories Morrison (2013); Squire et al. (2013); Morrison et al. (2013); Burby et al. (2015); Brizard et al. (2016). The development of geometric algorithms for these systems can be challenging. However, recently significant advances have been achieved in the development of structure preserving geometric algorithms for charged particle dynamics Qin and Guan (2008); Qin et al. (2009); Guan et al. (2010); Squire et al. (2012a); Qin et al. (2013); Liu et al. (2014); Zhang et al. (2014, 2015); Ellison et al. (2015a); He et al. (2015b, c); Ellison et al. (2015b); Liu et al. (2015); He et al. (2016), the Vlasov-Maxwell systems Squire et al. (2012b, c); Xiao et al. (2013); Kraus (2013); Evstatiev and Shadwick (2013); Shadwick et al. (2014); Xiao et al. (2015b, a); Crouseilles et al. (2015); Qin et al. (2015); He et al. (2015a); Qin et al. (2016); Webb (2016), compressible ideal MHD Zhou et al. (2014, 2016), and incompressible fluids Pavlov et al. (2011); Gawlik et al. (2011). All of these methods have demonstrated unparalleled long-term numerical accuracy and fidelity compared with conventional methods. As a side note, we point out that for infinite-dimensional Hamiltonian systems, an alternative viewpoint is to treat them as multi-symplectic systems Bridges (1997); Marsden et al. (1998), and corresponding multi-symplectic algorithms Reich (2000); Sun and Qin. (2000); Bridges and Reich (2001); Wang and Qin (2001); Hong and Qin (2002); Chen et al. (2002); Guo and Wu (2003) have also been developed.
In the present work, an explicit, high-order, noncanonical symplectic algorithm for integrating the compressible ideal two-fluid system is developed. We discretize the fluid as particles in the Lagrangian description, which naturally guarantees conservation of the density. The electromagnetic fields and internal energy are discretized over a cubic mesh by using the theory of discrete exterior calculus (DEC) Hirani (2003). High-order Whitney interpolating forms Xiao et al. (2015a) are used to ensure the gauge symmetry of Maxwell’s equations. The discrete Poisson bracket for the ideal two-fluid system is obtained by the similar technique that is used in obtaining the discrete Vlasov-Maxwell bracket Xiao et al. (2015a), and the final numerical scheme is constructed by the powerful Hamiltonian split method He et al. (2015a, c); Xiao et al. (2015a). We note that for the existing structure preserving method for the compressible fluid Zhou et al. (2014), all fields are discretized over a moving mesh, which does not apply to cases where the mesh deforms significantly during the evolution, such as in a rotating fluid. This difficulty is overcome by using a fixed mesh rather than a moving one for discretizing the electromagnetic and internal energy fields in our method. The conservation of symplectic structure guarantees that the numerical errors of all invariants such as the total energy and momentum are bounded within a small value during the simulations Hairer et al. (2006). Therefore, this method is most suitable for solving long-term multi-scale problems.
The paper is organized as follows. In Sec. II the Hamiltonian theory of the ideal two-fluid system is reviewed and the geometric structure preserving method is developed. Two numerical examples, the dispersion relation of waves in an ideal two-fluid system and the oscillating two-stream instability, are given in Sec. III. Finally, in Sec. IV we conclude.
Ii Structure preserving discretization for ideal two-fluid systems
The starting point of our development is the Lagrangian of the ideal two-fluid system, written in terms of Lagrangian variables, which is quite similar to the Lagrangian for the Vlasov-Maxwell system except for the addition of internal energy terms (see e.g. Charidakos et al. (2014)). This Lagrangian is given as follows:
where , , and are the mass, charge, and initial number density distribution of species , respectively, and are current position and velocity of fluid elements for species labeled by , which we take to be the initial value of in the configuration space, is the Jacobian of the coordinate transformation from the initial value to , is the internal energy per unit mass for species , and is the electromagnetic vector potential. In the arguments of the fields and we suppress the time variable. In this Lagrangian, we have ignored the entropy term in the internal energy, assuming barotropic fluids, and adopted the temporal gauge with . The permittivity and permeability are set to unity for simplicity.
Evolution equations are obtained upon variation of the action as in Hamilton’s principle
giving rise the equations of motion via
where and are the electromagnetic fields. These equations are exactly the ideal two-fluid equations in the Lagrangian variable description.
Now we discretize the Lagrangian using a method very similar to that for the discretization of the Vlasov-Maxwell system in Ref. Xiao et al. (2015a). The electromagnetic fields and internal energy are sampled over a cubic mesh, while the fluid is discretized into finite-sized smooth particle Monaghan (1992); Xiao et al. (2013, 2015a) moving between mesh grids. Modeling fluids using a set of Lagrangian particles is also the key idea of the smoothed-particle-hydrodynamics (SPH) method Monaghan (1992); Bonet and Lok (1999); Price and Monaghan (2004). However, the difference is that our internal energy fields are calculated on fixed mesh grids. Therefore, the method developed in the present study more closely resembles the structure-preserving symplectic particle-in-cell method of Ref. Xiao et al. (2015a) The resulting discrete Lagrangian is
Here, the subscript denotes the -th particle of species , and are Whitney interpolating maps for discrete 0-forms and 1-forms Hirani (2003); Whitney (1957); Desbrun et al. (2008); Xiao et al. (2015a), is discrete internal energy per unit volume for species , is the discrete curl operator that is defined in Eq. (15), are indices for the discrete 0-form, 1-form, 2-form, respectively. To simplify the notation, the grid size has been set to unity. The Whitney maps are defined as follows:
where the one-dimensional interpolation function is chosen in this paper to be
The equations of motion arising from the action with Lagrangian of (6) are the following:
Next we introduce two discrete fields and , which are discrete electromagnetic fields. We will make use of the following properties of the interpolating forms Xiao et al. (2015a); Hirani (2003); Whitney (1957),
The continuity equations for the densities are automatically satisfied, as can be shown by directly calculating the time derivative of ,
where can be viewed as the discrete momentum density over mesh grids, and is a discrete version of . So, Eq. (23) is essentially a kind of discrete continuity equation.
To construct the geometric structure preserving algorithm, the Hamiltonian theory for the discretized system is considered. Note that the only difference between the two-fluid Lagrangian and the Vlasov-Maxwell Lagrangian is the internal energy term, which can be written as a function of . Thus, the discrete Poisson structure of the ideal two-fluid system can be chosen to be the same as that for the Vlasov-Maxwell system Xiao et al. (2015a), which is
And the two-fluid Hamiltonian is
It turns out that the exact solutions for all sub-systems can be found and computed explicitly. The exact solutions for , , , and have been derived in Ref. Xiao et al. (2015a). They are
The solutions and are similar to . For , the exact evolution equations are
Using the property Eq. (11) of the Whitney interpolating forms, the exact solution can be written as
Here, can be interpreted as a discrete version of the continuous Newton’s second law,
At first look, it seems different from Newton’s law in the Lagrangian form derived in Ref. Morrison (1998), i.e.,
This is because the in Eq. (43) is defined as the internal energy per mass. Upon letting , the relation between and is
The final geometric structure-preserving scheme can be constructed from these exact solutions. For example, a first-order scheme can be chosen as
and a second-order scheme can be constructed as
The th-order scheme can be derived from the th-order scheme by using
Iii Numerical examples
To verify the practicability of our explicit high-order noncanonical symplectic algorithm for ideal two-fluid systems, we apply it to two physics problems. In the first problem we examine the dispersion relation of an electron-deuterium plasma, while the second concerns the oscillation two-stream instability.
For the electron-deuterium plasma, parameters of the unperturbed uniform plasma are chosen as follows,
where , is the speed of light in the vacuum, is the constant external magnetic field. This plasma supports both electron waves and ion waves, and their frequencies are very different since the deuterium ion is much heavier than the electron. The simulation is carried out in a mesh, and the periodical boundary condition is adopted in all directions. The grid size is choosen to be , the time step is set to . The simulation is initialized with stationary fluid particles being equally spaced with 4 particles per grid cell.
To numerically obtain the dispersion relation, the simulation is carried out with a small random perturbation. The space-time dependence of one field component is transformed into space, and a contour plot the field component in the space is used to make correspondence with the linear dispersion relation of the discrete system. In Fig. 1, such a contour plot is compared with the theoretical dispersion relation in both high frequency and low frequency ranges. We can see that the dispersion relation obtained by our geometric two-fluid algorithm agrees very well with the theory over the frequency range of the simulation. As expected, the total energy of the system is bounded to be within a interval of its intitial value of for all simulation time-steps, which is plotted in Fig. 2.
The second example is the well-known oscillating two-stream instability Nishikawa (1968); Morales et al. (1974). We consider the case of an unmagnetized cold two-fluid model, and compare with stability condition that was previously studied in Ref. Qin and Davidson (2014). We simulate an electron-positron plasma, with system parameters given as follows:
wiht the relative drift velocity between electrons and positrons chosen to be . The simulation domain is a mesh. Initial perturbations with two different wave numbers, and , are tested. According to the theoretical prediction of Ref. Qin and Davidson (2014), the mode with is stable while that with is unstable. Both of these predictions are confirmed by the simulation using our algorithm, as seen in Fig. 3. The evolution of the perturbed electrostatic field of the unstable mode is plotted in Fig. 4, which displays the space-time dependence of during the nonlinear evolution of the instability.