A collision method for kinetic simulation of Penning ion sources

Abstract

The present invention relates to the technical field of Penning ion sources and discloses a collision method for kinetic simulation of a Penning ion source. The method comprises integrating an ion distribution function of a current time step in a three-dimensional velocity space to obtain a first parameter, and performing a weighted integration of the first parameter and an axial velocity component in the ion distribution function to obtain a second parameter. The method has a simple expression for the numerical collision term of the entire algorithm, without any differential or integral operators, but only simple moment integrals, exponential function operations, and four arithmetic operations. However, the method is fully capable of simulating the transport of ions perpendicular to the magnetic field caused by the collision. Compared with existing technologies, the method avoids numerical instability and divergence problems in the simulation by omitting some physical steps, simplifying the collision calculation algorithm, and reducing the computational complexity of the simulation. However, the method can still efficiently and accurately describe the transport process of ions perpendicular to the magnetic field.

Description

Technical Field

[0001] The invention relates to the technical field of Penning ion sources, and in particular 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. Kinetic models include PIC (particle in cell) programs and continuum programs. In the field of kinetic simulation of Penning ion sources, the Lorentz collision model and Fokker-Planck collision term model are mainly used to describe the complex collision process between ions. The collision terms of these models contain rich physical meanings, but the 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 the Penning ion source (cylindrical structure, axial magnetic field direction, ion introduction port located on the side of the anode tube), the defects of the existing technology are particularly prominent. Among them, the differential operator and integral operator of the Lorentz and Fokker-Planck collision term models require the integral operator to traverse the entire spatial velocity distribution to calculate the statistics of the collision probability, while the differential operator is sensitive to the local gradient. The combined effect of the two requires a huge amount of calculation and is prone to numerical instability, resulting in numerical divergence of the distribution function during iteration, making the calculation impossible. This makes it difficult for the existing technology to efficiently and accurately describe the ion transport process in the direction perpendicular to the magnetic field. Summary of the Invention

[0004] In view of the deficiencies in 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, comprising:

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

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

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

[0010] S4: constructing a local Maxwell distribution function according to the first parameter, the second parameter, and the third parameter, applying a relaxation of the collision time to the local Maxwell distribution function, and obtaining an updated amount of the six-dimensional distribution function;

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

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

[0013] Perform three-dimensional velocity space integration on the ion distribution function of the current time step to obtain the initial integration result;

[0014] The initial integration result is mapped to a spatial position to obtain the first parameter.

[0015] Furthermore, 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:

[0016] Extract the first parameter to obtain the axial velocity component;

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

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

[0019] Furthermore, the energy term of the ion distribution function is integrated to obtain the total energy, and the total energy and the second parameter are subjected to energy conservation transformation to obtain the third parameter, including:

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

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

[0022] The energy of disordered thermal motion is converted to obtain the third parameter.

[0023] Furthermore, a local Maxwell distribution function is constructed according to the first parameter, the second parameter, and the third parameter, and a relaxation of the collision time is applied to the local Maxwell distribution function to obtain an update of the six-dimensional distribution function, including:

[0024] constructing a local Maxwell distribution function according to 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;

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

[0027] The distribution deviation and the time step scale factor are calculated to obtain the updated amount of the six-dimensional distribution function.

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

[0029] The updated amount of the six-dimensional distribution function is mapped in the gyrocentric coordinate system to establish an angle mapping relationship;

[0030] The angle mapping relationship is transformed into a center to obtain the updated value of the five-dimensional convolution center distribution function.

[0031] Furthermore, the collision time is converted to obtain the time step scale factor, including:

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

[0033] The time-step scale factor is calculated by comparing the collision time and the fixed simulation time step.

[0034] Furthermore, the construction of the local Maxwell distribution function needs to meet the following three constraints:

[0035] The first constraint condition: the first parameter remains consistent 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: the total energy remains constant before and after the distribution function is constructed.

[0038] Furthermore, the central transformation includes:

[0039] Establish a mapping table according to the coordinates of the gyration center and the gyration radius vector;

[0040] Converting the ion distribution function into a cyclotron angle function based on a mapping table;

[0041] The cyclotron angle function is integrated from 0 to 2π to achieve dimensionality reduction.

[0042] Furthermore, the calculation steps of the collision method do not include the coupling operation of differential operators and integral operators, but 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 and calculates the distribution function including the collision effect. The overall numerical collision term expression is simple, without any differential operators or integral operators, only simple moment integrals, exponential function operations and four arithmetic operations. However, it is fully capable of simulating the transport of ions perpendicular to the magnetic field caused by the collision. Compared with existing technologies, it avoids the numerical instability and divergence problems in the simulation by omitting some physical steps, simplifying the collision calculation algorithm, and reducing the computational complexity of the simulation. However, it can still efficiently and accurately describe the transport process of ions perpendicular to the magnetic field. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 It is a flow chart of the steps of the collision method of the Penning ion source kinetic simulation of the present invention. DETAILED DESCRIPTION

[0046] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0047] refer to Figure 1 , a collision method for kinetic simulation of a Penning ion source, comprising:

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

[0049] S2: Perform 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 energy conservation transformation on the total energy and the second parameter to obtain the third parameter;

[0051] S4: constructing a local Maxwell distribution function according to the first parameter, the second parameter, and the third parameter, applying a relaxation of the collision time to the local Maxwell distribution function, and obtaining an updated amount of the six-dimensional distribution function;

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

[0053] By integrating the ion distribution function in the three-dimensional velocity space, the complex collision terms are converted into energy conservation transformations based on integral parameters and local Maxwell distribution construction, avoiding the direct processing of a large number of coupled calculations of differential and integral operators, effectively reducing the computational complexity. At the same time, the local Maxwell distribution function is used to impose collision time relaxation to obtain the distribution function update, which alleviates the numerical instability problems caused by the full-space traversal of the integral operator and the local sensitive calculation of the differential operator, suppresses the numerical divergence in the distribution function iteration, and ensures the normal progress of the calculation. On this basis, by constructing six-dimensional and five-dimensional distribution function updates, the changes in ion distribution can be analyzed more clearly, thereby efficiently and accurately describing the ion transport process in the direction perpendicular to the magnetic field, providing a reliable calculation method for the kinetic simulation of the Penning ion source, and improving the simulation accuracy and efficiency.

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

[0055] Performing a three-dimensional velocity space integration on the ion distribution function of the current time step to obtain an initial integration result, specifically comprising: performing a three-dimensional velocity space numerical integration on the ion distribution function at each grid point in the three-dimensional physical space, wherein 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 of the physical position point, i.e., the particle number density;

[0056] The initial integration result is mapped to a spatial position to obtain a first parameter, specifically including: for each physical grid point, according to the magnetic moment, magnetic field strength and cyclotron phase angle of the ion located at this physical grid point, the corresponding cyclotron radius vector component is calculated, and then the cyclotron radius vector component is subtracted from the physical position to obtain the cyclotron center coordinate corresponding to the physical point, and at the same time, the particle number density value calculated at this physical grid point is associated with the cyclotron center coordinate, and finally, through a bilinear interpolation algorithm, the density value associated with the cyclotron center coordinate is matched to the standard cyclotron center grid point 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 ion located at this physical grid point, the corresponding cyclotron radius vector component is calculated as follows: for the ion at the physical space grid point, the ion magnetic moment, local magnetic field strength and cyclotron phase angle (the instantaneous angular position of the ion on the cyclotron circle, ranging from 0 to 360 degrees) of the point are obtained, and the constant factor is multiplied by the square root of the magnetic moment and then divided 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 ion cyclotron motion occurs in a plane perpendicular to the direction of the magnetic field, the coordinate axis direction in the device coordinate system is selected in the plane with the magnetic field direction as the normal as the reference, and the cyclotron radius is multiplied by The rotation phase angle is used to obtain the horizontal component of the device coordinate system. Multiply the rotation radius by The cyclotron phase angle is rotated to obtain the vertical component of the device coordinate system, and the direction of the device coordinate system parallel to the magnetic field is set to zero, thereby obtaining the three-dimensional displacement vector component (horizontal component, vertical component, 0) of the ion generated at a specific phase.

[0058] In the kinetic simulation of the Penning ion source, the three-dimensional velocity space integration of the ion distribution function of the current time step is performed. The complex full-space velocity distribution statistics are converted into structured numerical integration by dividing the grid units with fixed spacing and traversing and accumulating them, which reduces the traversal range of the integral operator. At the same time, when mapping the spatial position, the conversion between the physical position and the gyrocenter coordinates is accurately processed through the gyroradius vector calculation and bilinear interpolation algorithm, avoiding the redundant calculation caused by the differential operator's sensitivity to local gradients, reducing the overall computational complexity, and effectively avoiding the numerical instability problem caused by excessive computational effort, ensuring the efficient and stable operation of the computational process.

[0059] By calculating the cyclotron radius vector component based on the ion magnetic moment, magnetic field strength and cyclotron phase angle, the cyclotron motion characteristics of 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 can be used to match it to the standard grid points. The calculated particle number density in the cyclotron center coordinate system can truly reflect the distribution and transport laws of ions in the direction perpendicular to the magnetic field, overcoming the defects of traditional models in describing this direction and improving the accuracy of Penning ion source kinetic simulation.

[0060] In one 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 the second parameter, including:

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

[0062] The axial velocity component and ion distribution function are calculated to obtain the axial momentum integral value, which specifically includes: traversing all velocity space grid cells, extracting the axial velocity component value of the center point parallel to the magnetic field direction 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 unit volume, where the velocity unit volume is the product of the grid spacing in the three velocity directions), sums the axial momentum contributions of all velocity units to obtain the axial momentum integral value of the physical grid point;

[0063] The axial momentum integral value and the particle number density parameter are calculated to obtain a second parameter, specifically comprising: dividing the axial momentum integral value by the first parameter at the same position to obtain the second parameter, the second parameter representing 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 it, the huge amount of calculation required for the integral operator to traverse the entire spatial velocity distribution is avoided, and the complex and sensitive calculation of the differential operator to the local gradient is reduced. While simplifying the calculation process, the calculation efficiency is greatly improved, and large amounts of data can be processed quickly, saving computing resources and time costs.

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

[0066] In one case of this embodiment, the energy term of the ion distribution function is integrated to obtain the total energy, and the total energy and the second parameter are subjected to energy conservation transformation to obtain the third parameter, including:

[0067] The kinetic energy term and the magnetic moment energy term of the ion distribution function are integrated to obtain the total energy, specifically including: traversing all three-dimensional velocity space grid cells, and 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 the ion mass by half the square of the axial velocity, adding the magnetic moment value multiplied by the local magnetic field strength value, multiplying the sum of these two terms by the ion distribution function value of the cell, and then multiplying it by the volume of the velocity cell (that is, the product of the grid spacing in the three velocity directions), and then summing the calculation results of all velocity cells on the physical grid point. The sum result is the total energy of the physical grid point;

[0068] Perform kinetic energy analysis on the second parameter to obtain the directional motion energy component, and deduct the directional motion energy component from the total energy to obtain the disordered thermal motion energy. Specifically, the following steps are performed: multiply the ion mass value by half the square of the second parameter (axial velocity), and then multiply it by the first parameter (particle number density) to obtain the directional motion energy component of all ions at that 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. The disordered thermal motion energy represents the sum of the energy of the disordered thermal motion (random motion) of all ions at that physical grid point.

[0069] The disordered thermal motion energy is converted to obtain the third parameter, specifically including: multiplying the Boltzmann constant by the first parameter (particle number density) and then multiplying it by 1.5 to obtain the temperature coefficient, and dividing the disordered thermal motion energy by the temperature coefficient to obtain the third parameter of the physical grid point, which is the temperature.

[0070] By accurately integrating the energy term of the ion distribution function, the computational complexity is effectively reduced. In the process of obtaining the total energy, only the three-dimensional velocity space grid cells are traversed and calculated, avoiding the huge amount of calculation for traversing the entire space velocity distribution. Compared with the traditional model in which the integral operator needs to traverse the entire space velocity distribution, the computing resource consumption and time cost are greatly reduced. At the same time, the total energy and the second parameter are transformed by energy conservation 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 computational efficiency.

[0071] By separating the disordered thermal motion energy from the total energy and further converting it into temperature parameters, 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 the ions can be described more precisely. Compared with the defect of the existing technology that it is difficult to accurately describe the transport process of ions in the direction perpendicular to the magnetic field, the precise calculation and analysis of ion energy in the present invention enables the physical behavior of ions in a complex magnetic field environment to be more realistically reflected in the Penning ion source kinetic simulation, especially when describing the transport characteristics of ions perpendicular to the magnetic field, which can effectively avoid the problem of numerical divergence.

[0072] In one case of this embodiment, a local Maxwell distribution function is constructed based on the first parameter, the second parameter, and the third parameter, and a collision time relaxation is applied to the local Maxwell distribution function to obtain a six-dimensional distribution function update, including:

[0073] The local Maxwell distribution function is constructed based on the first parameter, the second parameter and the third parameter, specifically including: the particle number density and ion mass of Multiply the power to get the numerator , the azimuth integration result , Boltzmann constant and temperature Multiply, and multiply the result Power operation to get the denominator , divide the numerator by the denominator to get the leading coefficient term ; The velocity component and axial drift velocity The difference is squared and then multiplied by the ion mass , and obtain the weighted acceleration term , the base unit value 2, the Boltzmann constant and temperature Multiplying them together gives the thermal velocity term , divide the weighted velocity term by the thermal velocity term and take the negative sign to get the velocity exponent term , apply the natural exponential function to the speed exponential term, and get the first exponential term ; Magnetic moment Multiply by the axial magnetic field strength , we get the magnetic energy term , the Boltzmann constant and temperature Multiply them to get the thermal energy term , minus the result of dividing the magnetic energy term by the thermal energy term, and then applying the natural exponential function, we get the second exponential term , continuously multiply the leading coefficient term, the first exponential term, and the second exponential term to obtain the local Maxwell distribution function , when applied specifically, it can be achieved through the following calculation formula, for example:

[0074] ;

[0075] Where, represents the local Maxwell distribution function, represents the position-dependent particle number density, represents the ion mass, represents the Boltzmann constant, Indicates temperature, represents the natural exponential function, represents the velocity component parallel to the magnetic field, represents the overall axial velocity, represents the magnetic moment, represents the axial magnetic field strength, It represents the azimuthal integral result of the vertical magnetic field direction in velocity space;

[0076] Calculate the deviation between the current ion distribution function and the local Maxwell distribution function to obtain the distribution deviation. Specifically, the following steps are performed: at each grid point position in the three-dimensional physical space, for all discrete velocity units in the three-dimensional velocity space corresponding to each grid point position, for each such phase space unit, directly subtract the currently actually measured non-equilibrium distribution function value from the Maxwell distribution function value at that position. The difference obtained is the distribution deviation at that phase space unit. After traversing all physical positions and velocity units, a distribution function deviation data set covering the entire six-dimensional space is obtained. This set is the distribution deviation.

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

[0078] The distribution deviation and time step scale factor are calculated to obtain the six-dimensional distribution function update, which specifically includes: Divide by the collision time , and minus the result to get the negative time factor term , for the non-equilibrium ion distribution function and the local Maxwell distribution function Perform difference calculation 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 , when applied specifically, it can be achieved through the following calculation formula, for example:

[0079] ;

[0080] Where, represents the update amount of the six-dimensional distribution function, represents the simulation time step, represents the collision time, represents the nonequilibrium ion distribution function.

[0081] By directly calculating the deviation between the current distribution function and the local Maxwell distribution, and constructing a time step scaling factor based on the collision time, the evolution of the distribution function caused by the collision is converted into a quantifiable algebraic operation, avoiding the problem of strong coupling between differential and integral operators in traditional models. This approach not only retains the key physical mechanisms of thermal motion (velocity exponential term) and magnetic confinement (magnetic energy exponential term) in velocity space, but also simplifies the cross-region statistical correlation through localized approximation. It has a more delicate description capability for ion transport under axial magnetic field constraints in cylindrical structures, especially the vertical magnetic field transport process near the side outlet, thereby achieving improved numerical stability and optimized computational efficiency, ensuring the feasibility of long-term simulation, and more accurately capturing the relaxation dynamics from non-equilibrium to equilibrium.

[0082] By constructing a local Maxwell distribution including an exponential factor of the magnetic energy term, directly linking physical quantities such as magnetic moment, axial magnetic field strength and temperature, the cyclotron motion and energy exchange process of ions under magnetic field constraints are accurately described. 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 caused by global statistics, so that the velocity component, drift velocity and magnetic moment coupling in the direction perpendicular to the magnetic field can be clearly analyzed. By calculating the updated amount of collision time relaxation, the complex collision process is simplified to the 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 suitable for the strong gradient transport scenario at the side outlet of the anode tube, improving the ion extraction efficiency and beam quality.

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

[0084] The six-dimensional distribution function update amount is mapped in the gyrocenter coordinate system to establish an angle mapping relationship, specifically including: for the ions at each grid point in the three-dimensional physical space, the six-dimensional distribution function update amount of the physical grid point is distributed to the corresponding angular position of the standard gyrocenter grid point according to the gyrocenter coordinates and phase angle through a bilinear interpolation algorithm (that is, each gyrocenter coordinate corresponds to the update amount data of different phase angles), and the "gyrocenter coordinate" is established. Phase Angle The three-dimensional angle mapping relationship of the "distribution update amount";

[0085] Perform center transformation on the angle mapping relationship to obtain the updated value of the five-dimensional gyrocentric center distribution function, specifically including: Take the reciprocal to get the normalized coefficient , determine the rotation phase angle The lower limit of the integral is 0 and the upper limit is 2π, which constitutes a complete cyclotron period interval. and the current rotation phase angle Enter the radius of gyration function , get the radius of gyration vector , the coordinates of the rotation center Add the vector to the gyration radius vector term to get the offset position , from the offset position and the current rotation phase angle Extract the six-dimensional distribution function update amount at the position to obtain the position-related update amount term , set the position-related update term to the integrand of the phase angle integral, and in the integral interval On the integrand, the current rotation phase angle Perform definite integral calculation to obtain the integral term , multiply the integral term by the normalization coefficient to obtain the updated value of the five-dimensional gyrocentric distribution function , when applied specifically, it can be achieved through the following calculation formula, for example:

[0086] ;

[0087] Where, represents the update amount of the five-dimensional gyrocentric distribution function, represents the coordinates of the center of rotation, represents the gyration radius vector, represents the rotation phase angle, represents the integral of the cyclotron phase angle, represents the rotation phase angle period range.

[0088] By establishing a three-dimensional angle mapping relationship in the cyclotron center coordinate system, the refined allocation of the six-dimensional distribution function update amount is achieved. 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 position of the standard cyclotron center grid point according to the cyclotron center coordinates and phase angles, thus constructing the "cyclotron center coordinate system". This process avoids indiscriminate traversal of the entire space and only distributes data at specific angle positions. It not only retains the distribution characteristics under different phase angles, but also reduces redundant calculations through grid standardization. Compared with the extensive processing of the velocity distribution in the entire space in traditional methods, it significantly improves the spatial resolution of data distribution, provides clearly structured input data for subsequent integral calculations based on the cyclotron period, reduces the computational complexity caused by data coupling from the underlying architecture, and lays a stable data foundation for the iteration of distribution functions in complex magnetic field environments.

[0089] Through period normalization and definite integral processing of the cyclotron phase angle, the six-dimensional distribution function update is effectively reduced to a five-dimensional cyclotron center distribution function update. By designing the normalization coefficient and the integral of the complete cyclotron period interval, combined with vector operations on the cyclotron radius vector and the offset position, the position-dependent update is precisely extracted as the integrand. The statistical averaging of the cyclotron phase angle is achieved through definite integral calculations. This process cleverly decouples the strong coupling between the differential and integral operators in the Lorentz and Fokker-Planck models. The periodic averaging of the integral operator avoids the ergodic calculation of the full-space velocity distribution, while optimizing the constraints of the local gradient calculation of the differential operator, significantly reducing the computational effort. Furthermore, the mathematical treatment of the definite integral effectively suppresses the accumulation of numerical noise and avoids the numerical divergence of the distribution function caused by operator coupling in traditional methods. In particular, for the complex boundary conditions between the axial magnetic field and the side exit in the cylindrical Penning ion source structure, the method can more efficiently capture the transport characteristics of ions perpendicular to the magnetic field, providing a reliable technical solution for high-precision simulation of ion collision processes and transport behavior.

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

[0091] The first and third parameters are analyzed to obtain the collision time, specifically comprising: at each grid point in the three-dimensional physical space, the third parameter (temperature) is raised to the cube of two, multiplied by the square root of the ion mass, multiplied by the square of the vacuum dielectric constant and the constant coefficient, and the result of this calculation is divided by the fourth power of the elementary charge, and then divided by the local first parameter (particle number density), to obtain the average time interval between effective collisions of ions at that physical location, i.e., the collision time;

[0092] The collision time and the fixed simulation time step are calculated to obtain the time step scale factor, specifically including: dividing 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 proportional degree of the collision relaxation process experienced by the ions within a single simulation time step. The time step scale factor is calculated independently for each physical grid point to form a scale factor distribution in three-dimensional space.

[0093] By analyzing the particle number density (the first parameter) and temperature (the third parameter), the temperature is calculated at each grid point in the three-dimensional physical space to obtain the time step scale factor. This factor represents the proportion of ions undergoing the collision relaxation process within a single simulation time step and is calculated independently at each grid point to form a three-dimensional spatial distribution. This process not only avoids the complex coupled operations of differential and integral operators in traditional models and significantly reduces the statistical calculation amount of traversing the velocity distribution in the entire space, fundamentally reducing algorithm complexity and improving simulation efficiency, but also effectively alleviates the numerical divergence problems 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 operation of the calculation. At the same time, the three-dimensional spatial scale factor distribution can accurately capture the local collision characteristics of ions transported in the perpendicular magnetic field direction under the axial magnetic field in the cylindrical structure. It is particularly suitable for simulating transport processes in complex boundary regions such as the side outlet of the anode cylinder, providing a more accurate numerical description of the plasma dynamic behavior in the Penning ion source.

[0094] In one case of this embodiment, the construction of the local Maxwell distribution function needs to meet the following three constraints:

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

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

[0097] The third constraint: the total energy remains constant before and after the distribution function is constructed.

[0098] By constructing a localized Maxwell distribution function that satisfies three major constraints, the shortcomings of existing Penning ion source kinetic simulation techniques are effectively overcome. First, the first and second parameters, as well as the total energy, are strictly maintained constant before and after construction, ensuring the conservation of physical quantities and the physical consistency of the simulation model, fundamentally guaranteeing the accuracy of the description of the ion transport process. Second, to address the computational complexity and poor numerical stability of the Lorentz and Fokker-Planck collision term models, the localized Maxwell distribution function avoids the traversal statistics of the integral operator over the entire spatial velocity distribution through localization. This also reduces the complexity of the differential operator's sensitivity to local gradients, significantly reducing the computational effort and improving simulation efficiency. Furthermore, this scheme effectively suppresses numerical divergence and enhances the stability of the iterative process, enabling the distribution function to evolve reliably during complex collisions. This allows for an efficient description of ion transport characteristics, particularly in the description of ion transport processes perpendicular to the magnetic field.

[0099] In one aspect of this embodiment, the central transformation includes:

[0100] The mapping table is established based on the gyration center coordinates and gyration radius, specifically including: for each physical grid point, the distribution function value of the physical point is associated with the gyration center coordinates, and for each gyration center grid point, the physical position (gyration center coordinates) is calculated by reverse calculation. Radius of gyration), and query the data of the adjacent physical grid points, establish a coordinate two-way lookup table, and obtain a mapping table containing the correspondence between the physical position and the coordinates of the gyration center and the interpolation weight;

[0101] The ion distribution function is converted into a cyclotron angle function based on a mapping table, specifically including: reconstructing the distribution function into a form that explicitly depends on the cyclotron phase angle by querying the mapping table, specifically: for each cyclotron center grid point and phase angle, the corresponding physical position coordinates are calculated using the mapping relationship, and the original distribution function value at the position is read. It is necessary to traverse all cyclotron 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 under the phase angle to form a new function with the phase angle as the independent variable, thereby obtaining a five-dimensional cyclotron angle function including the cyclotron center position, velocity, magnetic moment and phase angle;

[0102] The cyclotron angle function is integrated circularly from 0 to 2π to complete dimensionality reduction, specifically including: dividing the phase angle range from 0 to 360 degrees into predetermined parts (36 minutes, and the step size is 10 degrees), traversing all discrete phase angle points for each cyclotron center point, reading the function values ​​corresponding to the two adjacent phase angles in turn, adding the two function values ​​and multiplying them by half of the step size, accumulating the calculation results of all adjacent points, and dividing the accumulated total by 2π to obtain the average distribution function value of the cyclotron center point, thereby eliminating the phase angle dimension and reducing the original function to a five-dimensional distribution function that only depends on the cyclotron center position, axial velocity and magnetic moment.

[0103] By establishing a mapping table and a bidirectional lookup table, the correspondence between the physical position and the coordinates of the cyclotron center and the interpolation weight are constructed, laying an efficient data foundation for subsequent calculations. When processing collision terms, there is no need to indiscriminately traverse the velocity distribution in the entire space. Instead, the mapping table is used to accurately locate the associated physical grid points, and the distribution function value under the corresponding phase angle is directly extracted, which greatly reduces the amount of calculation. The cyclotron angle function is circularly integrated to complete the dimensionality reduction. After eliminating the phase angle dimension, the function form is simplified and the coupling complexity of the differential operator and the integral operator is reduced. This not only alleviates the problem of computational explosion caused by the joint action of the two, but also significantly improves numerical stability, effectively avoiding the numerical divergence of the distribution function during iteration, so that dynamic simulation calculations can be carried out stably and efficiently.

[0104] By converting the ion distribution function into a five-dimensional cyclotron angle function including the phase angle, the influence of the phase angle in the ion cyclotron motion is fully taken into account, and the distribution characteristics of ions under different cyclotron phases can be captured in detail. By integrating the cyclotron angle function in equidistant portions, the function values ​​within the phase angle range are weighted averaged, eliminating the interference of local phase fluctuations, so that the obtained average distribution function value can more accurately reflect the overall transport trend of ions in the direction perpendicular to the magnetic field. Compared with the description difficulties caused by complex calculations and numerical instability in the existing technology, the present invention provides technical support for accurately characterizing the transport process of ions in the Penning ion source through reasonable coordinate transformation and dimensionality reduction, thereby improving the practicality of the kinetic model in this scenario.

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

[0106] The first operation mode: single moment integral operation;

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

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

[0109] In the kinetic simulation scenario of the Penning ion source, the calculation steps of the algorithm in this method do not include the coupling operations of differential operators and integral operators. Only a single moment integral operation and two exponential function operations are used. This fundamentally solves many problems caused by the coupling of the two operators in the existing technology, avoids the huge amount of calculation caused by the integral operator traversing the entire space velocity distribution to calculate the statistical collision probability and the differential operator sensitively calculating the local gradient, greatly improves the calculation efficiency, and at the same time eliminates the hidden dangers of numerical instability caused by the joint action of the two, 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.

[0110] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the 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: include: S1: Integrate the ion distribution function of the current time step in the three-dimensional velocity space to obtain the first parameter, which is the particle number density in the cyclotron center coordinate system; S2: Perform a weighted integration on the first parameter and the axial velocity component in the ion distribution function to obtain a second parameter, which is the axial velocity; S3: Integrate the energy term of the ion distribution function to obtain the total energy, perform energy conservation transformation on the total energy and the second parameter to obtain the third parameter, which is temperature; S4: Construct a local Maxwell distribution function based on the first parameter, the second parameter, and the third parameter, apply collision time relaxation to the local Maxwell distribution function, and obtain a six-dimensional distribution function update, including: The local Maxwell distribution function is constructed according to the first parameter, the second parameter and the third parameter, specifically including: ； Where, represents the local Maxwell distribution function, represents the particle number density in the cyclotron center coordinate system, represents the ion mass, represents the Boltzmann constant, Indicates temperature, represents the natural exponential function, represents the velocity component parallel to the magnetic field, represents the axial speed, represents the magnetic moment, represents the axial magnetic field strength, It represents the azimuthal integral result of the direction perpendicular to the magnetic field in velocity space; 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 scale factor; Calculate the distribution deviation and the time step scale factor to obtain the updated amount of the six-dimensional distribution function; S5: Analyze the updated amount of the six-dimensional distribution function to obtain the updated amount of the five-dimensional distribution function, including: The updated amount of the six-dimensional distribution function is mapped in the gyrocentric coordinate system to establish an angle mapping relationship; Perform center transformation on the angle mapping relationship to obtain the updated value of the five-dimensional convolution center distribution function; The central transformation includes: Establish a mapping table according to the coordinates of the gyration center and the gyration radius vector; Converting the ion distribution function into a cyclotron angle function based on a mapping table; The cyclotron angle function is integrated from 0 to 2π to achieve dimensionality reduction.

2. The collision method for kinetic simulation of a Penning ion source according to claim 1, characterized in that: Integrate the ion distribution function of the current time step in the three-dimensional velocity space to obtain the first parameter, including: Perform three-dimensional velocity space integration on the ion distribution function of the current time step to obtain the initial integration result; The initial integration result is mapped to a spatial position to obtain the first parameter.

3. The collision method for kinetic simulation of a Penning ion source 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 a second parameter, including: Extract the first parameter to obtain the axial velocity component; The axial velocity component and ion distribution function are calculated to obtain the axial momentum integral value; The axial momentum integral value and the particle number density parameter are calculated to obtain the second parameter.

4. The collision method for kinetic simulation of a Penning ion source according to claim 3, characterized in that: Integrate the energy term of the ion distribution function to obtain the total energy, and perform energy conservation transformation on the total energy and the second parameter to obtain the third parameter, including: Integrate the kinetic energy term and magnetic moment energy term of the ion distribution function to obtain the total energy; Perform kinetic energy analysis on the second parameter to obtain the directional motion energy component, and deduct the directional motion energy component from the total energy to obtain the disordered thermal motion energy; The energy of disordered thermal motion is converted to obtain the third parameter.

5. The collision method for kinetic simulation of a Penning ion source according to claim 1, characterized in that: The collision time is converted to obtain the time step scale factor, including: Analyze the first parameter and the third parameter to obtain the collision time; The time-step scale factor is calculated by comparing the collision time and the fixed simulation time step.

6. The collision method for kinetic simulation of a Penning ion source according to claim 1, characterized in that: The construction of the local Maxwell distribution function needs to meet the following three constraints: 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: the total energy remains constant before and after the distribution function is constructed.

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

Images