Collision method for Penning ion source kinetic simulation

Abstract

The invention relates to the technical field of Penning ion sources, and discloses a Penning ion source kinetic simulation collision method, which comprises the following steps: carrying out integration on an ion distribution function of a current time step in a three-dimensional velocity space to obtain a first parameter, carrying out weighted integration on the first parameter and an axial velocity component in the ion distribution function to obtain a second parameter, according to the method, the overall numerical value collision term expression of the algorithm is simple, no differential operator or integral operator exists, only simple moment integral, exponential function operation and arithmetic operation are needed, but the transportation, perpendicular to the magnetic field direction, of ions brought by collision can be completely simulated, and compared with the prior art, the method has the advantages that the simulation efficiency is high; according to the method, a part of physical steps are omitted, a collision calculation algorithm is simplified, the calculation amount of simulation is reduced, the problems of numerical instability and divergence in simulation are solved, and the transportation process of ions in the direction perpendicular to the magnetic field can still be efficiently and accurately described.

Description

Technical Field

[0001] The present invention relates to the technical field of Penning ion sources, and particularly to a collision method for kinetic simulation of a Penning ion source. Background Art

[0002] At present, plasma simulation is mainly divided into two categories: fluid models and kinetic models. In kinetic models, there are also PIC (particle in cell) programs and continuum programs. In the field of kinetic simulation of Penning ion sources, for the description of complex ion-ion collision processes, Lorentz collision models and Fokker-Planck collision term models are mainly used. The collision terms of these models contain rich physical meanings, but their expressions are complex, involving a large number of differential operators and integral operators, and there is a coupling relationship between the two.

[0003] In the kinetic simulation scenario of a Penning ion source (cylindrical structure, magnetic field direction is axial, and the ion extraction port is located on the side of the anode cylinder), the defects of the prior art are particularly prominent. The differential operators and integral operators of the Lorentz and Fokker-Planck collision term models, due to the fact that the integral operator needs to traverse the full-space velocity distribution for the statistical calculation of the collision probability, and the differential operator is sensitive to the local gradient calculation, the combined effect of the two results in a huge amount of calculation, which easily leads to numerical instability, causing the distribution function to diverge numerically during iteration, making the calculation unable to proceed normally. This makes it difficult for the prior art to efficiently and accurately describe the ion transport process perpendicular to the magnetic field direction. Summary of the Invention

[0004] Aiming at the deficiencies of the prior art, the present invention provides a collision method for kinetic simulation of a Penning ion source to solve the above problems.

[0005] The above technical objectives of the present invention are achieved through the following technical solutions:

[0006] A collision method for kinetic simulation of a Penning ion source includes:

[0007] S1: Integrate the ion distribution function at the current time step in three-dimensional velocity space to obtain a first parameter;

[0008] S2: Perform a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain a second parameter;

[0009] S3: Integrate the energy term of the ion distribution function to obtain the total energy, and perform an energy conservation transformation on the total energy and the second parameter to obtain a third parameter;

[0010] S4: Construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter, and apply relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update;

[0011] S5: Analyze the six-dimensional distribution function update to obtain the five-dimensional distribution function update.

[0012] Furthermore, integrate the ion distribution function at the current time step in the three-dimensional velocity space to obtain the first parameter, including:

[0013] Integrate the ion distribution function at the current time step in the three-dimensional velocity space to obtain the initial integration result;

[0014] Perform a spatial position mapping on the initial integration result to obtain the first parameter.

[0015] Furthermore, perform a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain the second parameter, including:

[0016] Extract the axial velocity component from the first parameter;

[0017] Calculate the axial velocity component and the ion distribution function to obtain the axial momentum integral value;

[0018] Calculate the axial momentum integral value and the particle number density parameter to obtain the second parameter.

[0019] Furthermore, integrate the energy term of the ion distribution function to obtain the total energy, and perform an energy conservation transformation on the total energy and the second parameter to obtain the third parameter, including:

[0020] Integrate the kinetic energy term and the magnetic moment energy term of the ion distribution function to obtain the total energy;

[0021] Perform a kinetic energy analysis on the second parameter to obtain the directed motion energy component, and subtract the directed motion energy component from the total energy to obtain the disordered thermal motion energy;

[0022] Transform the disordered thermal motion energy to obtain the third parameter.

[0023] Furthermore, construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter, and apply relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update, including:

[0024] Construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter;

[0025] Calculate the deviation between the current ion distribution function and the local Maxwell distribution function to obtain the distribution deviation degree;

[0026] Convert the collision time to obtain the time step ratio factor;

[0027] Perform operations on the distribution deviation degree and the time step ratio factor to obtain the six-dimensional distribution function update amount.

[0028] Furthermore, analyze the six-dimensional distribution function update amount to obtain the five-dimensional distribution function update amount, including:

[0029] Map the six-dimensional distribution function update amount in the gyrocenter coordinate system to establish an angular mapping relationship;

[0030] Perform a central transformation on the angular mapping relationship to obtain the five-dimensional gyrocenter distribution function update amount.

[0031] Furthermore, convert the collision time to obtain the time step ratio factor, including:

[0032] Analyze the first parameter and the third parameter to obtain the collision time;

[0033] Calculate the collision time and the fixed simulation time step to obtain the time step ratio factor.

[0034] Furthermore, the construction of the local Maxwell distribution function needs to satisfy the following three constraint conditions:

[0035] The first constraint condition: the first parameter remains the same before and after constructing the distribution function;

[0036] The second constraint condition: the second parameter remains unchanged before and after constructing the distribution function;

[0037] The third constraint condition: the total energy remains constant before and after constructing the distribution function.

[0038] Furthermore, the central transformation includes:

[0039] Establish a mapping table according to the gyrocenter coordinates and the gyro radius vector;

[0040] Based on the mapping table, transform the ion distribution function into a gyroangle function;

[0041] Perform a circular integration of the gyroangle function from 0 to 2π to complete the dimensionality reduction.

[0042] Furthermore, the calculation steps of the collision method do not include the coupled operation of differential operators and integral operators, and only include three operation methods.

[0043] In summary, the present invention mainly has the following beneficial effects:

[0044] This method is based on the Krook collision term model to calculate the distribution function including collision effects. The overall numerical collision term expression is simple, without any differential or integral operators, only simple moment integrals, exponential function operations, and arithmetic operations. However, it can fully simulate the ion transport perpendicular to the magnetic field direction caused by collisions. Compared with existing technologies, by omitting some physical steps, simplifying the collision calculation algorithm, and reducing the computational amount of the simulation, it can avoid the problems of numerical instability and divergence in the simulation, but still can efficiently and accurately describe the ion transport process perpendicular to the magnetic field direction. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 is the flowchart of the collision method for the Penning ion source kinetic simulation of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0046] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0047] Refer to Figure 1 , a collision method for the Penning ion source kinetic simulation, including:

[0048] S1: Integrate the ion distribution function at the current time step in the three-dimensional velocity space to obtain the first parameter;

[0049] S2: Perform a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain the second parameter;

[0050] S3: Integrate the energy term of the ion distribution function to obtain the total energy, and perform an energy conservation transformation on the total energy and the second parameter to obtain the third parameter;

[0051] S4: Construct a local Maxwell distribution function according to the first parameter, the second parameter, and the third parameter, and apply relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update amount;

[0052] S5: Analyze the six-dimensional distribution function update amount to obtain the five-dimensional distribution function update amount.

[0053] By integrating the ion distribution function in the three-dimensional velocity space, the complex collision term is transformed into an energy conservation transformation based on integral parameters and the construction of a local Maxwell distribution, avoiding the direct handling of the coupled calculation of a large number of differential and integral operators, effectively reducing the computational complexity. At the same time, the distribution function update is obtained by applying collision time relaxation using the local Maxwell distribution function, alleviating the numerical instability problem caused by the full-space traversal of the integral operator and the locally sensitive calculation of the differential operator, suppressing the numerical divergence in the distribution function iteration, and ensuring the normal progress of the calculation. On this basis, by constructing the six-dimensional and five-dimensional distribution function updates, the ion distribution changes are analyzed more clearly, so as to efficiently and accurately describe the transport process of ions in the direction perpendicular to the magnetic field, providing a reliable calculation method for the kinetic simulation of Penning ion sources and improving the simulation accuracy and efficiency.

[0054] In one case of this embodiment, the ion distribution function at the current time step is integrated in the three-dimensional velocity space to obtain a first parameter, including:

[0055] Integrating the ion distribution function at the current time step in the three-dimensional velocity space to obtain an initial integration result, specifically including: at each grid point in the three-dimensional physical space, performing numerical integration of the ion distribution function in the three-dimensional velocity space, where the velocity space is divided into grid cells with a fixed spacing, traversing all velocity cells, multiplying the ion distribution function value at the center point of each cell by the width of the cell in the three velocity directions, and accumulating the product results of all cells to obtain the initial integration value, i.e., the particle number density, at this physical position point;

[0056] Performing a spatial position mapping on the initial integration result to obtain the first parameter, specifically including: for each physical grid point, calculating the corresponding cyclotron radius vector component according to the magnetic moment, magnetic field strength, and cyclotron phase angle of the ions located at this physical grid point, then subtracting the cyclotron radius vector component from the physical position to obtain the cyclotron center coordinates corresponding to this physical point, and at the same time associating the particle number density value calculated at this physical grid point with the cyclotron center coordinates. Finally, through the bilinear interpolation algorithm, the density value associated with the cyclotron center coordinates is matched to the standard cyclotron center grid points to obtain the first parameter, i.e., the particle number density in the cyclotron center coordinate system;

[0057] Among them, according to the magnetic moment, magnetic field strength, and cyclotron phase angle of the ions at this physical grid point, the corresponding cyclotron radius vector components are calculated as follows: For the ions at the physical space grid point, obtain the ion magnetic moment, local magnetic field strength, and cyclotron phase angle (the instantaneous angular position of the ion on the cyclotron circumference, with a value range of 0 to 360 degrees) at this point. Multiply the square root of the magnetic moment by a constant factor and divide by the square root of the magnetic field strength to obtain the cyclotron radius. The constant factor is determined by basic physical constants such as the ion mass. When the cyclotron motion of the ions occurs in a plane perpendicular to the magnetic field direction, with the magnetic field direction as the normal, select the axis direction in the device coordinate system in this plane as the reference. Multiply the cyclotron radius by the cyclotron phase angle to obtain the horizontal component in the device coordinate system. Multiply the cyclotron radius by the cyclotron phase angle to obtain the vertical component in the device coordinate system. Set the direction parallel to the magnetic field in the device coordinate system to zero, and then obtain the three-dimensional displacement vector components (horizontal component, vertical component, 0) of the ion generated at a specific phase.

[0058] In the Penning ion source kinetic simulation, by performing a three-dimensional velocity space integration on the ion distribution function at the current time step, in the way of dividing into fixed-spacing grid cells and traversing and accumulating, the complex full-space velocity distribution statistics are transformed into a structured numerical integration, reducing the traversal range of the integration operator. At the same time, when mapping the spatial position, through the cyclotron radius vector calculation and bilinear interpolation algorithm, the conversion between the physical position and the cyclotron center coordinates is accurately processed, avoiding the redundant calculation caused by the local gradient sensitivity of the differential operator, reducing the overall computational complexity, effectively avoiding the numerical instability problem caused by excessive computational amount, and ensuring the efficient and stable operation of the calculation process.

[0059] By calculating the cyclotron radius vector components based on the ion magnetic moment, magnetic field strength, and cyclotron phase angle, the cyclotron motion characteristics of the ions in the magnetic field can be accurately captured, the physical position can be accurately mapped to the cyclotron center coordinates, and the bilinear interpolation algorithm is used to match to the standard grid points, so that the particle number density in the cyclotron center coordinate system obtained by the calculation can truly reflect the distribution and transport law of the ions in the direction perpendicular to the magnetic field, overcoming the defects of the traditional model in the description of this direction and improving the accuracy of the Penning ion source kinetic simulation.

[0060] In a case of this embodiment, a weighted integration is performed on the first parameter and the axial velocity component in the ion distribution function to obtain a second parameter, including:

[0061] Extract the first parameter to obtain the axial velocity component, specifically including: In the three-dimensional velocity space, for each grid cell in the three-dimensional velocity space, directly extract the velocity component value parallel to the magnetic field direction from the three-dimensional velocity coordinates of the center point of this cell as the axial velocity component of this cell;

[0062] Calculate the axial velocity component and the ion distribution function to obtain the integrated value of the axial momentum, which specifically includes: traversing all velocity space grid cells, extracting the axial velocity component value in the direction parallel to the magnetic field at the center point of each cell, obtaining the ion distribution function value corresponding to the center point of the cell, and calculating the axial momentum contribution of the cell (axial velocity component value distribution function value velocity cell volume, where the velocity cell volume is the product of the grid spacings in the three velocity directions), summing up the axial momentum contributions of all velocity cells to obtain the integrated value of the axial momentum at this physical grid point;

[0063] Calculate the axial momentum integrated value and the particle number density parameter to obtain the second parameter, which specifically includes: dividing the axial momentum integrated value by the first parameter at the same position to obtain the second parameter, and the second parameter represents the axial velocity.

[0064] By directly extracting the axial velocity component from the center point of the three-dimensional velocity space grid cell and performing weighted integral calculation based on this, the huge computational amount of the integral operator traversing the full-space velocity distribution is avoided, the complex and sensitive calculation of the differential operator for the local gradient is reduced, while simplifying the calculation process, the computational efficiency is greatly improved, a large amount of data can be processed quickly, and the computational resources and time costs are saved.

[0065] Through the targeted calculation of the axial velocity component and the ion distribution function, the second parameter is obtained through the reasonable operation of the axial momentum integrated value and the particle number density parameter. This calculation method optimizes the calculation logic, reduces the accumulation of numerical errors caused by complex operations, effectively controls the stability of the distribution function during the iterative calculation process, avoids the problem of numerical divergence, and ensures the accuracy of the calculation results.

[0066] In one case of this embodiment, integrate the energy term of the ion distribution function to obtain the total energy, and perform an energy conservation transformation on the total energy and the second parameter, including:

[0067] Integrate the kinetic energy term and the magnetic moment energy term of the ion distribution function to obtain the total energy, which specifically includes: traversing all three-dimensional velocity space grid cells, for each velocity cell, obtaining the velocity component value (axial velocity) parallel to the magnetic field direction and the magnetic moment value perpendicular to the magnetic field direction at the center point of the cell, and calculating the energy contribution of the cell, specifically: multiplying half of the product of the ion mass and the square of the axial velocity, adding the product of the magnetic moment value and the local magnetic field intensity value, multiplying the sum of these two terms by the ion distribution function value of the cell, then multiplying by the volume of the velocity cell (i.e., the product of the grid spacings in the three velocity directions), and then accumulating and summing up the calculation results of all velocity cells at this physical grid point, and the summation result is the total energy of this physical grid point;

[0068] Perform kinetic energy analysis on the second parameter to obtain the directional motion energy component, and subtract the directional motion energy component from the total energy to obtain the disordered thermal motion energy, specifically including: multiplying half of the square of the second parameter (axial velocity) by the first parameter (particle number density) and then multiplying by the ion mass value to obtain the directional motion energy component of all ions at this point due to the overall motion along the magnetic field direction. Subtract the directional motion energy component from the total energy to obtain the disordered thermal motion energy, which represents the total energy of all ions' disordered thermal motion (random motion) at this physical grid point;

[0069] Convert the disordered thermal motion energy to obtain the third parameter, specifically including: multiplying the Boltzmann constant by the first parameter (particle number density) and then multiplying by 1.5 to obtain the temperature coefficient. Divide the disordered thermal motion energy by the temperature coefficient to obtain the third parameter at this physical grid point, and the third parameter is the temperature.

[0070] By precisely integrating the energy term of the ion distribution function, the calculation complexity is effectively reduced. During the process of obtaining the total energy, only the three-dimensional velocity space grid cells are traversed for calculation, avoiding the huge computational amount of traversing the full-space velocity distribution. Compared with the traditional model where the integral operator needs to traverse the full-space velocity distribution, the consumption of computing resources and time cost are significantly reduced. At the same time, the total energy and the second parameter are subjected to an energy conservation transformation to obtain the third parameter. Through step-by-step calculation and parameter conversion, the calculation process is made clearer and more orderly, reducing the risk of numerical instability caused by complex calculations and significantly improving the calculation efficiency.

[0071] By separating the disordered thermal motion energy from the total energy and further converting it to obtain the temperature parameter, the accuracy of the simulation results is effectively improved. By distinguishing the directional motion energy component and the disordered thermal motion energy, the motion state and energy distribution of ions can be described more precisely. Compared with the defect of the prior art that it is difficult to accurately describe the ion transport process perpendicular to the magnetic field direction, the accurate calculation and analysis of ion energy in the present invention enable a more realistic reflection of the physical behavior of ions in a complex magnetic field environment during the kinetic simulation of a Penning ion source. Especially when describing the transport characteristics of ions perpendicular to the magnetic field direction, the problem of numerical divergence can be effectively avoided.

[0072] In one case of this embodiment, construct a local Maxwell distribution function according to the first parameter, the second parameter, and the third parameter, and apply relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update amount, including:

[0073] Construct a local Maxwell distribution function according to the first parameter, the second parameter, and the third parameter, specifically including: multiplying the particle number density by the power of the ion mass to obtain the numerator term , multiply the azimuthal integration result , the Boltzmann constant and the temperature , and perform an -th power operation on the product result to obtain the denominator term . Divide the numerator term by the denominator term to obtain the prefactor term ; Square the difference between the velocity component and the axial drift velocity , and then multiply the result by the ion mass to obtain the weighted velocity term . Multiply the base unit value 2, the Boltzmann constant and the temperature to obtain the thermal velocity term . Divide the weighted velocity term by the thermal velocity term and then take the negative sign to obtain the velocity exponent term . Apply the natural exponential function to the velocity exponent term to obtain the first exponential term ; Multiply the magnetic moment by the axial magnetic field strength to obtain the magnetic energy term . Multiply the Boltzmann constant and the temperature to obtain the thermal energy term . Take the negative sign of the result of dividing the magnetic energy term by the thermal energy term, and then apply the natural exponential function to obtain the second exponential term . Continuously multiply the prefactor term, the first exponential term, and the second exponential term to obtain the local Maxwell distribution function . When specifically applied, it can be implemented through the following calculation formula, for example: ;

[0074] In the formula, represents the local Maxwell distribution function, represents the particle number density related to the position, represents the ion mass, represents the Boltzmann constant, represents the temperature, represents the natural exponential function, represents the velocity component parallel to the magnetic field, represents the axial overall velocity, represents the magnetic moment, represents the axial magnetic field strength, represents the azimuthal integration result in the velocity space perpendicular to the magnetic field direction;

[0075] Calculate the deviation between the current ion distribution function and the local Maxwell distribution function to obtain the distribution deviation degree, which specifically includes: at each grid point position in the three-dimensional physical space, for all discrete velocity cells in the three-dimensional velocity space corresponding to each grid point position, for each such phase space cell, directly subtract the value of the Maxwell distribution function at this position from the actually measured non-equilibrium distribution function value at present. The obtained difference is the distribution deviation amount on this phase space cell. After traversing all physical positions and velocity cells, a set of distribution function deviation data covering the entire six-dimensional space is obtained, and this set is the distribution deviation degree;

[0076] Convert the collision time to obtain the time step scale factor;

[0077] Perform operations on the distribution deviation degree and the time step scale factor to obtain the six-dimensional distribution function update amount, which specifically includes: divide the simulation time step by the collision time and perform a negative sign operation on the calculation result to obtain the negative time factor term , calculate the difference between the non-equilibrium ion distribution function and the local Maxwell distribution function to obtain the distribution function difference term , multiply the negative time factor term by the distribution function difference term to obtain the six-dimensional distribution function update amount . When specifically applied, it can be implemented through the following calculation formula. For example: ;

[0078] In the formula, represents the six-dimensional distribution function update amount, represents the simulation time step, represents the collision time, represents the non-equilibrium ion distribution function.

[0079] By directly calculating the deviation between the current distribution function and the local Maxwell distribution and constructing a time step scale factor in combination with the collision time, the evolution of the distribution function caused by collisions is transformed into a quantifiable algebraic operation, avoiding the strong coupling problem of differential and integral operators in traditional models. This processing method not only retains the key physical mechanisms of thermal motion (velocity exponential term) and magnetic confinement (magnetic energy exponential term) in the velocity space, but also simplifies the cross-region statistical correlation through local approximation. For the ion transport under axial magnetic field confinement in a cylindrical structure, especially the vertical magnetic field transport process near the side outlet, it has a more delicate description ability, achieving an improvement in numerical stability and an optimization of calculation efficiency, ensuring the feasibility of long-time simulation and more accurately capturing the relaxation dynamics from the non-equilibrium state to the equilibrium state.

[0080] By constructing a local Maxwell distribution that includes the exponential factor of the magnetic energy term, directly correlating physical quantities such as magnetic moment, axial magnetic field intensity, and temperature, accurately describing the cyclotron motion and energy exchange process of ions under magnetic field confinement, the calculation of the distribution deviation focuses on the deviation between the measured value of each phase space unit and the local equilibrium state, avoiding the information ambiguity brought by global statistics, enabling the clear analysis of the coupling effect of the velocity component perpendicular to the magnetic field direction, drift velocity, and magnetic moment. Through the calculation of the update amount of collision time relaxation, the complex collision process is simplified to an exponential relaxation of the non-equilibrium distribution to the local equilibrium state, which not only retains the physical essence of the collision term but also avoids the calculation ambiguity caused by the coupling of differential-integral operators. It is especially applicable to the strong gradient transport scenario at the side extraction port of the anode cylinder, improving the ion extraction efficiency and beam quality.

[0081] In one case of this embodiment, analyzing the update amount of the six-dimensional distribution function to obtain the update amount of the five-dimensional distribution function, including:

[0082] Mapping the update amount of the six-dimensional distribution function in the cyclotron center coordinate system to establish an angular mapping relationship, specifically including: for the ions at each grid point in the three-dimensional physical space, distributing the update amount of the six-dimensional distribution function of the physical grid point to the corresponding angular positions of the standard cyclotron center grid points through the bilinear interpolation algorithm according to the cyclotron center coordinates and phase angles (that is, the update amount data corresponding to different phase angles for each cyclotron center coordinate), establishing a three-dimensional angular mapping relationship of "cyclotron center coordinate phase angle distribution update amount";

[0083] Performing a central transformation on the angular mapping relationship to obtain the update amount of the five-dimensional cyclotron center distribution function, specifically including: taking the reciprocal of the period range to obtain the normalization coefficient , determining that the lower limit of the integral of the cyclotron phase angle is 0 and the upper limit is 2π to form a complete cyclotron period interval, inputting the magnetic moment and the current cyclotron phase angle into the cyclotron radius function to obtain the cyclotron radius vector , adding the cyclotron center coordinates vectorially with the cyclotron radius vector term to obtain the offset position , extracting the update amount of the six-dimensional distribution function from the offset position and the current cyclotron phase angle to obtain the position-related update amount term , setting the position-related update amount term as the integrand of the phase angle integral, and performing a definite integral calculation of the integrand with respect to the current cyclotron phase angle on the integral interval to obtain the integral term , multiply the integral term by a normalization coefficient to obtain the update amount of the five-dimensional gyrocenter distribution function , when specifically applied, it can be achieved through the following calculation formula, for example: ;

[0084] In the formula, represents the update amount of the five-dimensional gyrocenter distribution function, represents the gyrocenter coordinate, represents the gyroradius vector, represents the gyro-phase angle, represents the integral of the gyro-phase angle, represents the gyro-phase angle of the period range.

[0085] By establishing a three-dimensional angular mapping relationship in the gyrocenter coordinate system, the refined distribution of the six-dimensional distribution function update amount is realized. For the ions at the three-dimensional physical space grid points, the bilinear interpolation algorithm is used to accurately map the update amount to the corresponding angular positions of the standard gyrocenter grid points according to the gyrocenter coordinates and phase angles, and the corresponding relationship of "gyrocenter coordinate phase angle distribution update amount" is constructed. This process avoids the non-differential traversal of the entire space and only performs data distribution at specific angular positions, which not only retains the distribution characteristics under different phase angles but also reduces redundant calculations through grid standardization. Compared with the rough processing of the velocity distribution in the entire space in the traditional method, the spatial resolution of data distribution is significantly improved, providing clearly structured input data for the subsequent integral calculation based on the gyro-period, reducing the computational complexity caused by data coupling from the underlying architecture, and laying a stable data foundation for the distribution function iteration in complex magnetic field environments.

[0086] Through periodic normalization and definite integral processing of the gyrophase angle, the update quantity of the six-dimensional distribution function is effectively reduced to the update quantity of the five-dimensional gyrocenter distribution function. By using the design of the normalization coefficient and the integral over the complete gyroperiod interval, combined with the vector operation of the gyroradius vector and the offset position, the position-related update quantity is accurately extracted as the integrand. Through definite integral calculation, the statistical average of the gyrophase angle is achieved. This process cleverly decouples the strong coupling relationship between the differential operator and the integral operator in the Lorentz and Fokker-Planck models. It not only avoids the traversal calculation of the full-space velocity distribution through the periodic average of the integral operator but also optimizes the constraint conditions through the local gradient calculation of the differential operator, significantly reducing the computational amount. At the same time, the mathematical processing of the definite integral effectively suppresses the accumulation of numerical noise and avoids the problem of numerical divergence of the distribution function caused by operator coupling in traditional methods. Especially for the complex boundary conditions of the axial magnetic field and the side extraction port in the Penning ion source cylinder structure, it can more efficiently capture the transport characteristics of ions in the direction perpendicular to the magnetic field, providing a reliable technical solution for high-precision simulation of ion collision processes and transport behaviors.

[0087] In one case of this embodiment, the collision time is converted to obtain the time step scale factor, including:

[0088] Analyze the first parameter and the third parameter to obtain the collision time, specifically including: at each grid point position in the three-dimensional physical space, raise the third parameter (temperature) to the three-halves power, multiply by the square root of the ion mass, then multiply by the square of the vacuum permittivity and the constant coefficient, and divide this calculation result by the fourth power of the elementary charge, and then divide by the local first parameter (particle number density) to obtain the average time interval for effective collisions of ions at this physical position point, that is, the collision time;

[0089] Calculate the collision time and the fixed simulation time step to obtain the time step scale factor, specifically including: divide the preset fixed simulation time step value by the collision time to obtain a dimensionless ratio, which is the time step scale factor. The time step scale factor represents the proportion of the ion experiencing the collision relaxation process within a single simulation time step. Among them, the time step scale factor is calculated independently for each physical grid point, and thus a scale factor distribution in the three-dimensional space can be formed.

[0090] By analyzing the particle number density (the first parameter) and temperature (the third parameter), at each grid point position in the three-dimensional physical space, the temperature is calculated to obtain the time-step scale factor, which represents the proportion of the collision relaxation process experienced by ions within a single simulation time step, and each grid point is calculated independently to form a three-dimensional spatial distribution. This process not only avoids the coupled operations of complex differential and integral operators in traditional models, greatly reduces the statistical calculation amount of traversing the full-space velocity distribution, fundamentally reduces the algorithm complexity and improves the simulation efficiency, but also effectively alleviates the numerical divergence problem caused by the gradient sensitivity of the differential operator and the wide statistical range of the integral operator through dimensionless processing, making the distribution function iteration more stable and ensuring the normal progress of the calculation. At the same time, relying on the fine distribution of the scale factor in the three-dimensional space, the local collision characteristics differences of ions during the transport in the direction perpendicular to the magnetic field under the axial magnetic field in the cylindrical structure are accurately captured, especially suitable for simulating the transport process in complex boundary regions such as the side outlet of the anode cylinder, providing a more accurate numerical description of the plasma dynamics behavior in the Penning ion source.

[0091] In one case of this embodiment, the construction of the local Maxwell distribution function needs to satisfy the following three constraint conditions:

[0092] The first constraint condition: The first parameter remains consistent before and after constructing the distribution function;

[0093] The second constraint condition: The second parameter remains unchanged before and after constructing the distribution function;

[0094] The third constraint condition: The total energy remains constant before and after constructing the distribution function.

[0095] By constructing the local Maxwell distribution function that satisfies the three constraint conditions, the defects of the existing Penning ion source kinetic simulation technology are effectively overcome. First, the first parameter, the second parameter, and the total energy are strictly kept constant before and after construction, ensuring the conservation of physical quantities and the physical self-consistency of the simulation model, fundamentally guaranteeing the accuracy of the description of the ion transport process. Secondly, for the problems of large calculation amount and poor numerical stability in the Lorentz and Fokker-Planck collision term models, the local Maxwell distribution function avoids the traversal statistics of the integral operator for the full-space velocity distribution through localization processing, and at the same time reduces the complexity of the differential operator's sensitive calculation of the local gradient, greatly reducing the calculation amount and improving the simulation efficiency. Moreover, this scheme effectively suppresses the numerical divergence problem, enhances the stability of the iterative process, enables the distribution function to maintain reliable evolution during complex collision processes, especially in the description of the ion transport process perpendicular to the magnetic field direction, and can achieve an efficient description of the ion transport characteristics.

[0096] In one case of this embodiment, the central transformation includes:

[0097] Establish a mapping table based on the guiding center coordinates and the gyroradius, specifically including: for each physical grid point, associate the distribution function value of this physical point with the guiding center coordinates. At the same time, for each guiding center grid point, calculate the physical position (guiding center coordinates gyroradius) through reverse calculation, and query the data of adjacent physical grid points to establish a two-way lookup table of coordinates, obtaining a mapping table containing the corresponding relationship between physical positions and guiding center coordinates and interpolation weights;

[0098] Convert the ion distribution function into a gyroangle function based on the mapping table, specifically including: by querying the mapping table, reconstruct the distribution function into a form that explicitly depends on the gyro-phase angle. Specifically, for each guiding center grid point and phase angle, use the mapping relationship to calculate the corresponding physical position coordinates, read the original distribution function value at this position. It is necessary to traverse all guiding center coordinates and discrete phase angles (such as 0°, 90°, 180°, and 270°), locate the associated physical grid points through the mapping table, and directly extract the distribution function value at this phase angle to form a new function with the phase angle as the independent variable, obtaining a five-dimensional gyroangle function containing the guiding center position, velocity, magnetic moment, and phase angle;

[0099] Perform a circular integral of the gyroangle function from 0 to 2π to complete dimensionality reduction, specifically including: equally divide the phase angle range from 0 to 360 degrees into a predetermined number of parts (36 parts, with a step size of 10 degrees). For each guiding center point, traverse all discrete phase angle points, sequentially read the function values corresponding to two adjacent phase angles, add the two function values and multiply by half of the step size, accumulate the calculation results of all adjacent points, and divide the accumulated sum by 2π to obtain the average distribution function value of this guiding center point, thereby eliminating the phase angle dimension and reducing the original function to a five-dimensional distribution function that only depends on the guiding center position, axial velocity, and magnetic moment.

[0100] By establishing a mapping table and a two-way lookup table, the corresponding relationship between physical positions and guiding center coordinates and interpolation weights is constructed, laying an efficient data foundation for subsequent calculations. When dealing with the collision term, there is no need to traverse the full-space velocity distribution indiscriminately. Instead, with the help of the mapping table, the associated physical grid points are accurately located, and the distribution function value at the corresponding phase angle is directly extracted, significantly reducing the amount of calculation. Performing a circular integral on the gyroangle function to complete dimensionality reduction and eliminating the phase angle dimension simplifies the function form and reduces the coupling complexity of differential operators and integral operators. This not only alleviates the problem of explosion in the amount of calculation caused by their combined action but also significantly improves the numerical stability, effectively avoiding numerical divergence in the distribution function during iteration, enabling the kinetic simulation calculation to proceed stably and efficiently.

[0101] By transforming the ion distribution function into a five-dimensional gyration angle function that includes the phase angle, the influence of the phase angle in the ion gyration motion is fully considered, and the distribution characteristics of ions at different gyration phases can be carefully captured. By performing equidistant segmentation integration on the gyration angle function and weighted averaging the function values within the phase angle range, the interference of local phase fluctuations is eliminated, enabling the obtained average distribution function value to more accurately reflect the overall transport trend of ions in the direction perpendicular to the magnetic field. Compared with the description dilemma caused by complex calculations and numerical instability in the prior art, the present invention provides technical support for accurately depicting the transport process of ions in a Penning ion source through reasonable coordinate transformation and dimensionality reduction, and improves the practicality of the kinetic model in this scenario.

[0102] In one case of this embodiment, the calculation steps of the collision method do not include the coupled operation of differential operators and integral operators, and only include three operation methods:

[0103] The first operation method: single-moment integral operation;

[0104] The second operation method: exponential function operation;

[0105] The third operation method: four arithmetic operations.

[0106] In the kinetic simulation scenario of a Penning ion source, the calculation steps of the algorithm in this method do not include the coupled operation of differential operators and integral operators, and only use single-moment integral operation and two exponential function operations, fundamentally solving many problems caused by the coupling of the two operators in the prior art, avoiding the huge computational amount caused by the integral operator traversing the full-space velocity distribution to statistically calculate the collision probability and the differential operator sensitively calculating the local gradient, greatly improving the computational efficiency. At the same time, the hidden danger of numerical instability caused by their combined action is eliminated, effectively preventing the numerical divergence of the distribution function during the iteration process, and ensuring that the calculation can be carried out stably and normally. This enables the algorithm to more efficiently and accurately describe the transport process of ions in the direction perpendicular to the magnetic field.

[0107] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A collision method for kinetic simulation of a Penning ion source, characterized in that including: S1: Integrate the ion distribution function at the current time step in the three-dimensional velocity space to obtain the first parameter; S2: Perform a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain the second parameter; S3: Integrate the energy term of the ion distribution function to obtain the total energy, and perform an energy conservation transformation on the total energy and the second parameter to obtain the third parameter; S4: Construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter, and apply relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update; S5: Analyze the six-dimensional distribution function update to obtain the five-dimensional distribution function update.

2. The collision method for Penning ion source kinetic simulation according to claim 1, wherein Integrating the ion distribution function at the current time step in the three-dimensional velocity space to obtain the first parameter includes: Integrate the ion distribution function at the current time step in the three-dimensional velocity space to obtain the initial integration result; Perform a spatial position mapping on the initial integration result to obtain the first parameter.

3. A collision method for Penning ion source kinetic simulation according to claim 2, characterized in that, Performing a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain the second parameter includes: Extract the axial velocity component from the first parameter; Calculate the axial momentum integral value for the axial velocity component and the ion distribution function; Calculate the second parameter for the axial momentum integral value and the particle number density parameter.

4. A collision method for Penning ion source kinetic simulation according to claim 3, characterized in that, Integrating the energy term of the ion distribution function to obtain the total energy, and performing an energy conservation transformation on the total energy and the second parameter to obtain the third parameter includes: Integrate the kinetic energy term and the magnetic moment energy term of the ion distribution function to obtain the total energy; Perform a kinetic energy analysis on the second parameter to obtain the directional motion energy component, and subtract the directional motion energy component from the total energy to obtain the disordered thermal motion energy; Transform the disordered thermal motion energy to obtain the third parameter.

5. A collision method for Penning ion source kinetic simulation according to claim 4, characterized in that, Constructing a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter, and applying relaxation of the collision time to the local Maxwell distribution function to obtain the six-dimensional distribution function update includes: Construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter; Calculate the deviation between the current ion distribution function and the local Maxwell distribution function to obtain the distribution deviation; Convert the collision time to obtain the time step ratio factor; Perform an operation on the distribution deviation and the time step ratio factor to obtain the six-dimensional distribution function update.

6. The collision method for Penning ion source kinetic simulation according to claim 5, wherein Analyzing the six-dimensional distribution function update to obtain the five-dimensional distribution function update includes: Perform a mapping on the six-dimensional distribution function update in the gyrocenter coordinate system to establish an angular mapping relationship; Perform a central transformation on the angular mapping relationship to obtain the five-dimensional gyrocenter distribution function update.

7. A collision method for Penning ion source kinetic simulation according to claim 5, characterized in that Converting the collision time to obtain the time step ratio factor includes: Analyze the first parameter and the third parameter to obtain the collision time; Calculate the time step ratio factor for the collision time and the fixed simulation time step.

8. A collision method for Penning ion source kinetic simulation according to claim 5, characterized in that, The construction of the local Maxwell distribution function needs to satisfy the following three constraint conditions: The first constraint condition: The first parameter remains consistent before and after constructing the distribution function; The second constraint condition: The second parameter remains unchanged before and after constructing the distribution function; The third constraint condition: The total energy remains constant before and after constructing the distribution function.

9. A collision method for Penning ion source kinetic simulation according to claim 6, characterized in that, The said central transformation includes: Establishing a mapping table based on the guiding center coordinates and the gyroradius vector; Converting the ion distribution function into a gyrophase function based on the mapping table; Performing a circular integral of the gyrophase function from 0 to 2π to complete the dimensionality reduction.

10. A collision method for Penning ion source kinetic simulation according to claim 1, characterized in that The calculation steps of the collision method do not include the coupled operation of differential operators and integral operators, and only include three operation methods.

Images