Tokamak anti-shear alfvén continuum mode damping calculation method

Abstract

The application discloses a Tokamak anti-shear Alfven eigenmode continuum damping calculation method. The application can effectively avoid the numerical divergence and instability problems caused by the traditional real axis integral across the Alfven resonance point, improve the extraction precision of small-scale continuum damping rate, and solve the complex global eigenvalue problem into the local integral along the radial direction and the matching problem by combining the shooting integral with the iterative solution, and improve the parameter scanning efficiency by constructing the Alfven continuum equation based on the equilibrium reconstruction data, automatically identifying the Alfven resonance point, constructing the integral path around the singularity on the complex plane, and combining the logarithmic derivative matching conditions of the left and right side solutions at the matching point to solve the complex eigenfrequency.

Description

Technical Field

[0001] This invention relates to the technical field of magnetically confined nuclear fusion plasma, specifically to a method for calculating the damping of the continuous spectrum of the anti-shear Alfvén eigenmode in a tokamak. Background Technology

[0002] Magnetic confinement fusion is a crucial pathway to achieving clean and sustainable energy. In tokamak devices, advanced operating modes (such as anti-shear configurations) are considered one of the core operating scenarios for future fusion reactors (such as ITER and CFETR) due to their superior confinement performance. The Reversed Shear Alfvén Eigenmode (RSAE), also known as the Alfvén cascade mode, is one of the most likely operating modes to occur in advanced tokamak devices. Its occurrence, conditions, and stability characteristics directly relate to the engineering implementation of fusion reactors, the diagnosis of plasma operating conditions during operation, and the evaluation of fusion confinement performance. Its clear engineering application value is mainly reflected in the following two aspects:

[0003] First, it serves as a diagnostic platform for advanced operating modes. The frequency of RSAE affects the minimum safety factor ( The RSAE is extremely sensitive, and its frequency evolution over time and plasma parameters is closely related to changes in the safety factor profile. By experimentally measuring the frequency evolution of RSAE, the safety factor profile inside the plasma can be retrieved with high precision. Compared to other active diagnostic methods that require a neutral beam, such as the Motional Stark Effect (MSE), RSAE not only provides non-disturbing diagnostic basis for the establishment, maintenance, and control of the anti-shear configuration, which is beneficial to improving the accuracy and real-time performance of device operation control, but also effectively avoids interference with the plasma state, making it a non-destructive auxiliary diagnostic method.

[0004] Second, it is an important reference indicator for evaluating the performance of combustion plasma. In future combustion plasma devices, the high energy generated by the fusion reaction... High-energy particles generated by assisted heating (such as neutral beam injection (NBI) and electron cyclotron resonance heating (ECRH)) can excite RSAEs through resonant interactions. Experimental measurements have observed shear Alfvén mode instabilities driven by high-energy particles in tokamak devices with high-power assisted heating (such as JET and JT-60). Excited RSAEs lead to radial transport and redistribution of high-energy particles, which can cause loss of high-energy particles, reduce fusion power, and even damage the first wall in severe cases. Research results indicate that the minimum safety factor is a sensitive parameter controlling mode stability: in deep antimagnetic shear configurations, the minimum surface can form a transport barrier, resulting in good confinement of high-energy particles; while in intermediate antimagnetic shear configurations, the confinement of high-energy particles degrades significantly, potentially leading to avalanche transport. Therefore, the assessment results of RSAE stability can be directly used for operating parameter selection, heating scheme optimization, configuration control, and wall load risk warning, which has significant practical value. Existing experiments and studies have shown that the minimum safety factor and antimagnetic shear strength are important parameters affecting RSAE stability. Under different antimagnetic shear configurations, the confinement state of high-energy particles varies significantly, leading to different RSAE driving conditions and transport consequences. Therefore, quantitatively assessing whether RSAEs will be excited by high-energy particles before or during device operation is a crucial step in optimizing the operation mode and controlling risks in advanced tokamak systems.

[0005] Whether high-energy particles can successfully excite RSAE depends on whether their driving growth rate can overcome various damping mechanisms from the thermal plasma. Among various damping mechanisms (such as Landau damping, electron collision damping, radiation damping, etc.), continuum-spectrum damping is recognized as one of the most important and difficult-to-calculate mechanisms.

[0006] The physical essence of continuous spectrum damping is resonant absorption. In tokamak plasmas, the frequency of sheared Alfvén waves satisfies the local dispersion relation. Parallel wavenumber Variation with radial position. Due to safety factor. Radial dependence, local Alfvén frequency As the radial position changes continuously, the so-called Alfvén continuous spectrum is formed—that is, the frequency band formed by the local Alfvén frequencies corresponding to different radial positions. When the frequency of a discrete eigenmode (such as RSAE) happens to fall within the frequency range of the Alfvén continuous spectrum, the eigenmode and the continuous spectrum resonate at the resonance layer (i.e., the Alfvén resonance point, satisfying...). Strong coupling occurs between the eigenmode frequency and the continuous spectrum (the position of the resonant layer), and the energy of the eigenmode is resonantly transferred to the sheared Alfvén continuous spectrum, resulting in irreversible loss of eigenmode energy. This process is known as continuous spectrum damping, and its strength depends on the relative position of the eigenmode frequency and the continuous spectrum, the width of the resonant layer, and the plasma parameters.

[0007] From an engineering application perspective, the continuum damping rate is a fundamental input for determining whether a RSAE can be excited, whether suppression measures are needed, and how to optimize operating parameters. Only by obtaining relatively accurate continuum damping assessment results can we further compare them with high-energy particle driving terms, thereby providing RSAE stability criteria, instability thresholds, and risk level assessment results, and providing a basis for the parameter design, experimental scheme formulation, and operational control of tokamak devices.

[0008] Currently, the calculation of Alfvén eigenmodes in tokamak mainly employs global MHD eigenvalue codes (such as NOVA). These methods discretize the magnetohydrodynamic equations in the global space, directly solve for the eigenvalues ​​and eigenfunctions, and then evaluate the stability and damping characteristics of the mode.

[0009] When conducting engineering and quantitative evaluation of RSAE continuum damping under advanced tokamak operating modes, existing solutions have the following key shortcomings:

[0010] First, it is difficult to stably handle the singularities near the Alfvén resonance point, which affects the reliability of the damping assessment results.

[0011] RSAE continuous spectrum damping originates from the resonant coupling between the eigenmodes and the Alfvén continuous spectrum. The differential equation describing RSAE contains coefficients. At the Alfvén resonance point, The equations exhibit physical singularities. Traditional MHD code directly discretizes and integrates across this singularity on the real-axis grid, leading to numerical divergence or the introduction of non-physical artificial dissipations. These issues directly affect the calculation results of the continuous-spectrum damping rate, making it difficult to meet the computational stability and reliability requirements for device operation evaluation and control decisions.

[0012] Second, the extraction accuracy of small damping ratio is insufficient, making it difficult to support the quantitative judgment of the high-energy particle driving threshold.

[0013] Under actual tokamak operating parameters, the continuous-spectral damping ratio of the RSAE is typically much smaller than its mode frequency, usually... Even smaller. In global MHD code, such small damping is easily masked by numerical discretization errors, insufficient mesh resolution, or improper boundary condition handling, leading to insufficient confidence in the calculation results. Therefore, existing methods are difficult to use directly as reliable quantitative tools for RSAE instability threshold assessment, high-energy particle confinement risk prediction, and runtime parameter optimization.

[0014] Third, the computational cost is high, making it difficult to meet the needs of rapid parameter scanning in engineering applications.

[0015] In practical applications, RSAE continuum damping assessment often requires multiple calculations for parameters such as minimum safety factor, antimagnetic shear strength, density profile, pressure gradient, magnetic field strength, and heating scheme to determine the safe operating window and optimize the operating scheme. Global MHD eigenvalue calculation requires high-resolution discretization of the entire plasma region and global matrix solving. For research scenarios that require systematic scanning of equilibrium parameters (such as minimum safety factor and pressure gradient), the computational cost is high and it is difficult to meet the needs of large-scale parameter scanning.

[0016] Fourth, it lacks adaptability to actual equipment operation optimization.

[0017] Existing methods are more geared towards general state analysis than specifically optimizing for the RSAE continuum damping under antimagnetic shear configurations. In practical applications, facility operators are more concerned with whether the RSAE can be excited, the magnitude of the damping rate, which operating parameters will significantly alter the stability boundary, and how to adjust the safety factor profile or heating scheme to reduce risk, given an equilibrium profile and high-energy particle conditions. Existing global MHD codes typically cannot directly provide high-precision, repeatable, and rapidly scannable damping assessment results for engineering decision-making, limiting their application in advanced operating mode optimization.

[0018] In summary, existing global MHD eigenvalue methods suffer from several drawbacks when used for RSAE continuum damping calculations. These include instability in handling resonance singularities, insufficient accuracy in extracting small damping rates, low parameter scanning efficiency, and inadequate adaptability to engineering applications. Consequently, these methods struggle to meet the future application requirements of fusion devices for quantitative assessment of RSAE stability, prediction of high-energy particle transport risks, and optimization of operating parameters. Therefore, it is necessary to develop a dedicated method for RSAE continuum damping calculations to improve computational accuracy, numerical stability, and parameter scanning efficiency. Summary of the Invention

[0019] The technical problem this invention aims to solve is that existing global MHD eigenvalue methods, when used for RSAE continuum damping calculations, suffer from instability in handling resonance singularities, insufficient accuracy in extracting small damping rates, low parameter scanning efficiency, and insufficient adaptability to engineering applications. These issues make it difficult to meet the application requirements of future fusion devices for quantitative assessment of RSAE stability, prediction of high-energy particle transport risks, and optimization of operating parameters. The purpose is to provide a tokamak antishear Alfvén eigenmode continuum damping calculation method, thereby solving the aforementioned technical problems.

[0020] This invention is achieved through the following technical solution:

[0021] In a first aspect, the present invention provides a method for calculating the continuous spectrum damping of the tokamak anti-shear Alfvén eigenmode, comprising:

[0022] Tokamak equilibrium reconstruction data and preset circumferential modulus are obtained, wherein the equilibrium reconstruction data includes safety factor profile, plasma density, ion temperature, electron temperature, magnetic field strength, large radius and small radius;

[0023] The matching point is determined based on the minimum value location of the safety factor profile, and the polar modulus is determined based on the preset circumferential modulus.

[0024] Based on the equilibrium reconstruction data, the preset circumferential and poloidal moduli, the Alfvén continuous spectrum equation characterizing the radial mode structure of the anti-shear Alfvén intrinsic modes is constructed, and the Alfvén resonance point position in the Alfvén continuous spectrum equation is identified.

[0025] Construct an integration path around the Alfvén resonance point in the complex plane, analytically extend the Alfvén continuous spectrum equation on the real axis to the complex plane, and perform target integration from the left and right boundaries of the radial computational domain along the integration path to the matching point to obtain the left-end solution and the right-end solution at the matching point.

[0026] The residual function is determined by the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point, and the complex eigenfrequency when the residual function is zero is solved iteratively.

[0027] The imaginary part of the complex intrinsic frequency is used as the output of the continuous spectrum damping rate.

[0028] Further, the step of determining the matching point based on the minimum position of the safety factor profile and determining the polar modulus based on the preset circumferential modulus includes:

[0029] A one-dimensional search is performed on the safety factor profile along the radial direction to determine the position where the first derivative of the safety factor profile with respect to the radial position is zero, and this position is determined as the matching point.

[0030] Calculate the product of the preset circumferential modulus and the minimum value of the safety factor profile, and round the product to obtain the polar modulus.

[0031] Furthermore, the step of constructing the Alfvén continuous spectrum equation characterizing the radial mode structure of the anti-shear Alfvén eigenmodes based on equilibrium reconstruction data, preset circumferential moduli, and poloidal moduli includes:

[0032] Based on the safety factor profile, large radius, preset circumferential modulus, and polar modulus, the radial distribution of parallel wavenumbers is determined; wherein, the parallel wavenumber is the product of the preset circumferential modulus and the safety factor profile, minus the polar modulus, and then divided by the product of the safety factor profile and the large radius.

[0033] The distribution of the poloidal wavenumber along the radial direction is determined based on the poloidal modulus and the radial position; wherein the poloidal wavenumber is the poloidal modulus divided by the current radial position;

[0034] The radial distribution of Alfvén velocity is determined based on magnetic field strength and plasma density; wherein, the Alfvén velocity is the square root of the magnetic field strength divided by the product of vacuum permeability, plasma density, and ion mass.

[0035] The radial distribution of the target resonance coefficient is determined based on the parallel wavenumber distribution, the Alfvén velocity distribution, and the geodesic frequency; wherein the target resonance coefficient is the square of the current frequency minus the product of the square of the parallel wavenumber distribution and the square of the Alfvén velocity distribution, and then minus the square of the geodesic frequency.

[0036] Based on the target resonance coefficient and poloidal wavenumber distribution, a second-order singular differential equation with electrostatic potential as the solution variable is established as the Alfvén continuous spectrum equation.

[0037] Furthermore, the step of identifying the location of the Alfvén resonance point in the Alfvén continuous spectrum equation includes:

[0038] For the current frequency, calculate the real part of the target resonance coefficient point by point along the radial grid, and scan the sign change of the real part along the radial direction;

[0039] When the product of the real values ​​of two adjacent radial grid points is negative, it is determined that there is an Alfvén resonance point between the two radial grid points. The radial position of the Alfvén resonance point is calculated by linear interpolation, and the positions of the left and right resonance layers are output.

[0040] When the product of the real values ​​at all adjacent radial grid points is non-negative, the current frequency is determined to be outside the continuous spectrum.

[0041] Furthermore, the step of constructing an integration path around the Alfvén resonance point in the complex plane, analytically extending the Alfvén continuity spectrum equation on the real axis to the complex plane, and performing target integration from the left and right boundaries of the radial computational domain along the integration path towards the matching point to obtain the left-hand and right-hand solutions at the matching point includes:

[0042] Determine the detour radius; wherein the detour radius is less than the radial distance between the Alfvén resonance point and the matching point, and such that the constructed integral path does not touch other Alfvén resonance points;

[0043] For the Alfvén resonance point located to the left of the matching point, construct the left integration path; wherein, integrate along the real axis from the left boundary to the left Alfvén resonance point position minus the radius of the circumference; from this position, enter the complex plane along the first complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the second complex direction coefficient, reaching the left Alfvén resonance point position plus the radius of the circumference; from this position, integrate along the real axis to the matching point;

[0044] For the Alfvén resonance point located to the right of the matching point, construct a right-side integration path; wherein, integrate along the real axis from the right boundary to the position of the right Alfvén resonance point plus the detour radius; from this position, enter the complex plane along the second complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the first complex direction coefficient, reaching the position of the right Alfvén resonance point minus the detour radius; from this position, integrate along the real axis to the matching point;

[0045] The Alfvén continuous spectrum equation is analytically extended along the left and right integration paths, and non-zero amplitude and zero derivative are used as boundary conditions. The equations are then integrated towards the matching point along the left and right integration paths, respectively, to obtain the left-hand and right-hand solutions.

[0046] Further, the step of determining the residual function based on the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point includes:

[0047] The left-hand logarithmic derivative is calculated based on the function value of the left-hand solution at the matching point and its first derivative with respect to the radial position; wherein, the left-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value;

[0048] The right-hand logarithmic derivative is calculated based on the function value of the right-hand solution at the matching point and its first derivative with respect to the radial position; wherein, the right-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value;

[0049] The difference between the logarithmic derivative on the left and the logarithmic derivative on the right is defined as the residual function; wherein the matching condition for the logarithmic derivative is that the residual function takes the value of zero.

[0050] Furthermore, the step of iteratively solving for the complex eigenfrequency when the residual function is zero includes:

[0051] The local Alfvén frequency at the matching point is determined based on the values ​​of the parallel wavenumber and the Alfvén velocity at the matching point.

[0052] Based on the local Alfvén frequency, the search intervals for the real part and the imaginary part of the complex frequency are determined. A two-dimensional grid is established within the search intervals, and the modulus of the residual function at each grid point is calculated. The complex frequency corresponding to the minimum value of the modulus is determined as the initial guessed complex frequency.

[0053] Using the initial guessed complex frequency as the initial value, the complex frequency is iteratively adjusted so that the magnitude of the residual function approaches zero until the preset convergence condition is met. The complex frequency that meets the preset convergence condition is determined as the complex eigenfrequency.

[0054] Furthermore, after outputting the imaginary part of the complex eigenfrequency as the continuous spectrum damping rate, the method further includes:

[0055] To obtain the driving growth rate of high-energy particles in a tokamak device;

[0056] The continuous spectrum damping rate is compared with the high-energy particle-driven growth rate. When the high-energy particle-driven growth rate is greater than the absolute value of the continuous spectrum damping rate, the anti-shear Alfvén eigenmode is determined to be unstable.

[0057] In a second aspect, the present invention provides an electronic device, comprising: a memory, and one or more processors communicatively connected to the memory; the memory stores instructions executable by the one or more processors, the instructions being executed by the one or more processors to cause the one or more processors to implement the method described above.

[0058] Thirdly, the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method.

[0059] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0060] This invention constructs the Alfvén continuous spectrum equation based on equilibrium reconstruction data, automatically identifies Alfvén resonance points, constructs an integral path around the singularity on the complex plane, and solves for the complex eigenfrequency by combining the matching condition of the logarithmic derivatives of the solutions on both sides at the matching point. This effectively avoids the numerical divergence and instability problems caused by traditional real-axis integration crossing the Alfvén resonance point, and improves the extraction accuracy of small-scale continuous spectrum damping rates. At the same time, by adopting a combination of target integration and iterative solution, the complex global eigenvalue problem is transformed into a local integration and matching problem along the radial direction, reducing computational costs and improving parameter scanning efficiency. Furthermore, after outputting the imaginary part of the complex eigenfrequency as the continuous spectrum damping rate, it can be compared with the high-energy particle-driven growth rate, thus providing a quantitative basis for RSAE stability judgment, operating window determination, and tokamak device operating parameter optimization. Attached Figure Description

[0061] To more clearly illustrate the technical solutions of the exemplary embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of the present invention and should not be considered as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort. In the drawings:

[0062] Figure 1 A flowchart illustrating an accurate calculation method for the continuous spectrum damping of the tokamak anti-shear Alfvén eigenmode based on the complex path shooting method, provided in the embodiments of this specification;

[0063] Figure 2 This is a schematic diagram of the Alfvén continuous spectrum and RSAE eigenfrequencys provided in the embodiments of this specification;

[0064] Figure 3 This is a schematic diagram of the complex plane path integral provided in the embodiments of this specification;

[0065] Figure 4 This is a schematic diagram of the module structure provided in the embodiments of this specification; wherein, Figure 4 The vertical axis is The horizontal axis is the radius r;

[0066] Figure 5 This is a schematic diagram showing the matching of the logarithmic derivatives of the solutions on both sides of the module structure provided in the embodiments of this specification at the matching point; wherein, Figure 5 The vertical axis is The horizontal axis is the radius r;

[0067] Figure 6 The real part of the frequency provided in the embodiments of this specification varies. A schematic diagram illustrating the changing trend; where, Figure 5 The vertical axis represents frequency;

[0068] Figure 7 The frequency imaginary part provided in the embodiments of this specification follows A schematic diagram illustrating the changing trend;

[0069] Figure 8 This is a schematic diagram of the mode structure and mode frequencies calculated by NOVA as provided in the embodiments of this specification. Detailed Implementation

[0070] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the embodiments and accompanying drawings. The illustrative embodiments and descriptions of the present invention are only used to explain the present invention and are not intended to limit the present invention.

[0071] This embodiment provides a method for calculating the continuous spectrum damping of the tokamak antishear Alfvén eigenmode, including:

[0072] Step S11: Obtain tokamak equilibrium reconstruction data and preset circumferential modulus, wherein the equilibrium reconstruction data includes safety factor profile, plasma density, ion temperature, electron temperature, magnetic field strength, large radius and small radius.

[0073] In this embodiment, the tokamak equilibrium reconstruction data can be plasma equilibrium configuration data obtained by inverting diagnostic data such as magnetic measurement signals and external coil currents during a tokamak device discharge experiment. This equilibrium reconstruction data carries radial distribution information of various physical quantities within the plasma.

[0074] Specifically, the equilibrium reconstruction data can be provided by a structured data file output by the equilibrium reconstruction program. The executing entity can read and extract the radial distribution of the safety factor profile, plasma density profile, ion temperature profile, and electron temperature profile along the normalized small radius from this file, and simultaneously extract the magnetic field strength, large radius, and small radius as magnetic field geometric parameters.

[0075] More specifically, the normalized small radius is the ratio of the radial position to the small radius, ranging from 0 to 1, corresponding to the plasma core and the outermost magnetic surface, respectively. The safety factor profile characterizes the helicity of the magnetic field lines. In a typical anti-shear configuration, the safety factor profile exhibits a non-monotonic shape that first decreases and then increases, with a minimum region.

[0076] The plasma density profile is expressed in terms of particle number density, while the ion temperature profile and electron temperature profile describe the radial distribution of thermal kinetic energy of ion and electron populations, respectively.

[0077] In this embodiment, the executing entity can receive a preset circumferential modulus as input. The preset circumferential modulus represents the number of periods of the anti-shear Alfvén eigenmode along the tokamak circumference and is determined by the most unstable circumferential modulus observed in experiments.

[0078] Step S12: Determine the matching point based on the minimum position of the safety factor profile, and determine the polar modulus based on the preset circumferential modulus.

[0079] In this embodiment, the matching point is the convergence position of the target integral, and also the spatial reference point for matching the logarithmic derivatives of the eigenfunctions on the left and right sides.

[0080] Specifically, the executing entity can scan the values ​​of the safety factor profile point by point along the radial grid, locate the grid interval where the minimum value is located, and solve the precise radial position where the first derivative of the safety factor with respect to the radial position is 0 through one-dimensional search, and determine the radial position as the matching point.

[0081] In this embodiment, the polar modulus represents the number of periods of the anti-shear Alfvén eigenmode along the tokamak polar direction. The execution entity can extract the minimum value of the safety factor profile at the matching point, calculate the product of the preset circumferential modulus and the minimum value, and take the integer value of the product to obtain the polar modulus.

[0082] Understandably, the above calculation method ensures that near the minimum safety factor surface, the difference between the product of the preset circumferential modulus and the minimum safety factor minus the polar modulus approaches zero, thus satisfying the physical resonance condition.

[0083] Step S13: Based on the equilibrium reconstruction data, preset circumferential modulus and poloidal modulus, construct the Alfvén continuous spectrum equation characterizing the radial mode structure of the anti-shear Alfvén intrinsic mode, and identify the Alfvén resonance point position in the Alfvén continuous spectrum equation.

[0084] In this embodiment, the Alfvén continuous spectrum equation is a second-order singular differential equation used to describe the radial electrostatic potential distribution structure of the anti-shear Alfvén eigenmodes.

[0085] Specifically, the executing entity can read the aforementioned parameters and perform the following coefficient quantization calculations at each grid point in the radial direction:

[0086] First, calculate the product of the preset circumferential modulus and the safety factor profile, subtract the polar modulus, and divide by the product of the safety factor profile and the large radius to obtain the radial distribution of the parallel wavenumber.

[0087] Second, the radial distribution of Alfvén velocity is obtained by dividing the magnetic field strength by the square root of the product of vacuum permeability, plasma density, and ion mass.

[0088] Third, divide the poloidal modulus by the current radial position to obtain the radial distribution of the poloidal wavenumber.

[0089] Then, for the currently input probe frequency, the executing entity can calculate the square of that frequency, subtract the product of the square of the parallel wavenumber and the square of the Alfvén velocity, and then subtract the square of the geodesic frequency. The final calculation result is used as the radial distribution of the target resonance coefficient. Based on the target resonance coefficient and the poloidal wavenumber, the executing entity can establish a second-order singular differential equation with electrostatic potential as the solution variable, as the Alfvén continuous spectrum equation.

[0090] In this embodiment, the executing entity can further identify the location of the Alfvén resonance point based on the target resonance coefficient. For the current frequency, the executing entity can extract the real part value of the target resonance coefficient point by point on the radial grid and scan the sign change of the real part value between adjacent grid points. When the product of the real part values ​​on two adjacent grid points is negative, it indicates that the real part has undergone a zero-crossing deflection in that interval, and the executing entity can determine that an Alfvén resonance point exists in that interval.

[0091] The execution unit can perform linear interpolation using the real part values ​​at two grid points as weights to obtain the precise radial position where the real part is strictly 0. The execution unit can extract the interpolated positions to the left and right of the matching point, outputting them as the left and right resonance layer positions, respectively. If, for the current frequency, the product of the real part values ​​at all adjacent grid points is non-negative, then the current frequency is determined to be outside the Alfvén continuous spectrum.

[0092] Step S14: Construct an integration path around the Alfvén resonance point on the complex plane, analytically extend the Alfvén continuous spectrum equation on the real axis to the complex plane, and perform target integration from the left and right boundaries of the radial computational domain along the integration path to the matching point to obtain the left-hand and right-hand solutions at the matching point.

[0093] In this embodiment, in order to avoid the singularity of the differential equation caused by the target resonance coefficient being 0 at the resonance point, the executing entity can analytically extend the Alfvén continuous spectrum equation on the real axis to the complex plane.

[0094] Specifically, the executing entity can construct a left-side integral path for the left-side resonance layer position with a preset detour radius as the interval.

[0095] More specifically, the integral is performed linearly along the real axis from the left boundary of the radial computational domain to the left resonance layer position minus the circumduction radius. Then, it deflects along a first complex direction (increasing real part and decreasing imaginary part) into the lower complex plane, traversing a circumduction path segment within the complex plane. After circumduction, it returns to the real axis along a second complex direction (increasing real part and increasing imaginary part), reaching the left resonance layer position plus the circumduction radius. Finally, the integral continues linearly along the real axis to the matching point.

[0096] For the right-side resonance layer location, the executing entity can construct a symmetrical right-side integration path.

[0097] More specifically, integrate along the real axis from the right boundary to the location of the right resonant layer plus the radius of travel. Deflect along the second complex direction into the complex plane for circumduction. Then return along the real axis along the first complex direction to the location of the right resonant layer minus the radius of travel. Finally, integrate along the real axis to the matching point.

[0098] The execution entity can use a variable-step-size ordinary differential equation numerical solver to perform the above integration. For example, it can use MATLAB's `ode15s` function to perform continuous integration of the complex path segment and the real-axis path segment. After the complex path segment integration is completed, the execution entity can use the obtained complex numerical result as the initial condition for the next real-axis integration segment. The left and right boundary conditions are set to have a preset non-zero amplitude value for the electrostatic potential and a first derivative of 0. The execution entity can integrate along the left and right integration paths to the matching point, extract the electrostatic potential value and first derivative when the left side reaches the matching point as the left-hand solution, and extract the electrostatic potential value and first derivative when the right side reaches the matching point as the right-hand solution.

[0099] Step S15: Determine the residual function based on the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point, and iteratively solve for the complex eigenfrequency when the residual function is zero.

[0100] In this embodiment, at the matching point, the executing entity can calculate the logarithmic derivatives of the left-hand and right-hand solutions respectively, that is, by dividing the first derivative by the corresponding electrostatic potential value, the left-hand and right-hand logarithmic derivatives are obtained respectively. The executing entity can subtract the right-hand logarithmic derivative from the left-hand logarithmic derivative and use the difference to construct a residual function.

[0101] In this embodiment, the executing entity can employ a two-stage strategy to iteratively solve for the complex eigenfrequency that makes the residual function zero. In the coarse scan stage, the local Alfvén frequency at the matching point is calculated, which is the product of the parallel wavenumber and the Alfvén velocity. The executing entity can then construct a two-dimensional search grid for complex frequencies based on this local Alfvén frequency, traverse the grid nodes, and complete all the aforementioned calculation steps to find the complex frequency that minimizes the magnitude of the residual function, using this as the initial guessed frequency.

[0102] In the fine solution stage, the execution entity can start from the initial guessed frequency and continuously update the complex frequency using a nonlinear numerical iterative algorithm. For example, the Levenberg-Marquardt algorithm can be used to iteratively solve the complex frequency with a residual function of zero, and the residual function can be recalculated in each iteration.

[0103] When the magnitude of the residual function is less than the preset function tolerance, and the frequency change in adjacent iterations is less than the preset step size tolerance, the iteration is terminated, and the complex frequency that converges at this time is determined as the complex eigenfrequency.

[0104] Step S16: The imaginary part of the complex intrinsic frequency is output as the continuous spectrum damping rate.

[0105] In this embodiment, the executing entity can extract the imaginary part data from the converged complex eigenfrequency and use the imaginary part data as the continuous spectrum damping rate of the tokamak antishear Alfvén eigenmode.

[0106] In this embodiment, the executing entity can output the continuous spectrum damping rate, along with related data such as the real part of the intrinsic frequency, the radial distribution of the intrinsic function, and the location of the resonant layer, to a storage medium or display terminal in a structured format.

[0107] In some implementations, the steps of determining the matching point based on the minimum position of the safety factor profile and determining the polar modulus based on a preset circumferential modulus include:

[0108] Step S121: Perform a one-dimensional search along the radial direction of the safety factor profile to determine the position where the first derivative of the safety factor profile with respect to the radial position is zero, and determine this position as the matching point.

[0109] In this embodiment, the location of the matching point utilizes the non-monotonic characteristic of the safety factor profile under the antishear configuration. In the antishear configuration, the safety factor profile first decreases and then increases from the plasma core outwards, reaching a minimum value at a certain radial position in the middle. The matching point is precisely the radial position where this minimum value is located, and this position is also the spatial center of the antishear Alfvén eigenmode structure.

[0110] In one possible and specific implementation, the implementing entity can perform a one-dimensional search to locate matching points, and the specific process may include:

[0111] First, the executing entity can read the discrete values ​​of the safety factor profile at each grid point from the acquired radial grid. The radial grid is a one-dimensional ordered sequence of points along a normalized small radius direction, with the number of grid points being a preset interpolation grid number. The discrete values ​​of the safety factor profile are arranged in ascending order of radial position, forming a one-dimensional array.

[0112] Then, the executing entity can advance point by point along the radial direction from smallest to largest, approximating the first derivative of the safety factor with respect to the radial position between every two adjacent grid points. Specifically, the executing entity can extract the difference in the safety factor between two adjacent grid points and extract the radial coordinate spacing between these two grid points. The ratio of this difference to the radial coordinate spacing can be calculated and used as an approximate value of the first derivative over this interval. The sign of this approximate value indicates the increasing or decreasing trend of the safety factor over this interval. A positive value indicates that the safety factor increases with increasing radial position, while a negative value indicates that the safety factor decreases with increasing radial position.

[0113] Next, the execution unit can scan the sign changes of the approximate values ​​of the first derivative in adjacent intervals. In the backshear configuration, the safety factor profile first decreases and then increases, so the sign of the first derivative changes from negative to positive radially. When the approximate value of the first derivative in the current grid interval is detected to be negative, while the approximate value of the first derivative in the next adjacent interval is positive, it is determined that the zero-crossing deflection point of the first derivative lies within the grid range corresponding to these two intervals. The execution unit can further narrow the search range, performing a secondary search within this region with finer step sizes, gradually approaching the radial position where the first derivative is precisely zero.

[0114] Finally, the executing entity can determine the precise radial position where the first derivative is 0 as the matching point. If the absolute value of the approximate value of the first derivative has fallen below the preset convergence threshold during the quadratic search, the current search position can be directly locked as the matching point. Simultaneously, the safety factor value corresponding to this matching point can be recorded and used as the minimum safety factor, which will serve as the input parameter for subsequent calculations of the polar modulus.

[0115] Step S122: Calculate the product of the preset circumferential modulus and the minimum value of the safety factor profile, and round the product to obtain the polar modulus.

[0116] In this embodiment, the determination of the polar modulus can rely on two input variables: a preset circumferential modulus and the minimum safety factor obtained during the one-dimensional search and positioning process.

[0117] In one possible and specific implementation, step S122 may include:

[0118] First, the executing entity can read the preset circumferential modulus and the minimum safety factor, and calculate their product. It can be understood that, near the magnetic surface where the minimum safety factor is located, the product of the preset circumferential modulus and the minimum safety factor is approximately equal to the poloidal modulus corresponding to the anti-shear Alfvén eigenmode at that location. This is a rule determined by the resonance constraint condition holding at that physical surface.

[0119] Then, the executing entity can perform a rounding operation on the product of the aforementioned preset circumferential modulus and the minimum safety factor. This rounding operation can specifically employ rounding rules. It is understood that the poloidal modulus, representing the complete number of spatial periods, must be an integer. Simultaneously, the integer closest to the theoretical product value minimizes resonance mismatch error, thereby ensuring the best mode structure match of the anti-shear Alfvén eigenmodes near the magnetic surface of the minimum safety factor.

[0120] Finally, the executing entity can output the rounded integer value as the polar modulus for use in the subsequent construction of the Alfvén continuous spectrum equation.

[0121] In some embodiments, the step of constructing the Alfvén continuum equation characterizing the radial mode structure of the anti-shear Alfvén eigenmodes based on equilibrium reconstruction data, preset circumferential moduli, and poloidal moduli includes:

[0122] Step S131: Determine the radial distribution of parallel wavenumbers based on the safety factor profile, large radius, preset circumferential modulus, and polar modulus; wherein, the parallel wavenumber is the product of the preset circumferential modulus and the safety factor profile minus the polar modulus, and then divided by the product of the safety factor profile and the large radius.

[0123] In this embodiment, the parallel wavenumber is the pivotal quantity connecting the equilibrium profile safety factor and the modulus, and its calculation involves four input quantities: the radial distribution of the safety factor profile, the large radius, the preset circumferential modulus, and the aforementioned polar modulus.

[0124] Specifically, the parallel wavenumber can be calculated independently at each grid point in the radial direction, and its physical meaning is the wavenumber component of the sheared Alfvén wave along the direction of the magnetic field lines.

[0125] In one possible and specific implementation, the process of performing parallel wavenumber calculation may specifically include:

[0126] The execution entity first obtains the values ​​of the safety factor profile at each radial grid point and the value of the large radius from step S11. It then obtains the preset circumferential modulus from step S11 and the polar modulus from step S12. For each grid point on the radial grid, the execution entity reads the safety factor value at that grid point and performs the following calculations.

[0127] The first step is to calculate the numerator. The executing entity can multiply the preset circumferential modulus by the safety factor value at the current grid point to obtain a first intermediate quantity. The polar modulus is then subtracted from the first intermediate quantity to obtain the numerator value. The sign of this numerator value reflects the degree of deviation of the combination of the circumferential and polar moduli from the safety factor at the current radial position.

[0128] The second step is to calculate the denominator. The executing entity can multiply the safety factor value at the current grid point by the large radius to obtain the denominator value. The large radius is a fixed value in a tokamak and does not change with radial position; therefore, the radial variation of the denominator value is determined by the safety factor profile.

[0129] The third step is to calculate the parallel wavenumber. The executing entity divides the numerator value obtained in the first step by the denominator value obtained in the second step to obtain the parallel wavenumber value at that grid point. The above three steps can be repeated for all radial grid points to obtain the complete radial distribution of the parallel wavenumber.

[0130] Step S132: Determine the radial distribution of the poloidal wavenumber based on the poloidal modulus and radial position; wherein the poloidal wavenumber is the poloidal modulus divided by the current radial position.

[0131] In this embodiment, the poloidal wavenumber is the wavenumber component that characterizes the sheared Alfvén wave along the poloidal direction.

[0132] In one possible and specific implementation, performing a specific calculation of the poloidal wavenumber may include:

[0133] The executing entity can obtain the poloidal modulus from step S12. For each grid point on the radial grid, the radial coordinate value of that grid point is read, and the poloidal modulus is divided by that radial coordinate value to obtain the poloidal wave number at that grid point. Since the poloidal modulus is a fixed integer value, the variation of the poloidal wave number along the radial direction is determined by the radial coordinate value; that is, the closer to the core (the smaller the radial position), the larger the poloidal wave number. The closer to the edge (the larger the radial position), the smaller the poloidal wave number.

[0134] The executing entity can repeat the above operation for all radial grid points to obtain the complete radial distribution of the poloidal wavenumber.

[0135] Step S133: Determine the radial distribution of Alfvén velocity based on magnetic field strength and plasma density; wherein, the Alfvén velocity is the square root of the magnetic field strength divided by the product of vacuum permeability, plasma density, and ion mass.

[0136] In this embodiment, the Alfvén velocity is the characteristic propagation velocity of a sheared Alfvén wave in a plasma medium, and the vacuum permeability is a physical constant.

[0137] In one possible and specific implementation, the executing entity can obtain the magnetic field strength and plasma density profile from step S11. For each grid point on the radial grid, the plasma density value at that grid point is read, and the following calculations are performed:

[0138] The first step is to calculate the product term in the denominator.

[0139] The intermediate product can be obtained by multiplying the vacuum permeability, the plasma density at the current grid point, and the ion mass. Here, the vacuum permeability is a fixed physical constant, and the ion mass represents the mass of the dominant ion species in the plasma.

[0140] The second step is to calculate the denominator.

[0141] Take the square root of the intermediate product obtained in the first step to get the denominator. This denominator value characterizes the inertial effect of the plasma medium; that is, the greater the density, the larger the denominator value.

[0142] The third step is to calculate the Alfvén velocity.

[0143] Dividing the magnetic field strength by the denominator obtained in the second step yields the Alfvén velocity value at that grid point. In this embodiment, the magnetic field strength can be treated as a given parameter; therefore, the radial variation of the Alfvén velocity is mainly determined by the relative relationship between the plasma density profile and the magnetic field strength, meaning that the higher the density region, the lower the Alfvén velocity.

[0144] The executing entity can repeat the above three steps for all radial grid points to obtain the complete radial distribution of Alfvén velocity.

[0145] Step S134: Determine the radial distribution of the target resonance coefficient based on the parallel wavenumber distribution, Alfvén velocity distribution, and geodesic frequency; wherein the target resonance coefficient is the square of the current frequency minus the product of the square of the parallel wavenumber distribution and the square of the Alfvén velocity distribution, and then minus the square of the geodesic frequency.

[0146] In one possible and specific implementation, the implementing entity may perform a target resonance coefficient calculation, which may include:

[0147] The execution entity can obtain the radial distribution of the parallel wavenumber from step S131, the radial distribution of the Alfvén velocity from step S133, and further determine or obtain the radial distribution of the geodesic frequency based on the input balance parameters. For each grid point on the radial grid, the execution entity reads the parallel wavenumber value, Alfvén velocity value, and geodesic frequency value at that grid point and performs the following calculations:

[0148] The first step is to calculate the frequency squared term. This can be done by squaring the current frequency. The current frequency is the trial frequency during the iteration process; its real part represents the oscillation frequency of the eigenmode, and its imaginary part represents damping or growth.

[0149] The second step is to calculate the contributions of the parallel wavenumber and the Alfvén velocity. The execution entity can square the parallel wavenumber value at the current grid point to obtain the squared parallel wavenumber. It can also square the Alfvén velocity value at the current grid point to obtain the squared Alfvén velocity. Multiplying the squared parallel wavenumber and the squared Alfvén velocity yields the wavenumber-velocity coupling term. The physical meaning of this term can be the square of the local Alfvén continuum frequency at the current radial position, determined by both the magnetic field geometry and the plasma density.

[0150] The third step is to calculate the squared term of the geodesic acoustic frequency. This can be done by squaring the geodesic acoustic frequency value at the current grid point to obtain the squared geodesic acoustic frequency. The geodesic acoustic frequency is related to the plasma pressure gradient and reflects the correction of the continuous spectrum frequency by the finite pressure effect.

[0151] The fourth step is to synthesize the target resonance coefficient. The executing entity can subtract the wavenumber-velocity coupling term obtained in the second step from the squared frequency value obtained in the first step, and then subtract the squared geodesic frequency value obtained in the third step to obtain the target resonance coefficient value at that grid point. The sign of the target resonance coefficient determines the relative relationship between the current frequency and the local Alfvén continuous spectrum frequency. Specifically, a positive value indicates that the current frequency is higher than the continuous spectrum frequency, a negative value indicates that the current frequency is lower than the continuous spectrum frequency, and a zero value indicates that the two resonate.

[0152] The executing entity can repeat the above four steps for all radial grid points to obtain the complete radial distribution of the target resonance coefficient.

[0153] In this embodiment, the physical condition for the target resonance coefficient to be zero can be: the square of the current frequency is exactly equal to the product of the square of the parallel wavenumber and the square of the Alfvén velocity, plus the sum of the squares of the geodesic frequency. This condition is the Alfvén resonance condition.

[0154] Step S135: Based on the target resonance coefficient and poloidal wavenumber distribution, establish a second-order singular differential equation with electrostatic potential as the solution variable, which serves as the Alfvén continuous spectrum equation.

[0155] In this embodiment, after the parallel wavenumber distribution, poloidal wavenumber distribution, Alfvén velocity distribution and target resonance coefficient distribution are all calculated, the executing entity can substitute the target resonance coefficient distribution, poloidal wavenumber distribution and geometric coupling terms that characterize the average curvature and circumferential coupling effect into the second-order differential equation with electrostatic potential as the solution variable to form the Alfvén continuous spectrum equation.

[0156] Specifically, the equation includes a term that radially differentiates the combination of the target resonance coefficient and the electrostatic potential derivative, as well as an electrostatic potential constraint term that includes a geometric coupling term, a poloidal wavenumber distribution, and a target resonance coefficient distribution. The former term reflects the radial propagation characteristics of the electrostatic potential under radial non-uniformity conditions, while the latter term reflects the constraint effect of geometric effects and poloidal wavenumbers on the electrostatic potential mode structure.

[0157] The executing entity can fill the radial distribution of the target resonance coefficient obtained in step S134, the radial distribution of the poloidal wavenumber obtained in step S132, and the geometric coupling term into the corresponding positions of the equation to form a complete radial control equation.

[0158] Within the radial region where the target resonance coefficient is not zero, the equation is a conventional second-order differential equation. Near the Alfvén resonance point where the target resonance coefficient approaches zero, the equation coefficients exhibit singular behavior, and the solution displays singularity at this point.

[0159] In some embodiments, the step of identifying the location of Alfvén resonance points in the Alfvén continuity spectrum equation includes:

[0160] Step S136: For the current frequency, calculate the real part of the target resonance coefficient point by point along the radial grid, and scan the sign change of the real part along the radial direction.

[0161] In this embodiment, during the aforementioned steps, the executing entity has obtained the complete radial distribution of the target resonance coefficients. Since the Alfvén resonance point corresponds to the radial position where the real part of the target resonance coefficients crosses zero, the executing entity can use the current probe frequency as input to extract the real part values ​​of the target resonance coefficients point by point on the radial grid and scan the radial trend of the variation of the real part values.

[0162] Specifically, the executing entity can sequentially read the target resonance coefficient values ​​at each grid point in ascending order of radial position, and extract the real part from each target resonance coefficient value to form a real part numerical sequence corresponding one-to-one with the radial grid. This real part numerical sequence reflects the radial distribution of the target resonance coefficient under the current frequency conditions. Then, the executing entity can sequentially compare the signs of the real part values ​​at adjacent grid points along the radial direction to determine whether the real part of the target resonance coefficient undergoes a zero-crossing change from positive to negative or from negative to positive within the adjacent grid interval. If such a zero-crossing change occurs, it indicates the existence of an Alfvén resonance point within the corresponding radial interval.

[0163] Step S137: When the product of the real values ​​of two adjacent radial grid points is negative, it is determined that there is an Alfvén resonance point between the two radial grid points. The radial position of the Alfvén resonance point is calculated by linear interpolation, and the positions of the left and right resonance layers are output.

[0164] In this embodiment, when the executing entity detects that the real part values ​​at two adjacent radial grid points have opposite signs, it indicates that there is a zero-crossing point of the real part between the two grid points. Therefore, the executing entity can further calculate the product of the real part values ​​at the two adjacent grid points. When the product is negative, it can be determined that there is an Alfvén resonance point between the two radial grid points.

[0165] After determining the interval, the precise radial position of the Alfvén resonance point can be determined using linear interpolation. Specifically, the radial coordinate values ​​and corresponding real part values ​​of two adjacent radial grid points can be read. Based on the linear relationship between the real part values ​​of the two endpoints and the radial position, the radial position corresponding to the point where the real part is exactly zero can be calculated, and this radial position is determined as the location of the Alfvén resonance point.

[0166] Since the mode structure of the anti-shear Alfvén eigenmodes often corresponds to a resonant layer on each side of the matching point, the execution entity classifies the identified Alfvén resonant points according to their positional relationship relative to the matching point after completing the entire radial interval scan. Alfvén resonant points located to the left of the matching point are output as the left resonant layer position, and Alfvén resonant points located to the right of the matching point are output as the right resonant layer position. These left and right resonant layer positions will serve as input parameters for the subsequent construction of the complex plane integration path.

[0167] Step S138: When the product of the real values ​​on all adjacent radial grid points is non-negative, it is determined that the current frequency is not within the range of the continuous spectrum.

[0168] In this embodiment, after completing a segment-by-segment scan along the entire radial grid, if the executing entity does not detect a negative product of the real values ​​at any pair of adjacent radial grid points, it indicates that the real part of the target resonance coefficient does not have a zero-crossing position in the entire radial range. At this time, it can be determined that the current probe frequency does not correspond to the Alfvén resonance point in the radial calculation domain, that is, the current frequency is not within the continuous spectrum range.

[0169] In this scenario, the executing entity can mark the current frequency as a frequency without a resonance point and determine that the continuous spectrum damping corresponding to the current frequency is zero, thus ending the resonance point localization process for that frequency. During subsequent frequency searches, the executing entity can continue to repeat the aforementioned real part extraction, sign scanning, interval determination, and linear interpolation localization steps for new trial frequencies until a frequency point with both a left and right resonance layer is identified.

[0170] In some embodiments, the step of constructing an integration path around the Alfvén resonance point on the complex plane, analytically extending the Alfvén continuity spectrum equation on the real axis to the complex plane, and performing target integration from the left and right boundaries of the radial computational domain along the integration path towards the matching point to obtain the left-hand and right-hand solutions at the matching point includes:

[0171] Step S141: Determine the detour radius; wherein the detour radius is less than the radial distance between the Alfvén resonance point and the matching point, and the constructed integral path does not touch other Alfvén resonance points.

[0172] In this embodiment, to avoid the singular behavior of the equation coefficients caused by passing through the Alfvén resonance point when directly integrating along the real axis, the executing entity can first determine the detour radius used to construct the complex plane detour path. The detour radius is used to limit the local range within which the integration path deviates from the real axis and bypasses the Alfvén resonance point, and its value simultaneously satisfies two constraints:

[0173] First, the detour radius is smaller than the radial distance between the corresponding Alfvén resonance point and the matching point, so as to ensure that the detour path does not cross the matching point.

[0174] Second, the radius of the circumference should be such that the constructed integration path does not touch other Alfvén resonance points when circumnavigating in the complex plane, so as to avoid new singular points entering the vicinity of the integration path.

[0175] Specifically, the executing entity reads the positions of the left Alfvén resonance point, the right Alfvén resonance point, and the matching point, calculates the radial distances between the left and right Alfvén resonance points and the matching point, and combines this with the interval information between other Alfvén resonance points and the current resonance point to examine each candidate detour radius. When a candidate detour radius simultaneously satisfies the conditions of not crossing the matching point and not touching other Alfvén resonance points, the executing entity determines this candidate detour radius as the detour radius used in the current integration process and uses this detour radius to construct both the left and right integration paths.

[0176] Step S142: For the Alfvén resonance point located to the left of the matching point, construct the left integration path; wherein, integrate along the real axis from the left boundary to the left Alfvén resonance point position minus the detour radius; from this position, enter the complex plane along the first complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the second complex direction coefficient, reaching the left Alfvén resonance point position plus the detour radius; from this position, integrate along the real axis to the matching point.

[0177] In this embodiment, for the Alfvén resonance point located to the left of the matching point, a left-side integration path can be constructed. The left-side integration path is formed by sequentially connecting the first left-side real axis path segment, the left-side complex plane circumduction path segment, and the second left-side real axis path segment.

[0178] Specifically, starting from the left boundary of the radial computational domain, the integral proceeds inward along the real axis until it reaches the left Alfvén resonance point minus the radius of the traversal. This path segment is then defined as the first real axis path segment on the left. Next, the execution entity uses this position as the starting point for traversing the complex plane and enters the complex plane along the first complex direction coefficient. This first complex direction coefficient corresponds to the direction in which the real part increases and the imaginary part decreases. The integration path continues along this direction, deviating from the real axis and entering the region below the complex plane. After entering the complex plane, the integration continues along the preset complex plane path, ensuring that the entire path bypasses the location of the left Alfvén resonance point.

[0179] After completing the circumduction, return to the real axis along the second complex direction coefficient. The second complex direction coefficient corresponds to the direction in which the real part increases and the imaginary part increases. Return to the real axis along this direction and reach the left Alfvén resonance point plus the circumduction radius.

[0180] Finally, the executing entity takes this position as the starting point and continues to integrate inward along the real axis until it reaches the matching point, and determines this path segment as the second real axis path segment on the left.

[0181] With the above construction, the left integration path bypasses the singular position corresponding to the left Alfvén resonance point on the real axis and forms a continuous and smooth bypass channel in the complex plane, so that subsequent numerical integration can continue to advance to the matching point without passing through the singular point.

[0182] Step S143: For the Alfvén resonance point located to the right of the matching point, construct the right integration path; wherein, integrate along the real axis from the right boundary to the position of the right Alfvén resonance point plus the detour radius; from this position, enter the complex plane along the second complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the first complex direction coefficient, reaching the position of the right Alfvén resonance point minus the detour radius; from this position, integrate along the real axis to the matching point.

[0183] In this embodiment, for the Alfvén resonance point located to the right of the matching point, a right-side integration path is constructed. The right-side integration path is formed by sequentially connecting the right-side first real axis path segment, the right-side complex plane bypass path segment, and the right-side second real axis path segment, and its construction method corresponds to the left-side integration path in the spatial direction.

[0184] Step S144: Analytically extend the Alfvén continuous spectrum equation along the left and right integration paths, and integrate along the left and right integration paths toward the matching point with non-zero amplitude and zero derivative as boundary conditions to obtain the left-hand solution and the right-hand solution.

[0185] In this embodiment, after the left and right integration paths are constructed, the executing entity can perform an analytical extension of the Alfvén continuous spectrum equation along the left and right integration paths. Specifically, the executing entity no longer limits the integration variables to real radial coordinates on the real axis, but expands the integration variables to complex position parameters that change continuously along the complex plane path, and calculates the values ​​of the equation coefficients at the corresponding complex positions point by point during the path advancement, so that the Alfvén continuous spectrum equation originally defined on the real axis can be continuously valid along the complex plane path.

[0186] Before performing the target integration, boundary conditions are set at the left and right boundaries of the radial computational domain. These boundary conditions are defined as follows: the electrostatic potential has a preset non-zero amplitude value, and the first derivative of the electrostatic potential is set to 0. Then, the execution entity uses the left and right boundaries as initial integration starting points and performs numerical integration along the left and right integration paths, respectively. During integration, the electrostatic potential and its derivative are updated segment by segment according to the path progression order. Integration is advanced according to the real axis coordinates on the real axis path segments and according to the complex path parameters on the complex plane path segments. The electrostatic potential value and derivative at the end of the previous path segment are used as the initial conditions for the starting point of the next path segment.

[0187] When the integration along the left integration path reaches the matching point, the electrostatic potential value and its first derivative at that location are extracted and determined as the left-hand solution. When the integration along the right integration path reaches the matching point, the electrostatic potential value and its first derivative at that location are extracted and determined as the right-hand solution. These left-hand and right-hand solutions will serve as input data for subsequent logarithmic derivative matching and complex eigenfrequency calculations.

[0188] In some implementations, the step of determining the residual function based on the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point includes:

[0189] Step S151: Calculate the left-hand logarithmic derivative based on the function value of the left-hand solution at the matching point and its first derivative with respect to the radial position; wherein the left-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value.

[0190] In this embodiment, the executing entity reads the function value of the left-hand side solution at the matching point and reads the corresponding first-order derivative value. Then, the first-order derivative value is divided by the function value to obtain the left-hand logarithmic derivative. The left-hand logarithmic derivative characterizes the local rate of change of the left-hand side solution relative to its own amplitude at the matching point, and is used to describe the local morphological characteristics of the left-hand eigenfunction near the matching point.

[0191] Step S152: Calculate the right-hand logarithmic derivative based on the function value of the right-hand solution at the matching point and its first derivative with respect to the radial position; wherein the right-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value.

[0192] In this embodiment, corresponding to the calculation method of the left-hand logarithmic derivative, after completing the target integration along the right-hand integration path, the function value of the right-hand solution and the first-order derivative of this function value with respect to the radial position have been obtained at the matching point. The right-hand logarithmic derivative can then be calculated based on the function value and the first-order derivative value.

[0193] Specifically, the function value of the right-hand side solution at the matching point is read, and the corresponding first-order derivative value is also read. Then, the execution entity divides the first-order derivative value by the function value to obtain the right-hand logarithmic derivative.

[0194] Step S153: Determine the difference between the logarithmic derivative on the left and the logarithmic derivative on the right as the residual function; wherein, the matching condition for the logarithmic derivative is that the residual function takes a value of zero.

[0195] In this embodiment, after calculating the left-hand logarithmic derivative and the right-hand logarithmic derivative, the executing entity further constructs a residual function based on the difference between the two. Specifically, the right-hand logarithmic derivative is subtracted from the left-hand logarithmic derivative, and the resulting difference is determined as the residual function value corresponding to the current frequency.

[0196] Specifically, the residual function is used to measure the consistency of local connectivity between the left-hand and right-hand solutions at the matching point. When the logarithmic derivatives of the left and right sides are exactly equal, the difference is 0, indicating that the local rates of change of the eigenfunctions obtained from the left-hand integration and the right-hand integration are consistent at the matching point. In this case, the left-hand and right-hand solutions satisfy the logarithmic derivative matching condition. Conversely, when the difference is not zero, it indicates that there is still a difference in the local rates of change between the left-hand and right-hand solutions at the matching point, and they have not yet reached a matching state.

[0197] In some implementations, the step of iteratively solving for the complex eigenfrequency when the residual function is zero includes:

[0198] Step S154: Determine the local Alfvén frequency at the matching point based on the value of the parallel wavenumber at the matching point and the value of the Alfvén velocity at the matching point.

[0199] In this embodiment, the executing entity can read the radial position corresponding to the matching point, extract the parallel wave number value at that position from the parallel wavenumber distribution, and extract the Alfvén velocity value at that position from the Alfvén velocity distribution. Then, the parallel wave number value and the Alfvén velocity value are combined to obtain the local Alfvén frequency at the matching point. The local Alfvén frequency characterizes the characteristic frequency scale of the local continuous spectrum near the matching point, and is used to provide a frequency reference for setting the subsequent complex frequency search interval.

[0200] Step S155: Based on the local Alfvén frequency, determine the search interval for the real part of the complex frequency and the search interval for the imaginary part of the complex frequency. Establish a two-dimensional grid within the search interval, calculate the modulus of the residual function at each grid point, and determine the complex frequency corresponding to the minimum value of the modulus as the initial guessed complex frequency.

[0201] In this embodiment, after obtaining the local Alfvén frequency at the matching point, the search interval for the real part of the complex frequency and the search interval for the imaginary part of the complex frequency can be determined based on the local Alfvén frequency. A two-dimensional grid is established within the search interval, and the initial guessed complex frequency is determined by comparing the modulus of the residual function corresponding to each grid point.

[0202] Specifically, the execution entity sets the search range for the real part of the complex frequency centered on the local Alfvén frequency, ensuring that the real part search interval covers the possible oscillation frequency range of the eigenmode. Simultaneously, considering the characteristic that the imaginary part of continuous spectrum damping is negative, the search range for the imaginary part of the complex frequency is set, ensuring that the imaginary part search interval covers the possible variation range of the damping amount. Then, a two-dimensional grid is divided within the real and imaginary part search intervals to obtain multiple candidate complex frequency grid points.

[0203] For each candidate complex frequency grid point, the execution entity uses the complex frequency corresponding to that grid point as the current test frequency and re-executes the aforementioned target resonance coefficient calculation, Alfvén resonance point identification, complex plane integral path construction, target integration on both sides, and residual function calculation processes to obtain the residual function value corresponding to the candidate complex frequency. The execution entity further calculates the modulus of the residual function value to obtain the residual function modulus value corresponding to the candidate complex frequency. After the residual function modulus values ​​of all grid points have been calculated, the residual function modulus values ​​corresponding to all grid points can be compared, and the complex frequency corresponding to the smallest modulus value is determined as the initial guessed complex frequency.

[0204] Through the coarse scanning process described above, the executing entity transforms the originally unknown complex eigenfrequency problem into a grid search problem within a finite search interval, thereby obtaining an initial frequency position that is closer to the true solution, providing a stable initial value for subsequent fine iterations.

[0205] Step S156: Using the initial guessed complex frequency as the initial value, iteratively adjust the complex frequency to make the magnitude of the residual function approach zero until the preset convergence condition is met. The complex frequency that meets the preset convergence condition is determined as the complex eigenfrequency.

[0206] In this embodiment, after obtaining the initial guessed complex frequency, the execution subject uses the initial guessed complex frequency as the initial value to further execute the iterative solution process, continuously adjusting the real and imaginary parts of the complex frequency so that the magnitude of the residual function gradually decreases and approaches zero until the preset convergence condition is met.

[0207] Specifically, the initial guessed complex frequency is first used as the input calculation process for the current iteration frequency. Based on this current iteration frequency, the aforementioned residual function calculation steps are re-executed to obtain the residual function value and its magnitude corresponding to the current iteration frequency. Then, based on the changing trend of the residual function magnitude at the current iteration frequency and adjacent candidate frequencies, the real and imaginary parts of the complex frequency are adjusted so that the new trial frequency moves in the direction of decreasing residual function magnitude. After each frequency update, the execution entity re-executes the aforementioned calculation processes and obtains the residual function magnitude corresponding to the updated frequency.

[0208] More specifically, during the iteration process, the implementing entity continuously monitors two convergence metrics.

[0209] The first convergence criterion is whether the magnitude of the residual function is less than the preset function tolerance.

[0210] The second convergence criterion is whether the change in the complex frequency between two adjacent iterations is less than the preset step size tolerance.

[0211] When both of the above convergence indices are satisfied simultaneously, the current frequency is determined to have reached the preset convergence condition, and this current frequency is designated as the complex intrinsic frequency. If either convergence index is not yet satisfied, the next round of frequency adjustment and residual function calculation continues until the convergence condition is met.

[0212] Finally, the executing entity can output the complex frequency obtained when the preset convergence condition is met as the complex eigenfrequency. The real part of the complex frequency corresponds to the eigenfrequency of the inverse shear Alfvén eigenmode, and the imaginary part corresponds to the continuous spectrum damping.

[0213] In some embodiments, after outputting the imaginary part of the complex eigenfrequency as the continuous spectrum damping rate, the method further includes:

[0214] Step S17: Obtain the driving growth rate of high-energy particles in the tokamak device.

[0215] In this embodiment, after the execution entity completes the solution of the complex eigenfrequency and outputs the continuous spectrum damping rate, the driving growth rate of the high-energy particles in the tokamak device for the antishear Alfvén eigenmode can be further obtained.

[0216] Specifically, the executing entity can retrieve high-energy particle-related information from experimental diagnostic data or preprocessing results. This information may include the distribution intensity, energy level, and interaction parameters between the high-energy particles and plasma waves. Based on this information, the energy input capability of high-energy particles to mode perturbations can be quantified using a pre-defined physical model or empirical calculation module to obtain the corresponding driving growth rate value.

[0217] The driving growth rate characterizes the ability of high-energy particles to provide energy to the anti-shear Alfvén eigenmode; the larger the value, the stronger the instability driving effect on the mode.

[0218] Step S18: Compare the continuous spectrum damping rate with the high-energy particle-driven growth rate. When the high-energy particle-driven growth rate is greater than the absolute value of the continuous spectrum damping rate, the anti-shear Alfvén eigenmode is determined to be unstable.

[0219] In this embodiment, after obtaining the continuous spectrum damping rate and the high-energy particle-driven growth rate, the execution subject can compare the two to determine the stability state of the anti-shear Alfvén eigenmode.

[0220] Specifically, the continuous spectrum damping rate is solved as the imaginary part of the complex eigenfrequency in the aforementioned steps. Its value is negative, corresponding to the damping effect, that is, the energy of the RSAE is resonantly absorbed and irreversibly dissipated. The driving growth rate is positive, corresponding to the driving effect, that is, high-energy particles inject energy into the RSAE.

[0221] More specifically, the absolute value of the continuous spectrum damping rate is first calculated, and the physical strength of the damping is quantified as a positive number so that it is comparable to the driving growth rate on the same numerical scale.

[0222] Then, the driving growth rate can be numerically compared with the absolute value of the continuous spectrum damping rate. When the driving growth rate is greater than the absolute value of the continuous spectrum damping rate, it indicates that the rate at which high-energy particles inject energy into the RSAE exceeds the rate at which the thermal plasma continuous spectrum damping dissipates energy, and the RSAE can be determined to be in an unstable state. When the driving growth rate is less than or equal to the absolute value of the continuous spectrum damping rate, it indicates that the continuous spectrum damping is sufficient to dissipate the energy injected by high-energy particles, and the RSAE can be determined to be in a stable state.

[0223] This determination can be used to inform operational decisions for the tokamak device. When RSAE is determined to be unstable, it is necessary to reduce the auxiliary heating power to decrease the injection of high-energy particles, or adjust the safety factor profile to enhance continuous spectrum damping, thereby reducing the drive growth rate below the damping rate. When RSAE is determined to be stable, the current operating parameters are within a safe operating window.

[0224] like Figure 1 As shown, in a specific implementation scheme, a method for accurately calculating the continuous spectrum damping of the tokamak back-shear Alfvén eigenmode is provided based on the complex path shooting method, which may specifically include:

[0225] S1: Input of equilibrium and profile data.

[0226] Input tokamak equilibrium reconstruction data (such as the gfile file output by the EFIT program) and obtain the following physical quantities along the normalized small radius. Radial distribution:

[0227] Safety factor profile Plasma density ion temperature Electronic temperature Magnetic field strength Large radius Small radius .

[0228] Preprocessing of input data: Cubic spline interpolation is used to map each physical quantity onto a uniform radial grid, with the number of grid points... This ensures sufficient resolution for the localized RSAE model structure. A radial computational domain is defined. ,in , (Covering a local area of ​​RSAE). Select radial matching point. Minimum safety factor Location, determined by one-dimensional search. point.

[0229] Input parameters determined:

[0230]

[0231] S2: Construction of the Alfvén continuous spectrum equation.

[0232] Based on the ideal MHD model, a description of the electrostatic potential is derived. The second-order singular differential equation for a radial module structure. Under the approximation of a circular cross-section and a large ring diameter ratio, the equation can be written in the following standard form:

[0233]

[0234] The key coefficient is:

[0235]

[0236] in, , , , To measure the frequency of ground sound, It incorporates effects from favorable average curvature and toroidal couplings.

[0237] Singularity identification: for the frequency of the current iteration Scan on radial grid The sign change. When adjacent grid points satisfy At that time, the singularity position is calculated by linear interpolation. :

[0238]

[0239] Output the positions of the singularities on the left and right sides. (Left resonant layer) and (Right resonant layer). If no singularity is detected, the current frequency is not within the continuous spectrum range, and the continuous spectrum damping is zero.

[0240] S3: Complex path construction and singularity avoidance integrals.

[0241] To avoid numerical divergence caused by directly crossing singularities on the real axis, a bypass is constructed on the complex plane. The integration path.

[0242] Complex path construction method (taking the left segment as an example, the detour radius is taken as) m):

[0243]

[0244] Right side segment (from) arrive The structure is symmetrical, and the integration directions are opposite. (Circumference radius) Must meet At the same time, it ensures that the path does not touch other singularities, with typical values... m.

[0245] Extension of differential equations in the complex plane: The differential equations on the real axis are analytically extended to the complex plane. The integral on the complex path uses the variable-step Runge-Kutta method (such as MATLAB's ode15s function), with a relative tolerance set to... The absolute tolerance is set to For complex path segments, the integration variable is a dimensionless parameter. The actual radial coordinate is ,in For complex direction coefficients (e.g.) or ).

[0246] S4: Logarithmic derivative matching and eigenvalue solving.

[0247] After determining the complex integration path from step S3, this implementation scheme performs target integration from the left and right boundaries of the radial computational domain toward the matching point:

[0248] 1. Target shooting from the left side: from the left boundary Starting from point S3, integrate along the path constructed to the matching point. The solution on the left side is obtained. and its logarithmic derivative The left boundary condition is set to... , (i.e., the physical boundary condition where the amplitude of the mode far from the resonance region tends to zero).

[0249] 2. Target shooting from the left side: from the left boundary Starting from the point, integrate along the symmetric path to the matching point. The solution on the right side is obtained. and its logarithmic derivative The boundary conditions on the right side are also set to have small non-zero values ​​and a derivative of zero.

[0250] Matching conditions: Define residual function The eigenfrequency of RSAE satisfy When the residual function approaches zero, it indicates that the eigenfunctions on both sides can be smoothly connected at the matching point, and the corresponding frequency is the complex eigenfrequency of RSAE.

[0251] Iterative solution strategy:

[0252] 1. Coarse Scan: Performs a two-dimensional mesh scan over the target region in the complex plane. Typical scan range: Real part imaginary part (Negative values ​​correspond to damping). Mesh density Calculate the value at each grid point. Take the minimum point as the initial guess. .

[0253] 2. Detailed solution: with Using the initial values, the Levenberg-Marquardt algorithm (such as MATLAB's fsolve function) is used iteratively to solve the problem. Iteration termination condition: function tolerance Step size tolerance .generally It can converge in one iteration.

[0254] Output high-precision complex intrinsic frequency: ,in, (Negative value) represents the continuous spectrum damping ratio. Normalized damping ratio. The typical value is approximately This corresponds to the typical theoretical prediction value of RSAE.

[0255] S5: Output the results.

[0256] The output includes: the real part of the RSAE intrinsic frequency. (rad / s, convertible to kHz); Continuous spectrum damping rate (rad / s);

[0257] Normalized frequency and damping ratio ; Eigenfunction Profile and Singularity position , ;

[0258] Convergence history (number of iterations, final residual); all results are saved in a structured format (such as a text file) with metadata annotations (input parameters, code version, timestamp, etc.).

[0259] First, this implementation plan addresses the numerical singularity problem and achieves stable computation.

[0260] Because the complex path integration strategy is used in S3, the differential equation is analytically extended to the complex plane, and the integration path completely bypasses the Alfvén resonance point. The location of the point (e.g., the position of the point), thus completely avoiding singularities in numerical calculations. This allows the scheme to be used even in scenarios with extremely small damping (e.g., the position of the point). This approach can also obtain stable and convergent solutions, while the traditional MHD code often exhibits numerical divergence or non-physical oscillations under the same conditions. Therefore, this implementation scheme can stably solve the radial eigenfunctions of RSAE in the presence of continuum resonance, thus improving the reliability of the continuum damping calculation results.

[0261] Second, this implementation plan improves the extraction accuracy of damping rates in small-scale continuous spectrum.

[0262] Because this implementation scheme uses the logarithmic derivative matching condition instead of global field discretization in S4, it avoids the accumulation of numerical errors caused by insufficient grid resolution in traditional methods. The continuity condition at the matching point directly gives the eigenvalue equation. High-precision complex frequencies can be obtained by iteratively solving the problem using the Levenberg-Marquardt algorithm. Numerical tests show that for typical RSAE parameters ( The computational accuracy of this method can reach [percentage missing]. (Relative residual), while the calculation results of traditional global MHD code at this order of magnitude of damping rate are often overwhelmed by numerical noise, making it difficult to provide reliable quantitative values.

[0263] Third, this implementation plan improves the efficiency of parameter scanning calculation.

[0264] Because this implementation scheme adopts a target-shooting framework, it transforms the two-dimensional complex frequency search problem into a one-dimensional radial integration problem, avoiding the solution of large sparse matrices in the global MHD code. A single eigenvalue solution requires only 15-30 complex path integrals (including coarse scans and fine iterations), and the computation time for a single parameter point in the MATLAB environment is approximately 10-30 seconds. Therefore, when it is necessary to scan parameters such as the minimum safety factor, density profile, temperature profile, and modulus, this invention can significantly reduce the computation cost of a single set of parameters, making it suitable for batch computation and sensitivity analysis.

[0265] Fourth, the calculation results of this implementation plan are highly adaptable to engineering applications.

[0266] Because this implementation scheme directly solves the singular differential equations describing continuous spectrum resonant absorption in S2, the physical process is clear, and the contribution of continuous spectrum damping can be isolated and quantitatively evaluated. By adjusting the coefficients in the equations (such as adding finite Larmor radius corrections, kinetic effects, etc.), this implementation scheme can be easily extended to more complex physical models. Currently, this method has been successfully applied to the calculation of RSAE damping rate for typical anti-shear configurations in the HL-3 device and can be used in advanced operating scenarios for different devices such as DIII-D and ITER. By comparing with the high-energy particle-driven growth rate, this implementation scheme can be used to assess whether RSAE will be excited, thereby assisting in determining heating power, optimizing safety factor profiles, reducing the risk of high-energy particle loss, and supporting the design and optimization of advanced operating modes.

[0267] This implementation provides a high-precision, high-stability, and high-efficiency method for calculating RSAE continuum damping. For a given equilibrium configuration and the most unstable circumferential modulus determined by the high-energy particle distribution... This implementation scheme can directly output the continuous spectrum damping rate. This result can be correlated with the driving growth rate of high-energy particles. Perform quantitative comparisons: when The instability of the time-determined mode provides a quantitative basis for the configuration of auxiliary heating power, assessment of high-energy ion transport, and optimization design of advanced operating modes in fusion devices. Therefore, this implementation scheme is not only applicable to RSAE continuum damping calculations in existing tokamak experimental devices, but also to stability prediction and parameter optimization in advanced operating scenarios of future fusion reactors.

[0268] like Figure 2 As shown, the horizontal axis represents the radial position. The left vertical axis represents frequency. The right vertical axis represents the safety factor q. The figure shows the continuous spectrum gaps, the position of the RSAE eigenfrequency relative to the Alfvén continuous spectrum, and the Alfvén resonance point (the intersection of the continuous spectrum and the eigenfrequency, as shown by the red arrow in the figure).

[0269] like Figure 3 The diagram illustrates the integration path around the Alfonchi point in the complex plane. The horizontal axis represents the real axis. The vertical axis is the imaginary axis. The diagram is labeled: solid axis path (solid green line), complex plane path (blue or magenta / green), and singularity location. (Yellow dot), and the radius of travel. The annotation.

[0270] In one specific embodiment, a method for calculating the damping rate of the RSAE continuum under a typical anti-shear configuration of an HL-3 device (fixed) is provided. ) Implementation plan.

[0271] (a) Input data and calculation parameters.

[0272] This embodiment uses equilibrium data reconstructed from the 550ms experiment of the HL-3 tokamak's 11973rd gun. The main parameters of the HL-3 device are: large radius... m, small radius m, central magnetic field T, Alfvén frequency rad / s.

[0273] Safety factor profile Reconstructed using the EFIT program, it exhibits typical backshear characteristics. In the calculation, Fourth-order polynomial fitting was used:

[0274]

[0275] The fitting coefficient is: ; ; ; ; . matching point m is located at the minimum safety factor Location. Circular Module (The most unstable modulus determined by the distribution of high-energy particles), poloidal modulus (satisfy ).

[0276] Other physical parameters:

[0277] (Dimagnetic drift frequency); (Geometry acoustic frequency); m⁻¹; (Characterizing the adiabatic / fluid response of EP); the radial computational domain is set as follows: m, m.

[0278] (II) Calculation process

[0279] Step S1: Coarse scan to obtain initial guesses.

[0280] A two-dimensional grid is scanned in the complex frequency plane to find the minimum residual value as the initial value for iteration. Scan range:

[0281] 20 grid points

[0282] 20 grid points

[0283] After the scan is complete, the optimal initial value is obtained. .

[0284] Step S2: Automatic identification of resonance points.

[0285] The program automatically scans around the initial frequency determined by the scan. The sign change is used to detect the location of the Alfvén resonance point. Program output:

[0286] The resonance point on the left is located at 0.1962 meters, and the resonance point on the right is located at 0.3340 meters.

[0287] The left and right resonance points are symmetrically distributed at the matching point. The resonant points on both sides of m are spaced approximately 0.1378m apart, which is consistent with the physical expectations of RSAE.

[0288] Step S3: Complex path construction and singularity avoidance integral.

[0289] Taking the left segment as an example, the program constructs the following integration path:

[0290] Real axis path: from m-integral to m;

[0291] Complex path 1: along Directional detour, integration interval , m;

[0292] Complex path 2: along Directional detour, integration interval ;

[0293] Real axis path: from m integration to the matching point m;

[0294] right side segment Starting from m, construct a symmetrical integration path, adjusting the direction accordingly (first along m). , and then along ).

[0295] The integration of the differential equation was performed using the ode15s function in MATLAB, with a relative tolerance set to... The absolute tolerance is set to .

[0296] Step S4: Logarithmic derivative matching and eigenvalue solving

[0297] Calculate the logarithmic derivatives of both sides at the matching point:

[0298]

[0299] residual function The Levenberg-Marquardt algorithm (fsolve) is used iteratively to solve the problem. .

[0300] Iteration termination information:

[0301] The equation solution is complete, and the iterative solver stops. The reason for stopping is that the relative size of the current iteration step has decreased below the step size tolerance, and the magnitude of the residual function vector has approached zero according to the function tolerance, thus satisfying the preset convergence condition.

[0302] This indicates that the solution has converged successfully.

[0303] (III) Calculation Results

[0304] The complex eigenfrequency is obtained by solving. (Normalized frequency). Eigenfunction radial distribution as Figure 4 As shown in the figure, the real part, imaginary part, and modulus of the eigenfunction are plotted. The eigenfunction at the matching point... It reaches its peak value near m and then decays to both sides. At the resonance point ( m and Near m), the eigenfunctions exhibit significant local structural changes. The complex plane integration region (the bypass path interval) is marked with a shaded area in the figure, corresponding to the location of the resonance point. Figure 5 The solution on both sides is shown at the matching point. The matching of the logarithmic derivatives at point, the difference of the logarithmic derivatives at... The magnitude of the result verifies the accuracy of the solution.

[0305] (iv) Results Analysis

[0306] This embodiment successfully verifies the effectiveness of the proposed method. The complex path integration strategy successfully avoids the Alfvén singularity (…). m and m), fsolve converged successfully with no numerical divergence; the eigenfunction profile was smooth with no numerical oscillations near the singularity; the calculated normalized damping ratio Approximately The magnitude matches the theoretical prediction.

[0307] In one specific embodiment, a method is provided. An implementation plan for sensitivity analysis of the effect of parameters on the damping rate of the RSAE continuum spectrum.

[0308] (a) Purpose of Analysis

[0309] These are key parameters in the equation that include good mean curvature and circumferential coupling effects. This embodiment demonstrates how to use the present invention for scanning. The study investigates the influence of changes in equilibrium configuration geometry (such as curvature strength and circumferential coupling strength) on the damping rate of the RSAE continuum.

[0310] (II) Parameter Settings

[0311] Taking Example 1 as the baseline case ( ),exist Scan within the range, step size There are a total of 10 parameter points.

[0312] For each parameter point: circumferential modulus polar modulus (Fixed); match point Resonance points are automatically identified. The calculation process in Example 1 is executed, and the intrinsic frequencies are recorded.

[0313] (III) Calculation Results

[0314] Figure 6 and Figure 7 It showed that the frequency varied The changing trend. Figure 6 real part of frequency Follow It increases monotonically. The continuous spectral frequencies are also plotted in the figure. For reference only. It always lies within the gap in the continuous spectrum, which conforms to the basic characteristics of RSAE. Figure 7 For the imaginary part of the frequency (damping ratio) (Negative value, corresponding to damping) varies with The absolute value increases with curvature. This indicates that a stronger curvature effect enhances continuous spectrum damping, making the mode more stable.

[0315] (iv) Analytical Conclusions

[0316] This embodiment shows that: (Good mean curvature and circumferential coupling effect) have a significant impact on the continuous spectrum damping ratio of RSAE; this method can quantitatively give the damping ratio as a function of... The relationship between these changes provides a quantitative basis for comparing the stability of RSAE under different equilibrium configurations; the stronger the curvature effect ( The larger the value, the stronger the continuous spectrum damping, and the more stable the mode; conversely, a weak curvature configuration is conducive to the excitation of RSAE.

[0317] In one specific embodiment, an implementation scheme for comparative verification with existing MHD code is provided.

[0318] To verify the calculation accuracy of this implementation scheme, the same balance parameters as in Example 1 were selected, and a comparative calculation was performed using traditional MHD code (such as NOVA).

[0319] Comparison results: Traditional MHD code NOVA at the Alfvén resonance point ( m and Significant numerical singularities appear near m), the eigenfunctions cannot smoothly pass through these two positions, and traditional methods cannot effectively calculate the continuous spectrum damping ratio, such as... Figure 8 As shown. In contrast, the method in this implementation successfully avoids singularities by constructing a complex plane integral path. Combined with the target-shooting method framework, it can stably and accurately calculate the damping ratio and exhibits good numerical convergence (e.g., Figure 5 (As shown in the figure). This comparative verification fully demonstrates the significant advantages of this implementation scheme in solving numerical singularity problems and extracting high-precision damping rates.

[0320] According to an embodiment of the present invention, an electronic device is provided. The electronic device in this embodiment may include one or more of the following components: a processor, a network interface, memory, non-volatile memory, and one or more application programs, wherein the one or more application programs may be stored in the non-volatile memory and configured to be executed by one or more processors, and the one or more programs are configured to perform the methods described in the foregoing method embodiments.

[0321] According to embodiments of the present invention, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a computer, causes the computer to perform the method described in any of the above embodiments.

[0322] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for calculating the continuous spectrum damping of the tokamak anti-shear Alfvén eigenmode, characterized in that, include: Tokamak equilibrium reconstruction data and preset circumferential modulus are obtained, wherein the equilibrium reconstruction data includes safety factor profile, plasma density, ion temperature, electron temperature, magnetic field strength, large radius and small radius; The matching point is determined based on the minimum value location of the safety factor profile, and the polar modulus is determined based on the preset circumferential modulus. Based on the equilibrium reconstruction data, the preset circumferential and poloidal moduli, the Alfvén continuous spectrum equation characterizing the radial mode structure of the anti-shear Alfvén intrinsic modes is constructed, and the Alfvén resonance point position in the Alfvén continuous spectrum equation is identified. Construct an integration path around the Alfvén resonance point in the complex plane, analytically extend the Alfvén continuous spectrum equation on the real axis to the complex plane, and perform target integration from the left and right boundaries of the radial computational domain along the integration path to the matching point to obtain the left-end solution and the right-end solution at the matching point. The residual function is determined by the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point, and the complex eigenfrequency when the residual function is zero is solved iteratively. The imaginary part of the complex intrinsic frequency is used as the output of the continuous spectrum damping rate.

2. The method according to claim 1, characterized in that, The steps of determining the matching point based on the minimum position of the safety factor profile and determining the polar modulus based on the preset circumferential modulus include: A one-dimensional search is performed on the safety factor profile along the radial direction to determine the position where the first derivative of the safety factor profile with respect to the radial position is zero, and this position is determined as the matching point. Calculate the product of the preset circumferential modulus and the minimum value of the safety factor profile, and round the product to obtain the polar modulus.

3. The method according to claim 1, characterized in that, The steps of constructing the Alfvén continuous spectrum equation characterizing the radial mode structure of the anti-shear Alfvén intrinsic modes based on equilibrium reconstruction data, preset circumferential moduli, and polematic moduli include: Based on the safety factor profile, large radius, preset circumferential modulus, and polar modulus, the radial distribution of parallel wavenumbers is determined; wherein, the parallel wavenumber is the product of the preset circumferential modulus and the safety factor profile, minus the polar modulus, and then divided by the product of the safety factor profile and the large radius. The distribution of the poloidal wavenumber along the radial direction is determined based on the poloidal modulus and the radial position; wherein the poloidal wavenumber is the poloidal modulus divided by the current radial position; The radial distribution of Alfvén velocity is determined based on magnetic field strength and plasma density; wherein, the Alfvén velocity is the square root of the magnetic field strength divided by the product of vacuum permeability, plasma density, and ion mass. The radial distribution of the target resonance coefficient is determined based on the parallel wavenumber distribution, the Alfvén velocity distribution, and the geodesic frequency; wherein the target resonance coefficient is the square of the current frequency minus the product of the square of the parallel wavenumber distribution and the square of the Alfvén velocity distribution, and then minus the square of the geodesic frequency. Based on the target resonance coefficient and poloidal wavenumber distribution, a second-order singular differential equation with electrostatic potential as the solution variable is established as the Alfvén continuous spectrum equation.

4. The method according to claim 3, characterized in that, The step of identifying the location of Alfvén resonance points in the Alfvén continuity spectrum equation includes: For the current frequency, calculate the real part of the target resonance coefficient point by point along the radial grid, and scan the sign change of the real part along the radial direction; When the product of the real values ​​of two adjacent radial grid points is negative, it is determined that there is an Alfvén resonance point between the two radial grid points. The radial position of the Alfvén resonance point is calculated by linear interpolation, and the positions of the left and right resonance layers are output. When the product of the real values ​​at all adjacent radial grid points is non-negative, the current frequency is determined to be outside the continuous spectrum.

5. The method according to claim 4, characterized in that, The steps of constructing an integration path around the Alfvén resonance point in the complex plane, analytically extending the Alfvén continuity spectrum equation on the real axis to the complex plane, and performing target integration along the integration path from the left and right boundaries of the radial computational domain towards the matching point to obtain the left-hand and right-hand solutions at the matching point include: Determine the detour radius; wherein the detour radius is less than the radial distance between the Alfvén resonance point and the matching point, and such that the constructed integral path does not touch other Alfvén resonance points; For the Alfvén resonance point located to the left of the matching point, construct the left integration path; wherein, integrate along the real axis from the left boundary to the left Alfvén resonance point position minus the radius of the circumference; from this position, enter the complex plane along the first complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the second complex direction coefficient, reaching the left Alfvén resonance point position plus the radius of the circumference; from this position, integrate along the real axis to the matching point; For the Alfvén resonance point located to the right of the matching point, construct a right-side integration path; wherein, integrate along the real axis from the right boundary to the position of the right Alfvén resonance point plus the detour radius; from this position, enter the complex plane along the second complex direction coefficient, integrate along the path in the complex plane and return to the real axis along the first complex direction coefficient, reaching the position of the right Alfvén resonance point minus the detour radius; from this position, integrate along the real axis to the matching point; The Alfvén continuous spectrum equation is analytically extended along the left and right integration paths, and non-zero amplitude and zero derivative are used as boundary conditions. The equations are then integrated towards the matching point along the left and right integration paths, respectively, to obtain the left-hand and right-hand solutions.

6. The method according to claim 5, characterized in that, The step of determining the residual function based on the matching condition of the logarithmic derivatives of the left-hand and right-hand solutions at the matching point includes: The left-hand logarithmic derivative is calculated based on the function value of the left-hand solution at the matching point and its first derivative with respect to the radial position; wherein, the left-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value; The right-hand logarithmic derivative is calculated based on the function value of the right-hand solution at the matching point and its first derivative with respect to the radial position; wherein, the right-hand logarithmic derivative is the ratio of the first derivative with respect to the radial position to the function value; The difference between the logarithmic derivative on the left and the logarithmic derivative on the right is defined as the residual function; wherein the matching condition for the logarithmic derivative is that the residual function takes the value of zero.

7. The method according to claim 6, characterized in that, The step of iteratively solving for the complex eigenfrequency when the residual function is zero includes: The local Alfvén frequency at the matching point is determined based on the values ​​of the parallel wavenumber and the Alfvén velocity at the matching point. Based on the local Alfvén frequency, the search intervals for the real part and the imaginary part of the complex frequency are determined. A two-dimensional grid is established within the search intervals, and the modulus of the residual function at each grid point is calculated. The complex frequency corresponding to the minimum value of the modulus is determined as the initial guessed complex frequency. Using the initial guessed complex frequency as the initial value, the complex frequency is iteratively adjusted so that the magnitude of the residual function approaches zero until the preset convergence condition is met. The complex frequency that meets the preset convergence condition is determined as the complex eigenfrequency.

8. The method according to claim 1, characterized in that, After outputting the imaginary part of the complex eigenfrequency as the continuous spectrum damping rate, the method further includes: To obtain the driving growth rate of high-energy particles in a tokamak device; The continuous spectrum damping rate is compared with the high-energy particle-driven growth rate. When the high-energy particle-driven growth rate is greater than the absolute value of the continuous spectrum damping rate, the anti-shear Alfvén eigenmode is determined to be unstable.

9. An electronic device, characterized in that, include: A memory, and one or more processors communicatively connected to the memory; The memory stores instructions that can be executed by the one or more processors to cause the one or more processors to implement the method as described in any one of claims 1 to 8.

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

Images