High-temperature plasma kinetic simulation method and system based on moment equation solving

Abstract

The application relates to the technical field of computer simulation and numerical simulation, and discloses a high-temperature plasma kinetic simulation method and system based on moment equation solving, which comprises the following steps: initializing particle distribution and an electromagnetic field; driving the kinetic evolution of ions and electrons and synchronously updating weights; distributing particle current density and pressure to a spatial grid to obtain a global moment distribution; iteratively solving and updating the electromagnetic field in combination with the moment distribution; implicitly correcting ion weights based on a new electric field; performing data analysis and outputting results. The application innovatively uses a generalized Ohm's law based on a moment equation to iteratively solve the electromagnetic field, breaks through the limitations of traditional gyrokinetic simulation in the frequency and gradient scales, realizes high-fidelity simulation of key processes such as high-frequency compressional Alfven waves, and through the introduction of a special processing mechanism, enhances the robustness and convergence of the algorithm, and realizes high-precision and high-efficiency numerical simulation of multi-scale physical phenomena of high-parameter plasmas in a fusion reactor.

Description

Technical Field

[0001] This application relates to the field of computer simulation and numerical simulation technology, specifically to a method and system for simulating the kinetics of high-temperature plasma based on solving moment equations. Background Technology

[0002] Thermonuclear fusion is the core energy process of celestial bodies such as the sun and stars. It involves the combination of light atomic nuclei (such as deuterium and tritium) into heavy atomic nuclei under extremely high temperatures and pressures, releasing enormous amounts of energy. With its potential advantages such as abundant fuel sources, zero greenhouse gas emissions during the reaction process, and inherent high safety, controlled thermonuclear fusion is considered one of the key ways to solve future energy crises and environmental pollution problems.

[0003] In magnetically confined nuclear fusion devices (such as tokamas or stellarators), high-temperature plasma is confined by a strong magnetic field to maintain fusion conditions. However, plasma behavior is highly complex, involving various physical processes such as instabilities, turbulence, and transport. Direct experimental research is costly and presents challenges in terms of high-precision diagnostic techniques. Therefore, numerical simulation has become an indispensable technical tool for analyzing and understanding magnetically confined nuclear fusion plasma. It can reproduce experimental conditions, predict the spatiotemporal evolution of plasma, and gain a deeper understanding of the underlying physical mechanisms, thereby optimizing the physical design of the device. Through numerical simulation, researchers can assess confinement performance, reduce experimental operational risks, and accelerate the development of fusion energy.

[0004] Currently, mainstream plasma numerical simulation methods mainly revolve around two major frameworks: magnetohydrodynamic simulation and kinetic simulation, with a significant trade-off between the two in terms of physical description accuracy and computational cost.

[0005] a) Magnetohydrodynamic simulation:

[0006] Magnetohydrodynamic (MHD) models treat plasma as a continuous fluid in a three-dimensional configuration space, describing its macroscopic evolution by coupling fluid equations with Maxwell's equations. This method is computationally efficient and suitable for large-scale simulations across the entire device, making it valuable in engineering design and macroscopic stability analysis. However, the MHD approximation does not consider plasma kinetic effects (such as the finite Larmor radius effect and key physical processes like wave-particle interactions), limiting its accuracy in fusion plasma evolution analysis.

[0007] b) Kinetic simulation:

[0008] Kinetic simulations, based on kinetic equations, directly track the evolution of particle distribution functions in high-dimensional phase space. They can fully or partially describe kinetic effects such as wave-particle resonance and finite Larmor radii, but their computational complexity increases exponentially with dimensionality. These simulation methods mainly include full kinetics, cyclotron kinetics, and drift kinetics. Full kinetics encompasses all kinetic effects and has the highest physical fidelity among the three, making it suitable for studying intricate physical processes such as low-frequency turbulent transport, high-frequency wave heating, and current-driven processes. However, its computational cost is extremely high, and even with high-performance computing resources, long-term simulations at the experimental scale are difficult to achieve, severely limiting its engineering practicality. Cyclotron kinetics simulations, by averaging rapid cyclotron motion, reduce the phase space dimension from six to five, preserving key kinetic effects such as wave-particle resonance and finite Larmor radii under low-frequency perturbations, becoming a widely adopted compromise in current magnetic confinement fusion research. However, this method is based on frequencies far lower than the ion cyclotron frequency. The gyro radius is much smaller than the background gradient scale. Assumptions such as these cannot describe high-frequency fluctuation processes, and in regions with steep gradients (such as the tokamak platform region), due to... And thus, it fails. Furthermore, as a perturbation theory, it struggles to consistently handle the non-perturbation case of the evolution of the fully distributed function, and its further development faces both theoretical and numerical challenges. Drift physics describes the particle's guiding center motion as parallel motion along magnetic field lines and drift motion perpendicular to them, thereby simplifying the six-dimensional phase space to five dimensions. This method neglects the finite Larmor radius effect but does not require the perturbation field to be much smaller than the equilibrium field, and can describe non-perturbation evolution.

[0009] To overcome the limitations of cyclotron kinetic simulations while avoiding the prohibitive computational costs of full-scale kinetic simulations, the development of a novel hybrid simulation algorithm that balances physical completeness and computational feasibility is imperative. To this end, we propose a hybrid simulation algorithm combining full-scale kinetic ions and drift-scale kinetic electrons. This algorithm achieves a balance between efficiency and accuracy, aiming for more precise descriptions at key physical scales: by employing full-scale kinetic treatment for ions, it can capture the non-equilibrium dynamics of ions in steep gradients and their interactions with high-frequency fluctuations without approximations; simultaneously, by employing drift-scale kinetic descriptions for electrons, it maintains computational feasibility while ensuring necessary physical accuracy. This framework is expected to fill the gap between traditional cyclotron kinetic simulations and full-scale kinetic simulations, providing a crucial tool for in-depth research on multi-scale physics problems in fusion devices. Summary of the Invention

[0010] This application provides a high-temperature plasma kinetics simulation method based on solving the moment equation, which solves the accuracy and applicability bottlenecks faced by existing simulation methods in high-parameter plasma simulation of fusion reactors.

[0011] This application is achieved through the following technical solution:

[0012] Firstly, this application provides a method for simulating the kinetics of high-temperature plasma based on solving the moment equation, comprising the following steps:

[0013] The simulation parameters are set according to the characteristics of the target physical problem, the particle phase space distribution and electromagnetic field are initialized, and the initial simulation state is obtained.

[0014] Based on the initial simulation state, the dynamic evolution of ions and electrons is driven separately and their weights are updated synchronously to obtain the updated particle dynamic state;

[0015] Based on the updated particle dynamics state, the particle current density and pressure are distributed to the spatial grid points through the shape factor, and the moment distribution on the global grid points is obtained through data synchronization.

[0016] Combining the moment distribution and the electromagnetic field data from the previous time step, the electric field is iteratively solved and updated, and the magnetic field evolution is driven according to Faraday's law to obtain a self-consistent new electromagnetic field distribution;

[0017] The ion weights are implicitly corrected based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field.

[0018] Based on the data generated during the simulation, Fourier analysis and linear growth rate calculation are performed to output the simulation results of high-temperature plasma kinetics.

[0019] A further optimized scheme is that the motion of the driving ions and electrons adopts a hybrid simulation framework of all-kinetic ions and drifting electrons, wherein the ions are described by the Vlasov equations, and the drifting motion of the electrons is described by the following equations of motion:

[0020] ;

[0021] In the formula, R is the position of the electronic conductor. Let be the parallel velocity of the electron, and b be the unit vector of the magnetic field direction. and These represent the electronic mass and the ionic mass, respectively, and μ is the electron magnetic moment. To balance the magnetic field, To balance the electric field, ∇ is the gradient operator.

[0022] A further optimization involves introducing an iterative parallel electric field term during the iterative solution of the parallel Ohm's law. Its discrete form is defined as:

[0023] ;

[0024] In the formula, For grid point positions, The initial number of simulated particles allocated to each grid. For the number of iterations, and Let be the position and parallel velocity of the j-th simulated electron at time step n+1, respectively, and S be the shape factor function.

[0025] A further optimized solution is that, in the step of updating particle weights, the explicit promotion of ion weights is achieved through the following formula:

[0026] ;

[0027] In the formula, The ion weights after explicit promotion. The weight of the j-th simulated ion, Let n be the time step size and n be the time step index. and These are the electron temperature and the ion temperature, respectively. , They are respectively Xianghe To the ion velocity component, , , , and y, respectively, represent the z-direction component and y-direction component of the perturbation electric field; , These represent the x-direction component and z-direction component of the perturbation magnetic field, respectively. For ion positions, This represents the ion density.

[0028] A further optimization scheme is that the implicit correction of the ion weights is achieved through the following formula:

[0029] ;

[0030] In the formula, For the (n+1)th time step after implicit correction, The weights of each simulated ion. The ion weights after explicit promotion. Let be the vertical velocity of the ion at time step n+1. Let be the vertical component of the perturbation electric field at time step n+1.

[0031] A further optimization scheme is to use a spatial domain decomposition strategy for parallel computing, dividing the simulation region into intervals along the magnetic field direction (z-axis), with each process managing one interval, and dynamically remapping based on the new position of the particle after it is pushed.

[0032] A further optimization scheme is that the data synchronization includes: merging and summing the current density and pressure of overlapping grids in adjacent processes, and processing the boundary grid point data to meet the physical boundary conditions.

[0033] A further optimization scheme is that, in the iterative solution of the field equation, the error norm of the electric field update value is calculated after each iteration, and the iteration terminates when the error norm is less than a set threshold.

[0034] Secondly, this application provides a high-temperature plasma kinetic simulation system based on solving moment equations, used to implement the high-temperature plasma kinetic simulation method based on solving moment equations as described above; the system includes:

[0035] The initialization module is used to set simulation parameters according to the characteristics of the target physical problem, initialize the particle phase space distribution and electromagnetic field, and obtain the initial simulation state.

[0036] The particle propulsion module is communicatively connected to the initialization module and is used to propel the dynamic evolution of ions and electrons respectively and synchronously update their weights based on the initial simulation state, thereby obtaining the updated particle dynamic state.

[0037] The moment calculation and synchronization module is communicatively connected to the particle propulsion module. It is used to distribute the particle current density and pressure to the spatial grid points through the shape factor according to the updated particle dynamics state, and obtain the moment distribution on the global grid points through data synchronization.

[0038] The electromagnetic field evolution module is communicatively connected to the moment calculation and synchronization module. It is used to combine the moment distribution and the electromagnetic field data of the previous time step to iteratively solve the parallel and perpendicular Ohm's laws to update the electric field, and to drive the magnetic field evolution according to Faraday's law to obtain a self-consistent new electromagnetic field distribution.

[0039] The implicit weight correction module is communicatively connected to the electromagnetic field evolution module and is used to implicitly correct the ion weights based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field.

[0040] The data analysis and output module is communicatively connected to the implicit weight correction module. It is used to perform Fourier analysis and linear growth rate calculation based on the data generated during the simulation process, and output the high-temperature plasma kinetic simulation results.

[0041] A further optimization scheme is that the system is based on a spatial domain decomposition parallel architecture, in which each process in the particle propulsion module and electromagnetic field evolution module manages a local grid interval, and realizes dynamic load balancing remapping of particle data and cross-process synchronous communication of grid data through the message passing interface MPI.

[0042] Compared with the prior art, this application has the following advantages and beneficial effects:

[0043] By constructing a hybrid simulation framework that is self-consistently coupled between fully kinetic ions and drifting kinetic electrons, and innovatively employing the generalized Ohm's law based on the moment equation for iterative solution of the electromagnetic field, this method breaks through the limitations of traditional cyclotronic simulation in terms of frequency and gradient scale in terms of physical model. It can simulate key physical processes such as high-frequency compressed Alfvén waves and ion acoustic waves with high fidelity. In terms of numerical calculation, this method effectively overcomes the numerical cancellation problem in solving the field equation by introducing special processing mechanisms such as iterative parallel electric field terms, which significantly enhances the robustness and convergence of the algorithm. Ultimately, it achieves high-precision and high-efficiency numerical simulation of multi-scale physical phenomena of high-parameter plasmas in fusion reactors, providing a more advanced and reliable computational tool for magnetic confinement fusion research. Attached Figure Description

[0044] To more clearly illustrate the technical solutions of the exemplary embodiments of this application, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of this application and should not be considered as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort. In the drawings:

[0045] Figure 1 A flowchart of a high-temperature plasma kinetic simulation method based on solving moment equations provided in this application embodiment;

[0046] Figure 2 Another flowchart of the high-temperature plasma kinetic simulation method based on solving the moment equation provided in this application embodiment;

[0047] Figure 3 The diagram shows the real-frequency simulation results of compressed Alfvén waves provided in the embodiments of this application;

[0048] Figure 4 The following are the actual frequency simulation results of left-handed and right-handed spiral waves provided in the embodiments of this application;

[0049] Figure 5 The diagram shows the actual frequency simulation results of ion acoustic waves provided in the embodiments of this application;

[0050] Figure 6 The diagram shows the simulation results of the damping rate of the ion acoustic wave provided in the embodiments of this application. Detailed Implementation

[0051] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the embodiments and accompanying drawings. The illustrative embodiments and descriptions of this application are only for explaining this application and are not intended to limit this application.

[0052] First, some of the technical terms used in this application will be explained to help those skilled in the art understand this application.

[0053] Boris: Boris algorithm;

[0054] IAW: Ion Acoustic Wave;

[0055] MPI: Message Passing Interface.

[0056] Firstly, such as Figure 1 and Figure 2 As shown, this application provides a method for simulating the kinetics of high-temperature plasma based on solving the moment equation, including the following steps:

[0057] Step S1: Set simulation parameters according to the characteristics of the target physical problem, initialize the particle phase space distribution and electromagnetic field, and obtain the initial simulation state;

[0058] Step S2: Based on the initial simulation state, drive the dynamic evolution of ions and electrons respectively and update their weights synchronously to obtain the updated particle dynamic state;

[0059] Step S3: Based on the updated particle dynamics state, the particle current density and pressure are distributed to the spatial grid points through the shape factor, and the moment distribution on the global grid points is obtained through data synchronization;

[0060] Step S4: Combining the moment distribution and the electromagnetic field data from the previous time step, iteratively solve and update the electric field, and drive the magnetic field evolution according to Faraday's law to obtain a self-consistent new electromagnetic field distribution;

[0061] Step S5: Implicitly correct the ion weights based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field.

[0062] Step S6: Based on the data generated during the simulation, perform Fourier analysis and linear growth rate calculation, and output the high-temperature plasma kinetic simulation results.

[0063] This embodiment provides a novel kinetic simulation method for ions and electrons based on a fully kinetic approach. The core idea is to employ differentiated kinetic descriptions for ions and electrons, and to achieve efficient solution of the field equations through self-consistent coupling of kinetic equations and fluid moment equations. Ions are described using a fully kinetic approach to capture their complete dynamic behavior without approximations, particularly the physical processes related to high-frequency cyclotron motion. Electrons are described using a drift-motion approach, significantly improving computational efficiency while maintaining key physical accuracy. Based on this, the algorithm solves for the electric field by self-consistently coupling Maxwell's equations with the generalized Ohm's law derived from electron drift-motion kinetics, and updates the magnetic field according to Faraday's law, thus constructing a physically more complete and numerically self-consistent simulation framework.

[0064] In one embodiment, step S1: setting simulation parameters according to the characteristics of the target physical problem, initializing the particle phase space distribution and electromagnetic field, and obtaining the initial simulation state, is specifically implemented as follows:

[0065] Step S11: Based on the characteristics of the target physical problem, set the following key parameters of the system:

[0066] Geometric parameters: Define the simulation region in three orthogonal directions Grid division (number of grids respectively) );

[0067] Particle parameters: Determine the initial number of simulated particles allocated within each grid. ;

[0068] Time parameters: Set the time step used in the simulation. and total number of simulation steps ;

[0069] Physical parameters: Based on the set equilibrium distribution, initialize the steady-state density gradients of electrons and ions. With temperature gradient To effectively suppress aliasing errors in numerical simulations, appropriate wavenumber filtering conditions should be set to filter out high wavenumber components that cause aliasing.

[0070] Step S12: Based on the equilibrium distribution, randomly generate electrons and ions in the phase space and initialize their weights. At the same time, assign values ​​to the initial perturbation electric field and perturbation magnetic field in the simulation region to obtain the initial simulation state.

[0071] This embodiment systematically sets geometric, particle, temporal, and physical parameters, establishing a physically consistent and numerically robust starting point for the entire simulation. Precise mesh generation and particle allocation ensure geometric fidelity and statistical reliability, while pre-setting density and temperature gradients and actively suppressing aliasing errors (wavenumber filtering) guarantee the accuracy of subsequent dynamic evolution from the outset. Finally, by generating particles based on equilibrium distribution and initializing their weights, while simultaneously setting the perturbation electromagnetic field, this step successfully constructs a complete initial phase space state that can be used to drive evolution, laying a solid foundation for high-fidelity reproduction of complex multi-scale physical phenomena in high-temperature plasmas.

[0072] In one embodiment, step S2, based on the initial simulation state, involves driving the dynamic evolution of ions and electrons respectively and updating their weights synchronously to obtain the updated particle dynamic state, specifically including the following steps:

[0073] Step S21: The ion is described using a total kinetic theory based on the Vlasov equation, and its position is explicitly driven by the Boris scheme. and speed Simultaneously, the explicit part of the ion weights is updated, as shown in the following equation:

[0074] ;

[0075] In the formula, The ion weights after explicit promotion. The weight of the j-th simulated ion, Let n be the time step size and n be the time step index. and These are the electron temperature and the ion temperature, respectively. , Ion velocities To components and To the component, , These are the z-direction component and y-direction component of the perturbation electric field, respectively. , These represent the x-direction component and z-direction component of the perturbation magnetic field, respectively. For ion positions, Ion density;

[0076] Step S22: The electron is described by a drift physics, and its centripetal motion is controlled by the following set of equations;

[0077]

[0078]

[0079] In the formula, R is the position of the electronic conductor. Let be the parallel velocity of the electron, and b be the unit vector of the magnetic field direction. and Here, represents the electron mass and ion mass, respectively, and μ is the electron magnetic moment, which is a conserved quantity. To balance the magnetic field, To balance the electric field, ∇ is the gradient operator.

[0080] The electron's centering position was determined using a second-order precision numerical scheme. parallel velocity Promote and update the electronic weights according to the following formula. This completes the advancement of the electronic dynamic state:

[0081]

[0082] In the formula, For the n+1th time step The weight of each simulated electron.

[0083] This embodiment characterizes the multi-scale dynamics of different particles in an electromagnetic field by employing differentiated descriptions of ion kinetics and electron drift kinetics. Ions are explicitly propelled by the Boris scheme to achieve position and velocity, while electrons are tracked by the guide center motion equations to track their parallel velocity and vertical drift. The weights of both are updated synchronously to accurately capture the nonlinear effects of particle-field interactions. While ensuring a complete description of ion dynamics, the electron drift-dynamics approximation significantly reduces computational complexity, and the synchronous evolution mechanism of the weights effectively maintains the charge conservation and numerical stability of the system. This provides a high-precision and efficient foundation for subsequent moment calculations and self-consistent evolution of the electromagnetic field to update the particle dynamics state.

[0084] In one embodiment, step S3: Based on the updated particle dynamics state, the particle current density and pressure are distributed to the spatial grid points using a shape factor, and the moment distribution on the global grid points is obtained through data synchronization. This specifically includes the following steps:

[0085] Step S31: Calculate the lattice particle current density and pressure; Based on the updated particle state, use the shape factor Particle information is assigned to spatial grid points.

[0086] For the ion, calculate its current density in the x and y directions. , and parallel pressure The calculation formulas are as follows:

[0087]

[0088] For electrons, calculate their vertical pressure. parallel pressure The calculation formulas are as follows:

[0089]

[0090] Step S32: After each process completes its local calculation, the current density and pressure of the overlapping grids of adjacent processes are first merged and summed to obtain a globally consistent value; then the boundary grid point data are processed to satisfy the physical boundary conditions.

[0091] This embodiment introduces a shape factor to smoothly and continuously distribute discrete particle information to a spatial grid, realizing the mapping of physical quantities from particle description to field description, effectively avoiding numerical noise that may exist in the direct point particle method. By calculating the current density and pressure components of ions and electrons in specific directions, the multi-scale contribution of different particle types to the electromagnetic field source terms is characterized. Furthermore, by using cross-process grid data merging and summation operations, the consistency of the global moment distribution under the parallel computing architecture is ensured. At the same time, the integrity of the physical model is guaranteed through boundary condition processing, thus providing a high-precision, artifact-free field source data foundation for the subsequent self-consistent evolution of the electromagnetic field.

[0092] In one embodiment, step S4: combining the moment distribution and the electromagnetic field data from the previous time step, iteratively solving the parallel and perpendicular Ohm's laws and updating the electric field, and driving the magnetic field evolution according to Faraday's law to obtain a self-consistent new electromagnetic field distribution, specifically includes the following steps:

[0093] Step S41: Calculate the x-direction component of the iterative current density y-direction component of iterative current density and iterative parallel electric field terms ;

[0094] For ions, the vertical iteration current density at the lattice point is:

[0095]

[0096] When the electronic response is approximately adiabatic, the parallel electric field term in the parallel Ohm's law... With electron pressure gradient term This can lead to strong numerical cancellation behavior, directly affecting the numerical stability of the equations. To effectively mitigate this cancellation problem, this application addresses... The term was discretized and reconstructed, and expressed as:

[0097] ;

[0098] In the formula, For grid point positions, The initial number of simulated particles allocated to each grid. For the number of iterations, and Let be the position and parallel velocity of the j-th simulated electron at time step n+1, respectively, and S be the shape factor function. This discrete form naturally includes discrete particle effects and finite grid effects, thus better matching the behavior of the electron pressure gradient at the numerical level. When the iteration number k=0, all the above iteration terms are initialized to 0. After each process completes its local calculation, data synchronization is required, and the process is similar to the data synchronization operation in step S3;

[0099] Step S42: Substitute the current density and pressure of ions and electrons, as well as the iterative terms, into the parallel and perpendicular Ohm's laws, and solve simultaneously to obtain the updated electric field value. The specific equation is as follows:

[0100] Parallel Ohm's Law:

[0101] ;

[0102] Vertical Ohm's Law:

[0103] ;

[0104] The electric field update value can be obtained by simultaneously solving the above field equations using the finite difference method or the spectral method. Subsequently, on Perform filtering and calculate the error norm between the solution and the previous iteration. .like If the value is less than a set threshold, the iteration is considered converged, and the final solution is set to... Exit the iteration loop and continue with step S43; otherwise, return to step S41 for the next iteration.

[0105] Step S43: Obtain a self-consistent new electric field distribution Subsequently, based on Faraday's law, the magnetic field evolves to obtain a newer magnetic field. :

[0106] .

[0107] This step, by introducing a discretized reconstruction of the parallel electric field term, effectively alleviates the rigid numerical behavior caused by the coupling between this equation term and the electron pressure gradient term under the adiabatic electron response, significantly improving the numerical stability of the algorithm. This reconstruction naturally incorporates discrete particle effects and finite mesh effects, enabling a more accurate match to the physical behavior of the electron pressure gradient. Simultaneously, by iteratively solving for parallel and perpendicular Ohm's laws and performing convergence checks and filtering on the updated electric field values, the self-consistency and numerical robustness of the electromagnetic field evolution process are ensured, laying a crucial foundation for ultimately obtaining a physically reasonable and numerically stable self-consistent electromagnetic field distribution.

[0108] In one embodiment, the data synchronization in steps S3 and S4 includes: summing the current density and pressure of overlapping grids of adjacent processes, and processing the boundary grid point data to meet physical boundary conditions.

[0109] In one embodiment, step S5 involves implicitly correcting the ion weights by solving the following equation regarding the implicitly corrected ion weights. The equation is implemented as follows:

[0110] ;

[0111] In the formula, For the (n+1)th time step after implicit correction, The weights of each simulated ion. For explicit promotion of the second The weights of each simulated ion. Let be the vertical velocity of the ion at time step n+1. Let be the vertical component of the perturbation electric field at time step n+1.

[0112] This embodiment constructs a discrete equation for the implicitly corrected ion weights, coupling the weight update with the ion vertical velocity and the vertical component of the perturbation electric field at time step n+1, effectively incorporating the real-time feedback effect of the electric field on the ion distribution evolution. This implicit correction scheme significantly enhances the numerical stability of the algorithm under strong electric field perturbations. Furthermore, by introducing an electron-to-ion temperature ratio factor, it accurately reflects the physical scaling relationship of the electric field's influence on ion weights under different temperature conditions. Thus, while maintaining the self-consistency of particle simulation, it ensures the coordination and consistency between weight evolution and macroscopic field equations, providing a key guarantee for obtaining physically accurate final particle distribution.

[0113] In one embodiment, step S6: Based on the data generated during the simulation, Fourier analysis and linear growth rate calculation are performed to output the high-temperature plasma kinetic simulation results. Specifically, this is implemented as follows:

[0114] Fourier analysis is performed on the disturbed electromagnetic field data to obtain its power spectrum distribution in the wavenumber-frequency domain. Simultaneously, the spatiotemporal evolution data of the electromagnetic field is written to a file for subsequent visualization. If the simulation objective is to study linear eigenmodes, the linear growth rate of instability is calculated by fitting the particle energy flux. Furthermore, the spatiotemporal distribution data of key physical quantities such as particle current density, pressure, and energy flux are output to provide a complete data foundation for comprehensive physical diagnosis and in-depth analysis of the simulation results.

[0115] The effectiveness of this method in simulating different high-frequency fluctuation phenomena is verified through three typical embodiments below.

[0116] Example 1: Compressed Alfvén Wave

[0117] This embodiment aims to verify the algorithm's ability to simulate high-frequency magnetohydrodynamic-kinetic hybrid waves propagating perpendicular to the magnetic field direction. We use compressed Alfvén waves as an example to demonstrate the ability of the high-temperature plasma dynamics simulation method based on moment equation solving provided in this application to simulate high-frequency waves. Simulation parameters are set... The real frequency of the compressed Alfvén wave obtained from particle simulation is... See attached diagram for relationship diagram. Figure 3 The red solid line represents the theoretical solution of the dispersion relation, and the blue cross indicates the particle simulation results.

[0118] from Figure 3 As can be seen, the real frequency obtained by particle simulation is consistent with the theoretical result, and it can capture the mode conversion phenomenon of compressed Alfvén waves transforming into ionic Bernstein waves.

[0119] Example 2: Left-handed and right-handed waves

[0120] This embodiment demonstrates the accuracy of the algorithm in simulating high-frequency waves (such as left-handed and right-handed waves) propagating parallel to the magnetic field direction. The algorithm described in this application has the ability to accurately simulate high-frequency waves propagating parallel to magnetic field lines; specifically, left-handed and right-handed waves are used as verification examples. Simulation parameters are set as follows: The real frequencies of the left-handed and right-handed waves obtained from particle simulations are... See attached diagram for relationship diagram. Figure 4 The red solid line represents the theoretical solution of the dispersion relation, and the blue crosses indicate the particle simulation results. This result verifies that the high-temperature plasma dynamics simulation method based on the moment equation provided in this application can accurately calculate left-handed and right-handed waves within the specified parameter range.

[0121] Example 3: Ionized Acoustic Wave (IAW)

[0122] This embodiment verifies the effectiveness of the proposed method in classic scenarios with significant numerical cancellation problems. For particle simulation algorithms considering kinetic electrons, simulating ion acoustic waves (IAWs) is quite challenging due to the cancellation problem, especially accurately calculating the damping rate of the ion acoustic waves. The ion acoustic wave example aims to verify whether the high-temperature plasma dynamics simulation method based on solving the moment equation provided in this application can overcome the cancellation problem. The IAW parameters are selected as follows: The simulation parameters are It is important to emphasize that in order to accurately calculate the real frequency and damping rate of the ion acoustic wave, the simulation time step is very small. This significantly increases the computational cost. See particle simulation results for details. Figure 5 and 6 This indicates that the high-temperature plasma dynamics simulation method based on solving the moment equation provided in this application can accurately simulate ion acoustic waves, but the computational cost is too high. In practical applications, it is necessary to determine whether there may be a cancellation problem in the parameter space of the physical problem, and then reasonably set the simulation parameters.

[0123] Compared with the prior art, this application has the following significant advantages and positive effects:

[0124] This application employs a hybrid simulation framework combining fully kinetic ions and drifting kinetic electrons, overcoming the limitations of traditional cyclotron kinetic simulations in terms of frequency range and gradient adaptability. This model can self-consistently describe the high-frequency wave-particle interactions and physical processes in high-parameter plasmas of fusion reactors, as well as in strong gradient regions such as the slab region, significantly improving the realism and reliability of simulations of key physical phenomena.

[0125] In solving the electromagnetic field, this application constructs a self-consistent generalized Ohm's law by coupling Maxwell's equations with the electron drift kinetic velocity torque. This scheme effectively balances a complete description of electron kinetic effects with efficient numerical solutions to the field equations, achieving a high degree of compatibility between the two at the algorithmic level. This method eliminates the need to introduce electric potential and magnetic vector potential as intermediate variables, fundamentally avoiding gauge choice ambiguity and the numerical cancellation problem of Ampere's law, thus enhancing the numerical robustness and physical consistency of the algorithm.

[0126] Secondly, this application provides a high-temperature plasma kinetic simulation system based on solving moment equations, used to implement the high-temperature plasma kinetic simulation method based on solving moment equations as described above; the system includes:

[0127] The initialization module is used to set simulation parameters according to the characteristics of the target physical problem, initialize the particle phase space distribution and electromagnetic field, and obtain the initial simulation state.

[0128] The particle propulsion module is communicatively connected to the initialization module and is used to propel the dynamic evolution of ions and electrons respectively and synchronously update their weights based on the initial simulation state, thereby obtaining the updated particle dynamic state.

[0129] The moment calculation and synchronization module is communicatively connected to the particle propulsion module. It is used to distribute the particle current density and pressure to the spatial grid points through the shape factor according to the updated particle dynamics state, and obtain the moment distribution on the global grid points through data synchronization.

[0130] The electromagnetic field evolution module is communicatively connected to the moment calculation and synchronization module. It is used to combine the moment distribution and the electromagnetic field data of the previous time step to iteratively solve the parallel and perpendicular Ohm's laws to update the electric field, and to drive the magnetic field evolution according to Faraday's law to obtain a self-consistent new electromagnetic field distribution.

[0131] The implicit weight correction module is communicatively connected to the electromagnetic field evolution module and is used to implicitly correct the ion weights based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field.

[0132] The data analysis and output module is communicatively connected to the implicit weight correction module. It is used to perform Fourier analysis and linear growth rate calculation based on the data generated during the simulation process, and output the high-temperature plasma kinetic simulation results.

[0133] In one embodiment, the system is based on a spatial domain decomposition parallel architecture, wherein in the particle propulsion module and the electromagnetic field evolution module, each process manages a local grid interval and realizes dynamic load balancing remapping of particle data and cross-process synchronous communication of grid data through the message passing interface MPI.

[0134] The functions of each module in the high-temperature plasma kinetics simulation system based on the moment equation solution correspond to the steps in the embodiment of the high-temperature plasma kinetics simulation method based on the moment equation solution. Their functions and implementation processes will not be described in detail here.

[0135] Thirdly, embodiments of this application provide a high-temperature plasma kinetic simulation device based on solving the moment equation. The high-temperature plasma kinetic simulation device based on solving the moment equation can be a personal computer (PC), a laptop computer, a server, or other device with data processing capabilities.

[0136] In this embodiment of the application, the high-temperature plasma kinetics simulation device based on solving the moment equation may include a processor, a memory, a communication interface, and a communication bus.

[0137] The communication bus can be of any type and is used to interconnect the processor, memory, and communication interface.

[0138] The communication interface includes input / output (I / O) interfaces, physical interfaces, and logical interfaces. These interfaces are used for interconnecting devices within the high-temperature plasma kinetics simulation equipment based on moment equation solving, as well as for interconnecting the high-temperature plasma kinetics simulation equipment with other devices (such as other computing devices or user equipment). Physical interfaces can be Ethernet interfaces, fiber optic interfaces, ATM interfaces, etc.; user equipment can be displays, keyboards, etc.

[0139] Memory can be various types of storage media, such as random access memory (RAM), read-only memory (ROM), non-volatile RAM (NVRAM), flash memory, optical storage, hard disk, programmable ROM (PROM), erasable PROM (EPROM), electrically erasable PROM (EEPROM), etc.

[0140] The processor can be a general-purpose processor, which can call a high-temperature plasma kinetics simulation program based on solving the moment equations stored in memory and execute the high-temperature plasma kinetics simulation method based on solving the moment equations provided in the embodiments of this application. For example, the general-purpose processor can be a central processing unit (CPU). The method executed when the high-temperature plasma kinetics simulation program based on solving the moment equations is called can refer to the various embodiments of the high-temperature plasma kinetics simulation method based on solving the moment equations in this application, and will not be repeated here.

[0141] Fourthly, embodiments of this application also provide a readable storage medium.

[0142] This application stores a high-temperature plasma kinetic simulation program based on solving the moment equation on a readable storage medium, wherein when the high-temperature plasma kinetic simulation program based on solving the moment equation is executed by a processor, it implements the steps of the high-temperature plasma kinetic simulation method based on solving the moment equation as described above.

[0143] The method implemented when the high-temperature plasma kinetics simulation program based on solving the moment equation is executed can be referred to in various embodiments of the high-temperature plasma kinetics simulation method based on solving the moment equation of this application, and will not be repeated here.

[0144] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of this application. It should be understood that the above description is only a specific embodiment of this application and is not intended to limit the scope of protection of this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of protection of this application.

Claims

1. A method for simulating the kinetics of high-temperature plasma based on solving moment equations, characterized in that, Includes the following steps: The simulation parameters are set according to the characteristics of the target physical problem, the particle phase space distribution and electromagnetic field are initialized, and the initial simulation state is obtained. Based on the initial simulation state, the dynamic evolution of ions and electrons is driven separately and the particle weights are updated synchronously to obtain the updated particle dynamic state; The motion of the propelling ions and electrons is simulated using a hybrid framework of all-kinetic ions and drift-motion electrons. Ions are described kinetically using the Vlasov equations, while the drift-motion electrons are described based on the following equations of motion: ； In the formula, R is the position of the electronic conductor. Let be the parallel velocity of the electron, and b be the unit vector of the magnetic field direction. and These represent the electronic mass and the ionic mass, respectively, and μ is the electron magnetic moment. To balance the magnetic field, To balance the electric field, ∇ is the gradient operator; Based on the updated particle dynamics state, the particle current density and pressure are distributed to the spatial grid points through the shape factor, and the moment distribution on the global grid points is obtained through data synchronization. Combining the moment distribution and the electromagnetic field data from the previous time step, the electric field is iteratively solved and updated, and the magnetic field evolution is driven according to Faraday's law to obtain a self-consistent new electromagnetic field distribution; In the iterative solution of the parallel Ohm's law, an iterative parallel electric field term is introduced. Its discrete form is defined as: ； In the formula, For grid point positions, The initial number of simulated particles allocated to each grid. For the number of iterations, and Let be the position and parallel velocity of the j-th simulated electron at time step n+1, respectively, and S be the shape factor function; The ion weights are implicitly corrected based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field. Based on the data generated during the simulation, Fourier analysis and linear growth rate calculation are performed to output the simulation results of high-temperature plasma kinetics.

2. The high-temperature plasma kinetic simulation method based on solving the moment equation as described in claim 1, characterized in that, In the step of updating particle weights, the explicit push of ion weights is achieved through the following formula: ； In the formula, The ion weights after explicit promotion. The weight of the j-th simulated ion. Where n is the time step size, and n is the time step index. and These are the electron temperature and the ion temperature, respectively. , They are respectively Xianghe To the ion velocity component, , , , and y, respectively, represent the z-direction component and y-direction component of the perturbation electric field; , These represent the x-direction component and z-direction component of the perturbation magnetic field, respectively. The position of the ion. This represents the ion density.

3. The high-temperature plasma kinetic simulation method based on solving the moment equation according to claim 2, characterized in that, The implicit correction of the ion weights is achieved through the following formula: ； In the formula, For the (n+1)th time step after implicit correction, The weights of each simulated ion. The ion weights after explicit promotion. Let be the vertical velocity of the ion at time step n+1. Let be the vertical component of the perturbation electric field at time step n+1.

4. The high-temperature plasma kinetic simulation method based on solving the moment equation as described in claim 1, characterized in that, Parallel computation is performed using a spatial domain decomposition strategy. The simulation region is divided into intervals along the magnetic field direction (z-axis), with each process managing one interval. After a particle is pushed, it is dynamically remapped according to its new position.

5. The high-temperature plasma kinetic simulation method based on solving the moment equation according to claim 1, characterized in that, The data synchronization includes: summing the current density and pressure of overlapping grids in adjacent processes, and processing the boundary grid point data to meet physical boundary conditions.

6. The high-temperature plasma kinetic simulation method based on solving the moment equation according to claim 1, characterized in that, The iterative solution of the field equations calculates the error norm of the updated electric field value after each iteration, and the iteration terminates when the error norm is less than a set threshold.

7. A high-temperature plasma kinetics simulation system based on solving moment equations, characterized in that, A system for implementing the high-temperature plasma kinetics simulation method based on solving moment equations as described in any one of claims 1-6; the system comprises: The initialization module is used to set simulation parameters according to the characteristics of the target physical problem, initialize the particle phase space distribution and electromagnetic field, and obtain the initial simulation state. The particle propulsion module is communicatively connected to the initialization module and is used to propel the dynamic evolution of ions and electrons respectively and synchronously update their weights based on the initial simulation state, thereby obtaining the updated particle dynamic state. The moment calculation and synchronization module is communicatively connected to the particle propulsion module. It is used to distribute the particle current density and pressure to the spatial grid points through the shape factor according to the updated particle dynamics state, and obtain the moment distribution on the global grid points through data synchronization. The electromagnetic field evolution module is communicatively connected to the moment calculation and synchronization module. It is used to combine the moment distribution and the electromagnetic field data of the previous time step to iteratively solve the parallel and perpendicular Ohm's laws to update the electric field, and to drive the magnetic field evolution according to Faraday's law to obtain a self-consistent new electromagnetic field distribution. The implicit weight correction module is communicatively connected to the electromagnetic field evolution module and is used to implicitly correct the ion weights based on the electric field in the new electromagnetic field distribution to obtain the final particle weights that take into account the influence of the electric field. The data analysis and output module is communicatively connected to the implicit weight correction module. It is used to perform Fourier analysis and linear growth rate calculation based on the data generated during the simulation process, and output the high-temperature plasma kinetic simulation results.

8. The high-temperature plasma kinetic simulation system based on solving the moment equation as described in claim 7, characterized in that, The system is based on a spatial domain decomposition parallel architecture. In the particle propulsion module and the electromagnetic field evolution module, each process manages a local grid interval and realizes dynamic load balancing remapping of particle data and cross-process synchronous communication of grid data through the message passing interface MPI.

Images