Method, apparatus and storage medium for kinetic simulation of fusion reactor plasma
By combining fully kinetic ions and drifting kinetic electrons in a hybrid simulation method, the problems of insufficient accuracy and limited applicability of plasma simulation in existing technologies have been solved. This method enables efficient and accurate simulation of high-parameter plasmas, and has significant advantages, especially in describing the physical processes of high-frequency waves and the slab region.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- SOUTHWESTERN INST OF PHYSICS
- Filing Date
- 2026-03-11
- Publication Date
- 2026-06-12
AI Technical Summary
Existing numerical simulation methods for plasma suffer from insufficient accuracy and limited applicability in high-parameter plasma simulations of fusion reactors. In particular, traditional cyclotron kinetic methods have limitations in simulating high-frequency wave physics and strong gradient physics in the slab region.
A hybrid simulation method combining fully kinetic ions and drifting kinetic electrons is adopted. The Vlasov equations for ions are solved using the semi-implicit δf method, while the motion behavior of electrons is described using the drifting kinetic equations. The electric field is solved self-consistently based on the vertical Ohm's law and the implicit parallel Ampere's law, and the magnetic field is updated.
It achieves a more comprehensive and self-consistent kinetic numerical simulation of high-parameter plasma physics processes, overcomes the limitations of traditional methods in simulating high-frequency waves and the sump region, improves computational efficiency and accuracy, and can accurately describe the interaction between high-frequency waves and plasma as well as the sump region physics processes.
Smart Images
Figure CN122197456A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of plasma physics research technology, specifically to a kinetic simulation method, device, and storage medium for fusion reactor plasma. Background Technology
[0002] Achieving controlled nuclear fusion is one of the key pathways to solving humanity's future energy problems. The principle of nuclear fusion mimics the mechanism of the sun's light and heat emission, causing the nuclei of light elements (such as the hydrogen isotopes deuterium and tritium) to undergo a fusion reaction under extreme high temperature and pressure conditions, thereby releasing enormous energy. Compared to nuclear fission, nuclear fusion has advantages such as abundant fuel sources, inherent high safety, and no long-lived high-level radioactive waste, and is considered an ideal, clean, and ultimate energy source. Currently, magnetic confinement fusion, especially tokamak devices, is one of the most promising technological routes to achieve controlled nuclear fusion. However, how to stably confine high-temperature plasma (hundreds of millions of degrees Celsius) within a limited magnetic field for extended periods, reduce the thermal load on the device walls, and suppress various instabilities remains a core scientific and technological challenge.
[0003] Due to the extremely high cost of nuclear fusion plasma experiments and the limited availability of high-precision diagnostic methods, numerical simulation has become an indispensable technical tool for studying and predicting the complex behavior of magnetically confined plasmas. By constructing physical models and utilizing supercomputers for large-scale calculations, numerical simulation can deeply reveal the intrinsic physical mechanisms of various phenomena in plasma (such as instabilities, turbulence, and transport), thus providing crucial theoretical basis and quantitative guidance for the design optimization of experimental devices, the prediction and selection of operating parameters, and the formulation of plasma control strategies. Therefore, developing more accurate and efficient plasma numerical simulation methods is essential for accelerating the research and development of fusion energy.
[0004] Currently, mainstream plasma numerical simulation methods are mainly based on two major frameworks: magnetohydrodynamic simulation and kinetic simulation.
[0005] a) Magnetohydrodynamic simulation
[0006] Magnetohydrodynamic (MHD) models abstract plasma as a continuous fluid in a three-dimensional configuration space, systematically describing its macroscopic evolution behavior by simultaneously solving the fluid equations and Maxwell's equations. The greatest advantage of this model is its high computational efficiency, enabling relatively rapid simulation of the macroscopic evolution of plasma across the entire device scale, making it highly suitable for engineering design and large-scale stability analysis. However, because this model is based on fluid approximations and fails to include key kinetic effects such as the finite Larmor radius effect and wave-particle resonance, this inherent limitation results in significant shortcomings in accurately describing the microscopic physics and turbulent transport of fusion plasma.
[0007] b) Kinetic simulation
[0008] Kinetic models start from the most basic kinetic equations and accurately describe plasma behavior by solving for the evolution of particle distribution functions in phase space. These simulation methods mainly include full kinetics, cyclotron kinetics, and drift kinetics. Full kinetics can completely encompass all kinetic effects, such as wave-particle resonance and the finite Larmor radius effect, thus offering the highest accuracy and physical fidelity in predicting turbulent transport, high-frequency wave heating, and current-driven processes. However, because full kinetics requires analyzing a six-dimensional phase space, even the most powerful supercomputers today struggle to simulate experimental-scale plasmas and their long-term evolution, severely limiting its practicality in engineering applications. Cyclotron kinetics, by utilizing the characteristic of charged particles rapidly cyclotronizing around magnetic field lines under strong magnetic fields, performs cyclotron averaging on the kinetic equations, reducing the phase space dimension from six to five dimensions while preserving key kinetic effects such as wave-particle resonance and the finite Larmor radius under low-frequency perturbations. However, this simulation method is based on several orders of magnitude assumptions, including that the electromagnetic perturbation frequency is much lower than the ion cyclotron frequency. And the ion cyclotron radius is much smaller than the characteristic scale length of the background plasma. These assumptions severely limit their application: on the one hand, they cannot effectively describe high-frequency wave physical processes (such as...). (wave-particle interaction); on the other hand, in regions with strong gradients, such as the tokamak pedestal region, due to The fundamental assumptions no longer hold, causing the simulation to fail. Furthermore, cyclotron kinetics is essentially a perturbation theory, requiring the perturbation to be much smaller than the equilibrium field. Developing a non-perturbation simulation that can consistently describe the full-f distribution function requires significant expansion of the existing theoretical framework, currently facing considerable theoretical and numerical challenges. In contrast, drift kinetics also simplifies the phase space to five dimensions, achieving dimensionality reduction by describing the parallel motion of particles along magnetic field lines and their drift in the perpendicular direction. While this model neglects the finite Larmor radius effect, it does not rely on the assumption that the perturbation is much smaller than the equilibrium field, thus possessing the ability to describe non-perturbation evolution.
[0009] Therefore, despite the great success of cyclotron kinetic simulations, developing a new simulation method that can both reduce computational costs and more accurately capture the physics of cyclotron motion, especially the physics related to ion cyclotron motion, remains an urgent technical need in this field. Summary of the Invention
[0010] This invention addresses the problems of insufficient accuracy and limited applicability of existing simulation methods in high-parameter plasma simulation of fusion reactors. It provides a kinetic simulation method, equipment, storage medium, and program product for fusion reactor plasma, solves the compatibility problem of efficient solution of kinetic electron and field equations, and realizes a more comprehensive and self-consistent kinetic numerical simulation of high-parameter plasma physical processes.
[0011] The present invention is achieved through the following technical solution.
[0012] Firstly, a kinetic simulation method for fusion reactor plasma is provided, the method comprising:
[0013] The simulation region is initialized, and the ion weights and electron weights are driven by a first-order implicit scheme.
[0014] Based on the simulation initialization and the grid nodes of the simulation region, the particle display and push format is used to obtain the ion position and velocity, ion weight, electron magnetic moment and velocity, and electron weight, and to determine the preset interval to which each ion or electron belongs based on the updated coordinates of each ion or electron.
[0015] Based on the ion weight, the ion position and velocity, the current at the ion position is distributed to adjacent grid nodes according to the shape factor and the preset interval to obtain the ion current density on the grid node; based on the electron weight, the electron magnetic moment and velocity, the pressure and current at the electron position are distributed to adjacent grid nodes according to the shape factor and the preset interval to obtain the electron pressure and electron parallel current density.
[0016] Using the ion current density, the electron pressure, and the electron parallel current density, the implicit parallel Ampere's law and the perpendicular Ohm's law are solved using an iterative method to obtain the electric field strength. The magnetic field is then updated based on Faraday's law and the electric field strength to obtain the magnetic induction intensity.
[0017] Using the updated electric field strength and magnetic induction intensity, the ion weights and electron weights are updated via implicit particle propulsion. If the simulation time step is less than a first threshold, the following steps are performed: Based on the simulation initialization and the grid nodes of the simulation region, the particle explicit propulsion format is used to obtain the ion position and velocity, ion weight, electron magnetic moment and velocity, and electron weight. The preset interval to which each ion or electron belongs is determined according to the updated coordinates of each ion or electron.
[0018] In some embodiments, the simulation initialization includes setting: simulation region parameters, steady-state density gradients and temperature gradients for electrons and ions, wavenumber filtering conditions, electron and ion weight initialization, and initial perturbation electric and magnetic fields.
[0019] In some embodiments, based on the simulation initialization and the grid nodes of the simulation region, a particle display push format is used to obtain ion positions and velocities, ion weights, electron magnetic moments and velocities, and electron weights. Then, based on the updated coordinates of each ion or electron, a preset interval is determined, including:
[0020] The Boris scheme is used to drive ion position and velocity, and the explicit scheme is used to drive half-steps to obtain ion weights.
[0021] The electron magnetic moment and velocity are driven by the drift motion physical equation, and the electron weight is obtained by driving half-steps using an explicit scheme;
[0022] Based on the updated z-axis coordinates of each ion or electron, determine its corresponding preset. The interval, wherein the simulation space is divided into multiple consecutive intervals along the z-axis. In a range, each of the multiple processes is assigned and responsible for tracking all particles located within one of the ranges.
[0023] In some embodiments, the electric field strength is obtained by iteratively solving the implicit parallel Ampere's law and the perpendicular Ohm's law using the ion current density, the electron pressure, and the electron parallel current density. The magnetic field strength is then updated based on Faraday's law and the electric field strength to obtain the magnetic induction intensity, including:
[0024] The magnetic flux density is obtained by pushing Faraday's law half a step based on the explicit scheme.
[0025] For ions, the ion current density is iteratively calculated at the grid nodes;
[0026] For electrons, the parallel current density and electron pressure of the electrons are iteratively calculated at the grid nodes;
[0027] If the iteration error is less than or equal to the second threshold, the ion current density, electron parallel current density and electron pressure obtained by the iteration calculation are substituted into the implicit scheme parallel Ampere's law and perpendicular Ohm's law respectively to obtain the updated electric field intensity.
[0028] Based on the updated electric field strength, Faraday's law is applied using an implicit scheme to obtain the updated magnetic flux density.
[0029] In some embodiments, the method further includes: distributing the current at the ion position to adjacent grid nodes according to the shape factor and the preset interval based on the ion weight, the ion position, and the velocity to obtain the ion current density on the grid nodes; distributing the pressure and current at the electron position to adjacent grid nodes according to the shape factor and the preset interval based on the electron weight, the electron magnetic moment, and the velocity to obtain the electron pressure and electron parallel current density; and then...
[0030] The contributions of adjacent processes on the overlapping grid are summed to ensure the continuity of the calculation of the electron parallel current density and electron pressure within the preset interval.
[0031] Boundary condition constraints are applied to the boundary grid point data of the simulation region to meet the boundary requirements of the physical model.
[0032] In some embodiments, the method further includes: when the simulation time step is greater than or equal to the first threshold, processing the updated electric field strength and magnetic induction intensity and outputting simulation results, wherein processing the updated electric field strength and magnetic induction intensity and outputting simulation results includes:
[0033] The updated electric field strength and magnetic induction intensity are subjected to Fourier transform to obtain the distribution spectrum in the spatial wavenumber and frequency domain.
[0034] The updated electric field strength and magnetic induction intensity are written to a specified output file for subsequent plotting of the spatiotemporal evolution diagram; or
[0035] If linear eigenmode research is conducted, the linear growth rate of instability can be calculated based on particle energy flux, and the ion current density, electron pressure, electron parallel current density, and energy flux can be output simultaneously.
[0036] Secondly, a kinetic simulation system for fusion reactor plasma is provided, the system comprising:
[0037] The simulation initialization module is used to: initialize the simulation region and drive the ion weights and electron weights using a first-order implicit scheme;
[0038] The particle display and propulsion module is used to: based on the simulation initialization and the grid nodes of the simulation region, adopt the particle display and propulsion format to obtain ion positions and velocities, ion weights, electron magnetic moments and velocities, and electron weights, and determine the preset interval to which each ion or electron belongs based on the updated coordinates of each ion or electron;
[0039] The particle current density and pressure calculation module is used to: based on the ion weight, the ion position and velocity, according to the shape factor and the preset interval, distribute the current at the ion position to adjacent grid nodes to obtain the ion current density on the grid node; based on the electron weight, the electron magnetic moment and velocity, according to the shape factor and the preset interval, distribute the pressure and current at the electron position to adjacent grid nodes to obtain the electron pressure and electron parallel current density.
[0040] The field equation solving module is used to: use the ion current density, the electron pressure and the electron parallel current density to solve the implicit parallel Ampere's law and the perpendicular Ohm's law using an iterative method to obtain the electric field strength; and update the magnetic field based on Faraday's law and the electric field strength to obtain the magnetic induction intensity.
[0041] The implicit particle propulsion module is used to: update the ion weights and electron weights by implicitly propulsing particles using updated electric field strength and magnetic induction intensity, and return to the explicit particle propulsion module for processing if the simulation time step is less than a first threshold.
[0042] Thirdly, a kinetic simulation device for fusion reactor plasma is provided, the device comprising:
[0043] At least one processor;
[0044] At least one memory coupled to the at least one processor and storing instructions for execution by the at least one processor, the instructions implementing the method described in any of the above when executed by the at least one processor.
[0045] Fourthly, a computer-readable storage medium is provided, the computer-readable storage medium storing instructions that, when executed by a computer, cause the computer to perform the method described in any of the preceding embodiments.
[0046] Fifthly, a computer program product is provided, the computer program product including instructions, which, when executed by a computer, cause the computer to perform the method described in any of the above embodiments.
[0047] Compared with the prior art, the present invention has the following advantages and beneficial effects.
[0048] (1. The hybrid simulation model combining fully kinetic ions and drifting kinetic electrons effectively overcomes the limitations of the traditional cyclotron kinetic method in simulating high-frequency wave physics and strong gradient physics in the base zone. It can quasi-simulate key physical processes such as the interaction between high-frequency waves and plasma and the physics of the base zone in the high-parameter plasma of fusion reactors.)
[0049] (2. To address the challenge of balancing electron kinetic effects and computational efficiency in self-consistent solutions to electromagnetic fields, an implicit parallel Ampere's law and a perpendicular Ohm's law are coupled for the solution. This approach accurately captures electron kinetic behavior while ensuring numerical efficiency in solving the field equations. Compared to existing models based on parallel Ohm's law, which require handling higher-order electron velocity moments, the implicit parallel Ampere's law only needs to handle the first-order electron velocity moment, significantly reducing the numerical complexity of the cancellation problem. In simulations of low-frequency electromagnetic turbulence and other cases with severe cancellation problems, more accurate simulation results can be obtained with the same computational resources.)
[0050] (3. A semi-implicit scheme is used to advance the ion and electron weights, and the electromagnetic field equations are uniformly processed into a second-order time-accurate scheme. While maintaining high accuracy, this algorithm effectively relaxes the restrictions on the simulation time step, thereby significantly reducing the amount of computation and improving the feasibility and efficiency of long-term, high-precision simulations.) Attached Figure Description
[0051] To more clearly illustrate the technical solutions of the exemplary embodiments of the present invention, 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 the present invention and should not be regarded 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.
[0052] Figure 1 This is a flowchart of a kinetic simulation method for fusion reactor plasma according to an embodiment of the present invention.
[0053] Figure 2 This is a diagram of a kinetic simulation algorithm for fusion reactor plasma according to an embodiment of the present invention.
[0054] Figure 3 This is a schematic diagram of the real-frequency simulation results of compressed Alfvén waves according to an embodiment of the present invention.
[0055] Figure 4 This is a schematic diagram of the actual frequency simulation results of ion acoustic waves according to an embodiment of the present invention.
[0056] Figure 5 This is a schematic diagram of the simulation results of the ion acoustic damping rate according to an embodiment of the present invention.
[0057] Figure 6 This is a schematic diagram of the ITG real-frequency simulation results according to an embodiment of the present invention.
[0058] Figure 7 This is a schematic diagram of the simulation results of the ITG growth rate according to an embodiment of the present invention.
[0059] Figure 8 This is a structural block diagram of a kinetic simulation system for fusion reactor plasma according to an embodiment of the present invention.
[0060] Figure 9 This is a schematic diagram of a kinetic simulation device for fusion reactor plasma according to an embodiment of the present invention. Detailed Implementation
[0061] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.
[0062] To address the shortcomings of existing simulation methods in simulating high-parameter plasmas in fusion reactors, such as insufficient accuracy and limited applicability, this invention proposes a hybrid particle simulation scheme based on fully kinetic ions and drifting kinetic electrons. This scheme employs a semi-implicit δf method, solving the Vlasov equations for ions and describing the motion of electrons using drifting kinetic equations. It also provides a self-consistent solution for the electric field based on the perpendicular Ohm's law and the implicit parallel Ampere's law, while updating the magnetic field using Faraday's law. This solves the compatibility problem between kinetic electrons and the efficient solution of field equations, achieving a more comprehensive and self-consistent kinetic numerical simulation of high-parameter plasma physics processes.
[0063] On the one hand, the present invention provides a kinetic simulation method for fusion reactor plasma. Figure 1 This is a flowchart illustrating a kinetic simulation method for fusion reactor plasma according to an embodiment of the present invention. (Reference) Figure 1 The method includes steps S10 to S50. Steps S10 to S50 are described in detail below with reference to the accompanying drawings.
[0064] In S10, the simulation region is initialized and the ion weights and electron weights are driven by a first-order implicit scheme.
[0065] For example, based on the characteristics of high-frequency wave physics problems, the simulation region parameters are set, including: three orthogonal directions. Number of grids , and Number of particles in a single grid Time step and total simulation time steps The steady-state density gradients of electrons and ions are set based on the equilibrium distribution. With temperature gradient To effectively suppress aliasing errors in numerical simulations, wavenumber filtering conditions need to be configured.
[0066] Based on this, electrons and ions are randomly generated according to the equilibrium distribution, and their weights are initialized. At the same time, the initial perturbation electric field and perturbation magnetic field are set.
[0067] To suppress the odd-even disconnection oscillations introduced by the second-order particle weighted driving scheme, a first-order implicit particle weighted driving scheme is adopted in the simulation initialization phase.
[0068] In S20, based on the simulation initialization and the grid nodes of the simulation region, a particle explicit propulsion format is used to obtain ion positions and velocities, ion weights, electron magnetic moments and velocities, and electron weights. The preset interval to which each ion or electron belongs is determined based on its updated coordinates. S20 involves explicit ion propulsion and explicit electron propulsion, including the following S21 to S23.
[0069] In S21, the ion position and velocity are driven using the Boris scheme, and the ion weights are obtained by driving half-steps using an explicit scheme. For example, the ion position and velocity are driven using the Boris scheme, with the specific format as follows:
[0070]
[0071] in, Indicates the steady-state electric field strength; Indicates the time step; and For Boris format intermediate step velocity auxiliary variables; Indicates steady-state magnetic flux density; Represents the normalized unit magnetic field vector; and Indicates the rotation factor; Indicates time step The corresponding speed; This represents the velocity corresponding to time step n; Indicates time step The corresponding position; This indicates the position corresponding to time step n.
[0072] Ion weights are first obtained by explicitly propelling the process half-step. :
[0073]
[0074] in, This represents the weight of the j-th simulated ion at time step n; This indicates the position of the j-th simulated ion at time step n; Indicates electron temperature; Indicates ion temperature; This represents the particle density of the i-th type of ion; , , These represent the ions' velocities in the X-axis, Y-axis, and Z-axis directions at time step n, respectively. , , These represent the components of the electric field intensity in the X, Y, and Z axes at time step n, respectively. , These represent the components of the magnetic field strength in the X-axis and Z-axis directions at time step n, respectively.
[0075] The normalization method used is:
[0076] ,
[0077] in, The subscripts represent ions respectively. and electronics speed, Characterizing the speed of sound For the electron magnetic moment, The subscripts represent the ions respectively. and electronics temperature, Indicates ion mass. Indicates the ion cyclotron frequency. This represents the ion cyclotron radius calculated based on the speed of sound. and These are the electric field strength and the magnetic induction intensity, respectively. Indicates current density, Indicates pressure. This represents the steady-state particle density. A physical quantity that has the dimension of length.
[0078] In S22, the electron magnetic moment and velocity are driven by the drift motion physics equation, and the electron weight is obtained by driving half-steps using an explicit scheme.
[0079] Electron propulsion employs the drift motion physics equation:
[0080]
[0081] in, Indicates the position of the electronic guide; Represents the parallel velocity components of the electron; Indicates electron mass; Indicates ion mass; This represents the gradient.
[0082] The actual driving scheme using second-order time precision is as follows:
[0083]
[0084] in, Indicates the first A simulated electron at time step The position of the guide center at that time; Indicates the first A simulated electron at time step The position of the guide center at that time; Representing time layer Time; Representing time layer The time.
[0085] The electron magnetic moment is a conserved quantity. The electron weight is first obtained by pushing half a step using an explicit scheme. :
[0086]
[0087] in, is the dissipation coefficient, used for the parity-disconnection sawtooth oscillations that may be encountered in the damped steady-wave example.
[0088] In S23, the preset coordinates of each ion or electron are determined based on their updated z-axis coordinates. The interval, wherein the simulation space is divided into multiple consecutive intervals along the z-axis. In a range, each of the multiple processes is assigned and responsible for tracking all particles located within one of the ranges.
[0089] After the particles are propelled, the position coordinates of all particles are updated. Then, the algorithm updates the coordinates of each particle based on its position. The axis coordinates determine the preset to which they belong. The algorithm simulates a spatial interval and sends all data of the particle to the corresponding process responsible for that interval. In the parallel architecture of this algorithm, the simulation space is along... The axis is divided into multiple consecutive Each process is assigned a range and is specifically responsible for tracking all particles located within one of these ranges, thereby achieving a static distribution of computational load among processes.
[0090] In S30, based on ion weight, ion position and velocity, the current at the ion position is distributed to adjacent grid nodes according to the shape factor and preset interval to obtain the ion current density on the grid node; based on electron weight, electron magnetic moment and velocity, the pressure and current at the electron position are distributed to adjacent grid nodes according to the shape factor and preset interval to obtain the electron pressure and electron parallel current density.
[0091] For ions, according to the shape factor The ion position current is distributed to adjacent grid points, and the specific format is as follows:
[0092]
[0093] in, This indicates the weight of the j-th simulated ion after the push; Represents the spatial coordinates of grid nodes; This represents the total number of simulated ions within each grid cell; This represents the velocity component in the d direction of the j-th simulated ion at time step n+1. ; Indicates that it is located at The component of the explicit current density of ions in the d direction; This indicates the position of the j-th simulated ion at time step n+1.
[0094] Similarly, the electron pressure and parallel current density can be obtained:
[0095]
[0096] in, This represents the explicit parallel current density of electrons. This represents the weight of the j-th simulated electron after explicit push; This represents the parallel velocity of the j-th simulated electron at time step n+1; This indicates the position of the j-th simulated electron at time step n+1; Indicates electronically displayed vertical pressure; Let represent the magnetic moment of the j-th electron.
[0097] To ensure the correctness of the field in parallel computing, the following synchronizations are required after the local particle moment allocation is completed: First, accumulate the contributions from the overlapping grids of adjacent processes to ensure the continuity of current density and pressure within the computational domain; Second, apply boundary condition constraints to the boundary grid point data of the simulation region to ensure that they meet the boundary requirements of the physical model.
[0098] In S40, the electric field strength is obtained by iteratively solving the implicit parallel Ampere's law and the perpendicular Ohm's law using the ion current density, electron pressure, and electron parallel current density. The magnetic field strength is then updated based on Faraday's law and the electric field strength to obtain the magnetic induction intensity. S40 includes S41 to S45.
[0099] In S41, the magnetic flux density is obtained by pushing Faraday's law half a step according to the explicit scheme. Faraday's law is obtained by pushing it half a step according to the explicit scheme. :
[0100] .
[0101] In S42, for ions, the ion current density is iteratively calculated at the grid nodes. For ions, the iterative current density at the grid points is:
[0102]
[0103] in, This represents the ion iteration current density in the x-direction; This represents the x-component of the electric field intensity at iteration step k.
[0104] For electrons, the iterative parallel current density and pressure at the lattice points are calculated using the following formula:
[0105]
[0106] in, This represents the parallel current density of the electron iteration; This represents the vertical pressure of electron iteration.
[0107] When the number of iterations At this point, both the iterative current and pressure are set to 0.
[0108] After each process finishes its calculation, a data synchronization operation similar to that of S30 is required.
[0109] In S43, for electrons, the electron parallel current density and electron pressure are iteratively calculated at the grid nodes.
[0110] The current density and pressure of ions and electrons and iterative current density and pressure Substituting into the field equation, the specific results are as follows.
[0111] Implicit parallel Ampere's law:
[0112]
[0113] in, and These represent the electric field strength at iteration steps k and k+1, respectively. This refers to the plasma specific pressure.
[0114] Implicit vertical Ohm's law:
[0115]
[0116] in, This represents the vertical current density during ion iteration; The gradient operator is perpendicular to the direction of the magnetic field; It is the unit vector in the x-direction.
[0117] The field equations are solved using the finite difference method or the spectral method to obtain the updated electric field. After filtering, the iteration error of the solution needs to be calculated. When the error meets the preset accuracy requirement, the solution is considered convergent. Let The iteration loop terminates, and the program flow enters S44; if the condition is not met, it returns to S43 and continues iterating.
[0118] In S44, if the iteration error is less than or equal to the second threshold, the ion current density, electron parallel current density and electron pressure obtained by the iteration calculation are substituted into the implicit scheme parallel Ampere's law and perpendicular Ohm's law respectively to obtain the updated electric field strength.
[0119] In S45, after obtaining Subsequently, the implicit driving force of Faraday's law was obtained. :
[0120] .
[0121] In S50, the ion weights and electron weights are updated by implicitly pushing the particles using the updated electric field strength and magnetic induction strength. If the simulation time step is less than the first threshold, the simulation returns to S20.
[0122] use as well as Update ion and electron weights according to the following formula:
[0123]
[0124]
[0125] in, and These represent the weights of the j-th simulated ion and electron at time step n+1, respectively. , , Let x, y, and z represent the electric field intensity at time step n+1, respectively. , These represent the x and z components of the magnetic field strength at time step n+1, respectively.
[0126] In some embodiments, the method further includes S60: when the simulation time step is greater than or equal to a first threshold, processing the updated electric field intensity and magnetic induction intensity and outputting the simulation results, wherein processing the updated electric field intensity and magnetic induction intensity and outputting the simulation results includes: performing Fourier transform on the updated electric field intensity and magnetic induction intensity to obtain the distribution spectrum of spatial wavenumber and frequency domain; writing the updated electric field intensity and magnetic induction intensity into a specified output file for subsequent plotting of spatiotemporal evolution diagrams; or, if linear eigenmode research is carried out, calculating the linear growth rate of instability based on particle energy flux, and simultaneously outputting ion current density, electron pressure, electron parallel current density, and energy flux.
[0127] The method of the present invention will be described below through specific application examples.
[0128] Example 1: Compressed Alfvén Wave
[0129] The algorithm's ability to simulate high-frequency waves is demonstrated by taking compressed Alfvén waves as an example.
[0130] Simulation parameter selection The real frequency of the compressed Alfvén wave obtained from particle simulation and See relationship diagram Figure 3 In this diagram, the red solid line represents the theoretical solution of the dispersion relation, and the black triangle represents the particle simulation results.
[0131] 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.
[0132] Example 2: Ionized Acoustic Wave (IAW)
[0133] 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 compare the ability of our algorithm to overcome the cancellation problem with other mainstream models. The IAW parameters are selected as follows: The simulation parameters are For detailed particle simulation results, please refer to [link / reference]. Figure 4 and 5 In this diagram, blue crosses represent simulation results of the field equation model composed of parallel and perpendicular Ohm's laws, while green pentagrams indicate the results of simulations using the parallel Ohm's law. Replace with The simulation results of the field equation model are shown in red. The red dots represent the simulation results of the particle simulation algorithm that uses a first-order implicit scheme to drive the particle weights and implicit parallel Ampere's law and vertical Ohm's law as the field equations. The black triangles represent the simulation results of the particle simulation algorithm proposed in this invention that uses a second-order semi-implicit scheme to drive the particle weights and implicit parallel Ampere's law and vertical Ohm's law as the field equations.
[0134] Figure 4 and Figure 5 This demonstrates that, under the same simulation parameters, this algorithm can accurately calculate both the IAW real frequency and the damping ratio simultaneously. The damping ratio calculation results, in particular, highlight the advantages of this algorithm, as the damping ratios obtained from other models are either too high or too low, with relative errors that cannot be ignored. These calculation results indicate that this algorithm has a stronger ability to overcome the cancellation problem.
[0135] Example 3: Ion Temperature Gradient Mode (ITG)
[0136] Ion temperature gradient mode ITG real frequency Particle simulation programs that incorporate kinetic electronic responses often encounter cancellation problems when calculating ITG, making accurate calculations difficult. Example 3 compares the algorithm's ability to simulate ITG with other mainstream models, using ITG parameters... The simulation parameters are The actual frequency and growth rate results of ITG obtained from particle simulation are shown in the appendix. Figure 6 and 7 x-axis The time step used in parallel Ohm's law and implicit parallel Ampere's law is... The algorithm proposed in this invention (Second order) Accuracy scheme) Simulation time step .Will Replace with The ITG obtained by applying parallel Ohm's law is a damped solution, which does not conform to physics. Therefore, in Figure 6 The simulation results for this model are not shown in the document.
[0137] Depend on Figure 6 and Figure 7 It can be seen that the algorithm proposed in this invention can obtain more accurate ITG real frequency and growth rate at a larger time step, indicating that the algorithm is easier to handle the cancellation problem and can relax the time step limitation.
[0138] In this invention, ion dynamics is described using the Vlasov equations; to reduce computational costs, electron dynamics is described using drift kinetic equations; in constructing the field equations, the vertical Ampere's law is transformed into the vertical Ohm's law to avoid direct calculation of the electron vertical current; to solve the cancellation problem, a semi-implicit scheme is used to advance the particle weights, and a parallel electric field equation is constructed based on the implicit parallel Ampere's law; the magnetic field is self-consistently updated using Faraday's law.
[0139] Compared with the prior art, the present invention has the following significant advantages and positive effects.
[0140] (1. The hybrid simulation model combining fully kinetic ions and drifting kinetic electrons effectively overcomes the limitations of the traditional cyclotron kinetic method in simulating high-frequency wave physics and strong gradient physics in the base zone. It can quasi-simulate key physical processes such as the interaction between high-frequency waves and plasma and the physics of the base zone in the high-parameter plasma of fusion reactors.)
[0141] (2. To address the challenge of balancing electron kinetic effects and computational efficiency in self-consistent solutions to electromagnetic fields, an implicit parallel Ampere's law and a perpendicular Ohm's law are coupled for the solution. This approach accurately captures electron kinetic behavior while ensuring numerical efficiency in solving the field equations. Compared to existing models based on parallel Ohm's law, which require handling higher-order electron velocity moments, the implicit parallel Ampere's law only needs to handle the first-order electron velocity moment, significantly reducing the numerical complexity of the cancellation problem. In simulations of low-frequency electromagnetic turbulence and other cases with severe cancellation problems, more accurate simulation results can be obtained with the same computational resources.)
[0142] (3. A semi-implicit scheme is used to advance the ion and electron weights, and the electromagnetic field equations are uniformly processed into a second-order time-accurate scheme. While maintaining high accuracy, this algorithm effectively relaxes the restrictions on the simulation time step, thereby significantly reducing the amount of computation and improving the feasibility and efficiency of long-term, high-precision simulations.)
[0143] On the other hand, the present invention provides a kinetic simulation system for fusion reactor plasma. Figure 8 This is a structural block diagram of a kinetic simulation system for fusion reactor plasma according to an embodiment of the present invention. (Reference) Figure 8 The system includes: a simulation initialization module, a particle display and propulsion module, a particle current density and pressure calculation module, a field equation solving module, a particle implicit propulsion module, and a data processing and output module.
[0144] The simulation initialization module is used to: perform simulation initialization on the simulation region and drive the ion weights and electron weights using a first-order implicit scheme.
[0145] The particle display propulsion module is used to: based on the simulation initialization and the grid nodes of the simulation region, adopt the particle display propulsion format to obtain the ion position and velocity, ion weight, electron magnetic moment and velocity, and electron weight, and determine the preset interval to which each ion or electron belongs based on the updated coordinates of each ion or electron.
[0146] The particle current density and pressure calculation module is used to: based on ion weight, ion position and velocity, according to the shape factor and preset interval, distribute the current at the ion position to the adjacent grid nodes to obtain the ion current density on the grid node; based on electron weight, electron magnetic moment and velocity, according to the shape factor and preset interval, distribute the pressure and current at the electron position to the adjacent grid nodes to obtain the electron pressure and electron parallel current density.
[0147] The field equation solving module is used to: use ion current density, electron pressure, and electron parallel current density to solve the implicit parallel Ampere's law and perpendicular Ohm's law using an iterative method to obtain the electric field strength; and update the magnetic field based on Faraday's law and the electric field strength to obtain the magnetic induction intensity.
[0148] The implicit particle propulsion module is used to: update the ion weights and electron weights by implicitly propulsing particles using updated electric field strength and magnetic induction intensity, and return to the explicit particle propulsion module for processing if the simulation time step is less than a first threshold.
[0149] In some embodiments, the system further includes: a data processing and output module, used to process and output simulation results of updated electric field intensity and magnetic induction intensity when the simulation time step is greater than or equal to a first threshold, wherein processing and outputting simulation results of updated electric field intensity and magnetic induction intensity includes: performing Fourier transform on updated electric field intensity and magnetic induction intensity to obtain the distribution spectrum of spatial wavenumber and frequency domain; writing updated electric field intensity and magnetic induction intensity into a specified output file for subsequent plotting of spatiotemporal evolution diagram; or, if linear eigenmode research is carried out, calculating the linear growth rate of instability based on particle energy flux, and simultaneously outputting ion current density, electron pressure, electron parallel current density and energy flux.
[0150] Other implementation details of this system are described in the section on kinetic simulation methods for fusion reactor plasmas above, and will not be repeated here.
[0151] In implementing the functions of the integrated modules described above in hardware, this embodiment of the invention provides a structure for a kinetic simulation device for fusion reactor plasma as described in the above embodiments. Figure 9 This is a schematic diagram of a kinetic simulation device for fusion reactor plasma according to an embodiment of the present invention. (Reference) Figure 9 The kinetic simulation device for fusion reactor plasma includes: at least one processor; and at least one memory. The at least one memory is coupled to the at least one processor and stores instructions for execution by the at least one processor, which, when executed by the at least one processor, implement the methods described above.
[0152] A processor can be a set of logic blocks, modules, and circuits that implement or execute the various exemplary logic blocks, modules, and circuits described in connection with embodiments of the present invention. The processor can be a central processing unit, a general-purpose processor, a digital signal processor, an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. It can implement or execute the various exemplary logic blocks, modules, and circuits described in connection with embodiments of the present invention. A processor can also be a combination that implements computational functions, such as a combination of one or more microprocessors, a combination of a digital signal processor and a microprocessor, etc.
[0153] The memory may be read-only memory (ROM) or other types of static storage devices that can store static information and instructions, random access memory (RAM) or other types of dynamic storage devices that can store information and instructions, electrically erasable programmable read-only memory (EEPROM), disk storage media or other magnetic storage devices, or any other medium that can be used to carry or store desired program code in the form of instructions or data structures and that can be accessed by a computer, but is not limited thereto.
[0154] In one implementation, the memory can exist independently of the processor. The memory can be connected to the processor via a bus and used to store instructions or program code. When the processor calls and executes the instructions or program code stored in the memory, it can implement the method provided in the embodiments of the present invention. In another implementation, the memory can also be integrated with the processor.
[0155] On the other hand, the present invention also provides a computer-readable storage medium (e.g., a non-transitory computer-readable storage medium) storing computer program instructions that, when executed on a computer, cause the computer to perform the method as described in any of the above embodiments.
[0156] Exemplary examples show that the aforementioned computer-readable storage media may include, but are not limited to: magnetic storage devices (e.g., hard disks, floppy disks, or magnetic tapes), optical discs (e.g., compact disks (CDs), digital versatile disks (DVDs), etc.), smart cards, and flash memory devices (e.g., erasable programmable read-only memory (EPROMs), cards, sticks, or key drives, etc.). The various computer-readable storage media described in this invention may represent one or more devices and / or other machine-readable storage media for storing information. The term "machine-readable storage medium" may include, but is not limited to, wireless channels and various other media capable of storing, containing, and / or carrying instructions and / or data.
[0157] This invention provides a computer program that, when run on a computer, causes the computer to perform the method of any of the above embodiments.
[0158] This invention provides a computer program product containing instructions that, when run on a computer, cause the computer to perform the method of any of the above embodiments.
[0159] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.
Claims
1. A kinetic simulation method for fusion reactor plasma, characterized in that, The method includes: The simulation region is initialized, and the ion weights and electron weights are driven by a first-order implicit scheme. Based on the simulation initialization and the grid nodes of the simulation region, the particle display and push format is used to obtain the ion position and velocity, ion weight, electron magnetic moment and velocity, and electron weight, and to determine the preset interval to which each ion or electron belongs based on the updated coordinates of each ion or electron. Based on the ion weight, the ion position and velocity, the current at the ion position is distributed to adjacent grid nodes according to the shape factor and the preset interval to obtain the ion current density on the grid node; based on the electron weight, the electron magnetic moment and velocity, the pressure and current at the electron position are distributed to adjacent grid nodes according to the shape factor and the preset interval to obtain the electron pressure and electron parallel current density. Using the ion current density, the electron pressure, and the electron parallel current density, the implicit parallel Ampere's law and the perpendicular Ohm's law are solved using an iterative method to obtain the electric field strength. The magnetic field is then updated based on Faraday's law and the electric field strength to obtain the magnetic induction intensity. Using the updated electric field strength and magnetic induction intensity, the ion weights and electron weights are updated via implicit particle propulsion. If the simulation time step is less than a first threshold, the following steps are performed: Based on the simulation initialization and the grid nodes of the simulation region, the particle explicit propulsion format is used to obtain the ion position and velocity, ion weight, electron magnetic moment and velocity, and electron weight. The preset interval to which each ion or electron belongs is determined according to the updated coordinates of each ion or electron.
2. The method according to claim 1, characterized in that, The simulation initialization includes setting: simulation region parameters, steady-state density gradients and temperature gradients for electrons and ions, wavenumber filtering conditions, electron and ion weight initialization, and initial perturbation electric and magnetic fields.
3. The method according to claim 1 or 2, characterized in that, Based on the simulation initialization and the grid nodes of the simulation region, a particle display-driven format is used to obtain ion positions and velocities, ion weights, electron magnetic moments and velocities, and electron weights. Then, based on the updated coordinates of each ion or electron, the preset interval to which it belongs is determined, including: Ion position and velocity are driven by the Boris scheme, and ion weights are obtained by driving half-steps using an explicit scheme. The electron magnetic moment and velocity are driven by the drift motion physical equation, and the electron weight is obtained by driving half-steps using an explicit scheme; Based on the updated z-axis coordinates of each ion or electron, determine its corresponding preset. The interval, wherein the simulation space is divided into multiple consecutive intervals along the z-axis. In a range, each of the multiple processes is assigned and responsible for tracking all particles located within one of the ranges.
4. The method according to claim 1 or 2, characterized in that, Using the ion current density, the electron pressure, and the electron parallel current density, the electric field strength is obtained by iteratively solving the implicit parallel Ampere's law and the perpendicular Ohm's law. The magnetic field strength is then updated based on Faraday's law and the electric field strength to obtain the magnetic induction intensity, including: The magnetic flux density is obtained by pushing Faraday's law half a step based on the explicit scheme. For ions, the ion current density is iteratively calculated at the grid nodes; For electrons, the parallel current density and electron pressure of the electrons are iteratively calculated at the grid nodes; If the iteration error is less than or equal to the second threshold, the ion current density, electron parallel current density and electron pressure obtained by the iteration calculation are substituted into the implicit scheme parallel Ampere's law and perpendicular Ohm's law respectively to obtain the updated electric field intensity. Based on the updated electric field strength, Faraday's law is applied using an implicit scheme to obtain the updated magnetic flux density.
5. The method according to claim 1 or 2, characterized in that, The method further includes: based on the ion weight, the ion position, and the velocity, and according to the shape factor and the preset interval, distributing the current at the ion position to adjacent grid nodes to obtain the ion current density at the grid nodes; based on the electron weight, the electron magnetic moment, and the velocity, and according to the shape factor and the preset interval, distributing the pressure and current at the electron position to adjacent grid nodes to obtain the electron pressure and electron parallel current density, and then... The contributions of adjacent processes on the overlapping grid are summed to ensure the continuity of the calculation of the electron parallel current density and electron pressure within the preset interval. Boundary condition constraints are applied to the boundary grid point data of the simulation region to meet the boundary requirements of the physical model.
6. The method according to claim 1 or 2, characterized in that, The method further includes: when the simulation time step is greater than or equal to the first threshold, processing the updated electric field strength and magnetic induction intensity and outputting the simulation results, wherein processing the updated electric field strength and magnetic induction intensity and outputting the simulation results includes: The updated electric field strength and magnetic induction intensity are subjected to Fourier transform to obtain the distribution spectrum in the spatial wavenumber and frequency domain. The updated electric field strength and magnetic induction intensity are written to a specified output file for subsequent plotting of the spatiotemporal evolution diagram; or If linear eigenmode research is conducted, the linear growth rate of instability can be calculated based on particle energy flux, and the ion current density, electron pressure, electron parallel current density, and energy flux can be output simultaneously.
7. A kinetic simulation system for fusion reactor plasma, characterized in that, The system includes: The simulation initialization module is used to: initialize the simulation region and drive the ion weights and electron weights using a first-order implicit scheme; The particle display and propulsion module is used to: based on the simulation initialization and the grid nodes of the simulation region, adopt the particle display and propulsion format to obtain ion positions and velocities, ion weights, electron magnetic moments and velocities, and electron weights, and determine the preset interval to which each ion or electron belongs based on the updated coordinates of each ion or electron; The particle current density and pressure calculation module is used to: based on the ion weight, the ion position and velocity, according to the shape factor and the preset interval, distribute the current at the ion position to adjacent grid nodes to obtain the ion current density on the grid node; based on the electron weight, the electron magnetic moment and velocity, according to the shape factor and the preset interval, distribute the pressure and current at the electron position to adjacent grid nodes to obtain the electron pressure and electron parallel current density. The field equation solving module is used to: use the ion current density, the electron pressure and the electron parallel current density to solve the implicit parallel Ampere's law and the perpendicular Ohm's law using an iterative method to obtain the electric field strength; and update the magnetic field based on Faraday's law and the electric field strength to obtain the magnetic induction intensity. The implicit particle propulsion module is used to: update the ion weights and electron weights by implicitly propulsing particles using updated electric field strength and magnetic induction intensity, and return to the explicit particle propulsion module for processing if the simulation time step is less than a first threshold.
8. A kinetic simulation device for fusion reactor plasma, characterized in that, The device includes: At least one processor; At least one memory coupled to the at least one processor and storing instructions for execution by the at least one processor, the instructions implementing the method of any one of claims 1 to 6 when executed by the at least one processor.
9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores instructions that, when executed by a computer, cause the computer to perform the method of any one of claims 1 to 6.
10. A computer program product, characterized in that, The computer program product includes instructions that, when executed by a computer, cause the computer to perform the method according to any one of claims 1 to 6.