Improved solution method for simulation of penning ion source cyclodynamics

Abstract

The application provides a simulation method for solving the improved Paning ion source cyclotron dynamics, and relates to the technical field of data processing.The method comprises the following steps: based on the motion equation, the original phase space coordinate system is transformed into a new phase space coordinate system through Lie transformation to obtain a generated vector field; based on the generated vector field, the transformation relationship between the original phase space coordinate system and the new phase space coordinate system is determined to obtain a simplified motion equation in the new phase space coordinate system; based on the simplified motion equation, the characteristic line method is adopted, and interpolation operation is performed on the axial position and parallel velocity dimensions to obtain a distribution function in the new phase space coordinate system; according to the distribution function in the new phase space coordinate system, the original phase space coordinate system is mapped back through the pullback transformation of Lie transformation, and finally, the original cyclotron center distribution function is obtained. The application simplifies the solving process of the charged particle cyclotron dynamics in the Paning ion source, and improves the simulation calculation efficiency.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to an improved solution method for cyclotron kinetics simulation of the Penning ion source. Background Technology

[0002] In the field of cyclotron kinetic simulation of Penning ion sources, the traditional approach to solve the Vlasov equation is to use the method of characteristics combined with high-dimensional interpolation techniques. Specifically, this method typically involves B-spline interpolation of the distribution function of the cyclotron center within a four-dimensional phase space (covering three spatial dimensions and one velocity dimension). However, since the computational complexity of interpolation tends to increase rapidly with increasing dimensionality, under some conventional grid settings, the number of addition and multiplication operations involved in each interpolation step may reach millions to tens of millions, and the memory usage may also be several GB to more than ten GB, which places high demands on computing equipment to a certain extent.

[0003] For example, in some common simulation scenarios, if the number of radial (x and y directions) grid points is set to 40 to 60, the number of axial (z direction) grid points to 150 to 250, and the number of parallel velocity dimension grid points to 50 to 80, in order to balance coverage of key areas of the discharge chamber and basic resolution requirements, the number of addition and multiplication operations involved in each step of four-dimensional B-spline interpolation may reach millions to tens of millions. At the same time, since the basis function data and intermediate calculation results required for four-dimensional interpolation need to be stored, the memory usage often reaches several GB to more than ten GB. Under such conditions, the simulation process of traditional methods usually needs to rely on servers with certain computing power to ensure running efficiency. If such simulations are carried out on ordinary office desktop computers, the calculation cycle may be long (e.g., the single simulation iteration takes a long time) or the running stability may be affected due to insufficient memory resources. This, to some extent, limits the flexible application of this method on more ordinary computing devices. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide an improved solution method for cyclotron kinetics simulation of Penning ion sources, which simplifies the solution process of cyclotron dynamics of charged particles in Penning ion sources and improves simulation calculation efficiency.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0006] Firstly, an improved solution method for cyclotron kinetics simulation of the Panning ion source, the method comprising:

[0007] Step 1: Obtain the electromagnetic field configuration of the Penning ion source, which includes a uniform magnetic field pointing along the axis and a three-dimensional spatially distributed potential.

[0008] Step 2: Based on the electromagnetic field configuration, establish the equation of motion of the gyroscope center in the original phase space coordinate system;

[0009] Step 3: Based on the equation of motion, the original phase space coordinate system is transformed into a new phase space coordinate system through Lie transformation to obtain the generated vector field;

[0010] Step 4: Based on the generated vector field, determine the transformation relationship between the original phase space coordinate system and the new phase space coordinate system, and obtain the simplified motion equations in the new phase space coordinate system;

[0011] Step 5: Based on the simplified equation of motion, the method of characteristics is used and interpolation is performed on the two dimensions of axial position and parallel velocity to obtain the distribution function in the new phase space coordinate system.

[0012] Step 6: Based on the distribution function in the new phase space coordinate system, the original phase space coordinate system is mapped back through the pullback transformation of the Lie transform, and finally the original cyclotron center distribution function is obtained.

[0013] In a second aspect, a computing device includes:

[0014] One or more processors;

[0015] A storage device for storing one or more programs that, when executed by one or more processors, cause the one or more processors to implement the method.

[0016] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0017] The above-described solution of the present invention has at least the following beneficial effects:

[0018] By using the Lie transformation to convert the original four-dimensional phase space coordinate system into a new phase space coordinate system, the motion equation of the cyclotron center is simplified to only needing to track the evolution of the axial position and parallel velocity. This simplifies the four-dimensional B-spline interpolation required to two-dimensional interpolation operations. Compared to four-dimensional interpolation, two-dimensional interpolation reduces the computational load and data storage requirements, effectively avoiding the problems of tens of millions of addition and multiplication operations and tens of gigabytes of memory usage caused by high-dimensional interpolation. This allows the simulation of the cyclotron kinetics of the Penning ion source to run stably without relying on high-performance computing servers or supercomputing centers, and ordinary office desktop computers and other conventional computing equipment, thus lowering the hardware threshold for simulation technology. Attached Figure Description

[0019] Figure 1 This is a flowchart illustrating the improved solution method for cyclotron kinetics simulation of the Penning ion source provided in an embodiment of the present invention.

[0020] Figure 2 This is a flowchart illustrating the process of establishing the equation of motion of the gyroscope center in the original phase space coordinate system based on the electromagnetic field configuration, according to an embodiment of the present invention. Detailed Implementation

[0021] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0022] like Figure 1 As shown, embodiments of the present invention propose an improved solution method for cyclotron kinetics simulation of the Penning ion source, the method comprising the following steps:

[0023] Step 1: Obtain the electromagnetic field configuration of the Penning ion source, which includes a uniform magnetic field pointing along the axis and a three-dimensional spatially distributed potential.

[0024] Step 2: Based on the electromagnetic field configuration, establish the equation of motion of the gyroscope center in the original phase space coordinate system;

[0025] Step 3: Based on the equation of motion, the original phase space coordinate system is transformed into a new phase space coordinate system through Lie transformation to obtain the generated vector field;

[0026] Step 4: Based on the generated vector field, determine the transformation relationship between the original phase space coordinate system and the new phase space coordinate system, and obtain the simplified motion equations in the new phase space coordinate system;

[0027] Step 5: Based on the simplified equation of motion, the method of characteristics is used and interpolation is performed on the two dimensions of axial position and parallel velocity to obtain the distribution function in the new phase space coordinate system.

[0028] Step 6: Based on the distribution function in the new phase space coordinate system, the original phase space coordinate system is mapped back through the pullback transformation of the Lie transform, and finally the original cyclotron center distribution function is obtained.

[0029] In this embodiment of the invention, the original four-dimensional phase space coordinate system is transformed into a new phase space coordinate system through Lie transformation, which simplifies the motion equation of the cyclotron center to only needing to track the evolution of the axial position and parallel velocity. This simplifies the four-dimensional B-spline interpolation required to two-dimensional interpolation operations. Compared with four-dimensional interpolation, the computational load and data storage requirements of two-dimensional interpolation are reduced, which can effectively avoid the problem of tens of millions of addition and multiplication operations and tens of GB of memory occupation caused by high-dimensional interpolation. This allows the simulation of the cyclotron kinetics of the Penning ion source to run stably without relying on high-performance computing servers or supercomputing centers, and ordinary office desktop computers and other conventional computing equipment can be used, thus lowering the hardware threshold of the simulation technology.

[0030] In a preferred embodiment of the present invention, step 1, obtaining the electromagnetic field configuration of the Penning ion source, the electromagnetic field configuration including a uniform magnetic field pointing along the axis and a three-dimensional spatially distributed potential, may include:

[0031] Step 101: Receive the discharge chamber geometric parameters of the Penning ion source. These parameters include the axial length, radial radius, and geometric positions of the cathode and anode. Based on these parameters, set a uniform magnetic field along the axial direction. Specifically, this includes receiving the discharge chamber geometric parameters via a data input interface. The axial length L refers to the vertical distance between the two end planes (planes perpendicular to the axial direction) of the discharge chamber, measured in millimeters, typically ranging from 10mm to 20mm. The radial radius r refers to the radius of the circular cross-section of the discharge chamber, i.e., the distance from the inner wall to the central axis, measured in millimeters, typically ranging from 4mm to 10mm. The cathode geometric position includes the axial coordinates. (The axial distance between the plane where the cathode is located and the plane at the left end of the discharge chamber, in millimeters, is usually...) =0, i.e., coinciding with the left-end plane) and radial range (The radial radius of the cathode, i.e., the distance from the central axis to the edge of the cathode, in millimeters); the anode geometry includes axial coordinates. (The axial distance between the plane where the anode is located and the plane at the left end of the discharge chamber, in millimeters, is usually...) =L, i.e., coinciding with the right-hand plane) and radial range (The radial distribution radius of the anode, i.e., the distance from the central axis to the inner wall of the anode, is usually expressed in meters.) =r, i.e., coinciding with the inner wall of the discharge chamber); after obtaining the above parameters, a uniform magnetic field along the axial direction (z direction) is set according to the radial radius r and axial length L. The direction of the magnetic field is consistent with the z-axis, and the magnetic field strength is... The value is determined by combining r and L (usually 0.1 to 1 Tesla).

[0032] Step 102: Input the discharge chamber geometric parameters and uniform magnetic field into the electric field solver, and calculate the three-dimensional spatial distribution of electric potential based on the electric field solver; combine the uniform magnetic field with the calculated three-dimensional spatial distribution of electric potential to generate the electromagnetic field configuration, specifically including: combining the above parameters (L, r, ... , , , ) and magnetic field strength Input electric field solver, the solver is based on the fundamental equations of electrostatics ,in, For the Laplace operator, combined with the boundary conditions (cathode potential) 1 = -V0, where V0 is 50~100 volts; anode potential Solving for 2 = 0 (i.e., grounded); when solving, consider the cylindrical boundary (r = Place = 2) Effect on the electric field: The three-dimensional space (x∈[-r, r], y∈[-r, r], z∈[0, L]) is discretized into a grid with mesh sizes Δx=Δy=Δz=0.02mm to 0.2mm. The electric potential at each grid point is obtained through iterative calculation. (x, y, z), i.e., the three-dimensional electric potential distribution; then, the integrated magnetic field distribution ( Along the z-direction, and (x, y, z) = ) and potential distribution (x, y, z), forming a spatial array containing all points ( , The electromagnetic field configuration of ).

[0033] In this embodiment, the magnetic field is set by geometric parameters to ensure that the magnetic field is precisely matched with the physical dimensions of the discharge chamber, providing reliable magnetic field parameters for subsequent force calculations.

[0034] like Figure 2 As shown, in another preferred embodiment of the present invention, step 2, establishing the equation of motion of the gyroscope center in the original phase space coordinate system according to the electromagnetic field configuration, may include:

[0035] Step 201: Based on the uniform magnetic field and the three-dimensional spatially distributed electric potential, calculate the electromagnetic force acting on the gyroscope center, specifically including:

[0036] First, calculate the electric field intensity based on the electric potential distributed in three-dimensional space, where the electric field component in the x-direction is... represents the electric field intensity at a point in space along the x-axis, with the SI unit being volts per meter (V / m). Its value is equal to a negative electric potential. The partial derivative of the electric potential at a point in three-dimensional space with respect to x (the x-coordinate of that point, in meters) is... ; electric field component in the y direction represents the electric field strength at a point in space along the y-axis, with the SI unit being volts per meter (V / m). Its value is equal to a negative electric potential. The partial derivative with respect to y (the y-coordinate of the point, in meters), i.e. z-direction electric field component The electric field intensity at a point in space along the z-axis is expressed in volts per meter (V / m), and its value is equal to the negative electric potential. The partial derivative with respect to z (the z-coordinate of the point, in meters), i.e. .

[0037] Step 202: Based on the electromagnetic force, establish the equation of motion of the cyclotron center along the x-direction in the original phase space coordinate system according to the electric field component in the y-direction and the uniform magnetic field; based on the electromagnetic force, establish the equation of motion of the cyclotron center along the y-direction in the original phase space coordinate system according to the electric field component in the x-direction and the uniform magnetic field; based on the electromagnetic force, express the velocity of the cyclotron center along the z-direction in the original phase space coordinate system as a parallel velocity, and establish the equation of motion of the cyclotron center along the parallel velocity direction in the original phase space coordinate system according to the electric field component in the z-direction, the ionic charge, and the ionic mass. Specifically, this includes: based on the Hamiltonian equation in non-canonical coordinates, the equation of motion of the cyclotron center is derived through Poisson brackets of coordinates and Hamiltonian, containing only the first derivative of time throughout, ensuring computational efficiency and accuracy; the velocity in the x-direction is defined as the first derivative of the x-coordinate with respect to time, i.e. This velocity is determined by x and the Hamiltonian. Poisson brackets Derivation, combined with uniform magnetic field and the electric field component in the y direction Ultimately The velocity in the y-direction is defined as the first derivative of the y-coordinate with respect to time, i.e. , through y and Poisson brackets Combined with a uniform magnetic field and electric field components in the x direction The derivation is as follows Parallel velocity in the z-direction Defined as the first derivative of the z-coordinate with respect to time. Directly from z and Poisson brackets Derivation; its rate of change over time is obtained through and Poisson brackets Derivation, combined with the z-direction electric field component The ionic charge q and ionic mass m were then used to derive the final result. .

[0038] In a preferred embodiment of the present invention, step 3, based on the equation of motion, transforming the original phase space coordinate system into a new phase space coordinate system through Lie transformation to obtain the generated vector field, may include:

[0039] Step 301: Based on the equations of motion, extract the Hamiltonian describing the energy state of the cyclotron center. Specifically, the Hamiltonian H of the cyclotron center is a physical quantity describing the total energy of the system, which consists of radial cyclotron energy, axial kinetic energy, and electric potential energy. Simultaneously, the potential distribution is constrained by the geometric parameters of the discharge chamber to ensure accurate energy calculation. The kinetic energy term reflects the energy possessed by the ions due to their motion, where the radial cyclotron energy is composed of the magnetic moment μ and the axial uniform magnetic field. It is determined that (the center of rotation does not participate in circular motion, and radial energy is quantized through magnetic moment), axial motion energy is proportional to the square of the parallel velocity, the ion mass is m, and the parallel velocity in the z-direction is... (Unit: m / s), therefore the kinetic energy term is The potential energy term reflects the energy an ion possesses due to its position in an electric field. It is calculated based on the principle that potential energy in an electric field equals the product of charge and potential. The ionic charge is q (for a single-charge ion, q = 1.6 × 10⁻⁶). -19 C), the potential distributed in three-dimensional space is (x, y, z) (unit V, determined by the electric field of the discharge chamber), therefore the potential energy term is q. (x, y, z); the complete Hamiltonian expression is: .

[0040] The discharge chamber has a cylindrical structure. Geometric parameters must first be calculated to determine the spatial range of the potential distribution. The volume V describes the three-dimensional size of the discharge chamber, and its calculation formula is V = πr. 2 L, where r is the radial radius of the discharge chamber (distance from the inner wall to the z-axis, in meters, e.g., r = 0.1 m), and L is the axial length (axial distance from the cathode to the anode, in meters, e.g., L = 0.5 m). The volume V determines the spatial boundary for potential calculation. The lateral surface area S describes the cylindrical surface range of the inner wall of the discharge chamber, calculated as S = 2πrL. The lateral surface area corresponds to the position of the anode (the anode is usually attached to the inner wall of the discharge chamber), therefore the potential must satisfy the boundary conditions at S. The constraint effect of geometric parameters limits the spatial range of potential calculation through V and sets the boundary conditions through S to ensure... The distribution of (x, y, z) inside the discharge chamber conforms to the actual electric field environment, thus ensuring the potential energy term q in the Hamiltonian. The calculation is accurate (if the boundary conditions are incorrect, the error in the potential energy term will lead to deviations in subsequent energy analysis).

[0041] Step 302a involves separating the Hamiltonian into zero-order and first-order Hamiltonian terms; and establishing an integral expression for the basic cyclotronic trajectory along the cyclotron center in a uniform magnetic field. Specifically, the electric field of the Penning ion source is predominantly axial, while the radial electric field is constrained by the magnetic field, causing only minor perturbations to the cyclotron center's motion. Therefore, the three-dimensional potential is... (x, y, z) can be decomposed into axial potentials (z), which varies only in the z-direction (independent of x and y), represents the dominant part of the electric field, such as from the cathode (z= , ( )=-V0, V0=50V~100V) to the anode (z= , ( The axial potential gradient at (0V) is the main energy source for axial acceleration of ions; the separation and physical meaning of the Hamiltonian term: According to the potential decomposition results, the Hamiltonian is separated into a zero-order Hamiltonian term. It is the dominant energy term of the system, and its calculation formula is: Among them, q (z) provides energy for axial acceleration of ions (such as when ions move from the cathode to the anode, potential energy is converted into kinetic energy). The energy corresponding to the radial cyclotron motion of the ion. The energy corresponding to the axial motion of the ions. The main laws of motion of the center of rotation (axial acceleration and radial rotation) are determined; the first-order Hamiltonian term. It is the system's disturbance energy term, calculated as follows: = (x, y, z).

[0042] The radial motion of the gyroscope center in the original phase space (old coordinates) is an electric drift motion, not a gyroscope motion at the particle level. The influence of this electric drift motion will be incorporated into the transformation relationship through subsequent coordinate transformations. In the transformed new phase space coordinate system, the radial (x, y directions) motion of the gyroscope center is completely eliminated, and it remains stationary. The core of establishing the integral expression is to quantify the cumulative effect of radial perturbation in the original phase space: the first-order Hamiltonian. The periodic variation of (x, y, z) with x and y originates from the radial electric drift motion of the cyclotron center, and the integral aims to capture the effect of this perturbation; therefore, the integral expression is set as the integrand, i.e., the first-order Hamiltonian. (x(t), y(t), z(t)), where (x(t) and y(t) are the coordinates of the radial electric drift trajectory of the center of rotation in the original phase space as a function of time), z(t) is the change of the axial position as a function of time; the integration variable is time t.

[0043] Step 302c involves integrating the integral expression along the trajectory of the cyclotron center in the electromagnetic field to obtain the gauge function for coordinate transformation. Specifically, the core of the integral calculation is obtaining the gauge function through trajectory solving and numerical integration. First, the first-order Hamiltonian at each time step needs to be determined. (t); (t)=q (x(t), y(t), z(t)), therefore, based on the non-regular coordinate Hamiltonian equation established in step 202 ( , , , The fourth-order Runge-Kutta method is used to calculate x(t), y(t), and z(t) within t∈[0,τ]. It is important to note that z(t) must be calculated simultaneously with x(t) and y(t). After subsequent coordinate transformations, the rates of change in the x and y directions are zero, and their trajectory evolution does not require additional tracking. However, the rate of change in the z direction is non-zero, making it a core dimension that needs to be focused on during trajectory calculation. The potential value is obtained by substituting the average values ​​of x(t), y(t), and z(t) into the three-dimensional potential distribution. (x,y,z), to obtain the corresponding value for each t. (t), and thus obtain (t)=q (t).

[0044] Numerical integration uses the trapezoidal rule, and the specific operation is as follows: time discretization divides the integration interval into [0,τ] (integrating once for each time step). The time interval is uniformly divided into N small time intervals Δt = (N is 8 or 16, matching the accuracy of trajectory solving within the time step), to obtain discrete time steps. =iΔt(i=0,1,2,...,N,where =0, =τ); the function takes values ​​based on the calculation of each discrete time step. corresponding ( )=q (x( ),y( ),z( The integral calculation yields the standard function using the trapezoidal rule formula:

[0045] , As a core parameter, the radial motion effect in the original phase space can be incorporated into the transformation relationship, simplifying the radial motion of the gyroscope center in the new phase space.

[0046] Step 303: Based on the gauge function, obtain the first-order generating vector through Poisson bracket operations; using the first-order generating vector, transform the coordinates in the original phase space coordinate system to the coordinates in the new phase space coordinate system to obtain the generating vector field. Specifically, this includes: Poisson brackets are core mathematical operators in analytical mechanics, primarily used to solve equations of motion. The Poisson brackets between coordinates and the Hamiltonian H represent the equation of motion for that coordinate, and the Poisson brackets between a physical quantity and H represent the rate of change of that physical quantity over time. Their definition is based on the fundamental form of the system. ,satisfy Z represents the original phase space coordinates. It is a basic form component.

[0047] In this invention, the original phase space coordinates are the three-dimensional real space coordinates (x, y, z) of the gyroscope center and the one-dimensional parallel velocity. and magnetic moment (The core dependent variable of the distribution function), the normalization function As the core physical quantity A in the calculation, it is represented by the five coordinates (x, y, z, ...) of the original phase space. μ) is used as the physical quantity B, through Poisson brackets The calculation yields four first-order generated vectors, each corresponding to a coordinate transformation deviation; the specific calculation takes the key vector as an example: when B = x (original x-coordinate), the Poisson brackets... Based on the basic form The variational relation derivation, combined with the correlation of partial derivatives with respect to y and z, ultimately yields a vector reflecting the deviation of the x-coordinate transformation; similarly, when B=y, the result is obtained through... The calculation yields the generated vector corresponding to the y-coordinate; when B=z, it is obtained through { The z-coordinate is obtained by generating the vector corresponding to the z-coordinate. The generation vector corresponding to the parallel velocity is obtained.

[0048] The above four first-order generating vectors are denoted as Together, they form a generating vector field; using this vector field, the original phase space coordinates (x, y, z, ...) are transformed. , The transformation to new phase space coordinates, retaining the transformation relationship to second-order quantities, is as follows: ;

[0049] ;

[0050] ;

[0051] ; (in Generate first-order vectors for each dimension. (The partial derivatives with respect to the corresponding coordinates are omitted here according to Einstein's summation convention, thus completing the coordinate mapping.)

[0052] In this embodiment, the generated vector field is obtained through Poisson bracket operations and coordinate transformation is completed, realizing a precise conversion from the original phase space to the new phase space, which creates conditions for the simplification of subsequent motion equations.

[0053] In a preferred embodiment of the present invention, step 4, constructing simplified motion equations in the new phase space coordinate system based on the generated vector, may include:

[0054] Step 401: Based on the generated vector field, establish the transformation relationship from the original phase space coordinate system to the new phase space coordinate system; according to the transformation relationship, determine the characteristic that the gyration center in the new phase space coordinate system remains stationary in the x and y directions, specifically including: based on the obtained generated vector field, corresponding to the original phase space coordinates (x, y, z, ... The trajectory deviation (magnetic moment μ remains unchanged in the cyclotron center coordinates and requires no transformation); the second-order transformation relationship between the new phase space coordinates and the original phase space coordinates includes the new phase space x-coordinate: (where x is the x-coordinate of the original phase space,) First-order direction generating vector, The second-order correction term for the direction; the y-coordinate of the new phase space: Where y is the y-coordinate of the original phase space; and z-coordinate of the new phase space. Where z is the z-coordinate of the original phase space; the parallel velocity in the new phase space is: ,in The original phase space parallel velocity; based on the above transformation relationship, the motion characteristics are analyzed: the first-order generated vector in the x-direction. Its core function is to counteract the radial electric field in the x-direction of the original phase space. Trajectory deviation caused by disturbance; first-order generated vector in the y-direction and second-order correction terms Used to counteract the radial electric field in the y-direction of the original phase space The trajectory deviation caused by the disturbance; therefore, the derivatives of the x and y coordinates with respect to time in the new phase space are both zero, that is, the gyrocenter remains stationary in the x and y directions of the new phase space.

[0055] Step 402: Based on the characteristics, the equations for the zero velocity of the gyroscope center along the x-direction in the new phase space coordinate system are obtained, as well as the equations for the zero velocity of the gyroscope center along the y-direction in the new phase space coordinate system. Specifically, the equations of motion in the x-direction are the equations for the x-coordinates of the new phase space coordinate system. It does not change with time t. According to the physical definition of velocity (velocity is the derivative of position with respect to time), the velocity in the x-direction of the new phase space is... The derivative with respect to t, i.e. ;because Since the derivative is constant, it is zero. Therefore, the equation of motion in the x-direction is: The equation of motion in the y-direction is the y-coordinate of the new phase space. The velocity in the y-direction of the new phase space does not change with time t. The derivative with respect to t, i.e. = ;because Since the derivative is constant, the equation of motion in the y-direction is zero. = =0.

[0056] Step 403: Based on the transformation relationship, transform the original equation of motion along the z-direction into the equation of motion of the gyroscope center along the z-direction in the new phase space coordinate system; based on the transformation relationship, transform the equation of motion parallel to the velocity direction in the original equation of motion into the dynamic equation of the parallel velocity change in the new phase space coordinate system. Specifically, this includes: clarifying that the physical definition of the velocity in the z-direction in the original phase space is the rate of change of position with respect to time; therefore, the equation of motion along the z-direction is... ,in is the velocity (in meters per second) of the original phase space gyroscope center along the axial direction (z-direction), z is the axial position coordinate of the gyroscope center in the original phase space (in meters), and t is the time variable (in seconds); according to the transformation relationship established in step 401, the z coordinates of the new and original phase spaces are transformed into... Parallel velocity transformation Based on the scalar invariance of the Lie transformation (ensuring consistency between the action and Hamiltonian transformations, i.e., coordinate transformations do not change the essence of physical quantities), the core idea of ​​the derivation is as follows: first, clarify the form of the target motion equations to be achieved in the new phase space, then determine the transformation relationship that satisfies the equations (i.e., first lock the motion equations before and after the transformation, then derive the transformation relationship in reverse); combining the generation vector compensation characteristics derived from scalar invariance, and the reverse mapping relationship between the original and new phase space physical quantities, the motion equations in the z-direction of the new phase space are finally obtained. This directly reflects the correspondence between the rate of change of axial position and axial velocity in the new phase space; in the original phase space, the equation of motion for parallel velocities is defined based on the Poisson brackets of analytical mechanics. ,in ( It is a zero-order Hamiltonian, determined by the axial electric field and the uniform magnetic field; The first-order Hamiltonian is obtained by subtracting the zeroth-order Hamiltonian from the total Hamiltonian; according to the Hamiltonian transformation relation of the Lie transform, the new phase space Hamiltonian is... ( (The equations of motion for the parallel velocities in the new phase space are also in Poisson bracket form, i.e., canceled out by the generated vectors and without radial perturbation terms). This form originates from the invariance of analytical mechanics laws under coordinate transformations, and is not derived by substituting coordinate transformation relationships; because Only with the axial electric field of the new phase space Relatedly, and combining the Poisson bracket operation rules (only the axial electric field affects the parallel velocity), the final simplified result is: Where q is the ionic charge and m is the ionic mass. Let be the electric field intensity component along the z-direction in the new phase space.

[0057] Step 404: Combine the equations for zero velocity in the x-direction, zero velocity in the y-direction, the equations for motion in the z-direction, and the dynamic equations for parallel velocity changes to form simplified equations of motion in the new phase space coordinate system. Specifically, this includes: (as in step 401) The core function of the generated vector is to synchronize x, y, z, The four dimensions are transformed to cancel out the motion in the x-direction of the original phase space caused by the radial electric field disturbance, so that the x-coordinate of the new phase space is... It does not change with time; based on the physical definition of velocity (velocity is the rate of change of position with respect to time), the equation of motion in the x-direction is derived as follows: This equation directly quantifies the characteristic of no motion in the x-direction of the new phase space; consistent with the logic in the x-direction, the transformation relationship in the y-direction... The generated vector in the middle is used to eliminate the motion in the y-direction of the original phase space caused by the radial electric field disturbance, so that the y-coordinate of the new phase space is... Remain stationary; based on the physical definition of velocity, the equation of motion in the y-direction is derived as follows: The quantification of the motionless characteristics in the y-direction of the new phase space is performed; the equation of motion in the z-direction derived in step 403 is used. This equation describes the axial position in the new phase space. The rate of change of time is equal to the axial velocity To clarify the dynamic relationship between the two ( For the axial coordinates of the new phase space, (This refers to the axial velocity of the new phase space).

[0058] The parallel velocity change dynamic equation derived in step 403 The equation is based on the Hamiltonian of the new phase space after the Lie transform. = (Contains only axial electric field and uniform magnetic field, without radial perturbation) Derivation: It is parallel velocity The rate of change over time, where q is the ionic charge and m is the ionic mass. The axial electric field intensity component in the new phase space reflects the change in axial velocity solely due to Driven by no additional disturbance; combining the above four equations according to the logic of radial stationary and axial motion, we obtain the simplified motion equations in the new phase space: , , , .

[0059] This embodiment transforms the complex motion in the original phase space into a simplified motion in the new phase space through a clear transformation relationship, providing a foundation for dimensionality reduction calculation and ensuring that the simplified logic of subsequent steps can be implemented.

[0060] In a preferred embodiment of the present invention, step 5, based on the simplified equation of motion, employs the method of characteristics and performs interpolation operations on the two dimensions of axial position and parallel velocity to obtain the distribution function in the new phase space coordinate system, may include:

[0061] Step 501: Based on the simplified equations of motion in the new phase space coordinate system, determine the trajectory of the gyroscope center in the new phase space. Specifically, this includes: determining the initial state parameters of the gyroscope center in the new phase space, and the initial time... When =0, axial position , parallel speed and the initial x' direction coordinates , Direction coordinates Because of the simplification of the equations of motion , Therefore , Keeping the initial values ​​constant, the trajectory only needs to be tracked. and The evolution of the time step Δt is defined, and its value must satisfy the Courant-Friedrichs-Lewy (CFL) condition. The core is to ensure numerical stability and prevent trajectory tracking from skipping grid points. For axial motion scenarios, the CFL condition is expressed as Δt ,in For the follow-up Dimensional grid spacing, The characteristic value of the parallel velocity in the new phase space.

[0062] The initial value for the iteration is determined by the initial time step. =0 (Initial axial position of the discharge chamber, such as 0mm) (Initial parallel velocity, such as 0) is the initial value. , Keeping the initial value (e.g., 3mm) unchanged, the number of iterations n increases from 0, and each iteration corresponds to a specific time step. = (Δt is a fixed time step that satisfies the CFL condition, such as...) s); The second-order Runge-Kutta method for intermediate quantity calculation includes the first step of intermediate quantity calculation ( , ):Will Moment , Substitute into the equation of motion and calculate = ( The first-order increment corresponds to ), = ( The first-order increment corresponds to Where q is the charge of a single ion. m is the mass of the proton. , For the axial electric field components of the new phase space); intermediate quantities in the second step ( , ):calculate Estimated value of time , Substituting into the equation of motion, we get , ( Since it is a constant, therefore ).

[0063] Iteratively updating trajectory parameters is the calculation Moment , Since the trajectory is the same at each step (no need for multiple iterations, only to determine a fixed trajectory shape), only 1-2 key iterations are required (e.g., from...). arrive ), then you can determine and With a fixed correspondence, the trajectory can be reused in each subsequent simulation step, only the distribution function of the interpolation is changed; the results of key iterations are recorded. and Corresponding values, combined with fixed values , The initial value, which constitutes the gyrocenter in the new phase space. - The complete motion trajectory in a two-dimensional plane requires no redundant calculation of multi-time data.

[0064] Step 502: Based on the motion trajectory, establish the solution path of the method of characteristics; based on the solution path, determine the grid points in the two dimensions of axial position and parallel velocity that need to be interpolated. Specifically, Step 502 establishes the solution path of the method of characteristics based on the motion trajectory obtained in Step 501. The fixed trajectory of each rotation center is directly used as a characteristic line, and all characteristic lines together constitute the solution path. The coordinates of the points on the characteristic lines are... It accurately reflects the relationship between the axial position of the gyroscope center and the parallel velocity in the new phase space, and divides the two-dimensional interpolation mesh based on the solution path, including clearly defining the new phase space. and The range of values ​​for (relative to the original phase space z, Consistent with the actual dimensions of the discharge chamber), where The range is from 0 to L (L is the classic axial length of the Penning ion source discharge chamber, which is 20 mm). for meters per second (covering the actual thermal velocity and axial motion velocity of the ions); in Divide the dimensions into P=200 uniform intervals (not exceeding 200), resulting in P+1=201 grid nodes. The node coordinates are... =0、 = , = … ,in =0.1mm ( for (Dimensional grid spacing); in The dimensions are divided into Q=64 grid nodes, corresponding to a grid spacing of Q−1=63, and the node coordinates are... = meters per second = … meters per second, of which ≈1587.3 m / s; each intersection of the two-dimensional grid is a grid point, with coordinates as ( , (i=1 to 201, j=1 to 64).

[0065] Determine the known points and the interpolation target points: The values ​​of the initial distribution function are preset on the grid points, which are the known points. The initial values ​​of the distribution function are directly recorded. When solving using the method of characteristics, it is necessary to trace back... The point in phase space where the center of rotation lies on the characteristic line ( , This point typically does not coincide with grid points, making it impossible to directly obtain the distribution function value. Such points are called interpolation target points. Interpolation is required to calculate the distribution function value of the target points, ensuring that all target points are within the grid. ∈[0,L]、 ∈[ , If the interpolation point exceeds the range of meters per second, the interpolation function value is set to zero, indicating that no ions are distributed to that phase space point.

[0066] Step 503 involves performing interpolation operations on the grid points in both axial position and parallel velocity dimensions to obtain the distribution function in the new phase space coordinate system. Specifically, this includes: first preparing the basic data for B-spline interpolation, i.e., quadratic B-spline interpolation (the basis functions are second-order polynomials, satisfying second-derivative continuity), and then... , The grid nodes of dimension are used to construct the interpolation basis functions and solve for the interpolation coefficients; where, Dimensional basis functions are denoted as ( (i is the index of the basis function, corresponding to) Dimensional grid). Dimensional basis functions are denoted as ( (j is the index of the basis function, corresponding to) Dimensional grid); basis functions use a sequence of heavy nodes ( = = = , = = = The intermediate nodes are evenly distributed in [ ]), and satisfy , Accuracy requirements; interpolation coefficients Using the initial values ​​of the distribution function of the known grid points, combined with the interpolation conditions Solving for the given values; collecting initial values ​​of known points: known grid points ( , Initial value of the distribution function ( , Preset completion (reflecting ion number density, such as axial inlet) =0、 =0 ( =0, =0)=1×10 18 m -3 s 3 , Increase the attenuation according to the Gaussian distribution).

[0067] Perform two-dimensional B-spline interpolation calculations: for any interpolation target point ( , Substituting into the two-dimensional B-spline interpolation function S( , )= (Summation range coverage) Dimension of all basis function indices i and (For all basis function indices j in dimension j), the distribution function value of the target point can be obtained directly. ( , )=S( , ); Integrate to obtain the complete distribution function: for all interpolation target points With known grid points Integration to form a complete distribution function in the new phase space This function can describe different (Axial position) and Distribution pattern of gyroscope centers at (parallel velocities).

[0068] This embodiment clearly defines the interpolation dimensions and grid point settings, avoids dimensional redundancy, provides structured data support for interpolation operations, and ensures interpolation accuracy.

[0069] In a preferred embodiment of the present invention, step 6, mapping back to the original phase space coordinate system based on the distribution function in the new phase space coordinate system through a pullback transformation of the Lie transform, and finally obtaining the original cyclotron center distribution function, may include:

[0070] Step 601: Based on the generated vector field, determine the pull-back transformation relationship from the new phase space coordinate system to the original phase space coordinate system; based on the distribution function in the new phase space coordinate system, and in conjunction with the pull-back transformation relationship, obtain the first-order and second-order correction terms of the distribution function in the original phase space coordinate system, specifically including: determining the pull-back transformation relationship: the generated vector field includes the generated vector in the x-direction. (Unit: meters, range from 0 to 1×10) -3 Meters, correcting radial perturbation in the x-direction, and generating vectors in the y-direction. (Unit: meters, range from 0 to 1×10) -3 Meters, correcting radial perturbation in the y-direction, and generating vector in the z-direction. (Unit: meters, range from 0 to 1×10) -3 Meters, correcting axial disturbance in the z-direction), parallel velocity generation vector. (Unit: meters per second, range from 0 to 1 × 10⁻⁶) 3 m / s (correcting for parallel velocity perturbations); the core of the pull-back transformation is to convert the physical quantities in the new phase space into the physical quantities in the original phase space. Specific relationships include the original x-coordinate:

[0071] x is the x-coordinate of the original phase space (unit: meters). The x-coordinate of the new phase space (unit: meters);

[0072] Original y-coordinate: y is the original phase space y-coordinate (unit: meters). The y-coordinate of the new phase space (unit: meters);

[0073] Original z-coordinate: z is the z-coordinate of the original phase space (unit: meters). The z-coordinate of the new phase space (unit: meters);

[0074] Original parallel velocity: , It is the original phase space parallel velocity (unit: meters per second). The parallel velocity in the new phase space (unit: m / s); then the first-order correction term of the distribution function in the original phase space is calculated. Based on the first-order term rule of scalar pullback transform, the original phase space distribution function With the new phase space distribution function The first-order association is: Higher-order terms, including first-order correction terms The first-order partial derivatives are calculated using the central difference method. The specific process is as follows: Dimensional division =200 intervals, maximum radial distance R=0.1 meters, grid spacing meters; non-boundary points (i=2 to 200, i is Dimension node index) uses central difference: Boundary point i=1 ( =0) Use forward differencing: Boundary point i=201 ( =R) uses backward differencing: Calculation method and Consistent Dimensional division =200 intervals, grid spacing rice; , The calculation logic is the same, only the grid spacing and node index of the corresponding dimension are replaced.

[0075] Recalculate the second-order correction term The formula is (where each second-order mixed partial derivative is calculated through two first-order central differences, with...) For example, for Taking the first-order partial derivative, we obtain all ; then Find the first-order partial derivatives, using central difference for non-boundary points (j=2 to 200): The boundary points are calculated using forward or backward differences; the calculation method for the remaining second-order mixed partial derivatives is similar, and the coefficient 21 is used to correct the superposition effect of second-order perturbations.

[0076] Step 602: Superimpose and combine the distribution function in the new phase space coordinate system with the first-order and second-order correction terms to obtain the complete distribution function in the original phase space coordinate system; based on the distribution function, obtain the final original cyclotron center distribution function, specifically including: superimposing and combining the distribution functions, and based on the pull-back transformation relationship, obtaining the new phase space distribution function. The original phase space coordinates are directly mapped (no additional mapping steps are required), and the complete distribution function calculation formula is as follows: ,in For first-order correction terms, The second-order correction term, after superposition, fully reflects the distribution characteristics of the original phase space cyclotron center; the superposition... It is the complete distribution function in the original phase space coordinate system, which can accurately describe the ion number density distribution at different radial positions, axial positions and parallel velocities in the original phase space, and can be directly used for plasma physics analysis or numerical simulation.

[0077] In this embodiment, a pullback transformation is established and a correction term is calculated to ensure the mapping accuracy from the new phase space to the original phase space and to avoid the transformation error affecting the reliability of the result.

[0078] Embodiments of the present invention also provide a computing device, including: a processor and a memory storing a computer program, wherein the computer program, when executed by the processor, performs the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0079] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0080] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. An improved solution method for cyclotron kinetics simulation of the Penning ion source, characterized in that, The method includes: Step 1: Obtain the electromagnetic field configuration of the Penning ion source, which includes a uniform magnetic field pointing along the axis and a three-dimensional spatially distributed potential. Step 2: Based on the electromagnetic field configuration, establish the equation of motion of the gyroscope center in the original phase space coordinate system; Step 3: Based on the equations of motion, the original phase space coordinate system is transformed into a new phase space coordinate system through Lie transformation to obtain the generated vector field. This includes: extracting the Hamiltonian describing the energy state of the cyclotron center based on the equations of motion; constructing a gauge function for coordinate transformation based on the Hamiltonian; obtaining a first-order generated vector through Poisson bracket operations based on the gauge function; and using the first-order generated vector to transform the coordinates in the original phase space coordinate system into the coordinates in the new phase space coordinate system to obtain the generated vector field. Step 4: Based on the generated vector field, determine the transformation relationship between the original phase space coordinate system and the new phase space coordinate system, and obtain the simplified equations of motion in the new phase space coordinate system. This includes: determining the transformation relationship from the original phase space coordinate system to the new phase space coordinate system based on the generated vector field; determining the characteristic that the gyroscope center remains stationary in the x and y directions in the new phase space coordinate system based on the transformation relationship; obtaining the equation for the gyroscope center having zero velocity in the x direction and the equation for the gyroscope center having zero velocity in the y direction in the new phase space coordinate system based on the transformation relationship; transforming the equation of motion along the z direction in the original equations of motion into the equation of motion along the z direction in the new phase space coordinate system based on the transformation relationship; transforming the equation of motion parallel to the velocity direction in the original equations of motion into the dynamic equation of parallel velocity change in the new phase space coordinate system based on the transformation relationship; and combining the equations for zero velocity in the x direction, zero velocity in the y direction, the equations of motion in the z direction, and the dynamic equation of parallel velocity change to form the simplified equations of motion in the new phase space coordinate system. Step 5: Based on the simplified equation of motion, the method of characteristics is used and interpolation is performed on the two dimensions of axial position and parallel velocity to obtain the distribution function in the new phase space coordinate system. Step 6: Based on the distribution function in the new phase space coordinate system, the original phase space coordinate system is mapped back through the pullback transformation of the Lie transform, and finally the original cyclotron center distribution function is obtained.

2. The improved solution method for cyclotron kinetics simulation of the Penning ion source according to claim 1, characterized in that, The electromagnetic field configuration of the Penning ion source is obtained, which includes a uniform magnetic field pointing along the axis and a three-dimensional spatially distributed electric potential, including: The geometric parameters of the discharge chamber of the Penning ion source are received, including the axial length and radial radius of the discharge chamber, as well as the geometric positions of the cathode and anode; based on the geometric parameters of the discharge chamber, a uniform magnetic field is set along the axial direction. The geometric parameters of the discharge chamber and the uniform magnetic field are input into the electric field solver, and the potential distributed in three-dimensional space is calculated based on the electric field solver. The uniform magnetic field is combined with the calculated potential distributed in three-dimensional space to generate the electromagnetic field configuration.

3. The improved solution method for cyclotron kinetics simulation of the Penning ion source according to claim 2, characterized in that, Based on the electromagnetic field configuration, the equations of motion for the gyroscope center in the original phase space coordinate system are established, including: Based on a uniform magnetic field and a three-dimensional spatially distributed electric potential, the electromagnetic force on the center of rotation is calculated. Based on electromagnetic force, the equations of motion of the gyrocenter along the x-direction in the original phase space coordinate system are established according to the electric field component in the y-direction and the uniform magnetic field. Based on electromagnetic force, the equations of motion of the gyrocenter along the y-direction in the original phase space coordinate system are established according to the electric field component in the x-direction and the uniform magnetic field. Based on electromagnetic force, the velocity of the gyrocenter along the z-direction in the original phase space coordinate system is expressed as a parallel velocity, and the equations of motion of the gyrocenter along the parallel velocity direction in the original phase space coordinate system are established according to the electric field component in the z-direction, the ionic charge, and the ionic mass.

4. The improved solution method for cyclotron kinetics simulation of the Penning ion source according to claim 3, characterized in that, Based on the Hamiltonian, construct the normalization functions for coordinate transformation, including: The Hamiltonian is divided into zero-order Hamiltonian terms and first-order Hamiltonian terms; and an integral expression for the basic cyclotronic motion trajectory along the center of cyclotron in a uniform magnetic field is established. Integrating the integral expression along the trajectory of the gyrocenter in the electromagnetic field yields the gauge function used for coordinate transformation.

5. The improved solution method for cyclotron kinetics simulation of the Penning ion source according to claim 4, characterized in that, Based on the simplified equations of motion, the method of characteristics is used, and interpolation is performed on the axial position and parallel velocity dimensions to obtain the distribution function in the new phase space coordinate system, including: Based on the simplified equations of motion in the new phase space coordinate system, the trajectory of the gyroscope center in the new phase space is determined. Based on the motion trajectory, establish the solution path of the method of characteristics; based on the solution path, determine the grid points in the two dimensions of axial position and parallel velocity that need to be interpolated. Interpolation is performed on the grid points in both axial position and parallel velocity dimensions to obtain the distribution function in the new phase space coordinate system.

6. The improved solution method for cyclotron kinetics simulation of the Penning ion source according to claim 5, characterized in that, Based on the distribution function in the new phase space coordinate system, the original phase space coordinate system is mapped back through the pullback transformation of the Lie transform, ultimately yielding the original cyclotron center distribution function, including: Based on the generated vector field, the pull-back transformation relationship from the new phase space coordinate system to the original phase space coordinate system was determined; based on the distribution function in the new phase space coordinate system, the first-order correction term and the second-order correction term of the distribution function in the original phase space coordinate system were obtained by combining the pull-back transformation relationship. The distribution function in the new phase space coordinate system is superimposed with the first-order and second-order correction terms to obtain the complete distribution function in the original phase space coordinate system; based on the distribution function, the final original cyclotron center distribution function is obtained.

7. A computing device, characterized in that, include: One or more processors; A storage device for storing one or more programs that, when executed by one or more processors, cause the one or more processors to implement the method as described in any one of claims 1 to 6.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 6.

Images