Temperature gradient monitoring method and system for bismuth target dissolution reactions

Abstract

The application provides a temperature gradient monitoring method and system for bismuth target dissolution reaction, relates to the technical field of data monitoring, and comprises the following steps: step 1, acquiring energy deposition density distribution data and nuclide activity data along the axial and radial directions inside the bismuth target, and converting the nuclide activity data into nuclide mass data at corresponding positions; step 2, based on the nuclide mass data, calculating the intensity of a photon source term generated by a spallation reaction inside the bismuth target and corresponding energy group division results, inputting the intensity of the photon source term and the corresponding energy group division results into a photon transport calculation model to perform photon transport calculation, and obtaining a photon radiation dose equivalent rate distribution outside the target body. The application eliminates radiation interference, accurately characterizes the three-dimensional non-uniform distribution of heat sources inside the bismuth target by dividing nested three-dimensional subspace domains, constructing a three-dimensional closed vortex loop and combining a comprehensive correction factor with a photoacoustic signal temperature effect correction, and realizes real and reliable monitoring of the axial and radial temperature gradients in the bismuth target dissolution reaction.

Description

Technical Field

[0001] This invention relates to the field of data monitoring technology, and in particular to a method and system for monitoring temperature gradients in bismuth target dissolution reactions. Background Technology

[0002] In accelerator-driven subcritical systems or spallation neutron sources, bismuth targets, as crucial spallation target materials, directly influence the thermal stress state and operational lifespan of the target under high-energy proton bombardment due to their temperature gradient distribution. Accurately obtaining axial and radial temperature gradient information within the bismuth target is valuable for assessing the target's thermal load and verifying thermal-hydraulic design. Existing temperature monitoring methods for spallation targets often employ indirect inversion based on externally measurable physical quantities. For example, measuring the photoacoustic spectral signals of characteristic gases surrounding the target and combining this with thermodynamic models to calculate the internal temperature distribution. However, this method requires further exploration in practical applications. The internal energy deposition distribution of bismuth targets under high-energy proton bombardment typically exhibits significant three-dimensional non-uniformity. Furthermore, the photons released during the decay of radioactive nuclides generated by the spallation reaction create a mixed radiation field around the target. When directly using the original photoacoustic spectral signals for temperature inversion without fully considering the non-uniformity of heat source distribution within the target and the interference of the radiation field on the photoacoustic signal, the obtained temperature gradient may deviate from the true value.

[0003] For example, in a spallation target thermal verification project, researchers attempted to monitor the target temperature gradient using helium photoacoustic sensors arranged around the target. During the actual measurements, they found that the axial temperature distribution at the target's front end, calculated based on a simplified one-dimensional heat source model and uncorrected photoacoustic signals, differed somewhat from the reference data obtained through subsequent three-dimensional neutron transport simulations combined with thermal analysis in local areas. This was particularly evident in local hotspot regions approximately 5 cm axially and 2 cm radially offset from the beam inrush point, where the discrepancies were relatively significant. Subsequent analysis suggested that this phenomenon might be related to the interference of the secondary photon field derived from the nuclide mass distribution within the target on the photoacoustic signal, and the fact that the three-dimensional non-uniform distribution characteristics of the heat source were not fully represented in the inversion model. Summary of the Invention

[0004] This invention provides a method and system for monitoring the temperature gradient of a bismuth target dissolution reaction, which fully restores the true spatial distribution of heat sources inside the target.

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

[0006] In a first aspect, a method for monitoring the temperature gradient during a bismuth target dissolution reaction, the method comprising:

[0007] Step 1: Obtain energy deposition density distribution data along the axial and radial directions inside the bismuth target and activity data of each nuclide, and convert the nuclide activity data into nuclide mass data at the corresponding positions;

[0008] Step 2: Based on the nuclide mass data, calculate the intensity of the photon source term generated by the spallation reaction inside the bismuth target and the corresponding energy group partitioning results. Input the photon source term intensity and the corresponding energy group partitioning results into the photon transport calculation model to perform photon transport calculation and obtain the photon radiation dose equivalent rate distribution outside the target.

[0009] Step 3: Extract multiple closed isosurface clusters from the energy deposition density distribution data and divide the target into nested three-dimensional subspace domains; using the point with the maximum energy deposition density as the base point, arrange exploration nodes radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point.

[0010] Step 4: Based on the energy deposition density distribution data of each three-dimensional subspace domain boundary and each spatial exploration node on the three-dimensional closed vortex loop, establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain, obtain the heat source distribution characteristic value, and normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target.

[0011] Step 5: Collect the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and use the photon radiation dose equivalent rate distribution to perform temperature effect correction on the original photoacoustic spectral signal to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

[0012] Secondly, a temperature gradient monitoring system for the bismuth target dissolution reaction includes:

[0013] The acquisition module is used to acquire energy deposition density distribution data along the axial and radial directions inside the bismuth target and the activity data of each nuclide, and convert the nuclide activity data into nuclide mass data at the corresponding location;

[0014] The partitioning module is used to calculate the intensity of photon source terms and the corresponding energy group partitioning results generated by the spallation reaction inside the bismuth target based on the nuclide mass data. The photon source term intensity and the corresponding energy group partitioning results are input into the photon transport calculation model to perform photon transport calculation and obtain the photon radiation dose equivalent rate distribution outside the target.

[0015] The module is used to extract multiple closed isosurface clusters from energy deposition density distribution data and divide the target into nested three-dimensional subspace domains. Taking the point with the maximum energy deposition density as the base point, exploration nodes are arranged radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point.

[0016] The solution module is used to establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain based on the energy deposition density distribution data of each spatial exploration node on the boundary of each three-dimensional subspace domain and the three-dimensional closed vorticity loop, to obtain the heat source distribution characteristic value, and to normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target.

[0017] The inversion module is used to acquire the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and to perform temperature effect correction on the original photoacoustic spectral signal using the photon radiation dose equivalent rate distribution to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

[0018] Thirdly, a computer-readable storage medium for storing a computer program for performing the method as described in the first aspect.

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

[0020] By converting nuclide activity data into nuclide mass data, and using this as a basis to calculate photon source term intensity and energy group partitioning results, combined with photon transport calculations to obtain the photon radiation dose equivalent rate distribution outside the target, the interference of the photon radiation field generated by the spallation reaction on the photoacoustic spectral signal can be effectively eliminated, and the temperature effect correction of the original photoacoustic signal can be achieved. Based on the energy deposition density isosurface cluster, the target is divided into nested three-dimensional subspace domains, and a three-dimensional closed vortex loop is constructed to arrange the exploration nodes, which is adapted to the three-dimensional non-uniform distribution characteristics of energy deposition inside the bismuth target. This overcomes the limitation that the simplified heat source model cannot characterize the local heat source differences, and fully restores the real spatial distribution law of heat sources inside the target.

[0021] By establishing and solving the characteristic equations of heat source contribution in each three-dimensional subspace domain, and then generating a comprehensive correction factor through normalized weighted fusion, the non-uniform heat source distribution inside the target can be quantitatively characterized and corrected, reducing the temperature gradient inversion bias caused by uneven heat source distribution and improving the accuracy of temperature inversion results. By combining the corrected photoacoustic signal with the comprehensive correction factor to invert the temperature gradient, the organic coupling of external non-invasive monitoring and internal three-dimensional thermal field characteristics can be achieved, which can truly reflect the thermal field distribution state of the bismuth target under dissolution reaction and high-energy proton bombardment. Attached Figure Description

[0022] Figure 1 This is a schematic flowchart of a temperature gradient monitoring method for bismuth target dissolution reaction provided in an embodiment of the present invention.

[0023] Figure 2 This is a schematic diagram of a temperature gradient monitoring system for bismuth target dissolution reaction provided in an embodiment of the present invention. Detailed Implementation

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

[0025] like Figure 1 As shown, embodiments of the present invention propose a method for monitoring the temperature gradient in a bismuth target dissolution reaction, the method comprising the following steps:

[0026] Step 1: Obtain energy deposition density distribution data along the axial and radial directions inside the bismuth target and activity data of each nuclide, and convert the nuclide activity data into nuclide mass data at the corresponding positions;

[0027] Step 2: Based on the nuclide mass data, calculate the intensity of the photon source term generated by the spallation reaction inside the bismuth target and the corresponding energy group partitioning results. Input the photon source term intensity and the corresponding energy group partitioning results into the photon transport calculation model to perform photon transport calculation and obtain the photon radiation dose equivalent rate distribution outside the target.

[0028] Step 3: Extract multiple closed isosurface clusters from the energy deposition density distribution data and divide the target into nested three-dimensional subspace domains; using the point with the maximum energy deposition density as the base point, arrange exploration nodes radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point.

[0029] Step 4: Based on the energy deposition density distribution data of each three-dimensional subspace domain boundary and each spatial exploration node on the three-dimensional closed vortex loop, establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain, obtain the heat source distribution characteristic value, and normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target.

[0030] Step 5: Collect the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and use the photon radiation dose equivalent rate distribution to perform temperature effect correction on the original photoacoustic spectral signal to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

[0031] In this embodiment of the invention, by converting nuclide activity data into nuclide mass data, and using this as a basis to calculate the photon source term intensity and energy group partitioning results, and combining photon transport calculations to obtain the photon radiation dose equivalent rate distribution outside the target, the interference of the photon radiation field generated by the spallation reaction on the photoacoustic spectral signal can be effectively eliminated, and the temperature effect correction of the original photoacoustic signal can be achieved. By relying on the energy deposition density isosurface cluster to divide the target into nested three-dimensional subspace domains, and constructing a three-dimensional closed vortex loop to arrange the exploration nodes, it can adapt to the three-dimensional non-uniform distribution characteristics of energy deposition inside the bismuth target, overcome the limitation that the simplified heat source model cannot characterize the local heat source differences, and completely restore the real heat source spatial distribution law inside the target.

[0032] By establishing and solving the characteristic equations of heat source contribution in each three-dimensional subspace domain, and then generating a comprehensive correction factor through normalized weighted fusion, the non-uniform heat source distribution inside the target can be quantitatively characterized and corrected, reducing the temperature gradient inversion bias caused by uneven heat source distribution and improving the accuracy of temperature inversion results. By combining the corrected photoacoustic signal with the comprehensive correction factor to invert the temperature gradient, the organic coupling of external non-invasive monitoring and internal three-dimensional thermal field characteristics can be achieved, which can truly reflect the thermal field distribution state of the bismuth target under dissolution reaction and high-energy proton bombardment.

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

[0034] Step 100: Construct a three-dimensional solid model containing the geometry, material composition, and density parameters of the bismuth target. Set the energy, flux, and spatial distribution parameters of the incident proton beam. Use the Monte Carlo method to simulate the nuclear reaction and energy loss process between protons and bismuth target nuclei. Step 101: During the simulation, divide the target space into grids and record proton transport trajectories and energy deposition events within each grid cell. Obtain energy deposition density distribution data at different axial depths and radial offset distances within the target. Simultaneously, by statistically analyzing the types of nuclear reactions and product nuclides occurring within each grid cell, obtain the nuclide activity data corresponding to each grid cell. Specifically, this includes:

[0035] Based on the actual operating conditions of the bismuth target in an accelerator-driven subcritical system or spallation neutron source, and considering the overall structural form, external assembly relationship, and working environment conditions of the target, a three-dimensional solid model is established that is completely consistent with the physical properties of the actual bismuth target, according to the geometric dimensions, material composition, and bulk density parameters of the bismuth target. This model can realistically reflect the spatial morphology, internal structure, material distribution, and physical properties of the bismuth target. After the three-dimensional solid model is constructed, the relevant parameters are precisely set according to the incident conditions under actual operating conditions. The incident conditions specifically include setting the incident proton beam to be perpendicular to the incident end face of the target along the axial direction of the bismuth target, setting the high-energy incident energy of a single proton to 1 GeV, and setting the beam current intensity to 1 mA. These values ​​represent the total number of protons incident to the end face of the target per unit time. At the same time, the proton beam is set to adopt a Gaussian spatial distribution on the incident cross section, the radial coverage radius of the proton beam is set to 50 mm, and the center of the proton beam is kept coincident with the axial center of the bismuth target to ensure that the simulated incident parameters are completely matched with the actual operating conditions.

[0036] Monte Carlo numerical simulation was used to perform a full-process simulation of the interaction between protons and a bismuth target. This simulation fully reproduced the complete physical process of the spallation reaction and intranuclear cascade reaction between the incident proton and the bismuth target nucleus. It also simulated the energy loss caused by continuous collisions between the proton and the nucleus and electrons during its transport within the target, as well as the gradual deposition and distribution of energy at different spatial locations within the target. During the Monte Carlo simulation, the entire spatial region occupied by the bismuth target was discretized into a grid, decomposing the continuous and undivided target space into several uniformly sized grid cells with uniquely defined spatial coordinates. Each independent grid cell served as the basic data statistical unit. Throughout the simulation, the transport path and trajectory of each proton within each grid cell, as well as the specific location and magnitude of energy deposition within the corresponding grid cell, were continuously tracked and recorded. After all simulations were completed, the results were analyzed. The energy deposition data of the grid cells are summarized, statistically analyzed, organized, and normalized to obtain the energy deposition density distribution data at different depths along the axial direction (0mm, incident end face), 10mm, 20mm, 30mm, 40mm, 50mm, 60mm, 70mm, 80mm, 90mm, and 100mm, and at different offset distances along the radial direction (0mm, axial center), 5mm, 10mm, 15mm, 20mm, 25mm, 30mm, 35mm, 40mm, 45mm, and 50mm (radial edge) inside the bismuth target. At the same time, the types of nuclear reactions occurring in each grid cell are classified and statistically analyzed. For each grid cell, the total number of spontaneous decays of various radionuclides in that cell per unit time is counted. The total number of decays per unit time is the activity data of the corresponding nuclide in that grid cell. Through the above statistical and calculation methods, the activity data of the corresponding nuclides in all grid cells are obtained.

[0037] Step 102: Based on the conversion relationship between nuclide activity and mass, for the nuclide activity data within each grid cell, unit conversion is performed according to the atomic mass number and decay constant of the corresponding nuclide, combined with Avogadro's constant, to convert the activity data into the corresponding nuclide mass value, forming nuclide mass data with a three-dimensional spatial coordinate index. Specifically, this includes: based on the quantitative physical conversion relationship between radioactive nuclide activity and nuclide mass, standardization unit conversion and numerical calculation are performed on the various types of nuclide activity data within each grid cell. During the calculation, the atomic mass number of the corresponding nuclide, the intrinsic physical constant characterizing the nuclide's decay rate (i.e., the nuclide decay constant), and Avogadro's constant are introduced as basic calculation parameters. The Avogadro constant uses a standard physical value of 6.022 × 10⁻⁶. The formula for calculating nuclide mass is: Nuclide mass = (Nuclide activity × Atomic mass number) ÷ (Decay constant × Avogadro constant). The activity of each type of nuclide in each grid cell is calculated according to this formula. The activity data is uniformly converted into nuclide mass values ​​corresponding to the spatial location. At the same time, the calculated nuclide mass values ​​are bound and associated with the three-dimensional spatial coordinates of the corresponding grid cell. Finally, nuclide mass data with three-dimensional spatial coordinate index and corresponding to any spatial location inside the bismuth target are formed.

[0038] This embodiment constructs a realistic three-dimensional solid model and uses the Monte Carlo method to simulate the interaction between protons and the bismuth target, which can realistically reproduce the interaction law of protons inside the target. By dividing the target space into a grid and recording and statistically analyzing multi-dimensional data, the energy deposition density distribution and nuclide activity data in the three-dimensional space inside the target can be obtained, which is adapted to the three-dimensional non-uniform distribution characteristics of energy deposition inside the bismuth target. The conversion of nuclide activity to nuclide mass is completed through a standardized formula, and the data is formed in the form of a three-dimensional coordinate index, which ensures the accuracy and spatial correspondence of the basic data and reduces the interference of data errors on the temperature gradient monitoring results from the source.

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

[0040] Step 200: Based on nuclide mass data with three-dimensional spatial coordinate index, combined with a pre-constructed nuclide decay characteristic database and energy photon emission characteristic database, the total intensity of energy photons released by all nuclides in each spatial grid cell is calculated according to the nuclide mass values ​​contained in each spatial grid cell, combined with the half-life parameters of the corresponding nuclides and the emission probability parameters of each energy photon. The energy of the energy photons is then merged and statistically analyzed according to eighteen preset energy group structures to obtain the energy photon source term intensity of each spatial grid cell and the corresponding energy group division results. Specifically, this includes:

[0041] Based on the established nuclide mass data with three-dimensional spatial coordinate indexing, the pre-built nuclide decay characteristic database and energy photon emission characteristic database were retrieved. The specific construction methods of the two databases are as follows: The nuclide decay characteristic database was first constructed by consulting nuclear structure and nuclear decay databases and nuclide physics handbooks to comprehensively collect all radioactive nuclides (mainly including bismuth isotopes, thallium isotopes, lead isotopes, polonium isotopes, etc.) that may be produced after the spallation reaction and intranuclear cascade reaction of the bismuth target under high-energy proton bombardment. Nuclides with a half-life greater than 1 microsecond and a significant contribution to photon radiation were screened out, while nuclides with no actual photon emission or a negligible half-life were excluded. Then, for the screened nuclides, their standard half-life parameters were collected one by one. After multiple cross-calibrations (compared with multiple authoritative databases to correct data deviations), the data were entered into the database according to nuclide type, and a database was established. A one-to-one correspondence between nuclide names and half-life parameters was established, and a data verification module was set up to ensure the accuracy and completeness of the parameters in the database. The energy photon emission characteristic database was constructed based on the nuclides in the nuclide decay characteristic database. The authoritative data mentioned above was consulted to collect the specific energy values ​​of all energy photons released during the decay process of each nuclide (including α decay and γ transitions accompanying β decay), as well as the emission probability corresponding to each energy photon (i.e., the probability of releasing a photon of that energy each time the nuclide decays). The collected photon energies and emission probabilities were sorted, and weak photons with energies below 0.01 MeV and emission probabilities below 0.1% (whose influence on the total radiation field can be ignored) were removed. The photons were classified and entered into the database in order of nuclide type and photon energy from low to high, thus constructing a complete correspondence system of nuclide, photon energy, and emission probability. The nuclide decay characteristics database specifically stores the half-life parameters of various radionuclides. Based on the actual products of the bismuth target spallation reaction, the half-lives and corresponding decay constants of major nuclides are selected, where the decay constant = ln² ÷ half-life. The energy photon emission characteristics database stores the emission probability parameters corresponding to photons of different energies released during the decay process of each nuclide. Specifically, for the aforementioned major nuclides, bismuth-210 decay primarily releases photons with an energy of 0.465 MeV, with an emission probability of 97.9%. The decay of thallium-208 primarily releases photons with an energy of 2.615 MELV, with an emission probability of 99.7%, and photons with an energy of 0.511 MELV, with an emission probability of 18.3%; the decay of lead-207 primarily releases photons with an energy of 0.803 MELV, with an emission probability of 84.5%, and photons with an energy of 0.609 MELV, with an emission probability of 4.9%; the decay of polonium-210 releases photons with an energy of 0.803 MELV, with an emission probability of 10.2%.

[0042] For each spatial grid cell, using the mass values ​​of various nuclides within that cell as the core calculation basis, and combining the half-life parameters of the corresponding nuclides and the emission probability parameters of photons at each energy level, the total intensity of energy photons released by all radioactive nuclides within a single spatial grid cell is calculated. The total energy photon intensity = (Mass of each nuclide × Photon emission probability of that nuclide × Decay constant of that nuclide); According to the preset eighteen energy group structures, all energy photons released in each spatial grid cell are classified, grouped, and statistically analyzed according to different energy ranges. The preset eighteen energy group structures are as follows: they are divided in order of energy from low to high, covering all photon energy ranges generated by the bismuth target spallation reaction. Energy group 1, 0 ≤ energy < 0.1 MeV, is a low-energy photon range, mainly consisting of weak decay photons; Energy group 2, 0.1 ≤ energy < 0.2 MeV, is also a low-energy photon range, containing a small amount of nuclide decay products; Energy group 3, 0.2 ≤ energy < 0.3 MeV... Group 1, 0.3 ≤ energy < 0.4 MELV, is a low-energy photon range, used to assist in monitoring weak radiation; Group 4, 0.3 ≤ energy < 0.4 MELV, is a medium-low energy photon range, containing some lead isotope decay photons; Group 5, 0.4 ≤ energy < 0.5 MELV, is a medium-low energy photon range, mainly bismuth-210 decay photons; Group 6, 0.5 ≤ energy < 0.6 MELV, is a medium-low energy photon range, containing 0.511 MELV photons from thallium-208 decay; Group 7, 0.6 ≤ energy < 0.8 MELV, is a medium-energy photon range, containing 0.609 MELV photons from lead-207 decay; Group 8, 0.8 ≤ energy < 1 MELV... The first energy group, 0.0 MeGV, is in the medium-energy photon range, containing 0.803 MeGV photons from the decay of lead-207 and polonium-210; energy group 9, 1.0 ≤ energy < 1.2 MeGV, is in the medium-energy photon range, containing a small number of secondary decay photons from nuclides; energy group 10, 1.2 ≤ energy < 1.5 MeGV, is in the medium-high energy photon range, used to characterize medium-energy radiation; energy group 11, 1.5 ≤ energy < 2.0 MeGV, is in the medium-high energy photon range, containing some secondary decay photons from heavy nuclei; energy group 12, 2.0 ≤ energy < 2.5 MeGV, is in the high-energy photon range, containing a small number of high-energy photons from the decay of nuclides; energy group 13, 2.5 ≤ energy < 1.5 MeGV, is in the medium-high energy photon range. The energy range is as follows: Group 14, 3.0 ≤ energy < 4.0 MELV, is a high-energy photon range, mainly consisting of 2.615 MELV photons from thallium-208 decay; Group 15, 4.0 ≤ energy < 5.0 MELV, is a high-energy photon range, containing a small number of spallation reaction-related photons; Group 16, 5.0 ≤ energy < 8.0 MELV, is a high-energy photon range, containing high-energy spallation reaction product photons; Group 17, 8.0 ≤ energy < 10.0 MELV, is a high-energy photon range, containing extremely high-energy decay photons; Group 18, energy ≥ 10 MELV.0 MeV represents the ultra-high energy photon range, which can be ignored but its integrity is preserved. During the aggregation and statistical analysis, all energy photons within each grid cell are classified according to the energy ranges of the aforementioned eighteen energy groups. The total intensity of all photons within each energy group is then statistically analyzed to determine the photon intensity value corresponding to each energy group. Finally, the energy photon source term intensity corresponding to each spatial grid cell is obtained, which is the sum of the photon intensities of all energy groups in that cell, along with the energy group partitioning result that precisely matches the energy photon source term intensity, i.e., the photon intensity distribution of each energy group in that cell.

[0043] Step 201: The energy photon source term intensity and energy group partitioning results of each grid cell are used as spatial distribution source terms input into the photon transport calculation model. This model includes the geometric parameters of the bismuth target, the geometric parameters of the cooling layer, and the material property parameters of each structural component. By simulating the photoelectric effect, Compton scattering, and pair production processes during the interaction between energy photons and matter, the photon radiation dose equivalent rate is calculated at four discrete spatial coordinate points corresponding to the radial coordinates of 0°, 90°, 180°, and 270° on five monitoring sections located 10 mm, 20 mm, 30 mm, 40 mm, and 50 mm axially from the incident surface around the bismuth target. This forms a photon radiation dose equivalent rate distribution corresponding to the spatial coordinates, specifically including:

[0044] The energy photon source term intensity and energy group partitioning results corresponding to all spatial grid cells are completely input into the pre-constructed photon transport calculation model as spatial distribution source terms. The construction of this photon transport calculation model matches the three-dimensional solid model in step 100. It contains complete bismuth target geometric parameters, which are completely consistent with the geometric dimensions and structural morphology of the three-dimensional solid model in step 100. It also includes the geometric parameters of the outer cooling layer, including specific parameters such as cooling layer thickness, inner diameter, outer diameter, and axial length. At the same time, it includes the material property parameters of the bismuth target body (bismuth material), cooling layer (usually helium or liquid metal), and various surrounding structural components, such as protective layers and supporting structures, including material density, atomic number, photoelectric absorption cross section, Compton scattering cross section, etc., to simulate the interaction between photons and matter, and can realistically reproduce the material environment in the photon transmission process.

[0045] By constructing a photon transport computational model, a comprehensive and detailed numerical simulation of the entire transmission process of energy photons from each spatial grid cell to a preset monitoring position outside the bismuth target is performed. This fully reproduces every physical step of the interaction between the energy photon and the surrounding matter during transmission, ensuring that the simulation process is highly consistent with the actual physical phenomenon. The specific simulation process is as follows: the model uses the energy photon source term intensity and energy group division results of each spatial grid cell obtained in step 200 as the initial input data, assigns a unique tracking identifier to each energy photon, and clarifies its initial energy (corresponding to the energy range of its energy group). The initial emission position (corresponding to the three-dimensional spatial coordinates of the grid cell) and initial transmission direction are used to establish a photon tracking log, which records the entire process data of each photon from emission to absorption, exiting the simulation area, or decaying. The model tracks each photon stepwise according to a preset spatial step size of 0.1 mm. After each step, it is determined whether the photon interacts with the matter in the current transmission path. The determination is based on the material properties of the current matter, calculating the probability of photoelectric effect, Compton scattering, and pair production interaction between the photon and the matter. If the probability is greater than a preset threshold, the interaction is considered complete. If the probability of interaction is less than a preset threshold, the photons are determined to have interacted, and the specific type of interaction is determined based on the probability of interaction. If the probability is less than a preset threshold, the photons are determined not to have interacted, and the process continues to advance along the current direction for the next step until the monitoring position is reached or the simulated area is exited.

[0046] The simulation process of the specific interaction is as follows. The photoelectric effect is simulated as follows: when a photon with low energy (primarily corresponding to energy groups 1 to 8) encounters an inner-shell electron of a material atom, the photoelectric effect is triggered. The model calculates the transfer of all the photon's energy to the inner-shell electron. According to the law of conservation of energy, the kinetic energy of the photoelectron is calculated: photoelectron kinetic energy = incident photon energy - binding energy of the inner-shell electron of the material atom. The binding energy of the inner-shell electron is determined by retrieving the energy level of the corresponding inner-shell electron (e.g., K-shell or L-shell electron) from the model's built-in atomic energy level parameter library based on the atom type along the current photon transmission path. This energy level is the electron binding energy, ensuring the binding energy... The values ​​are perfectly matched to the type of material and the inner shell of the electron. Simultaneously, the model simulates the trajectory of photoelectrons after they detach from the atom. The direction of photoelectron motion forms a certain angle with the incident photon direction, which is precisely calculated using the law of conservation of momentum. First, the momentum of the incident photon is calculated: photon momentum = photon energy ÷ speed of light. Then, the momentum of the photoelectron after detaching from the atom is calculated: photoelectron momentum = photoelectron mass × photoelectron velocity. Since momentum conservation must be followed in the photoelectric effect, the momentum vector of the incident photon is equal to the sum of the momentum vector of the photoelectron and the atomic recoil momentum vector. Therefore, the atomic recoil momentum is first calculated by vector subtraction: atomic recoil momentum = incident photon momentum - photoelectron momentum. The magnitude of the atomic recoil momentum = ...

[0047] ,in This is a temporary angle, i.e., the initial angle between the incident photon and the photoelectron, initially preset to 0°; then, based on the vector equilibrium relationship of momentum conservation, the angle between the photoelectron's motion direction and the photon's incident direction is calculated using the law of cosines. The cosine value of the angle = ( + - )÷(2×incident photon momentum×photoelectron momentum), and finally the specific value of the included angle is obtained by inverse cosine calculation, with the value ranging from 0° to 90°.

[0048] During its movement, photoelectrons continuously collide with atoms in the surrounding matter. Each collision transfers some of its energy to the colliding atom, gradually losing its own energy. The specific energy loss is calculated using the collision kinetic energy transfer formula: energy loss per collision = kinetic energy of photoelectron before collision × collision energy transfer coefficient. This coefficient is determined based on the type of colliding atom and the speed of the photoelectron, ranging from 0.01 to 0.1, until the photoelectron's kinetic energy is exhausted and completely absorbed by the surrounding matter. The model synchronously records the position, energy loss, and remaining kinetic energy of each collision in real time, forming a complete record of photoelectron energy loss. After the photoelectric effect occurs, the incident photon is completely absorbed because all its energy is transferred to the inner-shell electrons. The model immediately terminates the tracking of the photon and includes all of the photon's energy in the energy deposition of the currently interacting matter.

[0049] The Compton scattering physics process is simulated as follows: when a photon with moderate energy (mainly corresponding to photons in energy groups 9 to 14) undergoes an inelastic collision with the outer electrons of a material atom, Compton scattering is triggered. The model will use the Compton scattering formula (energy of the scattered photon = ...) to calculate the scattering process. The electron rest mass energy is taken as 0.511 megaelectron volts. The scattering angle (the angle between the scattering direction and the incident direction) is used to calculate the energy and scattering direction of the scattered photon. Simultaneously, the energy transferred to the outer electrons is calculated, simulating the trajectory and energy loss process of the recoil electrons (the recoil electrons collide with surrounding atoms, gradually losing energy until they are absorbed). The scattered photon continues to be tracked along a new direction until the next interaction occurs, it reaches the monitoring position, or passes out of the simulation area. The model records the photon's energy changes, direction changes, and energy loss data throughout the process.

[0050] The electron-electron pair effect is simulated as follows: when a photon reaches a high energy level (≥1.022 MeV), primarily corresponding to photons in energy groups 13 to 18, and moves near the atomic nucleus, the electron-electron pair effect is triggered. The model simulates the process of a high-energy photon converting into a pair of electrons and positrons. According to the law of conservation of energy, the kinetic energy of the electrons and positrons is calculated: total kinetic energy = photon energy - 2 × electron rest mass energy. The kinetic energy of the electrons and positrons is distributed in a 1:1 ratio. Simultaneously, the trajectories of the electrons and positrons are simulated. During their motion, the positrons annihilate with electrons in surrounding matter, converting into two photons with an energy of 0.511 MeV. These two newly generated photons are included in the tracking range and continue to be tracked step by step according to the above rules. The electrons collide with surrounding atoms, gradually losing energy until they are absorbed. The model simultaneously records the energy conversion, energy loss, and new photon generation data during this process.

[0051] Throughout the photon tracking process, the model records and stores data in real time, including the trajectory of each photon, the type of each interaction, the amount of energy transferred, the amount of energy lost, and changes in the direction of motion. Simultaneously, it statistically analyzes the number of photons arriving at each preset monitoring position (five monitoring sections at 10mm, 20mm, 30mm, 40mm, and 50mm axial distance from the incident surface along the periphery of the bismuth target), and the number of photons and the remaining energy of each photon at discrete coordinates at 0°, 90°, 180°, and 270° on the circumference of each section. Finally, the model, combining the calculation standard for photon radiation dose equivalent rate, calculates the dose equivalent rate through numerical integration based on the photon energy and quantity arriving at each monitoring point and the photon-matter interaction characteristics. That is, dose equivalent rate = (Single photon energy × number of photons of that energy arriving × radiation weighting factor) ÷ monitoring point area ÷ time, where the radiation weighting factor is determined according to the photon energy. For low-energy photons (photon energy < 1 MeV), the weighting factor is 1; for medium-high energy photons (photon energy 1 to 20 MeV), the weighting factor is 1; and for ultra-high energy photons (photon energy > 20 MeV), the weighting factor is 0.5. The photon radiation dose equivalent rate of each monitoring point is calculated to complete the entire simulation calculation process.

[0052] The specific monitoring setup involves setting up five monitoring sections at axial distances of 10mm, 20mm, 30mm, 40mm, and 50mm from the incident surface around the bismuth target. Each monitoring section is perpendicular to the bismuth target axis, and its diameter is consistent with the maximum radial dimension of the bismuth target, ensuring comprehensive coverage of the target's outer radiation field. On the circumference of each monitoring section, four discrete spatial coordinate points are extracted, corresponding to the radial coordinates 0° (right direction), 90° (upward direction), 180° (left direction), and 270° (downward direction). These four points are evenly distributed on the circumference of the monitoring section with consistent spacing, characterizing the radiation intensity in different directions of the section. The photon radiation dose equivalent rate is calculated for each discrete point, and the calculation results for all discrete spatial coordinate points are then bound to the corresponding three-dimensional spatial coordinates (axial distance, radial distance, and circumferential angle) to form a photon radiation dose equivalent rate distribution that precisely corresponds to the spatial coordinates.

[0053] This embodiment accurately calculates the photon source term intensity and completes energy group division by combining nuclide mass data with a professional database. Then, it simulates the interaction process between photons and matter through a photon transport model. This allows the acquisition of the photon radiation dose equivalent rate distribution at different spatial locations around the bismuth target, effectively eliminating the interference of spallation reaction photon radiation on the temperature monitoring signal and improving the reliability of the external monitoring signal.

[0054] In a preferred embodiment of the present invention, step 3 includes:

[0055] Step 300a: For each preset energy deposition density gradient threshold, traverse all grid cells of the target space grid, compare the energy deposition density value at each grid cell's eight vertices with the energy deposition density gradient threshold, and mark the grid cells crossed by the isosurface. Specifically, this includes: for the five preset energy deposition density gradient thresholds, the thresholds are set based on the energy deposition density distribution range inside the bismuth target, divided according to a uniform gradient to ensure complete coverage of the energy deposition changes from the target edge to the core. Considering the actual energy deposition characteristics of the bismuth target after high-energy proton bombardment, the target edge is less affected by spallation reactions, and the minimum energy deposition value is approximately 5× J / m³, therefore the first threshold is set to 5× J / m³ (low-energy deposition zone at the target edge); the core region of the target is where spallation reactions are concentrated, with the highest energy deposition value being approximately 2.5 × 10⁻⁶. J / m³, therefore the fifth threshold is set at 2.5× J / m³ (high-energy deposition zone at the core of the target); to fully capture the energy deposition gradient change from the edge to the core, a uniform gradient division method was adopted, with adjacent threshold intervals of 0.5× J / m³, thus setting the second threshold 1× J / m³, third threshold 1.5× J / m³, fourth threshold 2× J / m³, ensuring that the five thresholds can fully and uniformly cover all energy deposition intervals within the bismuth target, traversing all grid cells within the bismuth target spatial grid in ascending (or descending) order, strictly adhering to the three-dimensional coordinate order of the grid cells during the traversal to ensure no cell is missed; for each grid cell, reading the energy deposition density values ​​corresponding to the eight vertices of that cell (the grid cell is a cubic structure, and the eight vertices are defined as the eight vertices formed by extending along the six directions of axial +, axial -, radial +, radial -, circumferential +, and circumferential - from the geometric center of the grid cell to the grid boundary, namely (axial +, radial +, circumferential +), (axial +, radial +, circumferential -), (axial +, radial -, circumferential +), (axial +, radial -, circumferential -), (axial +, radial -, circumferential -), (axial -, radial +, circumferential -). (+), (axial-, radial-, circumferential-), (axial-, radial-, circumferential-), (axial-, radial-, circumferential-), (axial-, radial-, circumferential-), these values ​​come from the energy deposition density distribution data obtained in step 101, and correspond one-to-one with the three-dimensional coordinates of the grid cell vertices. The energy deposition density value of each vertex is compared one by one with the energy deposition density gradient threshold being calculated. If some of the eight vertices of the grid cell have values ​​greater than the threshold and some have values ​​less than the threshold, it means that there is a surface in the grid cell with an energy deposition density equal to the threshold, that is, the isosurface passes through the grid cell, and the grid cell is marked as the grid cell through which the isosurface passes. If all eight vertex values ​​are greater than the threshold or all are less than the threshold, it means that the isosurface does not pass through the grid cell, and no marking is made. The screening of all grid cells is completed, and the preliminary positioning of the isosurface passing through the grid is achieved.

[0056] Step 300b: For each marked grid cell, identify edges whose energy deposition density values ​​at both endpoints are greater than and less than the energy deposition density gradient threshold, respectively, as edges traversed by the isosurface; for each identified edge, calculate the intersection coordinates of the isosurface and the corresponding edge based on the ratio of the difference between the energy deposition density values ​​at the two endpoints and the energy deposition density gradient threshold, combined with the spatial coordinates of the two endpoints. Specifically, for each marked isosurface traversing the grid cell, iterate through the twelve edges corresponding to that grid cell (each cubic grid cell has twelve edges, each edge connecting two adjacent vertices), and judge the energy deposition density values ​​at the two endpoints of each edge one by one, identifying edges whose energy deposition density values ​​at both endpoints are greater than and less than the current energy deposition density gradient threshold, respectively. These edges are the edges actually traversed by the isosurface, and are marked separately, excluding edges whose values ​​at both endpoints are both greater than or less than the threshold (the isosurface does not traverse these edges); for each identified traversing edge, first calculate the intersection coordinates of the two endpoints of the edge. The difference ratio between the endpoint energy deposition density value and the current energy deposition density gradient threshold is calculated as follows: Difference Ratio = (Energy Deposition Density Gradient Threshold - Edge Start Energy Deposition Density Value) ÷ (Energy Deposition Density Value at Edge End - Edge Start Energy Deposition Density Value). The energy deposition density values ​​at the edge start and end are derived from the energy deposition density distribution data in step 101. If the edge end energy deposition density value is less than the start value, the numerator (Energy Deposition Density Gradient Threshold - Edge Start Energy Deposition Density Value) and the denominator (Energy Deposition Density Value at Edge End - Edge Start Energy Deposition Density Value) in the formula are simultaneously reversed to ensure that the difference ratio is within the range of 0 to 1, corresponding to the intersection point on the edge. Then, combining the three-dimensional spatial coordinates of the two endpoints of the edge (start coordinates and end coordinates, both of which are three-dimensional coordinates of the grid cell vertices), the intersection point coordinates of the isosurface and the edge are calculated through linear interpolation. The intersection point coordinates are calculated as follows: Intersection Point Coordinates = Edge Start Coordinates + (Energy End Coordinates - Edge Start Coordinates) × Difference Ratio. The calculation of the three-dimensional coordinates requires separate calculations for the axial, radial, and circumferential coordinates.

[0057] Step 300c: Within each marked grid cell, multiple intersection points are connected according to their spatial position to form triangular patches, generating isosurfaces corresponding to the energy deposition density gradient threshold. This process is repeated for all five thresholds to obtain five spatially closed isosurface clusters nested from the outside in, dividing the target's internal space into five nested three-dimensional subspace domains. Specifically, within each marked grid cell (a cubic grid cell with an axial length of 1mm, a radial width of 1mm, and a circumferential arc length of 1mm), the intersection points of all isosurfaces calculated through linear interpolation with the grid edges are arranged and connected sequentially according to their spatial adjacency within the grid cell, selecting spatial positions... Each set of three adjacent intersection points is used as a building unit. Independent triangular patches are formed by three-point coplanar constraints. The planar shape of the triangular patches can fit the smooth curved shape of the isosurface in the local space to the greatest extent. Within a single grid cell, all generated triangular patches are seamlessly spliced ​​and completely combined to fully cover the entire area traversed by the isosurface within the grid cell, eliminating gaps and overlaps between patches. This generates a closed isosurface that completely corresponds to the current energy deposition density gradient threshold, ensuring that the generated isosurface can accurately reflect the spatial distribution boundary of the energy deposition density inside the bismuth target under this threshold, and completely match the actual energy deposition distribution state obtained by simulation inside the target.

[0058] Following the complete calculation and construction process described above, including intersection sorting, triangular patch construction, seamless splicing, and isosurface generation, the operation is repeated for each of the five pre-set energy deposition density gradient thresholds, generating independent closed isosurfaces corresponding to the five different thresholds. The five generated closed isosurfaces are arranged sequentially according to energy deposition density values ​​from low to high (or from high to low), with each isosurface stably enclosing the exterior of an adjacent isosurface with a higher energy deposition density, ultimately forming five spatial closed isosurface clusters that are nested layer by layer from the outside in. Through the spatial segmentation effect of this nested isosurface cluster, the complete three-dimensional space inside the bismuth target is uniformly divided from the outside in into five sequentially nested three-dimensional sub-domains. The first sub-domain (outermost layer) is spatially located between the edge of the bismuth target and the isosurfaces corresponding to the first and second thresholds, with an energy deposition density range of 5 × 10⁻⁶. J / m³≤Energy deposition density<1× J / m³; The second subdomain, spatially located between the isosurfaces corresponding to the second and third thresholds, has an energy deposition density range of 1× J / m³ ≤ Energy deposition density < 1.5× J / m³; the third subdomain, spatially located between the isosurfaces corresponding to the third and fourth thresholds, has an energy deposition density range of 1.5 × J / m³≤Energy deposition density<2× J / m³; The fourth subdomain, spatially located between the isosurfaces corresponding to the fourth and fifth thresholds, has an energy deposition density range of 2× J / m³ ≤ Energy deposition density < 2.5× J / m³; The fifth subdomain (the core layer) is spatially located within the isosurface corresponding to the fifth threshold, with an energy deposition density range of ≥2.5× J / m³, the range of energy deposition density values ​​within each three-dimensional subdomain is relatively concentrated, and the energy deposition density gradient distribution remains uniform and stable.

[0059] Step 301: Using the spatial coordinate point corresponding to the maximum energy deposition density in the energy deposition density distribution data as the base point, spatial exploration nodes are set radially within each three-dimensional subspace domain to form a spatial exploration node distribution array within each three-dimensional subspace domain. Setting spatial exploration nodes radially within each three-dimensional subspace domain includes: uniformly setting five exploration nodes radially from the inner boundary to the outer boundary of the corresponding three-dimensional subspace domain; simultaneously, setting an exploration direction every 72 degrees circumferentially, so that the five exploration nodes in each exploration direction are located at different radial distances within the corresponding subspace domain. The location specifically includes: a comprehensive analysis of the acquired energy deposition density distribution data inside the bismuth target; a grid-by-grid traversal method is used to read the average energy deposition density of each cubic grid cell (1mm axial, 1mm radial, and 1mm circumferential) within the bismuth target spatial grid, while simultaneously reading the three-dimensional spatial coordinates (axial, radial, and circumferential) of the geometric center of each grid cell; by comparing the energy deposition density values ​​of all grid cells, the three-dimensional spatial coordinate point corresponding to the point with the maximum energy deposition density value is located, and this coordinate point is used as the base point for the entire target space exploration, clarifying the specific axial, radial, and circumferential coordinate values ​​of the base point. For each independent three-dimensional subspace domain, first determine the radial distance values ​​corresponding to the inner boundary (isosurface near the base point) and outer boundary (isosurface far from the base point) of the subspace domain (extracted from the three-dimensional coordinates of the isosurface). Starting from the inner boundary of the subspace domain and ending at the outer boundary, evenly deploy five exploration nodes along the radial direction (from the inner boundary to the outer boundary, perpendicular to the target axis). The node spacing is consistent. The spacing calculation formula is: node spacing = (radial distance of the outer boundary of the subspace domain - radial distance of the inner boundary) ÷ 4. All radial distances in the formula use a uniform unit (mm). Ensure that the five nodes are evenly distributed within the radial range of the subspace domain, the radial distance between adjacent nodes is completely equal, and the radial positions of the five nodes increase sequentially from the inner boundary to the outer boundary, covering the entire radial range of the subspace domain.

[0060] In the circumferential direction, using the circumferential coordinates of the base point as a reference, a detection direction is set at every 72° interval (360°÷5=72°, ensuring that the five detection directions are evenly distributed without angular deviation), for a total of five detection directions, corresponding to 0°, 72°, 144°, 216°, and 288° of the circumference, respectively. The five detection directions are evenly distributed within the 360° circumferential range, and the angle between two adjacent detection directions remains at 72°. In each detection direction, five detection nodes are simultaneously deployed according to the above radial layout rules, ensuring that the five detection nodes in each detection direction are evenly distributed. Distributed at different radial distances within the subspace domain, each node does not overlap with the others in either the circumferential or radial directions, forming a standardized layout of five radial points and five circumferential directions. This achieves comprehensive circumferential and radial coverage within the subspace domain, eliminating any blind spots. Within each three-dimensional subspace domain, a complete and uniform array of spatial exploration nodes is formed. Each node corresponds to a specific three-dimensional spatial coordinate. The axial coordinates are consistent with the axial range of the subspace domain, the radial coordinates are arranged sequentially according to the node spacing, and the circumferential coordinates correspond to five fixed exploration directions. The coordinates of each node are recorded and archived.

[0061] Step 302: Within each three-dimensional subspace domain, taking the base point as the starting origin, connect five exploration nodes in the same exploration direction radially outwards. Connect corresponding nodes at the same radial distance in adjacent exploration directions circumferentially. Simultaneously, connect corresponding nodes at the same radial distance in the same exploration direction along the axial direction between radial planes, forming a three-dimensional closed vortex loop encircling the base point and traversing the entire three-dimensional subspace domain. Specifically, this includes: first, defining the five preset circumferential exploration directions (0°, ..., ...) within the subspace domain. The 72°, 144°, 216°, and 288° directions are uniformly distributed based on the circumferential coordinates of the base point, and the angle between adjacent directions is strictly maintained at 72° to ensure circumferential coverage without blind spots. For each exploration direction, the three-dimensional coordinates of the projection point of the base point in that direction are first calculated (the axial direction is consistent with the base point, the radial direction is the radial coordinate of the base point, and the circumferential direction is the angle of the exploration direction). Using this projection point as the starting endpoint, five pre-deployed uniformly distributed exploration nodes are sequentially connected along the radial outward direction (from the base point to the outer boundary of the subspace domain, perpendicular to the target axis). Five nodes are evenly spaced radially, with the spacing calculated using the formula: Node Spacing = (Radial Distance from Outer Boundary of Subspace Domain - Radial Distance from Inner Boundary of Subspace Domain) ÷ 4. This results in five independent radial node lines, one for each exploration direction. Each radial line completely traverses the radial range of the subspace domain, starting at the base projection point and ending at the outermost exploration node in that direction. The five nodes along the line correspond sequentially to the radial positions of the subspace domain from the inside out. By collecting the energy deposition density values ​​of each node, the numerical difference and gradient change rate between adjacent nodes are calculated. That is, the radial energy deposition gradient change rate = energy deposition density difference between two adjacent nodes ÷ node spacing. The magnitude of this gradient change rate reflects the rate of change of energy deposition from the inside out in that exploration direction. A larger gradient change rate indicates faster radial energy deposition decay, while a smaller gradient change rate indicates more gradual decay. This quantifies the radial energy deposition gradient in that direction.

[0062] Using the same radial distance as a unified benchmark, that is, for the five radial distances corresponding to the five radial nodes, each radial distance is used as a benchmark to connect the corresponding nodes at that radial distance position in two adjacent detection directions, such as 0° and 72°, 72° and 144°, one by one. When connecting, it is ensured that the radial distances of adjacent nodes are completely consistent and the axial coordinates are aligned, ultimately forming five closed circumferential node connection lines, with each radial distance corresponding to one circumferential connection line; all circumferential connection lines are parallel to the target axis, and their axial range is completely consistent with the axial range of the three-dimensional subspace domain, covering the subspace. The entire axial length of the domain; by collecting the energy deposition density values ​​of all nodes along each circumferential line, comparing the numerical differences of different circumferential nodes at the same radial distance, the standard deviation and uniformity coefficient are calculated. The circumferential uniformity coefficient = circumferential energy deposition standard deviation ÷ mean energy deposition density of nodes along the same circumferential line; the magnitude of the standard deviation and uniformity coefficient reflects the circumferential energy deposition non-uniformity at the same radial position. The larger the coefficient, the more uneven the circumferential distribution, and the smaller the coefficient, the more uniform the distribution. At the same time, by directly comparing the values ​​of each node, the peak, valley and abnormal areas of circumferential energy deposition can be identified.

[0063] Between different radial planes, i.e., spatial planes corresponding to different radial distances, nodes with the same exploration direction and radial distance are selected along a direction completely consistent with the target axis. These nodes, located at different axial positions under the same circumferential angle and radial distance, are then sequentially connected to form multiple axial node lines. Each axial line corresponds to a fixed exploration direction and a fixed radial distance; that is, each combination of exploration direction and radial distance corresponds to one axial line. All axial lines are parallel to the target's radial direction. Data is collected from each axial line... The energy deposition density values ​​at each node are analyzed to determine their upward, stable, or declining trends. The axial gradient change rate is calculated, which is calculated as: axial energy deposition gradient change rate = energy deposition density difference between two adjacent nodes ÷ axial node spacing. The sign and magnitude of this gradient change rate reflect the axial energy deposition distribution characteristics at the same radial-circumferential position. A positive gradient indicates that energy deposition increases along the axial direction, a negative gradient indicates that energy deposition decreases along the axial direction, and a zero gradient indicates that the distribution is stable. This method captures the axial energy deposition distribution characteristics, fully covers the axial range of the subspace domain, and locates the concentrated and transitional sections of axial energy deposition.

[0064] The constructed radial, circumferential, and axial nodes are seamlessly connected. The connection process follows the node coordinate matching principle, that is, the outermost endpoint of the radial connection (the outermost exploration node) is connected to the corresponding node of the circumferential connection, the upper and lower axial endpoints of the circumferential connection are connected to the two end nodes of the corresponding axial connection, and the radial endpoint of the axial connection is then connected to the middle node of the corresponding radial connection to close the loop. This ensures that the three types of connections have no gaps or overlaps, forming a complete closed loop structure. Finally, a three-dimensional closed vorticity loop is constructed that is evenly distributed around the base point and runs through the entire three-dimensional subspace domain. This loop completely covers all axial, radial, and circumferential areas of the subspace domain without any exploration blind spots.

[0065] In this embodiment, the target space is divided into nested three-dimensional subspace domains by extracting multi-threshold isosurfaces, which can adapt to the three-dimensional non-uniform distribution characteristics of energy deposition inside the bismuth target. Then, by standardizing the layout of exploration nodes and constructing a three-dimensional closed vortex loop, a comprehensive spatial exploration of the heat source distribution inside the target can be achieved, effectively solving the problem that simplified models cannot characterize the local heat source differences inside the target.

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

[0067] Step 400: Extract the energy deposition density values ​​of each spatial exploration node on the three-dimensional closed vortex loop and the energy deposition density values ​​on the boundaries of each three-dimensional subspace domain; based on the general solution form of the heat conduction equation in spherical coordinates, use the energy deposition density values ​​on the boundaries of each three-dimensional subspace domain as boundary conditions, and use the energy deposition density values ​​of each spatial exploration node on the three-dimensional closed vortex loop as internal constraint conditions to establish the heat source contribution characteristic equation for each three-dimensional subspace domain; solve the heat source contribution characteristic equation for each three-dimensional subspace domain to obtain the heat source distribution characteristic values ​​corresponding to each three-dimensional subspace domain, specifically including:

[0068] The system iterates through each spatial exploration node of the three-dimensional closed vortex loop constructed in step 302, reading and extracting the energy deposition density value of each node. This value comes from the energy deposition density distribution data inside the bismuth target obtained in step 101, forming a node energy dataset with full three-dimensional coverage in the axial, radial, and circumferential directions. Simultaneously, it extracts the energy deposition density values ​​from the boundaries of each three-dimensional subspace domain segmented in step 300c, i.e., the isosurfaces corresponding to the energy deposition density gradient threshold. These values ​​also come from the energy deposition density distribution data in step 101. Based on the axisymmetric structural characteristics of the bismuth target, a spherical coordinate system is selected as the analysis benchmark, which better fits the spatial structure of the target. The general solution form of the heat conduction equation in the spherical coordinate system is used as the mathematical model for constructing the characteristic equation of the heat source contribution. This general solution form can describe the heat conduction law of a non-uniform heat source under spherically symmetric boundary conditions. The formula is as follows: ,in For the temperature distribution within the subspace domain, , , These represent the radial, polar, and circumferential coordinates of a spherical coordinate system. The internal heat source density constrained by the three-dimensional closed vorticity loop in step 302. The thermal conductivity of the bismuth target material is fixed (a known parameter). The extracted energy deposition density values ​​at the boundaries of each three-dimensional subspace domain are substituted into the boundary condition terms of the aforementioned spherical coordinate system heat conduction equation, defining them as the external heat flow boundaries of each subspace domain, thus clarifying the heat conduction constraints at the subspace domain boundaries. Simultaneously, the energy deposition density values ​​of each exploration node on the three-dimensional closed vorticity loop in step 302 are used as the heat source intensity constraints within the subspace domain and substituted into the internal heat source terms of the equation (…). By combining the thermal conductivity parameters of the bismuth target material, a heat source contribution characteristic equation is constructed for each three-dimensional subspace domain to ensure that the equation can reflect the heat source distribution and heat conduction relationship of the subspace domain.

[0069] The continuous heat source contribution characteristic equation is transformed into a discrete problem that can be quantified and solved using the finite element method. The specific process is as follows: The heat source contribution characteristic equation for each three-dimensional subspace domain is discretized. First, based on the spatial size and energy deposition distribution complexity of the subspace domain, each three-dimensional subspace domain is uniformly divided into several small tetrahedral finite element elements. The element size is controlled between 0.1 mm and 0.5 mm to ensure that the element division both conforms to the boundary shape of the subspace domain (i.e., the isosurface boundary determined in step 300c) and captures subtle changes in internal energy deposition. Then, using the finite element discretization principle, the continuous spherical coordinate system heat conduction equation is transformed into a discrete linear equation system. The specific form of this linear equation system is as follows: ,in, The stiffness matrix is ​​a finite element matrix whose elements are determined by the material properties, geometry, and node coordinates of each finite element. Specifically, based on the thermal conductivity of the bismuth target material (a fixed known parameter), the thermal conductivity coefficient of each finite element is first calculated according to the three-dimensional coordinates of the element, i.e., element thermal conductivity coefficient = bismuth target thermal conductivity × element cross-sectional area ÷ element length, where the element cross-sectional area and element length are calculated from the three-dimensional coordinates and geometry of the finite element. Then, through the element assembly method, the thermal conductivity matrices of all finite element elements are integrated one by one to finally form the overall stiffness matrix, and the matrix dimension is consistent with the total number of finite element nodes. For each finite element node, there is a temperature unknown vector containing the temperature values ​​to be determined for all discrete nodes. Each node corresponds to a temperature unknown, and these finite element nodes correspond one-to-one with the spatial exploration nodes in the subspace domain, that is, the nodes on the three-dimensional closed vorticity loop in step 302, to ensure that the temperature unknown can accurately match the actual exploration location. The load vector is generated step by step from the extracted external boundary constraint data and internal heat source intensity constraint data. The specific transformation process is as follows: First, the boundary energy deposition density values ​​of each three-dimensional subspace domain extracted in step 300c are converted into boundary heat flux density, which is calculated as boundary energy deposition density ÷ boundary element thickness. The calculated boundary heat flux density is then substituted into the boundary term of the load vector. Second, the node energy deposition density values ​​of the three-dimensional closed vorticity loop extracted in step 302 are converted into internal heat source loads, which are calculated as node energy deposition density × corresponding finite element volume. The calculated internal heat source loads are then substituted into the internal term of the load vector. After the above transformation and integration, the final load vector is formed.

[0070] Substitute all the aforementioned external boundary constraint data (step 300c) and internal heat source intensity constraint data (step 302) into the linear equation system and perform iterative calculations using the Gauss-Seidel iterative method. In the initial stage of the iteration, the initial temperature value for each node needs to be determined. This initial value is the estimated temperature value of the boundary of each subspace domain. The calculation process is based on the boundary energy deposition density value extracted in step 300c, combined with the thermal characteristic parameters of the bismuth target material, and derived through a simple thermal balance relationship. The estimated temperature value = boundary energy deposition density ÷ bismuth target volumetric specific heat capacity + ambient reference temperature. Here, the bismuth target volumetric specific heat capacity is a fixed known parameter reflecting the bismuth target material's ability to store heat; the ambient reference temperature is taken as 25℃, i.e., the bismuth... The initial ambient temperature of the target is a fixed value; the boundary energy deposition density is the energy deposition density value on the boundary of each three-dimensional subspace domain extracted in step 300c. The estimated temperature value of all nodes on the boundary of each subspace domain is calculated one by one using this formula, and this is used as the initial node temperature value in the initial stage of the iteration. This value is substituted into the linear equation system to complete the node temperature calculation for the first iteration. The node temperature values ​​obtained in each iteration are then substituted back into the equation system to continuously correct the temperature calculation results for each node. During the iteration process, the temperature change of all finite element nodes is monitored in real time, i.e., the temperature difference of the same node obtained in two adjacent iterations, until the temperature difference of all nodes is less than or equal to the preset accuracy requirement, i.e., the temperature difference ≤ 1× The iteration stops at ℃, and the final numerical solution contains the specific temperature value of each finite element node. These nodes are evenly distributed throughout the three-dimensional subspace domain, covering all axial, radial, and circumferential regions of the subspace domain, and correspond one-to-one with the spatial exploration nodes on the three-dimensional closed vortex loop constructed in step 302. The temperature value of each node corresponds to a specific three-dimensional spatial coordinate (axial, radial, and circumferential coordinates). Through the temperature values ​​of these nodes, the temperature at any spatial location within the subspace domain can be directly read. This not only presents the specific numerical differences in temperature of each node, but also reflects the continuous temperature distribution trend within the subspace domain through the temperature changes of adjacent nodes. For example, the temperature values ​​of nodes in the core region are higher, and the node temperature values ​​gradually decrease from the core to the edge, consistent with the attenuation law of the heat source from the core to the edge. At the same time, by comparing the temperature values ​​of nodes in different circumferential and axial directions, the differences in circumferential and axial temperature distribution within the subspace domain can be intuitively reflected, and the areas and ranges of high and low temperatures can be identified, comprehensively and accurately reflecting the temperature distribution status at each location within the entire three-dimensional subspace domain.

[0071] From the converged numerical solutions, key parameters that characterize the core features of heat source distribution within the subspace domain are extracted one by one. These parameters collectively constitute the characteristic values ​​of heat source distribution. Specifically, four core indicators are quantified: First, the maximum heat source intensity, which is the maximum value of the energy deposition density of all finite element nodes within the subspace domain. The specific value of this maximum value directly reflects the strongest region of the heat source within the subspace domain. The node coordinates corresponding to the maximum value are the spatial locations of the strongest heat source region. The magnitude of the maximum value directly reflects the degree of energy concentration of the heat source in that region; the larger the value, the higher the heat source intensity and the more concentrated the energy in that region. Second, the average heat source intensity, obtained by calculating the arithmetic mean of the energy deposition density of all nodes within the subspace domain. This average value directly reflects the overall intensity level of the heat source within the subspace domain. The higher the average value, the higher the overall energy level of the heat source within the entire subspace domain; a lower average value indicates a lower overall heat source intensity within the subspace domain. The first is relatively weak, which can intuitively reflect the overall distribution level of heat sources in the subspace domain; the second is the energy distribution variance, which quantifies the uniformity of heat source distribution within the subspace domain. The magnitude of the variance directly reflects the uniformity. The larger the variance, the greater the difference in energy deposition density among nodes, the more uneven the heat source distribution, and the existence of obvious energy concentration or low-energy areas; the smaller the variance, the smaller the difference in energy deposition density among nodes, and the more uniform the heat source distribution; the third is the thermal field contribution ratio, which is the ratio of the total heat source energy in the subspace domain (the sum of the sum of the energy deposition density of all nodes and the product of the unit volume) to the total heat source energy inside the bismuth target. This ratio directly reflects the degree of contribution of the heat source in the subspace domain to the overall thermal field of the target. The larger the ratio, the higher the proportion of the heat source energy in the subspace domain in the total heat source energy of the target, and the greater the impact on the overall thermal field of the target; the smaller the ratio, the smaller its contribution to the overall thermal field, which can clearly define the role weight of each subspace domain in the overall thermal field.

[0072] Based on the pre-defined weighting of each indicator—maximum intensity (0.3), average intensity (0.2), energy distribution variance (0.3), and thermal field contribution (0.2)—the four core indicators are weighted and integrated. The maximum intensity of the heat source directly determines the location and energy level of the core heat source within the subdomain, making it a crucial factor influencing the thermal distribution within the subdomain and thus highly important for characterizing heat source features; therefore, it is weighted at 0.3. The average intensity reflects the overall level of heat sources within the subdomain and serves as a supplementary indicator to the core intensity, reflecting the overall distribution; its importance is slightly lower than the maximum intensity and distribution variance, hence the weighting at 0.2. Variance directly quantifies the average distribution of heat sources. Uniformity, while the non-uniformity of heat sources is a key prerequisite for subsequent heat conduction analysis, and is as important as the maximum intensity. Therefore, its weight is consistent with the maximum intensity, so it is set to 0.3. Since it mainly reflects the degree of correlation between the subspace domain and the overall thermal field of the target, it is an indicator that connects the local and the whole. Its importance is lower than the core intensity and distribution uniformity, so it is set to the lowest weight of 0.2. The characteristic value of heat source distribution = (maximum intensity × 0.3) + (average intensity × 0.2) + (energy distribution variance × 0.3) + (thermal field contribution ratio × 0.2). The final characteristic value of heat source distribution quantifies the intensity, distribution pattern and contribution weight of heat sources in the subspace domain to the overall thermal field.

[0073] Step 401: Arrange the heat source distribution feature values ​​corresponding to each three-dimensional subspace domain according to the spatial topological order from the outside to the inside or from the inside to the outside. Perform normalized weighted fusion calculation with the volume proportion of each three-dimensional subspace domain as the weight to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target. Specifically, this includes: dividing each three-dimensional subspace domain, as formed in step 300c, into five corresponding heat source distribution feature values, arranging them according to the spatial topological order from the outside to the inside of the subspace domain to form an ordered feature value sequence. The arrangement is completely consistent with the spatial nesting structure of the heat sources inside the target, that is, the subspace domains corresponding to the energy deposition density gradient threshold from low to high are arranged sequentially from the outside to the inside to ensure the sorting of feature values. The order matches the distribution pattern of strong core heat sources and weak periphery heat sources inside the target. Then, weight normalization is performed, that is, the volume ratio of each three-dimensional subspace domain is calculated. This ratio is the initial analysis weight of each heat source distribution characteristic value. The volume ratio of the subspace domain = volume of a single three-dimensional subspace domain ÷ total volume of the bismuth target. The volume of a single three-dimensional subspace domain is calculated by the three-dimensional coordinates of its boundary isosurface, and the total volume of the bismuth target is calculated based on the actual geometric dimensions (axial length, radial radius) of the target. The volume ratio of all subspace domains is normalized. The normalized volume weight = volume ratio of a single subspace domain ÷ sum of the volume ratios of all subspace domains. This calculation ensures that the sum of the normalized weights of all subspace domains is 1.

[0074] Based on the normalized volume weights, the weighted summation and fusion calculations of the distribution characteristic values ​​of each sorted heat source are performed to generate a comprehensive correction factor. This factor is a key parameter that can comprehensively characterize the three-dimensional non-uniform distribution state of heat sources inside the target. Its value integrates the heat source intensity, distribution uniformity, and thermal field contribution weights of the five sub-space domains. The weighted fusion calculation formula is as follows: ,in As a comprehensive correction factor, For the first The characteristic value of heat source distribution in each subspace domain For the first The normalized volume weights of each subdomain are used to weight and fuse the characteristics of heat source intensity, distribution uniformity, and thermal field contribution of each subdomain with their corresponding volume weights. This not only preserves the independent heat source distribution characteristics of each subdomain but also takes into account the actual space ratio of each subdomain within the bismuth target. Furthermore, it conforms to the nested spatial structure of the target from the outside to the inside, and can quantify the distribution differences of heat sources in the axial, radial, and circumferential dimensions of the target from an overall perspective. This reflects the true distribution law of highly concentrated heat sources in the core area of ​​the target and gradually attenuating heat sources in the peripheral area, ensuring that the comprehensive correction factor can truly reflect the overall non-uniform distribution state of heat sources inside the target.

[0075] In this embodiment, by establishing a heat source contribution characteristic equation constrained by a three-dimensional vortex loop and combining it with the volume weight of the subspace domain, the generated comprehensive correction factor can fully reflect the three-dimensional non-uniform distribution characteristics of the heat source inside the target, overcoming the limitations of the uniform heat source assumption. The heat source distribution characteristic value obtained by solving and the final comprehensive correction factor are the key prerequisites for realizing high-precision simulation of the temperature field inside the bismuth target. Decomposing the complex three-dimensional heat source problem into the characteristic analysis of each subspace domain and then fusing them ensures the local accuracy of the analysis.

[0076] In a preferred embodiment of the present invention, step 5 includes:

[0077] Step 500: Extract the photon radiation dose equivalent rate values ​​corresponding to twenty discrete spatial coordinate points from the photon radiation dose equivalent rate distribution. Using the photon radiation dose equivalent rate values ​​at the twenty discrete spatial coordinate points, perform temperature effect correction processing on the original photoacoustic signal amplitude at the same discrete spatial coordinate point to obtain the corrected photoacoustic signal amplitude at each discrete spatial coordinate point. Specifically, this includes:

[0078] Using the photon radiation dose equivalent rate distribution results obtained in step 102, twenty discrete spatial coordinate points are located. These coordinate points are all located on the bismuth target surface, covering key axial and radial positions of the target, and correspond one-to-one with the spatial exploration node coordinates in step 302. The specific positions of the twenty discrete spatial coordinate points are as follows: with the center of the bismuth target as the origin, the circumferential direction is fixed at 0° (aligned with the circumferential exploration direction reference in step 302), five key positions in the core working area are selected along the axial direction (Z-axis), namely 2mm, 6mm, 10mm, 14mm, and 18mm, and four core surface regions are selected along the radial direction (R-axis). The key locations are 1mm, 2mm, 3mm, and 4mm. The five axial key locations and four radial key locations are interleaved to form 20 discrete spatial coordinate points. This ensures comprehensive coverage of the key heat source detection locations on the target surface, while also ensuring that each coordinate point is evenly distributed to reflect the radiation dose and temperature distribution characteristics of different areas on the target surface. The photon radiation dose equivalent rate value corresponding to each discrete coordinate point is extracted one by one, and the three-dimensional coordinate information (axial, radial, and circumferential) of each coordinate point is recorded to establish a one-to-one correspondence between discrete spatial coordinate points and photon radiation dose equivalent rates.

[0079] Using the extracted photon radiation dose equivalent rate values ​​at each discrete spatial coordinate point, temperature effect correction was applied to the pre-collected original photoacoustic signal amplitudes (collected on the target surface, corresponding one-to-one with the discrete coordinate points) at the same coordinate points. First, a temperature correlation model between the photon radiation dose equivalent rate and the photoacoustic signal amplitude was constructed. This involved building an experimental platform consistent with the actual working environment of the bismuth target, controlling the ambient temperature at 25℃, and ensuring the experimental environment was free from other interfering factors. Different gradients of photon radiation dose equivalent rates were set (covering the dose range that might occur in actual work; specifically, the values ​​were 0.1 mSv / h, 0.3 mSv / h, 0.5 mSv / h, 0.7 mSv / h, 0.9 mSv / h, and 1.1 mSv / h, a total of 6 gradients with a gradient interval of 0.2 mSv / h, covering both low-dose and high-dose interference scenarios). At each dose gradient, corresponding values ​​were collected. The photoacoustic signal amplitude at a location is recorded, along with the actual temperature at that location, obtained using a high-precision thermometer. Using the photon radiation dose equivalent rate as the independent variable and the photoacoustic signal amplitude interference as the dependent variable, a temperature correlation model is obtained by linear fitting based on the collected experimental data. The specific formula for this model is: Photoacoustic signal amplitude interference = Fitting coefficient × Photon radiation dose equivalent rate + Fitting intercept. The fitting coefficient and fitting intercept are both obtained through linear fitting of the aforementioned experimental data and are fixed constants. The fitting coefficient quantifies the change in photoacoustic signal amplitude interference when the photon radiation dose equivalent rate changes by 1 mSv / h, and the fitting intercept corrects for minor system deviations during the experiment. Using this linear formula, any photon radiation dose equivalent rate value can be directly substituted to calculate the corresponding photoacoustic signal amplitude interference, thereby quantifying the impact of radiation interference on the photoacoustic signal amplitude under different radiation doses.

[0080] Next, a correction calculation is performed: Corrected photoacoustic signal amplitude = Original photoacoustic signal amplitude - Photoacoustic signal amplitude interference. During the calculation, the photon radiation dose equivalent rate at each discrete coordinate point is substituted into the formula to calculate the photoacoustic signal amplitude interference corresponding to that coordinate point. Then, the original photoacoustic signal amplitude at that coordinate point is subtracted from the corresponding interference to complete the temperature effect correction for a single coordinate point. This process is iterated to complete the correction of all twenty discrete spatial coordinate points, ultimately obtaining the corrected photoacoustic signal amplitude at each discrete spatial coordinate point. This correction process eliminates the interference of photon radiation dose on the photoacoustic signal amplitude. Before correction, the original photoacoustic signal amplitude is affected by both temperature and radiation dose, resulting in large signal fluctuations and an inability to accurately reflect temperature changes. After correction, the signal interference caused by radiation dose is completely eliminated, and the change in the corrected photoacoustic signal amplitude is only related to the bismuth target surface temperature. That is, as the temperature increases, the corrected photoacoustic signal amplitude increases accordingly, and as the temperature decreases, the corrected photoacoustic signal amplitude decreases accordingly. The trend of signal value change is completely consistent with the trend of temperature change, ensuring that the corrected photoacoustic signal amplitude accurately reflects the temperature-related physical characteristics of the bismuth target surface.

[0081] Step 501: Based on the pre-established correspondence between the amplitude of the characteristic gas photoacoustic spectral signal and temperature, the calibrated photoacoustic signal amplitude at each discrete spatial coordinate point is converted into the surface temperature value at the corresponding discrete spatial coordinate point. Specifically, this includes: first, calling the pre-established database of the correspondence between the amplitude of the characteristic gas (helium) photoacoustic spectral signal and temperature. The construction process of this database is as follows: a constant temperature test platform is built, and the test environment temperature is controlled within the actual working temperature range of the bismuth target (25℃ to 150℃). Different constant temperature conditions are set in 1℃ increments. Under each constant temperature condition, the standard amplitude of the characteristic gas photoacoustic spectral signal is collected. Each temperature point is collected 10 times and the average value is taken to eliminate random errors. The temperature values ​​and standard amplitudes are stored one-to-one to form a basic database. The database is calibrated for accuracy, and the actual value of each temperature point is verified by a high-precision thermometer. Finally, a helium photoacoustic spectral signal amplitude-temperature correspondence database covering the entire working temperature range of the bismuth target with accurate mapping is formed. The amplitude of the corrected photoacoustic signal at each discrete spatial coordinate point obtained in step 500 is used as the input parameter and substituted into the database for matching. If the amplitude of the corrected photoacoustic signal is completely consistent with the standard amplitude in the database, the corresponding standard temperature value is directly read. If the amplitude of the corrected photoacoustic signal is between two standard amplitudes, linear interpolation is used for calculation. The interpolation calculation formula is: surface temperature value = low standard temperature + (corrected photoacoustic signal amplitude - low standard amplitude) × (high standard temperature - low standard temperature) ÷ (high standard amplitude - low standard amplitude). Wherein, the low standard temperature is the standard temperature value in the database that is closest to and slightly smaller than the amplitude of the corrected photoacoustic signal. For example, when the corrected amplitude is between the standard amplitudes corresponding to 25℃ and 26℃, the low standard temperature is 25℃, and the high standard temperature is the adjacent larger standard temperature value. Through this calculation, the amplitude of the corrected photoacoustic signal at each discrete coordinate point is converted into the surface temperature value at the corresponding coordinate point. This process realizes the conversion of the photoacoustic signal from a physical quantity (amplitude) to a temperature quantity (degrees Celsius) and records the surface temperature value at each coordinate point.

[0082] Step 502: Using the surface temperature value as a boundary condition, the energy deposition density distribution data is combined with a comprehensive correction factor to obtain the corrected internal heat source distribution. Based on the corrected internal heat source distribution and the surface temperature boundary condition, the heat conduction equation inside the bismuth target is solved to obtain the three-dimensional temperature distribution field along the axial and radial directions inside the bismuth target. Specifically, this includes: firstly, using the surface temperature values ​​at each discrete spatial coordinate point obtained in step 501 as the first type of boundary condition (i.e., given boundary temperature) for the heat conduction equation to clarify the temperature constraints at each key location on the bismuth target surface; simultaneously, combining the comprehensive correction factor generated in step 401 with the original energy deposition density distribution data obtained in step 101, and correcting the original energy deposition density using the comprehensive correction factor. The corrected energy deposition density equals the original energy deposition density. The degree × comprehensive correction factor is an overall correction coefficient obtained by weighting the heat source distribution characteristics (intensity, uniformity, and thermal field contribution) of each sub-space domain of the target with volume weights. Its value directly reflects the degree of deviation between the original energy deposition density and the actual thermal effect of the target. When the comprehensive correction factor > 1, it indicates that the original energy deposition density underestimates the heat source intensity, and the value increases after correction to match the actual thermal effect. When the comprehensive correction factor < 1, it indicates that the original energy deposition density overestimates the heat source intensity, and the value decreases after correction to return to the true level. This correction eliminates the deviation caused by the three-dimensional non-uniformity of the heat source in the original energy deposition density distribution, resulting in a corrected internal heat source distribution that better matches the actual thermal effect of the target, ensuring that the internal heat source data can accurately reflect the true distribution state of the heat source inside the target.

[0083] Based on the modified internal heat source distribution (internal heat source conditions) and surface temperature boundary conditions (external constraint conditions) generated above, a three-dimensional unsteady-state heat conduction equation for the interior of the bismuth target is constructed. This equation adopts the spherical coordinate system heat conduction equation form from step 400, and the equation parameters are supplemented and improved by combining the modified internal heat source distribution (added to the heat source terms of the equation) and surface temperature boundary conditions (added to the boundary constraint terms of the equation). The finite element analysis method is used to numerically solve the equation, and the iteration accuracy requirement of step 400 is followed during the solution process (temperature difference between adjacent iteration nodes ≤ 1× After multiple iterations and convergence, a complete three-dimensional temperature distribution field along the axial and radial directions is obtained inside the entire bismuth target. This distribution field can reflect the spatial temperature distribution state of the target under the action of the corrected internal heat source. Each spatial coordinate (axial, radial, circumferential) in the distribution field corresponds to a unique temperature value, and the temperature at any position of the target can be directly read. It shows the difference in temperature between the core area (axial 10mm, radial 0mm) and the outer area (axial 2mm, radial 4mm). By comparing the temperature values ​​at different axial positions, the temperature decrease / increase law from one end of the target to the other can be reflected. By comparing the temperature values ​​at different radial positions, the temperature decay law from the center to the edge can be reflected. The numerical change trend of the distribution field is highly matched with the distribution characteristics of the corrected internal heat source. The area with high heat source intensity corresponds to higher temperature value, and the area with uniform heat source distribution corresponds to a gradual temperature change.

[0084] Step 503: Extract temperature values ​​at 2 mm intervals along the axial direction and at 1 mm intervals along the radial direction from the three-dimensional temperature distribution field to form the temperature gradient distribution at different depths along the axial direction and at different offset distances along the radial direction inside the bismuth target. Specifically, this includes: calling the complete three-dimensional temperature distribution field data obtained in step 502, extracting data according to a preset spatial interval, that is, along the axial direction of the bismuth target, from one end of the target to the other end, extracting the temperature value at a position every 2 mm, and recording the axial coordinate and corresponding temperature value of the position; along the radial direction of the bismuth target, from the center (base point) of the target to the edge of the target, extracting the temperature value at a position every 1 mm, and recording the radial coordinate and corresponding temperature value of the position. The extracted temperature values ​​at various axial depths and radial offsets are categorized and integrated to create axial temperature sequences (axial coordinates and temperature values) and radial temperature sequences (radial coordinates and temperature values). The axial and radial temperature gradients are calculated as follows: Temperature gradient = Temperature difference between two adjacent extraction points ÷ Spatial interval between the two extraction points (2 mm axial interval, 1 mm radial interval). This formula is used to calculate the temperature gradients at each adjacent axial and radial extraction point, resulting in axial and radial temperature gradient sequences. Finally, these two sequences are integrated to form the temperature gradient sequence for the bismuth target at different depths along the axial direction and radial direction. The temperature gradient distribution at different offset distances clearly quantifies the rate (gradient magnitude) and direction (positive or negative) of temperature change within the target body with spatial position. The larger the absolute value of the gradient, the more drastic the temperature change within a unit spatial interval. For example, when the axial gradient is 5℃ / mm, the temperature increases by 5℃ for every 1mm of depth, which is much higher than in the region with a gradient of 0.5℃ / mm. A positive gradient value indicates that the temperature increases with increasing axial depth / radial offset distance, while a negative gradient value indicates that the temperature decreases with spatial position. The smaller the gradient value, the more gradual the temperature change, and the more uniform the heat source distribution and the more stable the heat conduction in the corresponding region.

[0085] In this embodiment, the original photoacoustic signal is corrected using the photon radiation dose rate, effectively eliminating the interference of the radiation field on the amplitude of the photoacoustic signal; the surface temperature is used as a boundary condition, and a corrected internal heat source is obtained by combining a comprehensive correction factor, making the solution conditions of the heat conduction equation closer to the actual physical conditions, thus ensuring the accuracy and reliability of the calculated three-dimensional temperature distribution field; the axial and radial interval temperature data are extracted and a gradient distribution is formed, realizing a refined description of the internal temperature field of the bismuth target.

[0086] like Figure 2 As shown, embodiments of the present invention also provide a temperature gradient monitoring system for bismuth target dissolution reactions, comprising:

[0087] The acquisition module is used to acquire energy deposition density distribution data along the axial and radial directions inside the bismuth target and the activity data of each nuclide, and convert the nuclide activity data into nuclide mass data at the corresponding location;

[0088] The partitioning module is used to calculate the intensity of photon source terms and the corresponding energy group partitioning results generated by the spallation reaction inside the bismuth target based on the nuclide mass data. The photon source term intensity and the corresponding energy group partitioning results are input into the photon transport calculation model to perform photon transport calculation and obtain the photon radiation dose equivalent rate distribution outside the target.

[0089] The module is used to extract multiple closed isosurface clusters from energy deposition density distribution data and divide the target into nested three-dimensional subspace domains. Taking the point with the maximum energy deposition density as the base point, exploration nodes are arranged radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point.

[0090] The solution module is used to establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain based on the energy deposition density distribution data of each spatial exploration node on the boundary of each three-dimensional subspace domain and the three-dimensional closed vorticity loop, to obtain the heat source distribution characteristic value, and to normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target.

[0091] The inversion module is used to acquire the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and to perform temperature effect correction on the original photoacoustic spectral signal using the photon radiation dose equivalent rate distribution to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

[0092] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

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

[0094] 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 for monitoring the temperature gradient in a bismuth target dissolution reaction, characterized in that, The method includes: Step 1: Obtain energy deposition density distribution data and nuclide activity data along the axial and radial directions inside the bismuth target, and convert the nuclide activity data into nuclide mass data at the corresponding locations. This includes: constructing a three-dimensional solid model containing the geometric dimensions, material composition, and density parameters of the bismuth target; setting the energy, flux, and spatial distribution parameters of the incident proton beam; simulating the nuclear reaction and energy loss process between protons and bismuth target nuclei using the Monte Carlo method; during the simulation, dividing the target space into grids and recording proton transport trajectories and energy deposition events in each grid cell to obtain energy deposition density distribution data at different depths along the axial direction and different offset distances in the radial direction inside the target; simultaneously, obtaining nuclide activity data corresponding to each grid cell by statistically analyzing the types of nuclear reactions and product nuclides occurring in each grid cell; based on the conversion relationship between nuclide activity and mass, for the nuclide activity data in each grid cell, performing unit conversion based on the atomic mass number and decay constant of the corresponding nuclide, combined with Avogadro's constant, to convert the activity data into the corresponding nuclide mass value, forming nuclide mass data with a three-dimensional spatial coordinate index; Step 2: Based on the nuclide mass data, calculate the photon source term intensity and corresponding energy group partitioning results generated by the spallation reaction inside the bismuth target. Input the photon source term intensity and corresponding energy group partitioning results into the photon transport calculation model to perform photon transport calculations and obtain the photon radiation dose equivalent rate distribution outside the target. This includes: based on the nuclide mass data with three-dimensional spatial coordinate index, combined with the pre-constructed nuclide decay characteristic database and energy photon emission characteristic database, according to the nuclide mass values ​​contained in each spatial grid cell in the nuclide mass data, combined with the half-life parameters of the corresponding nuclide and the emission probability parameters of each energy photon, calculate the total intensity of energy photons released by all nuclides in each spatial grid cell, and according to the pre-set eighteen energy group structures, merge and statistically analyze the energy photons to obtain the energy photon source term intensity and corresponding energy group partitioning results of each spatial grid cell. The energy group partitioning results are used as the energy photon source term intensity of each grid cell and the energy group partitioning results as the spatial distribution source term input into the photon transport calculation model. The photon transport calculation model includes the geometric structure parameters of the bismuth target, the geometric structure parameters of the cooling layer, and the material property parameters of each structural component. By simulating the photoelectric effect, Compton scattering, and electron-pair effect physical processes in the interaction between energy photons and matter, the photon radiation dose equivalent rate is calculated at four discrete spatial coordinate points on the circumference of each monitoring section, which are set at five monitoring sections at axial distances of 10 mm, 20 mm, 30 mm, 40 mm, and 50 mm from the incident surface along the periphery of the bismuth target, respectively, corresponding to the radial coordinates of 0 degrees, 90 degrees, 180 degrees, and 270 degrees. This forms a photon radiation dose equivalent rate distribution corresponding to the spatial coordinates. Step 3: Extract multiple closed isosurface clusters from the energy deposition density distribution data and divide the target into nested three-dimensional subspace domains; using the point with the maximum energy deposition density as the base point, arrange exploration nodes radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point. Step 4: Based on the energy deposition density distribution data of each three-dimensional subspace domain boundary and each spatial exploration node on the three-dimensional closed vortex loop, establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain, obtain the heat source distribution characteristic value, and normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target. Step 5: Collect the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and use the photon radiation dose equivalent rate distribution to perform temperature effect correction on the original photoacoustic spectral signal to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

2. The method for monitoring the temperature gradient in the bismuth target dissolution reaction according to claim 1, characterized in that, Step 3 includes: From the energy deposition density distribution data, based on five preset energy deposition density gradient thresholds, spatial closed isosurface clusters corresponding to different threshold levels are extracted. These spatial closed isosurface clusters divide the internal space of the target into five three-dimensional subspace domains nested sequentially from the outside to the inside. Using the spatial coordinate point corresponding to the maximum energy deposition density in the energy deposition density distribution data as the base point, spatial exploration nodes are set radially within each three-dimensional subspace domain to form a spatial exploration node distribution array within each three-dimensional subspace domain. Within each three-dimensional subspace domain, with the base point as the starting origin, five exploration nodes in the same exploration direction within the corresponding three-dimensional subspace domain are connected radially outward in sequence. In the circumferential direction, corresponding nodes at the same radial distance position in adjacent exploration directions are connected. At the same time, corresponding nodes at the same radial distance position in the same exploration direction are connected axially between each radial plane, forming a three-dimensional closed vortex loop that surrounds the base point and runs through the entire three-dimensional subspace domain.

3. The method for monitoring the temperature gradient in the bismuth target dissolution reaction according to claim 2, characterized in that, From the energy deposition density distribution data, based on five preset energy deposition density gradient thresholds, spatially closed isosurface clusters corresponding to different threshold levels are extracted. These spatially closed isosurface clusters divide the internal space of the target into five nested three-dimensional subspace domains, from the outside in, including: For each preset energy deposition density gradient threshold, traverse all grid cells of the target space grid, compare the energy deposition density value at the eight vertices of each grid cell with the magnitude of the energy deposition density gradient threshold, and mark the grid cells that the isosurface passes through. For each marked grid cell, identify the edges whose energy deposition density values ​​at both ends are greater than and less than the energy deposition density gradient threshold, respectively, and these edges are the edges through which the isosurface passes. For each identified edge, calculate the coordinates of the intersection point between the isosurface and the corresponding edge based on the ratio of the difference between the energy deposition density values ​​at the two endpoints of the corresponding edge and the energy deposition density gradient threshold, combined with the spatial coordinates of the two endpoints. Within each marked grid cell, multiple intersection points are connected according to their spatial position to form triangular patches, generating isosurfaces corresponding to the energy deposition density gradient threshold. The above process is repeated for all five thresholds to obtain five spatially closed isosurface clusters nested from the outside in, dividing the target's internal space into five nested three-dimensional subspace domains.

4. The method for monitoring the temperature gradient in the bismuth target dissolution reaction according to claim 3, characterized in that, Within each three-dimensional subspace domain, spatial exploration nodes are set radially, including: Starting from the inner boundary of the corresponding three-dimensional subspace domain and ending at the outer boundary, five exploration nodes are uniformly set radially. At the same time, an exploration direction is set every 72 degrees in the circumferential direction, so that the five exploration nodes in each exploration direction are located at different radial distances within the corresponding subspace domain.

5. The method for monitoring the temperature gradient in the bismuth target dissolution reaction according to claim 4, characterized in that, Step 4 includes: The energy deposition density values ​​of each spatial exploration node on the three-dimensional closed vortex loop and the energy deposition density values ​​on the boundaries of each three-dimensional subspace domain are extracted. Based on the general solution of the heat conduction equation in spherical coordinates, the energy deposition density values ​​on the boundaries of each three-dimensional subspace domain are used as boundary conditions, and the energy deposition density values ​​of each spatial exploration node on the three-dimensional closed vortex loop are used as internal constraints to establish the heat source contribution characteristic equation for each three-dimensional subspace domain. The heat source contribution characteristic equation for each three-dimensional subspace domain is solved to obtain the heat source distribution characteristic value corresponding to each three-dimensional subspace domain. The heat source distribution characteristic values ​​corresponding to each three-dimensional subspace domain are arranged in a spatial topological order from the outside to the inside or from the inside to the outside. The normalized weighted fusion calculation is performed with the volume ratio of each three-dimensional subspace domain as the weight to generate a comprehensive correction factor that characterizes the three-dimensional non-uniform distribution of heat sources inside the target.

6. The method for monitoring the temperature gradient in a bismuth target dissolution reaction according to claim 5, characterized in that, Step 5 includes: From the photon radiation dose equivalent rate distribution, the photon radiation dose equivalent rate values ​​corresponding to the positions of twenty discrete spatial coordinate points are extracted. Using the photon radiation dose equivalent rate values ​​at the positions of the twenty discrete spatial coordinate points, the original photoacoustic signal amplitude at the same discrete spatial coordinate point position is subjected to temperature effect correction processing to obtain the corrected photoacoustic signal amplitude at each discrete spatial coordinate point. Based on the pre-established correspondence between the amplitude of the characteristic gas photoacoustic spectral signal and the temperature, the amplitude of the corrected photoacoustic signal at each discrete spatial coordinate point is converted into the surface temperature value at the corresponding discrete spatial coordinate point. By using the surface temperature value as a boundary condition and combining the energy deposition density distribution data with a comprehensive correction factor, the corrected internal heat source distribution is obtained. Based on the corrected internal heat source distribution and the surface temperature boundary condition, the heat conduction equation inside the bismuth target is solved to obtain the three-dimensional temperature distribution field along the axial and radial directions inside the bismuth target. Temperature values ​​at 2 mm intervals along the axial direction and 1 mm intervals along the radial direction are extracted from the three-dimensional temperature distribution field to form the temperature gradient distribution at different depths along the axial direction and at different offset distances along the radial direction inside the bismuth target.

7. A temperature gradient monitoring system for bismuth target dissolution reaction, the system implementing the method as described in any one of claims 1 to 6, characterized in that, include: The acquisition module is used to acquire energy deposition density distribution data along the axial and radial directions inside the bismuth target and the activity data of each nuclide, and convert the nuclide activity data into nuclide mass data at the corresponding location; The partitioning module is used to calculate the intensity of photon source terms and the corresponding energy group partitioning results generated by the spallation reaction inside the bismuth target based on the nuclide mass data. The photon source term intensity and the corresponding energy group partitioning results are input into the photon transport calculation model to perform photon transport calculation and obtain the photon radiation dose equivalent rate distribution outside the target. The module is used to extract multiple closed isosurface clusters from energy deposition density distribution data and divide the target into nested three-dimensional subspace domains. Taking the point with the maximum energy deposition density as the base point, exploration nodes are arranged radially in each three-dimensional subspace domain to construct a three-dimensional closed vorticity loop around the base point. The solution module is used to establish and solve the heat source contribution characteristic equation of each three-dimensional subspace domain based on the energy deposition density distribution data of each spatial exploration node on the boundary of each three-dimensional subspace domain and the three-dimensional closed vorticity loop, to obtain the heat source distribution characteristic value, and to normalize and weight the heat source distribution characteristic value according to the spatial topology order to generate a comprehensive correction factor characterizing the three-dimensional non-uniform distribution of heat sources inside the target. The inversion module is used to acquire the original photoacoustic spectral signal of the characteristic gas around the bismuth target, and to perform temperature effect correction on the original photoacoustic spectral signal using the photon radiation dose equivalent rate distribution to obtain the corrected photoacoustic signal. Combined with the comprehensive correction factor, the temperature gradient distribution along the axial and radial directions inside the bismuth target is inverted.

8. A computer-readable storage medium, characterized in that, The computer-readable storage medium is used to store a computer program for performing the method as described in any one of claims 1 to 6.

Images