Method for reconstructing beam current loss heat source of electrostatic deflection plate of cyclotron

Abstract

The application provides a cyclotron electrostatic deflection plate beam loss heat source reconstruction method, and relates to the technical field of particle accelerator heat analysis and multi-physical field coupling calculation, and the method comprises the following steps: step 1, configuring beam dynamics simulation parameters, performing multi-particle tracking under a three-dimensional field diagram containing a central area electromagnetic field, a harmonic coil magnetic field and an electrostatic deflection electric field, and generating a particle loss list; step 2, converting the particle loss list into a beam source file of a Monte Carlo radiation transport program, configuring Monte Carlo parameters, including deflection plate materials, boundaries, particle types and energy spectra, performing particle transport simulation, and obtaining a three-dimensional energy deposition matrix on the electrostatic deflection plate element. The application realizes high-resolution heat source reconstruction, can reduce the calculation error of the peak temperature rise of the deflection plate, and can feed back the temperature field solving result to optimize the particle source, thereby providing a reliable basis for the thermal structure design and cooling optimization of the electrostatic deflection plate.

Description

Technical Field

[0001] This invention relates to the field of particle accelerator thermal analysis and multiphysics coupling calculation technology, and in particular to a method for reconstructing the heat source of beam loss in an electrostatic deflector plate of a cyclotron. Background Technology

[0002] In superconducting cyclotron accelerators, electrostatic deflection plates are key components of the beam extraction system. Their thermal performance has a significant impact on the stability and lifespan of the accelerator. During the deflection process, beam losses caused by various factors can lead to localized heat deposition on and inside the deflection plate. If the heat source distribution is not accurately assessed, it may pose a challenge to the rationality of the cooling system design.

[0003] Currently, in engineering, the thermal analysis of electrostatic deflection plates sometimes uses simplified heat source models for approximate estimation. For example, some methods approximate beam loss as a Gaussian distributed heat source. This method is relatively simple to model, but it may not accurately capture the details of local heat deposition in terms of reflecting the non-uniform distribution characteristics and energy spectrum differences of the actual beam loss. This limits the precise optimization of thermal management schemes. Existing heat source modeling methods still have room for improvement in terms of spatial resolution and energy distribution accuracy, especially under high-power beam conditions, where the need for refined modeling of heat source distribution has not been fully met. Summary of the Invention

[0004] The technical problem to be solved by this invention is to provide a method for reconstructing the heat source of beam loss in an electrostatic deflector plate of a cyclotron, thereby solving the problem of inaccurate heat source distribution affecting the thermal design and operational safety of the deflector plate.

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

[0006] A first aspect is a method for reconstructing the heat source of beam loss in a cyclotron electrostatic deflector plate, the method comprising:

[0007] Step 1: Configure beam dynamics simulation parameters, perform multi-particle tracking under a three-dimensional field diagram including the electromagnetic field of the central region, the magnetic field of the harmonic coil, and the electrostatic deflection electric field, and generate a particle loss list.

[0008] Step 2: Translate the particle loss list into a beam source file for the Monte Carlo radiation transport program, configure the Monte Carlo parameters, including deflection plate material, boundary, particle type and energy spectrum, execute particle transport simulation, and obtain the three-dimensional energy deposition matrix on the electrostatic deflection plate element.

[0009] Step 3: The weighted octree spatial partitioning algorithm is used to map the three-dimensional energy deposition matrix to the computational fluid dynamics grid to generate the volumetric heat source file;

[0010] Step 4: Load the volumetric heat source file into the fluid dynamics calculation software to solve the temperature field, obtain the solution results, and feed the solution results back to optimize the particle source.

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

[0012] One or more processors;

[0013] 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.

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

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

[0016] Multi-particle tracking is performed under a three-dimensional composite field map containing the electromagnetic field of the central region, the magnetic field of the harmonic coil, and the electrostatic deflection electric field. Using the coordinates and kinetic energy of the actual lost particles as input, the spot-like spatial distribution and energy spectrum differences of the beam loss are completely preserved, avoiding heat source distribution distortion from the outset. By employing a fourth-order Runge-Kutta numerical integration algorithm, the trajectory of particles in the composite field can be calculated more accurately, effectively ensuring the accuracy of the determination that the particle trajectory intersects with the deflection plate boundary, reducing the number of missed or misjudged lost particles, and thus improving the data reliability of the generated particle loss list. This particle loss list directly records the spatial coordinates and kinetic energy of each lost particle. When real information is transcribed into the beam source file of the Monte Carlo radiation transport program, no additional idealization assumptions are needed. It can directly reflect the location differences and energy distribution characteristics of actual beam loss, reducing the calculation error of the subsequent three-dimensional energy deposition matrix from the input end. By tracking loss particles and generating a high-confidence loss list, it can provide a reliable data foundation for the construction of the three-dimensional energy deposition matrix, volume heat source mapping and temperature field solution. Ultimately, it can effectively reduce the calculation error of peak temperature rise, and provide accurate heat source basis for the layout of cooling channels, material selection and structural strength optimization of electrostatic deflection plates, thereby avoiding the shortened device life or operational failure caused by inaccurate heat source assessment. Attached Figure Description

[0017] Figure 1 This is a schematic flowchart of the method for reconstructing the heat source of beam loss in a cyclotron electrostatic deflector plate provided in an embodiment of the present invention.

[0018] Figure 2 This is a flowchart illustrating step 3 provided in an embodiment of the present invention. Detailed Implementation

[0019] 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.

[0020] like Figure 1 As shown, an embodiment of the present invention proposes a method for reconstructing the heat source of beam loss in a cyclotron electrostatic deflector plate, the method comprising the following steps:

[0021] Step 1: Configure beam dynamics simulation parameters, perform multi-particle tracking under a three-dimensional field diagram including the electromagnetic field of the central region, the magnetic field of the harmonic coil, and the electrostatic deflection electric field, and generate a particle loss list.

[0022] Step 2: Translate the particle loss list into a beam source file for the Monte Carlo radiation transport program, configure the Monte Carlo parameters, including deflection plate material, boundary, particle type and energy spectrum, execute particle transport simulation, and obtain the three-dimensional energy deposition matrix on the electrostatic deflection plate element.

[0023] Step 3: The weighted octree spatial partitioning algorithm is used to map the three-dimensional energy deposition matrix to the computational fluid dynamics grid to generate the volumetric heat source file;

[0024] Step 4: Load the volumetric heat source file into the fluid dynamics calculation software to solve the temperature field, obtain the solution results, and feed the solution results back to optimize the particle source.

[0025] In this embodiment of the invention, multi-particle tracking is performed under a three-dimensional composite field map containing the central region electromagnetic field, the harmonic coil magnetic field, and the electrostatic deflection electric field. Simultaneously, using the coordinates and kinetic energy of the actual lost particles as input, the spot-like spatial distribution and energy spectrum differences of the beam loss can be completely preserved, avoiding the problem of heat source distribution distortion from the source. By employing a fourth-order Runge-Kutta numerical integration algorithm, the trajectory of particles in the composite field can be calculated more accurately, effectively ensuring the accuracy of the determination that the particle trajectory intersects with the deflection plate boundary, reducing the number of missed or misjudged lost particles, and thus improving the data reliability of the generated particle loss list. This particle loss list directly records the spatial coordinates of each lost particle. When real information such as beam density and kinetic energy is transcribed into the beam source file of the Monte Carlo radiation transport program, no additional idealization assumptions are required. This directly reflects the location differences and energy distribution characteristics of actual beam loss, reducing the calculation error of the subsequent three-dimensional energy deposition matrix from the input end. By accurately tracking loss particles and generating a high-confidence loss list, a reliable data foundation can be provided for the construction of the three-dimensional energy deposition matrix, volume heat source mapping, and temperature field solution. Ultimately, this can effectively reduce the calculation error of peak temperature rise and provide accurate heat source basis for the layout of cooling channels, material selection, and structural strength optimization of electrostatic deflection plates, thereby avoiding shortened device life or operational failures caused by inaccurate heat source assessment.

[0026] In a preferred embodiment of the present invention, step 1 includes:

[0027] Step 100: Configure beam dynamics simulation parameters. Under a three-dimensional field diagram containing the central region electromagnetic field, harmonic coil magnetic field, and electrostatic deflection electric field, use a fourth-order Runge-Kutta numerical integration algorithm for multi-particle tracking. When a particle trajectory intersects the geometric boundary of the electrostatic deflection plate, it is determined to be a lost particle, and the information of the lost particle is recorded, including the spatial coordinates and kinetic energy information of each lost particle. The use of the fourth-order Runge-Kutta numerical integration algorithm for multi-particle tracking includes: integrating the equations of motion for each particle in the three-dimensional field diagram containing the central region electromagnetic field, harmonic coil magnetic field, and electrostatic deflection electric field to calculate the particle trajectory. Specifically, this includes: clarifying the complete set of basic parameters for beam dynamics simulation, and all parameters are... The parameters must be strictly matched to the actual operating scenario of the electrostatic deflection plate of the superconducting cyclotron accelerator to provide accurate input conditions for trajectory calculation. Firstly, the spatial grid parameters use a polar coordinate grid to cover the entire region of particle motion. The radial range is set to 0 to 750 mm, which completely covers the inner and outer diameters of the electrostatic deflection plate; the axial range is set to -50 to 50 mm, which covers the thickness of the deflection plate and the surrounding field. The grid cell size is set to 0.1 mm both radially and axially. This size ensures the spatial resolution of the field map data, meets the accuracy requirements of field distribution for particle trajectory calculation, and avoids distortion of field information due to excessively coarse grids. Secondly, the macroparticle parameters are set, with the total number of macroparticles set to not less than 1 × 10⁻⁶. 5A sufficient number of macroparticles is required to meet statistical accuracy requirements and avoid deviations in beam loss distribution due to insufficient sample size; the statistical weight of each macroparticle is set according to the actual beam intensity, with a value range of 1×10. 2 Up to 1×10 4 One actual particle corresponds to one macroparticle; thirdly, the initial state parameters of the particles are: the initial energy range is set from 40keV to 230MeV, which covers the beam energy range when the deflector is working, ensuring that the simulated beam energy is consistent with the actual one; the initial phase is set from 115° to 155° and is uniformly distributed within this range, which conforms to the phase characteristics of accelerator beam injection; the initial position is concentrated at the beam injection point in the central region of the accelerator, specifically at radial r=5mm and axial z=0mm, matching the injection position of the actual beam; the initial velocity direction is radially outward along the accelerator, consistent with the tangent direction of the magnetic field in the central region, which fits the initial motion trend of the beam in the accelerator.

[0028] The system loads a pre-calculated three-dimensional field map containing the electromagnetic field of the central region, the magnetic field of the harmonic coil, and the electrostatic deflection electric field. The field map data is stored in Cartesian coordinates (x, y, z), where the electric field strength E(x, y, z) is in V / m and the magnetic field strength B(x, y, z) is in T. This field map can accurately reflect the superposition and distribution of the three fields at any point in space, and can directly provide a field environment that conforms to reality for the calculation of particle motion trajectory, avoiding trajectory deviation caused by the simplification of the field model.

[0029] After configuring the beam dynamics simulation parameters, based on the above parameters and the three-dimensional field diagram, the fourth-order Runge-Kutta numerical integration algorithm is used to discretize and solve the motion equations for each macroparticle. By iteratively calculating the position and velocity of the particles step by step, the complete motion trajectory of each macroparticle is finally obtained. The specific process is as follows: When a particle moves in the composite field composed of the electromagnetic field in the central region, the magnetic field of the harmonic coil, and the electrostatic deflection electric field, the change in its motion state is directly determined by the forces acting on it. Therefore, it is necessary to first derive the differential equations describing the particle's motion through a complete force analysis. Specifically, the electrostatic force acting on the particle is analyzed first. The electrostatic force is generated only by the electrostatic deflection electric field in the composite field, and its expression is: ,in, Let $\mathbf$ be the charge of a particle. Taking a proton as an example, its charge is $\mathbf$. ; Let be the electrostatic deflection electric field intensity vector at any point in space. The magnitude and direction of this vector are determined by the design parameters of the electrostatic deflection plate. For positively charged particles, such as protons, the direction of the electrostatic force is exactly the same as the direction of the electrostatic deflection electric field intensity vector. Next, we analyze the Lorentz force acting on the particle. The Lorentz force is generated by the magnetic field in the composite field, including the superposition effect of the magnetic field in the central region and the magnetic field of the harmonic coil. Its expression is: ,in, The velocity vector of the particle has three components. Corresponding to particles in Instantaneous velocities in three directions, and satisfying (That is, the velocity component is the first derivative of the position coordinate with respect to time). Let be the composite magnetic field strength vector at any point in space. It is the vector superposition of the magnetic field in the central region and the magnetic field of the harmonic coil at that point. The direction of the Lorentz force has a special characteristic: it is always perpendicular to the plane formed by the particle velocity vector and the composite magnetic field strength vector. This directional characteristic directly affects the trajectory of the particle. After clarifying the expressions and physical meanings of the electrostatic force and the Lorentz force, we derive the particle acceleration and the total force acting on the particle in the composite field using Newton's second law. It is the vector sum of the electrostatic force and the Lorentz force, i.e. = + According to Newton's second law ,in The mass of a particle, taking the proton as an example, is 1.673 × 10⁻⁶. -27 kg, Given the particle's acceleration vector, the expression for the particle's acceleration can be further derived. = This equation shows that the particle's acceleration is determined by both the total force and its own mass, serving as a crucial bridge connecting the force and the change in its state of motion. Finally, based on the derivative relationship between acceleration, velocity, and position, a system of first-order ordinary differential equations is obtained. Since acceleration is the first derivative of velocity with respect to time, i.e. And velocity is the first derivative of the position vector with respect to time, that is... ,in The spatial position vector of the particle directly reflects its instantaneous position in three-dimensional space. Substituting the acceleration expression into the derivative relationship, the particle's equation of motion can be ultimately decomposed into two interrelated sets of first-order ordinary differential equations: the position differential equation... This equation directly reflects the change of particle position with time; the rate of change of the position vector with respect to time is equal to the particle's instantaneous velocity vector. It is the core equation describing how the particle's position moves with time; the velocity differential equation... This equation reflects the change of particle velocity with time. The rate of change of the velocity vector with respect to time, i.e., the acceleration, is determined by the particle's charge, mass, and the cross product of the electric field strength, velocity, and magnetic field at a spatial point. It is the core equation describing how particle velocity changes with time. These two sets of first-order ordinary differential equations completely describe the motion of particles in a composite field.

[0030] After deriving the first-order ordinary differential equations (position and velocity equations) of particle motion, a numerical algorithm is needed to solve them to obtain the particle trajectory. Considering the subtle changes in the particle's motion state in the combined field, a high-precision algorithm capable of capturing these subtle changes is required. Therefore, a fourth-order Runge-Kutta numerical integration algorithm is adopted. The specific implementation process is as follows: First, the time step is determined, with Δt ≤ 5 ns. This step size is chosen because the particle in the combined field is subject to the combined effects of Lorentz force and electrostatic force, causing frequent changes in its motion direction and velocity. A step size that is too small will increase the computational load, while a step size that is too large may miss subtle motion changes, leading to trajectory deviations. Δt ≤ 5 ns achieves a balance between computational efficiency and trajectory accuracy, ensuring accurate capture of the particle's motion details. Next, the particle position and velocity within each time step t to t+Δt are iteratively calculated. The entire iterative process is divided into three closely connected stages. The calculation result of the previous stage is directly used as the input of the next stage. The purpose of this stage is to calculate the fourth-order iteration coefficients required for subsequent position updates based on the particle velocity at the current moment. Each coefficient corresponds to the velocity state at different stages, where the first-order coefficients... The initial velocity at time t is calculated using the following formula: ,in This is the particle velocity vector at time t, and this coefficient reflects the trend of the particle's position change within Δt under the initial velocity; second-order coefficient. The formula is based on the half-value of the velocity at time t and the first-order velocity coefficient. ,in This is the first-order velocity coefficient, used to correct the linear assumption of the initial velocity and introduce the initial effect of velocity changes; the third-order coefficient... The formula is based on the half-value of the velocity at time t and the second-order velocity coefficient. ,in It is a second-order velocity coefficient, which further corrects for the influence of velocity changes on position and improves accuracy; fourth-order coefficient The formula for calculating the full value of the velocity at time t and the third-order velocity coefficient is as follows: ,in It is a third-order velocity coefficient, which fully incorporates the influence of velocity changes to ensure high accuracy in coefficient calculation.

[0031] By weighted summing of the fourth-order iteration coefficients, the particle position and velocity at time t+Δt are obtained, completing one time-step iteration update. The position update uses the formula... ,in It's a weighting factor, 2 times. and This formula highlights the crucial role of intermediate stage coefficients in accuracy, ensuring high precision in position calculations and accurately reflecting particle position changes within Δt. Velocity updates utilize the following formula... Consistent with the position update logic, fourth-order velocity coefficients are integrated using weighted coefficients to ensure that velocity calculations accurately match the dynamic changes in acceleration; following the above iterative process, the position of each macroparticle is calculated step-by-step. With speed The preset maximum number of time steps is 1×10. 4 This value is derived based on the longest time it takes for a particle to travel from the injection point (radial r = 5 mm, axial z = 0 mm) to the deflection plate (radial 0 to 750 mm, axial -50 to 50 mm), using the particle's maximum velocity (corresponding to a proton velocity of approximately 2.1 × 10⁻⁶ m / s² at 230 MeV kinetic energy). 7 Based on calculations (m / s), the time it takes for a particle to travel from the injection point to the farthest point of the deflection plate is approximately 3.5 × 10⁻⁶ m / s. -5 Based on a time step of Δt = 5ns, it would require 7 × 10 ns. 3 Step, set 1×10 4 The number of steps can fully cover the entire motion process of the particle from injection to impact or exit, avoiding trajectory truncation due to insufficient steps. During the time-step iteration process, trajectory integrity checkpoints need to be executed simultaneously to monitor in real time whether the particle position exceeds the spatial range of the deflection plate (radial > 750mm or axial < -50mm or axial > 50mm), or whether the number of iteration steps reaches 1×10. 4 If any condition is met, the iteration stops. At the same time, the particle trajectory is checked for abnormal fluctuations. For example, if the position change between adjacent time steps is greater than 1 mm, it exceeds the normal motion range. If it exists, it is marked as an invalid trajectory and recalculated to ensure that the final complete motion trajectory of each macroparticle truly reflects the particle's motion state in the composite field.

[0032] After calculating the complete trajectory of each macroparticle, the next step is to determine whether the trajectory intersects with the electrostatic deflection plate to identify lost particles. A ray-convex polyhedron intersection determination algorithm (separation axis theorem) is used for accurate determination, while simultaneously recording key information about the lost particles. The specific implementation process is as follows: The separation axis theorem only applies to convex polyhedra, so the electrostatic deflection plate must first be abstracted into a high-precision convex polyhedron to ensure its geometry is completely consistent with the actual deflection plate. This convex polyhedron is enclosed by six rectangular faces, and the coordinates of all vertices are determined based on the actual design parameters of the deflection plate. With the accelerator center as the origin, the radial distance corresponding to the inner edge of the deflection plate is 100mm, and the outer edge... The corresponding radial distance is 700mm. Therefore, the x and y coordinates of the radial vertices of the convex polyhedron must satisfy a radial distance from the origin between 100mm and 700mm. The axial thickness of the deflection plate is 0.1mm, the axial coordinate of the lower surface is -0.05mm, and the axial coordinate of the upper surface is 0.05mm, corresponding to the z-coordinate range of the axial vertices of the convex polyhedron. The polar angle covered by the circumferential direction of the deflection plate ranges from 30° to 60° (the polar angle is the arctangent of y and x in the vertex coordinates), corresponding to the polar angle range of the circumferential vertices of the convex polyhedron. Combining the above parameters, the specific coordinates of the 12 edges and 8 vertices of the convex polyhedron are determined through geometric coordinate conversion, with the coordinate accuracy retained to 1×10⁻⁶. -6 For example, the coordinates of one vertex are obtained by multiplying the radial distance of the inner edge by the cosine of 30° to get the x-coordinate, multiplying by the sine of 30° to get the y-coordinate, and then taking the axial coordinate of the surface to get the z-coordinate. For the other vertex, the coordinates are obtained by multiplying the radial distance of the outer edge by the cosine and sine of 60° to get the x and y coordinates respectively, and taking the axial coordinate of the upper surface to get the z-coordinate. In this way, the spatial shape of the convex polyhedron is ensured to be completely matched with the actual deflection plate. The core logic of the separating axis theorem is to determine the spatial relationship between the ray and the convex polyhedron by the overlap of the projections. If there is an axis (called the separating axis) on which the projection range of the ray on the axis does not overlap with the projection range of the convex polyhedron on the axis, it means that the ray and the polyhedron do not intersect. If the projection ranges of the two overlap on all possible separating axes, it means that the ray and the polyhedron intersect, and the intersection position can be calculated by solving the simultaneous equations.

[0033] To ensure no scenario that could lead to the separation of the ray from the convex polyhedron is overlooked, potential separation axes are divided into two categories, totaling 18. The specific determination method is as follows: The first category consists of the face normal axes of the convex polyhedron. A convex polyhedron has six rectangular faces, each with a normal direction perpendicular to itself and pointing outwards. These six normal directions serve as the six separation axes: along the positive x-axis (corresponding to the normal of the face with the largest x-coordinate); along the negative x-axis (corresponding to the normal of the face with the smallest x-coordinate); along the y-axis; and along the x-axis. The positive y-axis (corresponding to the normal of the face with the largest y-coordinate of the polyhedron); the negative y-axis (corresponding to the normal of the face with the smallest y-coordinate of the polyhedron); the positive z-axis (corresponding to the normal of the face with the largest z-coordinate of the polyhedron, which is the upper surface of the deflection plate); and the negative z-axis (corresponding to the normal of the face with the smallest z-coordinate of the polyhedron, which is the lower surface of the deflection plate). The purpose of these axes is to detect the potential separation between the ray and the faces of the polyhedron. When a ray is parallel to a face, simply looking at the coordinate range is insufficient to determine whether they intersect. By projecting the ray onto the surface normal axis, the spatial distance between the ray and the surface can be visually observed, avoiding missed detections. The second type of trajectory—the edge cross product axis—is used. A convex polyhedron has 12 edges. First, the direction vector of each edge is determined (by subtracting the coordinates of the two endpoints of each edge; for example, subtracting the endpoint coordinates from the starting point coordinates of the edge). Then, the direction vector of the particle trajectory is determined (by the particle's velocity direction at the current time step). The particle trajectory direction vector is then cross-multiplied with the direction vector of each edge (a vector operation that results in a new vector perpendicular to the plane formed by the two original vectors). The resulting 12 new vector directions serve as 12 separation axes. These axes detect the separation of the ray as it moves along the edge of the polyhedron. When the ray extends along the edge direction, its coordinate range may overlap with the polyhedron, but they do not actually intersect. The projection of the cross product axis accurately distinguishes this situation, avoiding misjudgments. Together, these two types of axes form 18 potential separation axes, completely covering all directions in which the ray may separate from the convex polyhedron.

[0034] For each of the 18 potential separation axes, projection calculations and overlap checks are performed to ensure that the separation probability of each axis is verified. The specific logic is as follows: the projection interval calculation of the convex polyhedron involves traversing the 8 vertices of the convex polyhedron. For each vertex, its coordinates are multiplied by the direction vector of the current separation axis, i.e., the vertex's x-coordinate is multiplied by the x-component of the axis direction, the y-coordinate by the y-component of the axis direction, and the z-coordinate by the z-component of the axis direction. The three results are then added together to obtain the projection value of each vertex on the current axis. The precision of the projection value is retained to 1×10. -6m; After traversing all vertices, find the minimum and maximum values ​​among these projection values. The minimum value is taken as the left endpoint of the projection interval, and the maximum value is taken as the right endpoint. This gives the projection interval of the convex polyhedron on the current axis. The ray projection interval is calculated, where the ray is determined by the particle position at the current time step (i.e., the ray's starting point) and the particle velocity direction (i.e., the ray's direction). First, calculate the projection value of the ray's starting point on the current separation axis, using the same method as the vertex projection value, multiplying the starting point coordinates with the axis direction vector, retaining a precision of 1×10. -6 Since the ray extends infinitely along the velocity direction, we only need to consider extending it to the possible range of the convex polyhedron (to avoid meaningless calculations). Therefore, the ray projection interval is set to start from the initial projection value and end at the initial projection value plus 5mm. Because the maximum projection length of the convex polyhedron is about 3mm, the 5mm length can completely cover the projection range of the polyhedron without causing calculation redundancy due to an excessively long interval. If the left endpoint of the ray projection interval is greater than the right endpoint of the polyhedron projection interval (indicating that the ray is on the right side of the polyhedron projection), or the right endpoint of the ray projection interval is less than the left endpoint of the polyhedron projection interval (indicating that the ray is on the left side of the polyhedron projection), it means that the two intervals do not overlap. The current axis is the separation axis, and it can be directly determined that the ray does not intersect the polyhedron without verifying other axes. If the projection intervals of all 18 potential separation axes overlap, it is preliminarily determined that the ray intersects the polyhedron.

[0035] After preliminary determination of intersection, a validity check of the intersection determination is still required. This involves calculating the overlap length of the two projected intervals, which is calculated by subtracting the maximum value of the left endpoint from the minimum value of the right endpoints of the two intervals. Only when the overlap length is greater than 1 × 10⁻⁶ can the validity be verified. -6Only by matching the accuracy of the coordinates (m) can minor overlap misjudgments caused by calculation errors be eliminated; secondly, a reverse check is performed to see if the number of axes not judged as separation axes is 18, to avoid judgment deviations due to the omission of a certain axis, and to ensure the reliability of the intersection judgment results; after confirming the intersection of the ray and the convex polyhedron through the separation axis theorem and validity verification, the coordinates of the intersection point need to be further calculated to eliminate false intersection points, such as intersection points that exceed the actual range of the deflection plate. Finally, the lost particles are determined and key information is recorded. The specific implementation process is as follows: first, the faces in the convex polyhedron that may intersect with the ray are screened, and only the faces that overlap with the projection interval of the ray are retained in the projection overlap judgment, so as to reduce the amount of invalid calculation. For each screened face, its Given three non-collinear vertices, denoted as vertex 1 (x1, y1, z1), vertex 2 (x2, y2, z2), and vertex 3 (x3, y3, z3), derive the plane equation based on the coordinates of these three points. The standard form of the plane equation is Ax + By + Cz + D = 0, where the coefficients are calculated as follows: A = (y2 - y1)(z3 - z1) - (z2 - z1)(y3 - y1); B = (z2 - z1)(x3 - x1) - (x2 - x1)(z3 - z1); C = (x2 - x1)(y3 - y1) - (y2 - y1)(x3 - x1); D = -(Ax1 + By1 + Cz1). All coefficients are retained to a precision of 1 × 10⁻¹⁰. -6 This ensures that the plane's position in space perfectly matches the actual deflection plate surface; the particle position at the current time step is used as the ray starting point (x0, y0, z0), and the instantaneous velocity component of the particle ( , , Assuming the ray direction is 0, a time parameter t (t≥0, since the ray only extends along the velocity direction, t<0 has no physical meaning) is introduced to construct the ray equation. This equation can completely describe the trajectory of a ray in space; by substituting the ray equation into the plane equation and solving for the time parameter t through algebraic operations, the accuracy of the t value is preserved to 1×10⁻⁶. -9Since the ray extends only along the direction t≥0, the minimum value is selected from all positive t values. This value corresponds to the first intersection point of the ray and the convex polyhedron, i.e., the actual position where the particle impacts the deflection plate. Substituting this t value into the ray equation, the spatial coordinates (x, y, z) of the intersection point can be calculated. These coordinates are the instantaneous position where the particle intersects the boundary of the deflection plate. To eliminate false intersection points caused by calculation errors, if the coordinates exceed the actual range of the deflection plate, the intersection position must simultaneously meet the following three conditions to be considered a lost particle. Specifically, the radial range verification involves calculating the radial distance from the intersection point to the accelerator center, and the result must be between 100mm and 700m. The range of the intersection point is between m and the radial vertex range of the convex polyhedron. The axial range verification requires the z-coordinate of the intersection point to be between -0.05mm and 0.05mm, consistent with the thickness range of the deflection plate (0.1mm). The circumferential range verification is performed by calculating the polar angle of the intersection point using arctan2(y,x), and the result must be between 30° and 60°, consistent with the circumferential coverage range of the deflection plate. If the intersection position exceeds any of the above ranges, such as a z-coordinate of 1mm, which is significantly beyond the thickness of the deflection plate, the particle is re-determined as a non-loss particle to ensure that particles that have not collided with the deflection plate are not misjudged as losses, and particles that have collided with the plate are not missed as non-losses.

[0036] For particles identified as lost, core information must be recorded in real time. Specific recording content and requirements are as follows: First, record the instantaneous position (x, y, z) of the particle intersecting the deflection plate boundary; second, record the kinetic energy, based on the instantaneous velocity of the particle at the point of intersection (x, y, z). , , ) Calculation, first using the formula (Joule) = Calculate the kinetic energy in Joules, where m is the particle mass, such as 1.673 × 10⁻⁶ for a proton. -27 kg, and according to 1 eV = 1.602 × 10 -19 J is converted to electron volts, with the calculation accuracy retained to 1 eV, to accurately reflect the energy level when a particle hits the deflector plate; third, the particle type is labeled according to the actual beam type of the accelerator, such as proton or carbon ion. This type needs to match the particle-material interaction model in the subsequent Monte Carlo simulation, such as the interaction parameters of proton-molybdenum material corresponding to proton, to avoid errors in energy loss calculation due to mismatch of particle types; fourth, the statistical weight records the number of actual particles corresponding to the macroparticle, and the value is consistent with the macroparticle weight set in the beam dynamics simulation (100 to 10,000 actual particles correspond to 1 macroparticle).

[0037] Before storing the information of lost particles in the temporary database, an integrity check must be performed on each item to ensure data reliability. This check includes verifying whether the spatial coordinates meet the three range conditions for determining lost particles, eliminating false intersection data; verifying whether the kinetic energy value is within the range of 40keV to 230MeV (matching the initial energy range of the particles), avoiding abnormal kinetic energy caused by velocity calculation errors, such as negative kinetic energy or kinetic energy exceeding 230MeV; verifying whether the particle type is consistent with the preset beam type, avoiding incorrect type labeling; verifying whether the statistical weight is within the range of 100 to 10000, conforming to the setting logic of macroparticles; and verifying whether there is any missing information, such as missing kinetic energy values ​​or particle type labels. After all the verification items meet the requirements, the information is stored in the temporary database.

[0038] Step 101: Based on the lost particle information, generate a particle loss list containing particle spatial coordinates and kinetic energy. Specifically, Step 100 has stored all key information of lost particles (spatial coordinates, kinetic energy, particle type, statistical weight, time step identifier) ​​in a temporary database. However, this data is stored in an unstructured format, such as individual information fields being recorded independently without unified association. It cannot be directly called by subsequent Monte Carlo radiation transport programs, such as FLUKA and MCNP. It is necessary to convert the unstructured data into a structured particle loss list according to a standardized process. The specific implementation process is as follows: First, the data format is selected. Combining the core requirements of subsequent Monte Carlo radiation transport programs for data reading efficiency, multi-field collaborative storage, and structured association, it is finally determined to use the HDF5 format to store the particle loss list. The advantage of this format is that it supports efficient compression and fast reading of multi-dimensional and multi-field data, can completely preserve the structured association between data, such as the correspondence between datasets and field attributes, and can be adapted to the data reading interface of different Monte Carlo programs through the general HDF5 library, effectively avoiding data loss or distortion caused by format incompatibility, and providing reliable technical support for cross-program data transmission.

[0039] After determining to use HDF5 format to store the particle loss list, the next step is to create a structured dataset and define fields in the HDF5 file to achieve standardized storage of loss particle information. Specifically, a dataset named "Particle Loss Table" is created in the HDF5 file. This dataset is the core structure for storing loss particle information. Each loss particle corresponds to an independent data record, and each record must contain five fixed fields. The definition, data type, and data source of all fields must strictly match those in step 100 to ensure the integrity and traceability of the information. The specific field definitions are as follows: The first field is the particle type, with a string data type. The data in this field directly uses the type of the loss particle label from step 100, such as proton labeled as Proton and carbon ion labeled as Carbon-Ion. Its core function is to ensure that particle type information is consistently transmitted throughout the entire heat source reconstruction process; the second field is spatial coordinates, with a data type of three-dimensional floating-point array. This field stores the instantaneous coordinates (x, y, z) when the particle intersects the boundary of the deflection plate, with the coordinate unit being meters. The data precision is consistent with the coordinate precision recorded in step 100; the third field is kinetic energy, with a data type of floating-point. The data directly calls the kinetic energy result recorded in step 100, eliminating the need for repeated calculations. This avoids errors introduced by repeated calculations and ensures that the stored kinetic energy data is completely consistent with the energy of the particle when it actually hits the deflection plate; the fourth field is statistical weight, with a data type of integer. This field stores the actual number of particles corresponding to the current macroparticle, and the value is consistent with the macroparticle weight set in step 100 (1×10). 2 Up to 1×10 4 (Each actual particle corresponds to one macro particle); the fifth field is the time step identifier, which is an integer. This field stores the time step number when the particle is determined to be lost. The data comes from the time step record of the particle trajectory iteration calculation in step 100.

[0040] After writing all the data for the lost particles, a comprehensive data verification of the HDF5 file is required to ensure data reliability. This involves three steps: First, verifying field completeness to ensure no missing fields in each record, such as particle type, kinetic energy value, or spatial coordinates. Second, verifying the reasonableness of the data range, such as kinetic energy being within the range of 40keV to 230MeV and spatial coordinates being within the geometric boundaries of the deflection plate (radial 100mm to 700mm, axial -0.05mm to 0.05mm, circumferential 30° to 60°). Third, verifying the correctness of the data format, ensuring that the data types of each field meet the preset requirements, such as kinetic energy as floating-point, statistical weights as integers, and spatial coordinates as three-dimensional floating-point arrays. After successful verification, the HDF5 file is named the Particle Loss Table.

[0041] This embodiment performs multi-particle tracking under a three-dimensional composite field map, which can realistically reproduce the motion trajectory of particles in the actual field environment. Combined with the high-precision calculation of the fourth-order Runge-Kutta numerical integration algorithm, it can accurately determine whether a particle intersects with the deflection plate boundary, avoiding misjudgment or omission of lost particles due to trajectory calculation deviations. The particle loss list directly records the real spatial coordinates and kinetic energy of the lost particles, avoiding the problem of heat source distribution distortion caused by idealized assumptions from the source. The lost particles are determined through rigorous trajectory calculation and boundary comparison, and the core information of the particles is recorded completely. This information is structured into a list to ensure that each entry can reflect the actual state of the lost particles. Compared with simplified methods of losing particle statistics, such as only recording the number of lost particles while ignoring position and energy differences, this process can retain the detailed characteristics of beam loss, such as spot-like spatial distribution and energy spectrum differences, making the particle loss list more reliable.

[0042] In a preferred embodiment of the present invention, step 2 includes:

[0043] Step 200: Convert the particle loss list into a beam source file recognized by the Monte Carlo radiation transport program, and set the Monte Carlo parameters, including deflection plate material parameters, boundary conditions, particle type, and energy spectrum parameters. Specifically, this involves using a self-written Python script to complete data format conversion and preprocessing. The core task is to extract key particle information and adapt it to the input requirements of the Monte Carlo program. The specific operation is as follows: the script calls the HDF5 reading library to read the dataset storing the particle loss list in the HDF5 file, extracting the spatial coordinates, velocity, kinetic energy, particle type, and statistical weight of each lost particle line by line to ensure that no key information is omitted. The Monte Carlo program defaults to centimeters as the coordinate unit, while the HDF5 file uses meters. Each coordinate value needs to be multiplied by 100 to achieve unit conversion. For example, an original x-coordinate of 0.02 meters will be converted to 2 centimeters to avoid unit mismatch leading to incident position deviation. The Monte Carlo program requires the particle momentum direction to be a unit vector (magnitude of 1). Therefore, the velocity magnitude is first calculated based on the particle's instantaneous velocity components (u, v, w), and then each velocity component is divided by the magnitude to obtain the momentum direction unit vector. The Monte Carlo program defaults to gigaelectronvolts (GV) as the kinetic energy unit, while the HDF5 file uses electron volts (eV). The kinetic energy value needs to be multiplied by 10. -9The system performs unit conversion and generates a beam source file according to the beam source card format of the Monte Carlo program. The file contains two types of core cards: one is a basic particle attribute card, which records the particle type, the converted gigaelectronvolt-level kinetic energy, and the statistical weight (the statistical weight is consistent with the macroparticle weight set in step 100), ensuring that the simulated particle attributes match the actual beam; the other is a spatial position and momentum direction card, which records the converted centimeter-level x, y, and z coordinates and the momentum direction x, y, and z components, ensuring that the particle incident position and motion direction are consistent with the actual impact state. Each lost particle corresponds to a set of the above cards. The script automatically generates a complete beam source file by iterating through all particle data, without any manual intervention, thus avoiding information transmission deviations.

[0044] In the input file of the Monte Carlo program, a complete set of simulation parameters needs to be set based on the actual design parameters and operating conditions of the electrostatic deflection plate. All parameters are strictly matched with the actual scenario. The specific implementation process is as follows: the deflection plate material parameters are determined by creating a material definition card, specifying that the deflection plate material is a molybdenum alloy with a pure molybdenum elemental composition, where the atomic number of molybdenum is 42 and the mass number is 95.95; the physical properties are set as density of 10.28 g / cm³, thermal conductivity of 138 W / m / Kelvin, and average excitation energy of 286 eV. These parameters are used to accurately calculate the relationship between particles and material atoms. The interaction cross sections, such as ionization and scattering cross sections, are used to ensure that the simulation results of particle-material interactions conform to the actual physical processes. The geometric and boundary conditions must be set to fit the spatial shape and working environment of the deflection plate. All dimensional parameters are derived from the actual design data of the deflection plate in step 100. A cuboid is used to simulate the deflection plate, and its dimensions are completely consistent with the geometric boundaries of the deflection plate defined in step 100. The radial length is the difference between the outer radius and the inner radius of the deflection plate. For example, if the inner radius is set to 100mm and the outer radius to 700mm in step 100, then the radial length = 700mm - 10. 0mm = 600mm; the axial thickness is fixed at 0.1mm, which is the axial design thickness of the deflection plate in step 100; the circumferential length is the arc length corresponding to the circumferential polar angle range of the deflection plate, calculated as arc length = polar angle radian value × average radius, where the polar angle radian value is the radian value after conversion of the circumferential polar angle range of the deflection plate in step 100 (e.g., 30° to 60°), calculated as radian value = angle value × π / 180 (e.g., 60° converted to π / 3 radians, 30° converted to π / 6 radians, the polar angle radian value is the range difference, i.e., π / 3 - π / 6 = π / 6 radians), and the average radius is... The average value of the inner and outer radii of the deflection plate is calculated as: Average radius = (Inner radius + Outer radius) / 2. For example, if the inner radius is 100mm and the outer radius is 700mm, then the average radius = (100mm + 700mm) / 2 = 400mm. The geometric coordinate system must be consistent with the coordinate system of the beam source file to avoid particle incident position deviations due to coordinate system mismatch. Regarding boundary conditions, the area inside the deflection plate is defined as a material region where particles interact with molybdenum and participate in energy deposition calculations. The area outside the deflection plate is defined as a vacuum region with a vacuum level set to 1 × 10⁻⁶. -6 The millibar (consistent with the vacuum working environment inside the accelerator) is used to measure the energy deposition of particles after they pass through the deflection plate. The particles no longer participate in the energy deposition calculation (because there is no matter in the vacuum environment that can interact with the particles). The lower surface of the deflection plate (corresponding to the minimum axial coordinate) is defined as the convective boundary, which is used to correlate the cooling conditions of computational fluid dynamics. The other surfaces are defined as vacuum boundaries, which perfectly match the actual cooling structure design of the deflection plate.

[0045] The particle types and physical processes must be set to cover all particles and mechanisms involved in energy deposition, ensuring that the calculation process closely matches actual physical laws. The particle types must be clearly defined to include incident particles, such as protons, and secondary radiation particles (photons, electrons, positrons) that may be generated during the interaction between particles and materials. This ensures that all particles potentially involved in energy deposition are included in the simulation, avoiding incomplete energy deposition calculations due to the omission of secondary particles. The physical processes must be the key processes consistent with reality. The calculation logic for each process is as follows: the ionization loss process uses Bethe's formula to calculate the ionization energy loss of charged particles in the material. The complete expression of Bethe's formula is: ,in is a constant (0.307 MeV·cm² / gram), and Z is the atomic number of the material (42). The number of particles, such as protons. =1, A is the material mass number (95.95). It is the speed of light in a vacuum. It is the ratio of particle velocity to the speed of light (calculated as actual particle velocity / speed of light). The relativistic factor (calculated as follows) ), It is the electron rest energy (0.511 megaelectron volts). Given the average excitation energy of the material (286 electron volts), substituting the above parameters into the formula, the ionization energy loss rate of the charged particle per unit distance can be directly calculated. This process accounts for over 90% of the total energy loss of charged particles. To ensure accurate calculation of charged particle energy loss, the bremsstrahlung process simulates photon radiation generated by the interaction of high-energy charged particles with kinetic energies greater than 1 MeV with the atomic nuclei of materials. When calculating photon energy, the formula is: Photon Energy = Kinetic Energy Before Interaction - Kinetic Energy After Interaction (following the law of conservation of energy). The photon direction is derived through the conservation of momentum. This is the core mechanism for the generation of secondary photons by high-energy charged particles, ensuring that the simulation of the secondary photon generation process is consistent with reality. The elastic scattering process uses the Rutherford scattering method to calculate the particle scattering angle. Specific calculations must be based on the Rutherford scattering angle formula, the expression of which is... ,in The particle scattering angle (to be calculated). is the Coulomb vacuum permittivity, with a value of 8.854 × 10⁻⁶. -12 Farads per meter (F / m) The calculation and determination of the incident particle's kinetic energy follows a two-step logic. First, in step 101, the raw kinetic energy data of the lost particles is obtained through preliminary beam loss monitoring. This data, in electron volts (eV), is directly recorded in the particle loss list dataset. Second, during the beam source file conversion process in step 200, a self-written Python script is used to perform unit conversion on the raw kinetic energy data. Since the Monte Carlo program defaults to gigaelectron volts (GeV) as the unit of kinetic energy, the raw eV-level kinetic energy is multiplied by 10. -9 (The conversion is 1 gigaelectronvolt = 1 × 10⁻⁶) 9 The energy level is calculated from electron volts to gigaelectron volts, and the converted data is ultimately stored in the beam source file as the incident particle kinetic energy in elastic scattering calculations. The value of is ensured The kinetic energy of the particles lost in the actual beam is strictly consistent; b is the collision parameter (the shortest distance between the particle trajectory and the center of the atomic nucleus), generated through random sampling, and its value range is determined by the atomic density of molybdenum material. Based on the density of molybdenum of 10.28 grams per cubic centimeter and the mass number of 95.95, the atomic number density is calculated, and then the average distance between atoms is determined as the maximum value of b (10). -10 Rice to 10 -9 (within a range of meters), the sample is randomly generated between 0 and this maximum value. The closer the distance, the smaller b is, and the stronger the corresponding Coulomb interaction. The charge number of the incident particle, such as a proton. =1; Z is the target nucleus (molybdenum atom) charge number (Z=42, consistent with the atomic number of molybdenum in the deflection plate material parameters); e is the electron charge (value is 1.602×10). -19 C), the calculation first obtains the collision parameter b through random sampling, and then combines b with the incident particle kinetic energy E and charge number. Substituting parameters such as the target nucleus charge number Z into the above formula, the scattering angle can be obtained by solving. This process determines the direction of particle motion after scattering. This process changes the particle's trajectory, which in turn affects the distribution of particle energy loss (the spatial location of energy loss changes as the particle's path changes after scattering), ensuring that the simulation of the particle's trajectory and energy loss distribution closely matches reality.

[0046] The photoelectric effect and Compton scattering process are used to simulate the interaction between photons and electrons in materials. Specifically, low-energy photons (energy less than 1 MeV) interact with electrons in the molybdenum material primarily through the photoelectric effect. During this process, the photon energy is absorbed by the electrons, which then break free from atomic bonds to become photoelectrons. The kinetic energy of the photoelectrons is calculated as the incident photon energy minus the material ionization energy. The material ionization energy is directly taken as the average excitation energy of molybdenum, 286 EV (consistent with the average excitation energy set in the deflection plate material parameters). In other words, the kinetic energy gained by the photoelectrons is the energy remaining after the incident photon energy overcomes the binding energy of the molybdenum atoms to the electrons. High-energy photons (energy greater than 1 MeV) interact with electrons in the material primarily through Compton scattering. During this process, photons and electrons undergo elastic collisions. Upon impact, the photon energy decreases and its direction of motion changes. The energy of the scattered photon is calculated using the formula: the energy of the scattered photon equals the incident photon energy divided by [1 plus the incident photon energy divided by the electron rest energy, then multiplied by (1 minus the cosine of the scattering angle)]. The electron rest energy is 0.511 MeV (consistent with the electron rest energy parameter used in the ionization loss process), and the scattering angle is the deflection angle of the photon after the collision. Its value is generated through random sampling (the sampling range is 0 to π radians, which conforms to the physical range of scattering angle values). The above two processes correspond to the main interaction mechanisms between photons of different energies and electrons in the material. Through a clear division of energy ranges and calculation methods, the energy loss path of photons in molybdenum materials is fully covered, ensuring the accuracy of photon energy loss calculation.

[0047] Step 201: Based on the beam source file and Monte Carlo parameters, perform particle transport simulation. Calculate the particle-matter interaction using the Monte Carlo method to obtain the three-dimensional energy deposition matrix on the electrostatic deflection plate element. Step 201a simulates the transport process of incident particles in the deflection plate material using the Monte Carlo method, simulating the process of secondary radiation particles generated by the interaction between the incident particles and the deflection plate material. Specifically, the incident particle transport simulation needs to simultaneously complete trajectory calculation and energy loss calculation. During trajectory calculation, the initial position and momentum direction of the incident particles are extracted from the beam source file in step 200. Combined with the cuboid geometry of the deflection plate, the motion step size is determined based on the particle kinetic energy and the molybdenum material density. The higher the energy and the lower the material density, the larger the step size. In each step, the scattering angle is calculated in real time using the Rutherford scattering angle formula (consistent with step 200), and the particle motion direction is updated according to the scattering angle, ultimately forming a complete particle motion trajectory. The energy loss calculation is based on the above trajectory, using the Beth formula to calculate the ionization energy loss rate in each step, combined with the motion step size of that step, to obtain the energy loss amount of each step (loss amount = ionization energy loss rate × motion step size of that step), and the remaining kinetic energy of the particle is updated in real time. When the remaining kinetic energy of the particle is less than 1 kiloelectron volt, or when the particle passes through the deflection plate and enters the vacuum region, the transport simulation of the particle is stopped, and its complete trajectory and energy loss data of each step are recorded.

[0048] During the transport of the incident particle, secondary radiation particles are also generated. The specific simulation is as follows: Bremsstrahlung photons are generated in scenarios where the kinetic energy of the incident particle is greater than 1 MeV. At this time, the high-energy charged particle decelerates due to a Coulomb interaction with the molybdenum nucleus, generating a photon in the process. The photon's energy is the difference between the particle's kinetic energy before and after the interaction (photon energy = kinetic energy before interaction - kinetic energy after interaction, following the law of conservation of energy). The direction of the photon's motion is derived using the momentum conservation relationship (momentum conservation formula is...). = + ,in , For particle energy, For rest mass, The speed of light in a vacuum is used to generate secondary electrons. After generation, the particle type, energy, location of origin, and direction of motion are recorded. Secondary electrons are generated through inelastic collisions between the incident particle and electrons in the molybdenum material. After gaining energy, the electrons escape atomic bonds. Their kinetic energy is the difference between the kinetic energy of the incident particle before and after the interaction, minus the ionization energy of molybdenum (secondary electron kinetic energy = (kinetic energy before interaction - kinetic energy after interaction) - 286 electron volts, following the law of conservation of energy). The direction of motion is calculated using the Mott scattering formula (the deflection angle of the secondary electron relative to the direction of the incident particle). satisfy = ,in wave number = , Let be Planck's constant. The momentum of the secondary electron. (where is the scattering angle). After generation, the particle type, energy, generation location, and direction of motion are recorded. All generated secondary radiation particles are added to the simulation queue.

[0049] Step 201b involves tracking the motion of secondary radiation particles within the electrostatic deflection plate during the simulation, statistically analyzing the energy deposition of incident and secondary radiation particles within the electrostatic deflection plate, and obtaining a spatially resolved three-dimensional energy deposition matrix within the electrostatic deflection plate element. Specifically, secondary radiation particle tracking needs to be conducted separately for photons and secondary electrons. During photon tracking, the interaction mechanism is categorized by energy. Low-energy photons (energy less than 1 MeV) primarily interact with the molybdenum material through the photoelectric effect, their energy being absorbed by electrons in the material. These energy-absorbing electrons break free from atomic bonds and become photoelectrons. The kinetic energy of the photoelectrons is the incident photon energy minus... The ionization energy of molybdenum (286 electron volts) generates photoelectrons that become new particles and are added to the tracking queue. Energy loss is calculated according to the secondary electron tracking rule. High-energy photons (energy greater than 1 megaelectron volts) primarily undergo Compton scattering. The energy of the scattered photon is calculated using the formula: Scattered photon energy = Incident photon energy / [1 + Incident photon energy / (0.511 megaelectron volts) × (1 - Cosine of scattering angle)], where 0.511 megaelectron volts is the electron rest energy, the scattering angle is the deflection angle of the photon after collision (ranging from 0 to π radians, generated through random sampling), and the direction of motion of the scattered photon is... Through the derivation of momentum conservation, the momentum conservation relationship is: incident photon momentum = scattered photon momentum + recoil electron momentum, where the photon momentum equals the photon energy divided by the speed of light in a vacuum. When the photon energy is less than 1 kiloelectron volts, or when it passes through the deflection plate and enters the vacuum region, the tracking of the photon is stopped, and its total energy loss is recorded (total energy loss = initial energy - remaining energy). Secondary electron tracking uses the Beth formula, consistent with the calculation of incident particle energy loss, to calculate the ionization energy loss rate. Combined with the step size (the step size increases with increasing electron kinetic energy and decreasing molybdenum material density), the energy loss per step is obtained (loss = ionization energy loss). The remaining kinetic energy is updated in real time (remaining kinetic energy = current kinetic energy - loss). When the remaining kinetic energy of the secondary electron is less than 1 kiloelectron volt, or when it passes through the deflection plate, tracking stops, and its total energy loss is recorded (total energy loss = initial energy - remaining energy). The generation of the three-dimensional energy deposition matrix requires three steps: voxel mesh generation, energy deposition statistics, and matrix construction. The voxel mesh generation is based on the geometric boundary of the deflection plate. The space of the deflection plate is discretized into a uniform cubic mesh with a voxel volume of 0.01 cubic millimeters. Each voxel is a cube, and its side length is the cube root of the volume. Each voxel corresponds to a unique three-dimensional index. Here, i is the radial index, calculated as (the radial coordinate of the particle minus the minimum radial coordinate of the deflector plate) divided by 0.01 mm and rounded to the nearest integer; j is the axial index, calculated as (the axial coordinate of the particle minus the minimum axial coordinate of the deflector plate) divided by 0.01 mm and rounded to the nearest integer; k is the circumferential index, calculated as (the circumferential coordinate of the particle minus the minimum circumferential coordinate of the deflector plate) divided by 0.01 mm and rounded to the nearest integer. The actual spatial coordinates of the voxels can be deduced from these indices. For example, the radial coordinate equals the minimum radial coordinate of the deflector plate plus i multiplied by 0.01 mm, and the axial and circumferential coordinates are calculated similarly. During energy deposition statistics, based on the motion data of the incident particle and all secondary radiation particles at each step, the starting and ending coordinates of the particle's motion at that step are first used to determine all voxels traversed by the trajectory. Then, the indices of each voxel traversed are matched using the above index calculation method. The energy loss of this step is accumulated into the energy value of the corresponding voxel. After traversing all motion steps of all particles, the cumulative value of each voxel is the total energy deposition of that voxel. When constructing the matrix, the total energy deposition of each voxel is arranged in the order of voxel index (first fix j and k and increment i, then fix i and k and increment j, and finally fix i and j and increment k) to form a three-dimensional energy deposition matrix. The matrix is ​​stored in binary format, in which the total energy deposition of each voxel is recorded in the form of double-precision floating-point number (occupying 8 bytes). At the same time, a mapping table between index and actual coordinate is generated, which contains the radial coordinate, axial coordinate and circumferential coordinate of the voxel center corresponding to each index.

[0050] This embodiment converts the actual particle loss list into a beam source file without adding additional idealized assumptions. It directly reflects the positional differences and energy spectrum distribution of the actual beam loss. Simultaneously, by precisely setting parameters such as the deflector plate material and boundaries, it ensures that the input conditions of the Monte Carlo simulation are consistent with the actual deflector plate, reducing simulation errors from the input end. Simulating the transport and energy loss of incident particles, it also tracks the motion process of secondary radiation particles and statistically analyzes their energy deposition, ensuring that the energy contributions of both incident and secondary radiation particles are included in the calculation. This fully recreates the energy deposition process of beam loss within the deflector plate, avoiding distortion of the energy deposition matrix due to missing energy contributions. By dividing the space into voxels and statistically analyzing the total energy deposition of each voxel, the resulting three-dimensional energy deposition matrix possesses spatial resolution, accurately reflecting the spatial distribution characteristics of energy within the deflector plate. Compared to energy deposition calculations with limited resolution, this matrix provides accurate energy distribution data for subsequent volumetric heat source mapping and temperature field solutions. This supports the optimization of the electrostatic deflector plate cooling channel layout, material selection, and structural strength, avoiding design defects in the deflector plate caused by insufficient accuracy in energy deposition calculations, and reducing the risk of shortened device lifespan or operational failures.

[0051] like Figure 2As shown, in another preferred embodiment of the present invention, step 3 includes:

[0052] Step 300: Based on the three-dimensional energy deposition matrix, obtain the discrete scattered energy deposition field. Specifically, this includes: based on the three-dimensional energy deposition matrix output from the Monte Carlo radiation transport simulation in Step 200, this matrix is ​​generated by the Monte Carlo program and records the energy deposition data of each voxel within the deflection plate. A self-written parsing script reads the complete data from the matrix. The three-dimensional energy deposition matrix uses voxels as the basic unit, and each voxel corresponds to a set of data: the coordinates of the voxel center in three-dimensional space (i.e., the position values ​​along the three directions), and the total energy value deposited by all particles within that voxel. During the parsing process, the script extracts the center coordinates and corresponding total energy deposition value of each voxel one by one, forming a discrete point set containing all voxel information. Each point consists of three-dimensional coordinates and energy deposition values. This set is the discrete scattered energy deposition field. For example, if the center coordinates of a voxel in the three-dimensional energy deposition matrix are (2.1mm, 3.5mm, 1.8mm), the total energy deposition value within that voxel is 500. If the electron volts are 1, the discrete points formed after analysis are (2.1 mm, 3.5 mm, 1.8 mm, 500 electron volts). The set of all such points is the discrete scattered energy deposition field.

[0053] Step 301 involves partitioning the discrete scattered energy deposition field using a weighted octree spatial partitioning algorithm to obtain the partitioned discrete scattered energy deposition field. Step 301a involves dividing the three-dimensional space into multiple sub-regions based on the spatial coordinates of the energy deposition points in the discrete scattered energy deposition field using a weighted octree spatial partitioning algorithm. Specifically, this includes: determining the spatial boundary of the discrete scattered energy deposition field; firstly, traversing all energy deposition points and calculating their maximum and minimum coordinate values ​​in three directions; using the cuboid formed by these two sets of coordinates as the root node of the weighted octree, which represents the initial spatial range; and then recursively partitioning the root node along the midpoints of the three directions using the weighted octree spatial partitioning algorithm. Each partition divides the current region into 8 cubic sub-regions of equal volume (called sub-nodes). The initial partitioning stops when either sub-region contains no more than 50 energy deposition points or the side length of the sub-region is not less than 0.01 mm (consistent with the voxel side length of the three-dimensional energy deposition matrix). Further partitioning of the sub-region stops when either condition is met.

[0054] Step 301b involves adjusting the fineness of the spatial division based on the distribution density of energy deposition points within each sub-region to obtain the spatial division result. Specifically, this includes: calculating the distribution density of energy deposition points for each initially divided sub-region, where the density value is equal to the number of energy deposition points contained in the sub-region divided by the volume of the sub-region (the volume of the sub-region is the cube of its side length); setting two density thresholds, with a high density threshold of 10 points per cubic millimeter and a low density threshold of 2 points per cubic millimeter. If the density of a sub-region is higher than the high density threshold, it needs to be further divided into 8 smaller sub-regions in the same way, repeating the division until its density does not exceed the high density threshold; if the density of a sub-region is lower than the low density threshold, the division of the sub-region is stopped to avoid over-dividing in areas with sparse energy deposition points and wasting computational resources; if the density of a sub-region is between the low density threshold and the high density threshold, the current division state remains unchanged. Through the above adjustments, a spatial division result that meets the density requirements is finally obtained.

[0055] Step 301c: Based on the spatial partitioning results, organize the energy deposition points in the discrete scattered energy deposition field into corresponding sub-regions to form a partitioned and organized discrete scattered energy deposition field. Specifically, this includes: traversing all energy deposition points in the discrete scattered energy deposition field, matching each point with a corresponding sub-region, and for each point's three-dimensional coordinates, determining whether its coordinate values ​​along the three directions are within the corresponding directional boundary range of a certain sub-region, i.e., the point's coordinates are greater than the minimum boundary and less than the maximum boundary of the sub-region in that direction. If all three directions satisfy this condition, then the point is assigned to that sub-region. After completing the matching of all points, each sub-region contains a set of energy deposition points belonging to its range, forming a partitioned and organized discrete scattered energy deposition field.

[0056] Step 302: Based on the discrete scattered energy deposition field after partitioning and sorting, construct a local weighted interpolation model. Specifically, this includes: obtaining the cell edge lengths of the computational fluid dynamics (CFD) grid, i.e., the edge lengths of a single cell in the grid along three directions, which are known parameters; setting the interpolation influence radius of the center of each CFD grid cell to 1.5 times the cell edge length; selecting only energy deposition points within this radius to participate in the interpolation calculation of that grid cell. For example, if the cell edge length of the CFD grid is 0.2 mm, then the interpolation influence radius is 1.5 × 0.2 mm = 0.3 mm, and only energy deposition points within 0.3 mm of the grid cell center participate in the interpolation; traversing each CFD grid cell, first calculating the three-dimensional coordinates of its center (determined by the midpoint of the cell boundary coordinates); then, according to step 301... The obtained spatial partitioning results are used to find all sub-regions whose distance from the center is less than or equal to the influence radius. All energy deposition points are extracted from these sub-regions to form the local deposition point set of the grid cell, i.e., the energy deposition points that affect the interpolation of the grid cell. For each point in the local deposition point set, its Euclidean distance to the center of the computational fluid dynamics grid cell is calculated and denoted as d. The weight of this point is calculated using a Gaussian weighting function, as shown in the formula: , where exp is the natural exponential function and σ is the weighting coefficient (taken as half of the interpolation influence radius). According to this formula, the closer the point d is, the larger its weight, and the farther the point is, the larger its weight, and the smaller its weight.

[0057] When constructing a local weighted interpolation model, the interpolated energy deposition value at the center of the hydrodynamic grid cell is calculated based on the energy deposition value of the local deposition points and their corresponding weights. That is, the sum of the products of the weights and energy deposition values ​​of all local deposition points is first calculated, and then the sum of the weights of all local deposition points is calculated. The result of dividing the two is the interpolated energy deposition value at the center of the grid cell, thereby realizing the mapping of energy deposition values ​​from discrete points to grid cells.

[0058] Step 303: Based on the local weighted interpolation model, the energy deposition values ​​in the discrete scattered energy deposition field are mapped onto the computational fluid dynamics (CFD) grid. Based on the energy deposition values ​​mapped onto the CFD grid, the unit volume heat source value for each grid cell is generated, forming a volume heat source file. Specifically, this includes: The core of calculating the volume of a CFD grid cell is to solve the volume problem according to the cell's geometric parameters and the grid type. First, the geometric parameters of the CFD grid cell are defined, including the cell boundary coordinates and side lengths. Based on these, the volume of each cell is calculated. For structured grids (cells with regular shapes, mostly cubes or cuboids, and boundary coordinates with a regular distribution), since the lengths of the cell along the three directions, i.e., the side lengths in the corresponding directions, can be directly extracted from the grid parameters, the volume is calculated as the product of the side lengths in the three directions. For example, if the side lengths of a structured grid cell along the radial, axial, and circumferential directions are all 0.2 mm, its volume is radial side length × axial side length × circumferential side length, i.e., 0.2 mm × 0.2 mm × 0.2 mm = 0.008. Cubic millimeters; due to the need to use the International System of Units (SI) for subsequent calculations, cubic millimeters need to be converted to cubic meters. The conversion relationship is 1 cubic millimeter = 10 cubic meters. -9 Cubic meters, therefore 0.008 cubic millimeters = 0.008 × 10 -9 cubic meters = 8 × 10 -12 For unstructured meshes (where the element shape is irregular, commonly tetrahedron, pentahedron, or polyhedron, and the boundary coordinates have no fixed pattern and the edge lengths cannot be directly extracted), their volume cannot be calculated by simply multiplying the edge lengths. It is necessary to directly call the built-in volume calculation module in the computational fluid dynamics solver (such as ANSYS Fluent or OpenFOAM). The working logic of this module is as follows: it automatically reads all boundary coordinates of each element, including vertex coordinates, boundary face coordinates, and other geometric information. Based on the element's topology (e.g., a tetrahedron consists of 4 vertices, and a pentahedron consists of 5 faces), it selects the corresponding geometric operation method. For tetrahedral elements, the volume is obtained by calculating 1 / 6 of the scalar triple product of the three non-coplanar edge vectors. For polyhedral elements, it first decomposes them into multiple simple sub-elements such as tetrahedrons or prisms, then sums the volumes of each sub-element, and finally automatically outputs the volume value of each element. The entire process requires no manual intervention. The unit of the output volume is consistent with the solver's mesh coordinate unit by default, such as millimeters, and can be converted to cubic meters as needed.

[0059] The calculation of the heat source value per unit volume is based on the interpolated energy deposition value of the grid cells obtained in step 302, combined with the beam running time. First, two key parameters required for the calculation are determined: one is the interpolated energy deposition value of the grid cells output in step 302, i.e., the total energy deposited by all particles within that cell; the other is the beam running time (1 second for steady-state simulation, as steady-state thermal analysis requires calculating the energy input per unit time, i.e., power). Before calculation, units must be standardized; the interpolated energy deposition value needs to be converted from electron volts to joules (the conversion relationship is 1 electron volt = 1.602 × 10⁻⁶). -19 The unit is joules, and the grid cell volume needs to be converted to cubic meters. The running time is in seconds. The calculation method is: heat source value per unit volume = interpolated energy deposition value (joules) ÷ [grid cell volume (cubic meters) × beam running time (seconds)], with the unit being watts per cubic meter.

[0060] The generation of volumetric heat source files must be done in a format recognizable by the computational fluid dynamics (CFD) solver, organizing the correspondence between mesh element indices and unit volume heat source values. First, it must be ensured that the file format matches the target solver, such as ANSYS Fluent or OpenFOAM. Common formats include HDF5 format files and user-defined function source files. If generating an HDF5 format file, a mapping dataset between mesh element indices and unit volume heat source values ​​needs to be constructed. First, a unique index is assigned to each CFD mesh element (indexes increment sequentially by element number). Then, each index is associated with its corresponding unit volume heat source value, ensuring that one index corresponds to only one element's heat source value to avoid data confusion. For example, if a CFD mesh contains 1000 elements, when generating the HDF5 file, index 0 in the dataset corresponds to element 0 with a unit volume heat source value of 1.03 × 10⁻⁶. -5 Watts per cubic meter, index 1 corresponds to the heat source value per unit volume of unit 1, such as 1.21 × 10⁻⁶. -5 Watts per cubic meter, and so on, until all 1000 elements are mapped. If a user-defined function source file is generated, a data reading function needs to be written. The function needs to associate the mesh element identification information with the heat source value per unit volume, that is, call the corresponding heat source value per unit volume through the element number (or element ID). When writing the function, it is necessary to ensure that the function syntax conforms to the solver requirements (e.g., ANSYS Fluent's UDFs must follow C language syntax) and that the data reading logic is correct to avoid calling errors. The final generated volume heat source file can be directly imported into the computational fluid dynamics solver as the heat source input data for thermal field analysis.

[0061] In this embodiment, the spatial partitioning is dynamically adjusted according to the energy deposition point density using a weighted octree. Combined with local weighted interpolation, the speckled details of high-energy deposition regions can be preserved, avoiding the distortion of heat source distribution caused by interpolation methods. This provides a heat source input that is close to reality for solving the thermal field. The spatial partitioning level is reduced in low-energy deposition regions, while the partitioning is refined in high-energy deposition regions. This ensures accuracy while avoiding unnecessary consumption of computational resources, balancing the accuracy and efficiency of heat source reconstruction. The energy deposition simulated in Monte Carlo simulation is accurately mapped to the CFD mesh, achieving efficient coupling of beam dynamics, radiation transport, and thermal fluid analysis. This provides reliable data support for optimizing the cooling scheme of electrostatic deflection plates, thermal stress analysis, and structural design.

[0062] In a preferred embodiment of the present invention, step 4 includes:

[0063] Step 400: Load the volumetric heat source file into the computational fluid dynamics calculation software. Based on the loaded volumetric heat source file, set the thermal conductivity of the deflection plate material and the boundary heat transfer conditions. Specifically, this includes: First, loading the volumetric heat source file. The loading method depends on the file format. If it is a user-defined function source file, it must first be compiled according to the syntax requirements of the computational fluid dynamics solver (e.g., user-defined functions in ANSYS Fluent must be compiled according to C language syntax to ensure that the solver can recognize the code logic). After compilation, import it through the solver's user-defined module. After import, the module automatically reads the correlation between the mesh element identification information and the heat source value per unit volume in the file, and accurately matches the heat source value corresponding to each identification to the computational mesh element in the solver. If it is HDF5... The format file can be directly imported into the solver's data import module by selecting the volumetric heat source dataset. The solver will automatically read the mapping data between the mesh element indices and the heat source values ​​per unit volume, thus associating the heat source values ​​with the corresponding mesh elements. After the volumetric heat source file is loaded, two core parameters are set based on the actual characteristics of the deflection plate: one is the thermal conductivity of the deflection plate material. In the solver's material properties module, the actual material of the deflection plate, molybdenum alloy, is selected, and its thermal conductivity is set to 138 W·m. -1 ·K -1 The first value is the standard thermal conductivity of molybdenum at room temperature, used to accurately calculate the efficiency of heat conduction within the material. The second is the boundary heat transfer condition. Considering the cooling structure and vacuum environment during the operation of the deflector plate, boundary heat transfer conditions are defined for each surface. Specifically, the back of the deflector plate, which is in contact with the cooling system, is designated as the convective heat transfer boundary, with a convective heat transfer coefficient set at 10000 W·m⁻¹. -2 · K -1 This value represents the typical convective heat transfer coefficient of a water-cooled structure, closely reflecting actual cooling efficiency. The inlet water temperature of the cooling water source is set to 20℃ (the normal temperature of a room-temperature cooling water source). The remaining surfaces of the deflector plate that do not contact the cooling system are protected because the deflector plate operates in a vacuum environment (vacuum degree 1×10⁻⁶).-6 (mbar), with no heat exchange with the outside air, so it is set as an adiabatic boundary to avoid affecting the accuracy of the calculation due to false heat loss. After the parameters are set, the spatial range of all boundaries must be checked to ensure that the boundary range completely coincides with the corresponding surface of the deflection plate geometric model to prevent errors in heat conduction calculation due to boundary misalignment.

[0064] Step 401: Solve the steady-state equation for heat conduction using the thermal conductivity of the deflection plate material and the boundary heat transfer conditions to obtain the three-dimensional temperature distribution within the electrostatic deflection plate. Specifically, this includes: first, constructing the steady-state equation for heat conduction, in the form of... ,in For gradient operators, The thermal conductivity of the deflection plate material is 138 W·m. -1 ·K -1 ), The desired temperature value is... The value of the heat source per unit volume (from the volume heat source file loaded in step 400) represents the energy balance between heat conduction and volume heat source within the material, conforming to the physical laws of steady-state thermal analysis. Next, the equation is discretized using the finite volume method to each computational fluid dynamics grid element. During discretization, the gradient operator is used... This is converted into the ratio of temperature difference to distance between mesh cell nodes, such as the gradient along the x-direction. ,in , Temperature of adjacent nodes (for node spacing); substitute and Then, the discrete equation for each unit is: ( Temperature of adjacent units (The current unit temperature), then perform iterative solution and convergence determination, setting the solution residual convergence criterion to 10. -6 That is, the calculation error of the discrete equations of each unit must be less than 10. -6 The iterative calculation is initiated. During the iteration process, the solver automatically updates the temperature value of each grid cell until the residuals of all cells meet the convergence criteria. The calculation is then considered converged, and the three-dimensional temperature distribution is output. After convergence, the temperature values ​​of all grid cells are extracted from the solver post-processing module, and the data is organized according to the three-dimensional spatial coordinates (x, y, z) of the cells to form the three-dimensional temperature distribution within the electrostatic deflection plate.

[0065] Step 402 involves feeding back the three-dimensional temperature distribution to the beam dynamics simulation parameter configuration to optimize the particle source and electrostatic deflector plate structure design. Specifically, this includes: first, performing a temperature distribution analysis, statistically analyzing the highest temperature, temperature gradient, and spatial location of the hotspot region in the three-dimensional temperature distribution, and comparing it to the safe operating temperature of the molybdenum alloy (approximately 500℃). If the highest temperature is lower than the safe temperature, it indicates that the current particle source parameters and deflector plate structure meet thermal safety requirements. If the highest temperature exceeds the safe temperature, the optimization process needs to be initiated. If optimization is required, the particle source parameters are first optimized. Based on the spatial location of the hotspot region, the source of beam loss particles in that region is determined by backtracking the particle trajectory from the beam dynamics simulation. The particle source parameters are then adjusted, such as narrowing the particle incident phase range (reducing particle impacts on the hotspot region under the corresponding phase) and optimizing the initial particle energy distribution (reducing the proportion of high-energy particles incident in the hotspot region). After adjustment, steps 100 to 401 are re-executed. The simulation calculations continue until the hotspot temperature drops below the safe temperature. If the hotspot temperature still does not meet the standard after adjusting the particle source parameters, the deflection plate structure is optimized. This involves optimizing the deflection plate structure based on the hotspot location, such as adding cooling channels on the back side corresponding to the hotspot area (to increase the convection heat transfer area and improve heat dissipation efficiency) or thinning the deflection plate thickness in the hotspot area (to reduce heat accumulation). After updating the deflection plate geometric model, the calculations from steps 200 to 401 are re-executed to verify the effect of the structural optimization. Finally, the optimization closed-loop verification is performed, and the temperature analysis, parameter adjustment, and simulation calculation process is repeated until the temperature in all areas of the three-dimensional temperature distribution meets the safety requirements, thus achieving the coupled optimization of beam dynamics and thermal analysis.

[0066] This embodiment abandons the idealized Gaussian heat source assumption and solves the heat conduction equation based on the actual beam loss volume heat source and real materials and boundary parameters. This avoids temperature distribution distortion caused by parameter simplification and provides an accurate basis for the thermal safety assessment of the deflection plate. The temperature distribution results are fed back to the beam dynamics simulation to establish a closed loop of beam loss, energy deposition, temperature distribution, and parameter optimization, breaking through the limitation of separating beam analysis and thermal analysis. It can specifically solve the thermal problems in hot spots. The temperature distribution intuitively identifies the weak hot areas of the deflection plate, providing a clear direction for cooling channel layout, material selection, and structural size adjustment, avoiding resource waste caused by blind design, and improving the working stability and lifespan of the deflection plate.

[0067] 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.

[0068] 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.

[0069] 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. A method of beam current loss heat source reconstruction for electrostatic deflector of a cyclotron, characterized in that, The method comprises: Step 1, configuring beam dynamics simulation parameters, under the three-dimensional field map containing the central zone electromagnetic field, the harmonic coil magnetic field and the electrostatic deflection electric field, performing multi-particle tracking by using the fourth-order Runge-Kutta numerical integration algorithm, when the particle trajectory intersects with the geometric boundary of the electrostatic deflection plate, determining that the particle is lost, and recording the lost particle information, including the spatial coordinates and kinetic energy information of each lost particle; based on the lost particle information, generating a particle loss list containing particle spatial coordinates and kinetic energy; Step 2, converting the particle loss list into a beam source file of a Monte Carlo radiation transport program, configuring Monte Carlo parameters, including deflection plate material, boundary, particle type and energy spectrum, performing particle transport simulation, and obtaining a three-dimensional energy deposition matrix on the electrostatic deflection plate element, including: simulating the transport process of incident particles in the deflection plate material by the Monte Carlo method, and simulating the process of generating secondary radiation particles by the interaction between incident particles and the deflection plate material; tracking the motion of secondary radiation particles in the electrostatic deflection plate during the simulation, and statistically analyzing the energy deposition of incident particles and secondary radiation particles in the electrostatic deflection plate to obtain a three-dimensional energy deposition matrix of the internal space of the electrostatic deflection plate element; Step 3, based on the three-dimensional energy deposition matrix, obtaining a discrete scattered point energy deposition field; performing partitioning and arrangement on the discrete scattered point energy deposition field by using a weighted octree space division algorithm to obtain a partitioned and arranged discrete scattered point energy deposition field; based on the partitioned and arranged discrete scattered point energy deposition field, constructing a local weighted interpolation model; based on the local weighted interpolation model, mapping the energy deposition values in the discrete scattered point energy deposition field to the computational fluid dynamics grid, and based on the energy deposition values mapped to the computational fluid dynamics grid, generating the unit volume heat source value of each grid element to form a volume heat source file; Step 4, loading the volume heat source file into a fluid mechanics calculation software to solve the temperature field, obtaining a solving result, and feeding back the solving result to optimize the particle source.

2. The method for reconstructing the beam loss heat source of a cyclotron accelerator electrostatic deflector plate according to claim 1, characterized in that, The fourth-order Runge-Kutta numerical integration algorithm is used for multi-particle tracking, which comprises: Integrating the motion equation of each particle in the three-dimensional field map containing the central zone electromagnetic field, the harmonic coil magnetic field and the electrostatic deflection electric field to calculate the particle trajectory.

3. The method of claim 2, wherein the beam current loss heat source reconstruction is performed by a computer system. The step 2 comprises: Converting the particle loss list into a beam source file recognized by the Monte Carlo radiation transport program, and setting the Monte Carlo parameters, including the deflection plate material parameters, the boundary conditions, the particle type and the energy spectrum parameters; Based on the beam source file and the Monte Carlo parameters, performing particle transport simulation, calculating the interaction between particles and matter by the Monte Carlo method, and obtaining a three-dimensional energy deposition matrix on the electrostatic deflection plate element.

4. The method of claim 3, wherein the beam current loss heat source reconstruction is performed by a computer system. The weighted octree space division algorithm is used for partitioning and arranging the discrete scattered point energy deposition field to obtain a partitioned and arranged discrete scattered point energy deposition field, which comprises: Based on the spatial coordinates of the energy deposition points in the discrete scattered point energy deposition field, the three-dimensional space is divided into multiple sub-regions by using the weighted octree space division algorithm; According to the distribution density of the energy deposition points in each sub-region, the fineness of the space division is adjusted to obtain the space division result; Based on the space division result, the energy deposition points in the discrete scatter energy deposition field are organized into corresponding sub-regions to form a partitioned discrete scatter energy deposition field.

5. The method of claim 4, wherein the beam current loss heat source reconstruction is performed by a computer system. The step 4 comprises: loading the volumetric heat source file to the fluid mechanics calculation software, setting the deflection plate material thermal conductivity and boundary heat exchange condition based on the loaded volumetric heat source file; solving the heat conduction steady-state equation through the deflection plate material thermal conductivity and boundary heat exchange condition to obtain the three-dimensional temperature distribution in the electrostatic deflection plate; feeding back the three-dimensional temperature distribution to the beam dynamics simulation parameter configuration for optimizing the particle source and electrostatic deflection plate structure design.

6. A computing device, comprising: comprise: one or more processors; a storage device for storing one or more programs, when the one or more programs are executed by the one or more processors, so that the one or more processors implement the method as claimed in any one of claims 1 to 5.

7. A computer readable storage medium characterized in that, The computer readable storage medium stores a program which, when executed by a processor, implements the method as claimed in any one of claims 1 to 5.

Images