Fast solving and reduced-order simulation method of GIS substation TGPR considering frequency-dependent characteristics

Abstract

The application discloses a GIS substation TGPR fast solving and order reduction simulation method considering frequency variation characteristics, and comprises the following steps: S1, constructing a multi-physical field full-order initial model considering high-resistivity soil frequency variation characteristics; S2, performing adaptive mesh division based on an electromagnetic gradient posterior error indicator to generate optimized large-scale full-order system matrices C and G; S3, performing frequency variation parameter dimension reduction processing based on wide-frequency vector matching and recursive convolution; S4, performing system-level model order reduction based on a multi-point wide-frequency Krylov subspace; S5, performing super-fast time domain solving of transient ground potential rise and inverse projection reconstruction to restore the real transient ground potential rise waveform in three-dimensional space; the method cooperates adaptive mesh, Krylov subspace order reduction and recursive convolution, shortens the calculation time by several orders of magnitude under the premise of ensuring the three-dimensional full-wave calculation precision, improves the calculation speed of high-frequency transient solving, optimizes the calculation resource allocation, and completes the high-precision fast evaluation of the TGPR at a very low time cost.

Description

Technical Field

[0001] This invention relates to the field of substation simulation calculation technology, and in particular to a method for rapid solution and reduced-order simulation of TGPR for GIS substations that takes into account frequency-varying characteristics. Background Technology

[0002] In GIS substations, rapid transient overvoltage (VFTO) is an overvoltage phenomenon caused by the operation of disconnecting switches in the GIS. It poses a threat to the insulation safety of GIS equipment, and long-term operation accumulation accelerates insulation aging. Its radiated field may interfere with secondary equipment through electromagnetic coupling. Transient ground potential rise (TGPR) refers to a phenomenon in power systems or lightning protection systems where, due to transient processes such as lightning strikes, short-circuit faults, or switching operations, the potential of the grounding grid or grounding device rises significantly relative to the distant "zero potential" ground for a very short time. It can trigger ground potential backflash, damaging electronic equipment; generate step voltage and touch voltage, threatening personal safety; and interfere with secondary systems (such as protection and communication equipment), affecting power system stability. Therefore, in actual GIS substation construction, it is crucial to perform VFTO and TGPR simulation calculations in areas with high soil resistivity. However, existing simulation calculations mainly suffer from the following three significant defects:

[0003] 1. The "curse of dimensionality" results in an extremely large computational burden for high-frequency full-wave simulations. The transient signals generated by the operation of GIS disconnect switches contain high-frequency components reaching tens or even hundreds of megahertz (MHz). To accurately capture the propagation of high-frequency electromagnetic waves and the "skin effect" of metallic conductors, existing simulation methods must perform extremely high-density spatial meshing on the GIS shell and grounding grid. This causes the degrees of freedom (DoF) of the three-dimensional electromagnetic field solution matrix to easily reach millions or even tens of millions. Traditional full-scale solutions face a severe "curse of dimensionality," with single simulations often taking several days and even leading to computer memory overflow, making it impossible to meet the needs of rapid comparison of multiple solutions in engineering design.

[0004] 2. Static / manual mesh generation mechanisms are inefficient and wasteful of resources. Existing simulation software often uses globally uniform meshes or manually refined local meshes based on engineer experience. During transient waveform propagation, the region of concentrated electromagnetic energy changes dynamically over time. Static meshes cannot simultaneously meet the "accuracy requirements of high-frequency skin regions" and the "computational efficiency of low-frequency regions," wasting significant computational resources on non-critical regions with minimal electromagnetic field gradient changes.

[0005] 3. The "history convolution" of frequency-varying characteristics in time-domain integration is extremely time-consuming. The resistivity and dielectric constant of high-resistivity soils, as well as the internal impedance of grounding conductors, all exhibit strong "frequency-varying nonlinearity." Current techniques for considering these frequency-varying parameters in time-domain transient simulations typically require extremely complex inverse Fourier transforms or massive full-history time-domain convolution integrals. As the computational step size increases, the storage and retrieval of historical convolution data grows exponentially, severely slowing down the solution speed. Summary of the Invention

[0006] To overcome the above problems, the purpose of this invention is to provide a fast solution and order reduction simulation method for TGPR of GIS substations that takes into account frequency-varying characteristics. This method, through the synergy of adaptive mesh, Krylov subspace order reduction and recursive convolution, reduces the calculation time by several orders of magnitude while ensuring the accuracy of three-dimensional full-wave calculation, improves the calculation speed of high-frequency transient solution, optimizes the allocation of computing resources, and completes high-precision and fast evaluation of TGPR with extremely low time cost.

[0007] The technical solution adopted in this invention is:

[0008] A rapid solution and reduced-order simulation method for TGPR of GIS substations considering frequency-varying characteristics includes the following steps:

[0009] S1: Construct a multi-physics full-order initial model that considers the frequency-varying properties of high resistivity soil.

[0010] S2: For the full-order initial model in S1, adaptive mesh partitioning is performed based on the electromagnetic gradient posterior error indicator to generate the optimized large-scale full-order system matrices C and G.

[0011] S3: Based on wideband vector matching and recursive convolution, the frequency-varying parameter dimensionality reduction process transforms the time-consuming historical integral terms introduced by soil frequency variation in the C and G matrices of the large-scale full-order system matrix generated in S2 into local first-order recursive state variables.

[0012] S4: Based on the system-level model order reduction of multi-point broadband Krylov subspace, the large-scale matrices C and G containing recursive variables are compressed into micro matrices Cr and Gr of several hundred orders by using the orthogonal projection matrix V.

[0013] S5: Ultrafast time-domain solution and inverse projection reconstruction of transient ground potential rise, using a micromatrix. The solution is obtained by rapid stepping, and then multiplied by V to restore the true transient ground potential rise waveform in three-dimensional space.

[0014] As a further description of the present invention, the specific process of S1 is as follows:

[0015] S11: Based on the actual physical dimensions of the target GIS substation, establish a three-dimensional geometric topology model of the GIS substation; the components of the three-dimensional geometric topology model include the GIS shell, internal central conductive rod, SF6 insulating gas layer, insulating basin, grounding down conductor, and layered non-uniform high-resistivity soil.

[0016] S12: Establish a broadband model of frequency-varying media based on the Debye or Cole-Cole broadband dispersion model, and define the complex permittivity of high resistivity soil as a function of frequency.

[0017] S13: Apply a high-frequency transient excitation source to treat the arc breakdown between the contacts of the GIS disconnector as a time-varying nonlinear resistance model, and inject transient voltage or current excitation with a rise time of nanoseconds and a main frequency component of tens of megahertz into the central conductive rod.

[0018] As a further description of the present invention, when constructing the three-dimensional geometric topology model in S11, a fully matched layer or a second-order Mur absorbing boundary condition is set at the outer truncation boundary of the computational domain to simulate the extension of the earth to infinity and absorb the transient electromagnetic waves radiated outward.

[0019] As a further description of the present invention, the complex permittivity defined as a function of frequency in S12 is specifically in the following form:

[0020] .

[0021] in: It is the optical frequency dielectric constant; Let be the polarization amplitude of the k-th polarization relaxation process; Let k be the relaxation time corresponding to the kth relaxation process; The DC conductivity of the soil; It is the vacuum permittivity; The number of relaxation processes.

[0022] As a further description of the present invention, the specific process of S2 is as follows:

[0023] S21: Coarse trial calculation. Perform a full-order electromagnetic field trial calculation on the initial coarse grid for one excitation cycle to obtain the trial solution. and .

[0024] S22: Error assessment, calculate the posterior error indication value of the Maxwell equation residuals for each tetrahedral element K. .

[0025] S23: Adaptive adjustment yields a spatiotemporally optimized discrete spatial grid topology, and large-scale state-space stiffness / mass matrices C and G of millions of orders, but whose quantity is already optimized, derived from this grid.

[0026] As a further description of the present invention, the error indication value in S22 The calculation formula is:

[0027] .

[0028] in: Let K be the local posterior error indicator value; The characteristic size of unit K; and These are the numerical solutions of the electric and magnetic fields obtained by finite element time-domain calculations under the initial sparse mesh. Permeability; The curl of the electric field; This is the time partial derivative of the magnetic field; Denotes the L2 norm on cell K; Let K be any face on the boundary of element K; For this element Feature dimensions; Represents the set of all boundary surfaces of element K; This represents the jump in the tangential electric field at the common interface e between adjacent grids; This represents the L2 norm on interface e.

[0029] As a further description of the present invention, the specific process of adaptive adjustment in S23 is as follows:

[0030] S231: Set global error threshold , The value range is set to 0.01 to 0.1 according to the required calculation precision.

[0031] S232: Within each adaptive cycle, calculate the posterior error indication value for all grid cells. And calculate the global maximum error. and average error .

[0032] S233: For The element triggers the local subdivision and refinement of the Delaunay mesh, connecting the midpoints of each edge of element K to subdivide a tetrahedral element into eight sub-tetrahedral elements.

[0033] for The element triggers a mesh coarsening operation, which inversely transforms adjacent element groups that meet the conditions back to their parent element before refinement, and deletes redundant nodes and elements.

[0034] As a further description of the present invention, the specific process of S3 is as follows:

[0035] S31: Using the vector matching method, poles and residues are extracted, and frequency domain data of high-resistivity soil frequency-varying polarization parameters and skin impedance of metallic conductors are collected. A vector matching algorithm based on Sanathanan-Koerner iteration is then used to analyze the complex frequency domain response. The approximation is an Nth-order rational fraction, and the calculation formula is:

[0036] .

[0037] In the formula: It is a complex frequency variable; Let m be the m-th complex pole; For corresponding to the poles The residues; d and h are real coefficients.

[0038] S32: Based on the spatial mapping of recursive convolution states, rational fractions are transformed into temporal recursive convolution state space equations, and global history integrals are transformed into first-order recursive equations for local state variables. The discrete-time form of the first-order recursive equations is as follows:

[0039] .

[0040] in, Let be the time-domain local state variable corresponding to the m-th pole; This represents the discrete time step for time-domain simulation. For exponential recursion factors; This represents the contribution of the electric field intensity to the state variables within the current time step.

[0041] As a further description of the present invention, the specific process of S4 is as follows:

[0042] S41: Extraction of the full-order differential equation system. The full-order mass matrix C and full-order stiffness matrix G of the system are extracted from the large-scale state-space matrix equation of the three-dimensional full-wave electromagnetic field system. The large-scale state-space matrix equation is expressed as follows.

[0043] .

[0044] in, For unknowns at the nodes, Let be the spatial excitation vector.

[0045] S42: Select multiple key frequency points covering the main frequency range of the VFTO signal. At each frequency point Calculate the transition matrix at this point. and initial vector The calculation formula is:

[0046] .

[0047] .

[0048] S43: Using the modified Arnoldi orthogonalization iterative algorithm, generate the orthogonal projection basis matrix for multi-frequency joint matching. ,in In the millions, It has several hundred levels.

[0049] S44: By applying a consistent Galerkin projection mapping to the high-dimensional original system using the orthogonal projection matrix V, a compact reduced-order macromodel is obtained. The calculation formula is as follows:

[0050] .

[0051] As a further description of the present invention, the specific process of S5 is as follows:

[0052] S51: Embed the recursive convolution state equation from S3 into the reduced-order equation generated in S4, as follows:

[0053] .

[0054] in, For recursive convolution state variables, This represents the contribution to the reduced-order system after projection.

[0055] S52: The unconditionally stable Newmark-β method is used to solve the discrete-time step problem, obtaining the time series solution in the reduced-order space. .

[0056] S53: Utilizing the orthogonal projection matrix in S4 Inverse spatial dilation mapping is performed to losslessly restore low-dimensional data to the real electromagnetic fields on millions of grid nodes in three-dimensional physical space. Transient values, the mapping process is as follows:

[0057] .

[0058] S54: Extract the potential distribution of each observation point along the GIS shell and the high resistivity grounding network, and output high-precision transient ground potential rise and transient shell potential time-domain waveform diagram and spatial electromagnetic energy cloud map.

[0059] The beneficial effects of this invention are:

[0060] This invention presents a rapid solution and order-reduction simulation method for TGPR (Transient Response Rate) of GIS substations that considers frequency-varying characteristics. This method employs a multi-point Krylov subspace projection algorithm to transform the high-dimensional full-wave space equation matrix into a low-dimensional compact matrix. While ensuring the accuracy of broadband transient response calculation (error controlled within a specified threshold), it significantly reduces the system's degrees of freedom. It overcomes the limitations of traditional physical geometric grids on computational scale, significantly improves the computational speed of high-frequency transient solutions, and achieves a system-level improvement in the efficiency of complex three-dimensional simulations.

[0061] This invention presents a rapid solution and reduced-order simulation method for TGPR (Time-Domain Electromagnetic Ratio) in GIS substations, taking into account frequency-varying characteristics. This method addresses the computational waste caused by traditional global static mesh generation and the massive historical data storage burden generated by the frequency-varying nonlinear characteristics of high-resistivity soils during time-domain integration. It employs a spatial adaptive mesh generation technique based on a posteriori error indicator, combined with vector matching and recursive convolution algorithms to process frequency-varying parameters. This achieves on-demand dynamic densification of the spatial mesh and effective removal of the historical data burden from the time-domain integration. It optimizes computational resource allocation in both spatial and temporal dimensions, avoids computational redundancy, and further improves the efficiency of time-domain electromagnetic simulation under complex nonlinear media conditions.

[0062] This invention presents a rapid solution and order-reduction simulation method for TGPR (Transient Potential Rise) in GIS substations, taking into account frequency-varying characteristics. Addressing the limitation of existing electromagnetic transient simulation methods that rely on equivalent circuits to accurately characterize high-frequency and complex three-dimensional structural parameters, this method employs an order-reduction algorithm to extract and construct an equivalent "transient mathematical macro-model" from a high-precision three-dimensional full-physics model. This macro-model boasts computational speed equivalent to circuit-level node scale while preserving the physical field accuracy of three-dimensional full-wave electromagnetic simulation. This enables engineers to perform high-precision and rapid evaluation of transient potential rise (TGPR) at extremely low time costs when comparing and optimizing GIS grounding system design schemes in areas with high soil resistivity. Attached Figure Description

[0063] Figure 1 This is a flowchart of the rapid solution and order-reduced simulation method for TGPR of GIS substations that takes into account frequency-varying characteristics, as proposed in this invention.

[0064] Figure 2 The curve of frequency-variable grounding impedance versus frequency is shown in Example 7 of the proposed method for rapid solution and reduced-order simulation of TGPR in GIS substations that takes into account frequency-variable characteristics.

[0065] Figure 3 This is the DoF path reduction method in Implementation Example 7 of the GIS substation TGPR fast solution and order reduction simulation method considering frequency variation characteristics proposed in this invention.

[0066] Figure 4This is a comparison diagram of TGPR waveforms in Implementation Example 7 of the TGPR rapid solution and order reduction simulation method for GIS substations that takes into account frequency variation characteristics proposed in this invention.

[0067] Figure 5 This is a comparison of calculation times in Implementation Example 7 of the TGPR fast solution and order reduction simulation method for GIS substations that takes into account frequency variation characteristics proposed in this invention. Detailed Implementation

[0068] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the protection scope of the present invention.

[0069] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0070] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0071] This invention is described in detail with reference to the schematic diagrams. When detailing the embodiments of this invention, for ease of explanation, the cross-sectional views illustrating the device structure may be partially enlarged, not adhering to the usual scale. Furthermore, the schematic diagrams are merely examples and should not be construed as limiting the scope of protection of this invention. In actual fabrication, the three-dimensional spatial dimensions of length, width, and depth should be included.

[0072] Furthermore, in the description of this invention, it should be noted that the terms "upper," "lower," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. These terms are used solely for the convenience of describing the invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. In addition, the terms "first," "second," or "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0073] Unless otherwise explicitly specified and limited, the terms "installation," "connection," and "joining" in this invention should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; similarly, they can refer to mechanical connections, electrical connections, or direct connections, or indirect connections through an intermediate medium, or internal connections between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0074] like Figure 1 As shown, it illustrates a specific embodiment of the present invention:

[0075] Example 1:

[0076] A rapid solution and order-reduced simulation method for TGPR of GIS substations considering frequency-varying characteristics includes the following steps:

[0077] S1: Construct a multi-physics full-order initial model that takes into account the frequency-varying properties of high-resistivity soils; the purpose of this step is to establish a baseline physical and geometric framework for the propagation of high-frequency electromagnetic waves in complex GIS structures and non-uniform soil media.

[0078] S2: For the full-order initial model in S1, adaptive mesh partitioning is performed based on the electromagnetic gradient posterior error indicator to generate the optimized large-scale full-order system matrices C and G; this step is to overcome the "curse of dimensionality" of the global mesh caused by the extremely short wavelength of high frequency.

[0079] S3: Dimensionality reduction of frequency-varying parameters based on wideband vector matching and recursive convolution. For the large-scale full-order system matrix generated in S2, the time-consuming historical integral terms introduced by soil frequency variation in the C and G matrices are transformed into local first-order recursive state variables. The purpose of this step is to solve the problem of memory overflow and extremely low computational efficiency caused by full-history convolution when dealing with frequency-varying media in traditional time-domain full-wave simulation.

[0080] S4: System-level model order reduction based on multi-point broadband Krylov subspaces. Using the orthogonal projection matrix V, the large-scale matrices C and G containing recursive variables are compressed into miniature matrices of several hundred orders. and .

[0081] S5: Ultrafast time-domain solution and inverse projection reconstruction of transient ground potential rise, using a micromatrix. The solution is obtained by rapid stepping, and then multiplied by V to restore the actual transient ground potential rise waveform in three-dimensional space. This step uses a macro model with extremely low dimensions to complete the rapid time-domain stepping calculation and restores the data to the physical space to output the final result.

[0082] In this embodiment, the simulation method transforms the high-dimensional full-wave space equation matrix into a low-dimensional compact matrix using the multi-point Krylov subspace projection algorithm. This significantly reduces the system's degrees of freedom while maintaining the accuracy of broadband transient response calculations (error controlled within a specified threshold). It overcomes the limitations of traditional physical geometric grids on computational scale, significantly improving the computational speed of high-frequency transient solutions and achieving a system-level improvement in the efficiency of complex 3D simulations. Employing a spatial adaptive mesh generation technique based on a posteriori error indicators, combined with vector matching and recursive convolution algorithms to handle frequency-varying parameters, it achieves on-demand dynamic densification of the spatial mesh and effective removal of the burden of historical data integration in the time dimension. It optimizes computational resource allocation in both spatial and temporal dimensions, avoiding computational redundancy and further improving the efficiency of time-domain electromagnetic simulations under complex nonlinear media conditions. Through a reduction-order algorithm, it utilizes an extremely low-dimensional macro-model to complete ultra-fast time-domain step calculations, preserving the physical field accuracy of the 3D full-wave electromagnetic simulation, and achieving high-precision and rapid evaluation of transient potential rise (TGPR) at extremely low time cost.

[0083] Example 2:

[0084] Specifically, the process of S1 is as follows:

[0085] S11: Based on the actual physical dimensions of the target GIS substation, establish a three-dimensional geometric topology model of the GIS substation; the components of the three-dimensional geometric topology model include the GIS shell, internal central conductive rod, SF6 insulating gas layer, insulating basin, grounding down conductor, and layered non-uniform high-resistivity soil.

[0086] In this embodiment, when establishing the three-dimensional geometric model of the GIS substation, it is necessary to set truncation boundary conditions (ABC).

[0087] This condition is set to simulate the earth extending to infinity and absorbing outwardly radiated transient electromagnetic waves. At the truncated boundary of the computational domain, a perfectly matched layer (PML) or a second-order Mur absorption boundary condition is introduced. The PML layer ensures that electromagnetic waves achieve reflection-free attenuation across the entire frequency band and at arbitrary incident angles by constructing anisotropic artificial complex permittivity and complex permeability tensors.

[0088] S12: Establish a broadband model for frequency-varying media based on the Debye or Cole-Cole broadband dispersion model, defining the complex permittivity of high-resistivity soil as a function of frequency; specifically:

[0089] .

[0090] in: It is the optical frequency dielectric constant, and its value ranges from 2 to 10 depending on the soil type; The polarization amplitude of the k-th polarization relaxation process represents the contribution intensity of each relaxation process in the dielectric spectrum. The relaxation time, in seconds, corresponds to the k-th relaxation process and is used to describe how quickly the polarization process responds to frequency changes. The DC conductivity of the soil is expressed in Siemens units per meter, and the value ranges from 0.0001 to 0.1 S / m depending on the soil resistivity. Let be the vacuum permittivity, and take the value of . ; The number of relaxation processes is 1 to 5.

[0091] In this embodiment, the model can accurately characterize the dispersion characteristics of high resistivity soils in the 100kHz to 100MHz frequency band, where the dielectric constant and conductivity change strongly with frequency, thus providing an accurate frequency domain data basis for subsequent dimensionality reduction processing of frequency-varying parameters.

[0092] S13: Apply a high-frequency transient excitation source to simulate the fast transient overvoltage (VFTO) generated by the operation of the GIS disconnector (DS). The arc breakdown between the contacts of the GIS disconnector is equivalent to a time-varying nonlinear resistance model, such as the Toepler arc model or the Cassie-Mayr dynamic arc model. Inject transient voltage or current excitation with a rise time in the nanosecond range and a main frequency component as high as tens of megahertz into the central conductive rod.

[0093] In this embodiment, three-dimensional geometric modeling determines the boundary conditions and topological structure of the spatial electromagnetic field, while frequency-varying medium modeling assigns dynamically changing physical properties to specific regions (such as soil regions or conductor surfaces) within this structure. In the frequency domain, the full-order initial model is characterized by the three-dimensional full-wave Maxwell's equations that take into account frequency-varying constitutive relations.

[0094] .

[0095] in, The vector of spatial electric field intensity; Permeability; The VFTO excitation current source is injected externally; is the complex permittivity of the soil region in the geometric model.

[0096] Example 3:

[0097] Specifically, the process of S2 is as follows:

[0098] S21: Coarse trial calculation. Perform a full-order electromagnetic field trial calculation on the initial coarse grid for one excitation cycle to obtain the trial solution. and .

[0099] S22: Error assessment, calculate the posterior error indication value of the Maxwell equation residuals for each tetrahedral element K. Error indication value The calculation formula is:

[0100] .

[0101] in: Let be the local posterior error indicator value of element K, which is a dimensionless positive real number; The characteristic dimension of element K is defined as the diameter of the inscribed sphere or the length of the longest side of the element. and These are the numerical solutions of the electric and magnetic fields obtained by finite element time-domain calculations under an initial sparse grid. For non-ferromagnetic materials, the permeability is taken as the vacuum permeability. ; The curl of the electric field characterizes the vortex intensity of the electromagnetic field within the unit cell. This is the time partial derivative of the magnetic field; The L2 norm on unit K is represented by the square root of the square integral. Let K be any face on the boundary of element K; For this element Feature dimensions; Represents the set of all boundary surfaces of element K; This represents the jump in the tangential electric field at the common interface e between adjacent grids, i.e., the difference in the tangential electric field calculated from two adjacent elements. Let be the unit normal vector of interface e; This represents the L2 norm on interface e.

[0102] In this embodiment, the physical meaning of the posterior error indicator is that the first term measures the degree to which the electromagnetic field inside the unit satisfies Maxwell's equations, and the second term measures the continuity of the tangential electric field at the unit boundary; when A large value indicates a large electromagnetic field gradient or significant skin effect near the cell, requiring mesh refinement. This indicator can accurately quantify the energy concentration region and skin depth range during the propagation of high-frequency transient electromagnetic waves, providing a rigorous mathematical basis for adaptive mesh generation.

[0103] S23: Adaptive adjustment yields a spatiotemporally optimized discrete spatial grid topology, and large-scale state-space stiffness / mass matrices C and G of millions of orders, but with optimized quantity, derived from this grid. The specific process of adaptive adjustment is as follows:

[0104] S231: Set global error threshold , The value range is set to 0.01 to 0.1 according to the required calculation precision.

[0105] S232: Within each adaptive cycle, calculate the posterior error indication value for all grid cells. And calculate the global maximum error. and average error .

[0106] S233: For The element triggers the local subdivision and refinement of the Delaunay mesh, connecting the midpoints of each edge of element K to subdivide a tetrahedral element into eight sub-tetrahedral elements.

[0107] for The element triggers a mesh coarsening operation, which inversely transforms adjacent element groups that meet the conditions back to the parent element before refinement, and deletes redundant nodes and elements.

[0108] In this embodiment, the above-mentioned densification and coarsening operations are performed simultaneously within an adaptive cycle, and the time-domain trial solution on the new grid is recalculated after completion. This process is iterated until the global error converges to the preset tolerance range. Through this adaptive mechanism, areas with highly concentrated high-frequency transient currents—including the GIS shell connection flange, the shallow 0.5-meter range of high-resistivity soil, the area around the grounding grid intersection node, and the conductor surface with significant skin effect—are automatically densified to a feature size of 1 / 3 to 1 / 5 of the skin depth. For deep soil areas with minimal electric field gradient changes or areas far from the electromagnetic wave injection point, the grid is automatically merged to a larger feature size. Ultimately, the degrees of freedom of the full-order model are controlled to the absolute lower limit while meeting the high-frequency calculation accuracy requirements, reducing the number of grid cells by 80% to 95% compared to a globally uniform grid.

[0109] Example 4:

[0110] Specifically, the process of S3 is as follows:

[0111] S31: Using the vector matching method, poles and residues are extracted, and frequency domain data of high-resistivity soil frequency-varying polarization parameters and skin impedance of metallic conductors are collected. A vector matching algorithm based on Sanathanan-Koerner iteration is then used to analyze the complex frequency domain response. The approximation is an Nth-order rational fraction, and the calculation formula is:

[0112] .

[0113] in: It is a complex frequency variable; Let be the m-th complex pole, and it must be constrained by an algorithm to lie in the left half of the complex plane, i.e. This is to ensure the system causality and absolute stability of the established temporal recursive convolution model; For corresponding to the poles The residues can be real numbers or pairs of conjugate complex numbers; d and h are real coefficients, used to characterize the constant term and the linear term proportional to the frequency response, respectively; the order N of the rational fraction is determined according to the fitting accuracy requirements, and usually ranges from 10 to 50; the vector matching algorithm gradually corrects the pole positions through the Sanathanan-Koerner iterative process, so that the rational fraction best approximates the original frequency domain sampling data in the least squares sense, thereby providing an accurate and compact exponential function expansion form for subsequent recursive convolution.

[0114] S32: Based on the spatial mapping of recursive convolution states, the rational fraction is transformed into a time-domain recursive convolution state-space equation, and the global history integral is transformed into a first-order recursive equation for local state variables, so that the medium polarization response at the current time depends only on the local state variables of the previous time step; the discrete-time form of the first-order recursive equation is:

[0115] .

[0116] in, Let be the time-domain local state variable corresponding to the m-th pole, whose dimension is the same as that of the electric field vector; This is the discrete time step for time-domain simulation, and its value range is determined based on the highest frequency component, typically from 0.1 nanoseconds to 1 nanosecond. It is an exponential recursive factor that can be pre-calculated and stored during the simulation process; This represents the contribution of the electric field intensity to the state variables within the current time step, and can be approximated using trapezoidal integrals or linear interpolation. The core feature of this recursive convolution algorithm is:

[0117] Current moment Local state variables Dependent only on local state variables of the previous time t and the electric field strength within this time step This eliminates the need to store historical electric field data across the entire time axis; thus reducing the memory overhead of full-history convolution integrals required for processing frequency-varying media in traditional time-domain full-wave simulations. Reduce to ,in This represents the total number of time steps, thus completely eliminating the "curse of historical convolution" in the processing of frequency-varying parameters.

[0118] Example 5:

[0119] Specifically, the process of S4 is as follows:

[0120] S41: Extraction of the full-order differential equations: The full-order mass matrix C and full-order stiffness matrix G of the system are extracted from the large-scale state-space matrix equations characterizing the three-dimensional full-wave electromagnetic field system. The large-scale state-space matrix equations are expressed as follows:

[0121] .

[0122] in, For unknowns at nodes, Let be the spatial excitation vector.

[0123] S42: Select multiple key frequency points covering the main frequency range of the VFTO signal. At each frequency point Calculate the transition matrix at this point. and initial vector The calculation formula is:

[0124] .

[0125] .

[0126] S43: Using the modified Arnoldi orthogonalization iterative algorithm, generate the orthogonal projection basis matrix for multi-frequency joint matching. The transfer matrix at all frequencies and initial vector Calculate the matrix as input. The orthogonal basis vectors are calculated, and the basis vectors calculated for all frequency points are horizontally concatenated to generate a multi-frequency joint matching orthogonal projection matrix, where... In the millions, It has several hundred levels.

[0127] In this embodiment, the modified Arnoldi algorithm performs the following operations in each iteration step:

[0128] Calculate new basis vectors, perform Gram-Schmidt orthogonalization on all previously generated basis vectors, and remove linearly dependent components; the iteration terminates when: a preset order reduction q is reached, or the contribution of the newly added basis vectors to the system response is less than a set threshold (e.g., ...). The final projection basis matrix V has orthogonality, i.e. Furthermore, the subspace spanned by its columns can simultaneously approximate the transfer function behavior of the original system near multiple frequency points; the typical value of the reduced order q is 50 to 300, which is much smaller than the dimension n of the full-order model, thus laying the foundation for subsequent Galerkin projection.

[0129] S44: By applying a consistent Galerkin projection mapping to the high-dimensional original system using the orthogonal projection matrix V, a compact reduced-order macromodel is obtained. The calculation formula is as follows:

[0130] .

[0131] In this embodiment, for the reduced-order state variables and output variables ,satisfy in, These are the mass matrix and stiffness matrix after order reduction, respectively. This is the reduced-order input matrix; the dimension n of the original system is usually in the range of... The order of magnitude is reduced, while after projection, the dimension q of the reduced-order macromodel is compressed to 50 to 300; the projection process preserves the key dynamic characteristics of the original system: for selected multiple expansion frequencies The transfer function of the reduced-order system is consistent with that of the original system in the moment-matching sense, i.e., the former... The matrices are equal within the numerical error range; in this way, the large-scale three-dimensional full-wave electromagnetic field model, which was originally impossible to solve in engineering time, is transformed into a compact, reduced-order macro-model that can be solved in milliseconds on a regular computer; this reduced-order macro-model is mathematically rigorously equivalent to the input-output behavior of the original full-order system over a wide bandwidth, thus reducing the computational complexity from [previous method] to [new method] while ensuring the accuracy of high-frequency transient response calculation. Reduce to This achieved an acceleration of more than three orders of magnitude.

[0132] Example 6:

[0133] Specifically, the process of S5 is as follows:

[0134] S51: Embed the recursive convolution state equation in S3 into the reduced-order equation generated in S4 to form the following augmented system:

[0135] .

[0136] in, For recursive convolution state variables, This represents the contribution to the reduced-order system after projection.

[0137] S52: The unconditionally stable Newmark-β method is used to solve the discrete-time step problem, obtaining the time series solution in the reduced-order space. .

[0138] In this embodiment, the parameters of the Newmark-β method are set as follows: and These two parameters guarantee that the algorithm has second-order accuracy and unconditional numerical stability, that is, regardless of the time step. Regardless of the selection method, the numerical solution will not diverge; at each discrete time step The solution is obtained by following the prediction-correction process:

[0139] The first step is prediction. Initial values ​​of the state variables at time 1.

[0140] The second step is to update each local state variable using the recursive formula in S4. .

[0141] The third step is to solve the reduced-order linear equations to obtain the corrected system of equations. .

[0142] The fourth step is to update the velocity and acceleration terms for the next calculation.

[0143] In this embodiment, due to the reduced order matrix The scale is only a few hundred orders, and the size of the coefficient matrix of the linear equation system to be solved at each time step is . Its computational complexity for inversion or LU decomposition is O(n). Typically, the time taken is less than 1 millisecond; the entire time-domain simulation process usually requires tens of thousands to hundreds of thousands of time steps, with the total time taken controlled to tens of seconds to several minutes; in comparison, when directly solving a full-order system, each time step requires... Solving the matrix ( Each step takes several seconds or more, and the total time can take several days or even weeks. Therefore, by combining order reduction and recursive convolution, this invention achieves a computational speedup of more than three orders of magnitude.

[0144] S53: Utilizing the orthogonal projection matrix in S4 Inverse spatial dilation mapping is performed to losslessly restore low-dimensional data to the real electromagnetic fields on millions of grid nodes in three-dimensional physical space. Transient values, the mapping process is as follows:

[0145] .

[0146] S54: Extract the potential distribution of each observation point along the GIS shell and the high resistivity grounding network, and output high-precision transient ground potential rise and transient shell potential time-domain waveform diagram and spatial electromagnetic energy cloud map.

[0147] In this embodiment, the mapping process will have a dimension of The low-dimensional state sequence is precisely dilated back to dimension . The higher-dimensional physical space, in which , Since the projection basis matrix V is orthogonal and the subspace spanned by its columns contains the main dynamic characteristic modes of the original system over a wide bandwidth, the above dilation mapping can losslessly reconstruct the transient electromagnetic field values ​​at all grid nodes in the three-dimensional physical space within a specified error tolerance (usually a relative error of less than 1%); specifically including:

[0148] Restore the electric field vector at each finite element node The three components and the magnetic field vector The three components, and the Poynting vector, current density distribution, and energy density distribution calculated from them; based on this, the potential distribution along the GIS shell is further extracted, and the transient ground potential rise of each observation point relative to the infinitely distant reference ground is calculated. Simultaneously, the potential distribution on the surface of the GIS shell is extracted, and the transient shell potential is output. Time-domain waveforms; the final output includes: waveforms of TGPR and TEV at each observation point over time, a color cloud map of potential distribution on the entire grounding network, and a cloud map of spatial electromagnetic energy density distribution on a specified profile; the above output results are exported in standard data formats (such as VTK, MAT or CSV) and can be directly used for subsequent insulation coordination assessment, grounding system optimization design and determination of personal safety protection distance.

[0149] Example 7:

[0150] To verify the actual effect and speedup of the method of the present invention, this embodiment takes a 252kV GIS substation in Northwest China as an example and constructs a complete set of full-order and reduced-order comparison tests to quantify the speedup and accuracy loss of steps S1-S5 of the present invention.

[0151] I. Test Conditions

[0152] GIS substation main parameters:

[0153] (1) Voltage level: 252 kV, reference phase voltage .

[0154] (2) Three-phase common busbar length: 10m, outer shell diameter: 0.5m, central conductive rod diameter: 0.1m.

[0155] (3) Insulating gas (0.4MPa), 8 pot-type insulators.

[0156] (4) Four grounding leads, each 5m long, are connected from the outer shell to the horizontal grounding grid.

[0157] (5) Horizontal grounding grid: 5m×5m grid, 4 vertical grounding electrodes with a length of 5m, buried at a depth of 0.6m.

[0158] Soil parameters (determined from field survey data):

[0159] (1) Three-layer soil structure: (Top layer 0-8m) (Middle level 8-28m) (Below 28 m depth).

[0160] (2) Frequency-varying model: Debye broadband dispersion, , , , (correspond ).

[0161] VFTO Incentives:

[0162] (1) Amplitude (Typical level of contact breakdown in disconnector switches).

[0163] (2) Rising time 1ns, decay constant 50ns, corresponding to main frequency 1-100MHz.

[0164] (3) Excitation source model: Toepler time-varying nonlinear resistive arc model.

[0165] Computational domain and boundary:

[0166] (1) Total calculation area: 100m×100m×30m, with a 5m thick PML layer on the periphery.

[0167] (2) Initial grid: Tetrahedron with an average side length of 2 m (based on the highest frequency of 100 MHz wavelength) (Division).

[0168] (3) Boundary conditions: PML anisotropic complex medium absorption + second-order Mur first-order ABC.

[0169] II. Construction of the initial model of the full order in S1: According to the specific method of step S1 of this invention, a three-dimensional geometric topology model is constructed. The components include: GIS shell (aluminum), central conductive rod (copper), SF6 gas, basin insulator (epoxy resin), grounding down conductor, horizontal grounding grid (galvanized flat steel), and layered non-uniform soil.

[0170] The frequency-varying medium modeling adopts the Debye wideband dispersion model: where K=1 (a single-order Debye can characterize the dominant relaxation process in the 100kHz-100MHz frequency band).

[0171] In the frequency domain, the full-order initial model is characterized by Maxwell's equations. Frequency-variable grounding impedance. The calculation results are as follows Figure 2 As shown.

[0172] 3. S2 Adaptive Mesh Generation: Adaptive mesh generation is performed according to the specific method of S2 described in this invention.

[0173] (1) S21 coarse calculation: Perform full-order electromagnetic field calculation for one excitation cycle (200 ns) on the initial 800,000 DoF grid; the coarse calculation adopts the low-order approximation of the S1 frequency-varying model (only K=1 order Debye is taken) to balance optimization quality and computational efficiency.

[0174] (2) S22 Error Assessment: For each tetrahedral element K, calculate the posterior error indication value of Maxwell's equations.

[0175] (3) S23 Adaptive adjustment: Set global error threshold .

[0176] ① Encryption principles: The cells trigger Delaunay local subdivision (midpoint connection method, tetrahedron subdivided into 8 sub-tetrahedrons), which densifies the cells by 1.5 times to 1.2 million DoF.

[0177] ② Coarsening Criteria: The unit triggers coarsening (removing redundant nodes), coarsening deep soil units far from the electromagnetic energy concentration zone to 1 / 5 of their original size.

[0178] ③ Convergence: After 5 adaptive iterations, the global error converges to Below, the final DoF is 240,000 (70% less than the initial uniform grid and 80% less than the peak of 1.2 million).

[0179] The key areas for adaptive mesh densification are: the GIS shell connection flange and the area within 0.5m of the shallow soil layer (skin depth). ), around the grounding grid intersection nodes, and on the conductor surface (skin depth 1 / 5). DoF reduction path as follows: Figure 3 As shown.

[0180] IV. S3 Frequency-Variation Parameter Dimensionality Reduction (Vector Matching + Recursive Convolution): Following the specific method of S3 described in this invention, vector matching is performed on the aforementioned frequency-variable soil and skin conductor frequency domain data. The Sanathanan-Koerner iterative algorithm is used to progressively correct the pole positions, forcibly constraining all poles. To ensure causality and stability, in this embodiment, N=30 (range 10-50), with a fitting accuracy of 0.5%.

[0181] The rational fraction expansion obtained by vector matching is transformed into a time-domain recursive convolution state-space equation. This recursive convolution reduces the memory overhead of the full-history convolution integral in traditional time-domain full-wave simulation for handling frequency-varying media from O(Nt) to O(1).

[0182] V. S4 Krylov subspace system-level order reduction: According to the specific method of S4 described in this invention, the aforementioned full-order system is subjected to Krylov order reduction.

[0183] (1) Extraction of the full-order S41 system, where It is a sparse matrix. This is the activation vector.

[0184] (2) S42 multi-frequency selection, covering the VFTO main frequency band 1-100MHz, selects 7 frequency points at logarithmic intervals. , Calculated at each frequency point.

[0185] (3) The S43 modified Arnoldi iteration adopts the implicit restart Arnoldi method (IRAM), and the iteration termination condition is the contribution of the newly added basis vector to the system response. ; final generation orthogonal projection basis matrix of order (range 50-300) ,satisfy .

[0186] (4) S44 Galerkin projection, reduced order matrix As a dense matrix, LU decomposition can be completed in milliseconds. The reduced-order system is identical to the original system in the moment-matching sense for the first 2q moments, meaning that the transfer function is strictly equal around the 7 expanded frequency points.

[0187] VI. S5 Time-Domain Solution and Inverse Projection: Following the specific method of S5 described in this invention, the aforementioned recursive convolution state equation is embedded into a reduced-order equation. The parameters in this embodiment are:

[0188] (1) Time step: .

[0189] (2) Total simulation time: .

[0190] (3) Total number of time steps: .

[0191] (4) The coefficient matrix size for solving the linear equation system in one step is 200×200 (dense), and the LU decomposition takes 0.5ms.

[0192] In the inverse projection reconstruction, the 200-dimensional reduced-order solution is losslessly restored to the 240,000-dimensional full-order solution. Since V satisfies orthogonality and its column-spread subspace contains the main dynamic modes of the original system over a wide bandwidth, the relative error of the inverse projection is controlled within 0.5%. Observation points are extracted: along the GIS outer shell, at the four corners of the horizontal grounding grid, and at the end of the vertical grounding electrode. The output includes TGPR and TEV time-domain waveforms and a spatial electromagnetic energy cloud map, such as... Figure 4 As shown.

[0193] VII. Performance and Accuracy Comparison.

[0194] (a) Comparison of calculation time

[0195] On the same hardware platform, the computation time of the three methods is as follows: Figure 5 As shown in the figure. Detailed performance parameters are shown in Table 1.

[0196] Table 1 Comparison of performance parameters of the three methods

[0197]

[0198] Key acceleration ratio:

[0199] (1) The relative full-order FEM-EMTP of this method is 5000×(69.4h→50s).

[0200] (2) This method is relative to CDEGS HIFREQ: 216 × (3.0h→50s).

[0201] (3) DoF reduction: 1200000→200=6000×.

[0202] (4) Memory reduction: approximately 30×.

[0203] (II) Accuracy Comparison

[0204] Table 2 Comparison of TGPR calculation accuracy among three methods

[0205]

[0206] As can be seen from Table 2:

[0207] (1) The peak value of TGPR of this method is 43.98 kV, the peak error of the full-order reference is only 0.06%, the RMS error is 0.07%, and the waveform correlation coefficient is 0.998, which meets the engineering accuracy requirements.

[0208] (2) CDEGS HIFREQ has a typical error of 15-20% in the MHz range and cannot capture the nanosecond leading edge.

[0209] (III) Error Budget Decomposition

[0210] Table 3. Error budget decomposition of this method (independent contribution of each stage)

[0211]

[0212] In summary, in this embodiment, regarding speed, the four steps of adaptive meshing (S2, DoF reduced by 80%), recursive convolution (S3, avoiding full-history convolution), Krylov order reduction (S4, DoF reduced from 240,000 to 200), and Newmark stepping (S4 / S5, time step reduced from 0.1ns to 0.5ns) work together to reduce the engineering time of the 252kV GIS TGPR simulation from 69.4h (2.9 days) to 50s (0.8min), a speedup of 5000×. Compared with CDEGS HIFREQ, the speedup is approximately 216×.

[0213] In terms of accuracy, the relative error of the TGPR peak value of this method is 0.06%, the RMS error is 0.07%, and the waveform correlation coefficient is 0.998.

[0214] In terms of frequency-varying modeling, S3 maps the Debye frequency domain response to the time domain, avoiding the 30%+ error of the static conductivity method on the MHz scale. S4's multi-point wideband Krylov order reduction achieves moment matching with the original system at 7 expanded frequency points, ensuring faithful wideband transient response.

[0215] These three features combined transform GIS TGPR simulation from an offline task that takes several days into a multi-scheme comparison tool that can be run on-site, making it directly applicable for the optimized design of grounding systems in areas with high soil resistivity.

[0216] The preferred embodiments of the present invention have been described in detail above with reference to the accompanying drawings. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention.

[0217] Many other changes and modifications can be made without departing from the concept and scope of this invention. It should be understood that this invention is not limited to the specific embodiments, and the scope of this invention is defined by the appended claims.

Claims

1. A rapid solution and order-reduced simulation method for TGPR of GIS substations considering frequency-varying characteristics, characterized in that, Includes the following steps: S1: Construct a multi-physics full-order initial model that takes into account the frequency-varying properties of high resistivity soil; S2: For the full-order initial model in S1, adaptive mesh partitioning is performed based on the electromagnetic gradient posterior error indicator to generate the optimized large-scale full-order system matrices C and G. S3: Based on wideband vector matching and recursive convolution, the frequency-varying parameter dimensionality reduction process transforms the time-consuming historical integral terms introduced by soil frequency variation in the C and G matrices of the large-scale full-order system matrix generated in S2 into local first-order recursive state variables. S4: Based on the system-level model order reduction of multi-point broadband Krylov subspace, the large-scale matrices C and G containing recursive variables are compressed into micro matrices Cr and Gr of several hundred orders by using the orthogonal projection matrix V. S5: Ultrafast time-domain solution and inverse projection reconstruction of transient ground potential rise, using a micromatrix. The solution is obtained by rapid stepping, and then multiplied by V to restore the true transient ground potential rise waveform in three-dimensional space.

2. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 1, characterized in that, The specific process of S1 is as follows: S11: Based on the actual physical dimensions of the target GIS substation, establish a three-dimensional geometric topology model of the GIS substation; the components of the three-dimensional geometric topology model include the GIS shell, internal central conductive rod, SF6 insulating gas layer, insulating basin, grounding down conductor, and layered non-uniform high-resistivity soil. S12: Establish a broadband model of frequency-varying media based on the Debye or Cole-Cole broadband dispersion model, and define the complex permittivity of high resistivity soil as a function of frequency. S13: Apply a high-frequency transient excitation source to treat the arc breakdown between the contacts of the GIS disconnector as a time-varying nonlinear resistance model, and inject transient voltage or current excitation with a rise time of nanoseconds and a main frequency component of tens of megahertz into the central conductive rod.

3. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 2, characterized in that, When constructing the three-dimensional geometric topology model in S11, a fully matched layer or a second-order Mur absorbing boundary condition is set at the outer truncation boundary of the computational domain to simulate the extension of the earth to infinity and absorb the transient electromagnetic waves radiated outward.

4. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 2, characterized in that, In S12, the complex permittivity is defined as a function of frequency in the following specific form: ； in: It is the optical frequency dielectric constant; Let be the polarization amplitude of the k-th polarization relaxation process; Let k be the relaxation time corresponding to the kth relaxation process; The DC conductivity of the soil; It is the vacuum permittivity; The number of relaxation processes.

5. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 1, characterized in that, The specific process of S2 is as follows: S21: Coarse trial calculation. Perform a full-order electromagnetic field trial calculation on the initial coarse grid for one excitation cycle to obtain the trial solution. and ; S22: Error assessment, calculate the posterior error indication value of the Maxwell equation residuals for each tetrahedral element K. ; S23: Adaptive adjustment yields a spatiotemporally optimized discrete spatial grid topology, and large-scale state-space stiffness / mass matrices C and G of millions of orders, but whose quantity is already optimized, derived from this grid.

6. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 5, is characterized in that, The error indication value in S22 The calculation formula is: ； in: Let K be the local posterior error indicator value; The characteristic size of cell K; and These are the numerical solutions of the electric and magnetic fields obtained by finite element time-domain calculations under an initial sparse grid. Permeability; The curl of the electric field; This is the time partial derivative of the magnetic field; Denotes the L2 norm on cell K; Let K be any face on the boundary of element K; For this element Feature dimensions; Represents the set of all boundary surfaces of element K; This represents the jump in the tangential electric field at the common interface e between adjacent grids; This represents the L2 norm on interface e.

7. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 5, characterized in that, The specific process of adaptive adjustment in S23 is as follows: S231: Set global error threshold , The value range is set to 0.01 to 0.1 according to the required calculation precision; S232: Within each adaptive cycle, calculate the posterior error indication value for all grid cells. And calculate the global maximum error. and average error ; S233: For The element triggers the local subdivision and refinement of the Delaunay mesh, connecting the midpoints of each edge of element K, thus subdividing a tetrahedral element into eight sub-tetrahedral elements. for The element triggers a mesh coarsening operation, which inversely transforms adjacent element groups that meet the conditions back to the parent element before refinement, and deletes redundant nodes and elements.

8. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 1, characterized in that, The specific process of S3 is as follows: S31: Using the vector matching method, poles and residues are extracted, and frequency domain data of high-resistivity soil frequency-varying polarization parameters and skin impedance of metallic conductors are collected. A vector matching algorithm based on Sanathanan-Koerner iteration is then used to analyze the complex frequency domain response. The approximation is an Nth-order rational fraction, and the calculation formula is: ； In the formula: It is a complex frequency variable; Let m be the m-th complex pole; For corresponding to the poles The residues; d and h are real coefficients; S32: Based on the spatial mapping of recursive convolution states, rational fractions are transformed into temporal recursive convolution state space equations, and global history integrals are transformed into first-order recursive equations for local state variables. The discrete-time form of the first-order recursive equations is as follows: ； in, Let be the time-domain local state variable corresponding to the m-th pole; This represents the discrete time step for time-domain simulation. For exponential recursion factors; This represents the contribution of the electric field intensity to the state variables within the current time step.

9. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 1, characterized in that, The specific process of S4 is as follows: S41: Extraction of the full-order differential equations. The full-order mass matrix C and full-order stiffness matrix G of the system are extracted from the large-scale state-space matrix equations of the three-dimensional full-wave electromagnetic field system. The large-scale state-space matrix equations are expressed as follows: ； in, For unknowns at the nodes, Let it be the spatial excitation vector; S42: Select multiple key frequency points covering the main frequency range of the VFTO signal. At each frequency point Calculate the transition matrix at this point. and initial vector The calculation formula is: ； ； S43: Using the modified Arnoldi orthogonalization iterative algorithm, generate the orthogonal projection basis matrix for multi-frequency joint matching. ,in In the millions, It has several hundred levels; S44: By applying a consistent Galerkin projection mapping to the high-dimensional original system using the orthogonal projection matrix V, a compact reduced-order macromodel is obtained. The calculation formula is as follows: 。 10. The method for rapid solution and order-reduced simulation of TGPR in GIS substations considering frequency-varying characteristics according to claim 1, characterized in that, The specific process of S5 is as follows: S51: Embed the recursive convolution state equation from S3 into the reduced-order equation generated in S4, as follows: ； in, For recursive convolution state variables, This represents the contribution to the reduced-order system after projection. S52: The unconditionally stable Newmark-β method is used to solve the discrete-time step problem, obtaining the time series solution in the reduced-order space. ; S53: Utilizing the orthogonal projection matrix in S4 Inverse spatial dilation mapping is performed to losslessly restore low-dimensional data to the real electromagnetic fields on millions of grid nodes in three-dimensional physical space. Transient values, the mapping process is as follows: ； S54: Extract the potential distribution of each observation point along the GIS shell and the high resistivity grounding network, and output high-precision transient ground potential rise and transient shell potential time-domain waveform diagram and spatial electromagnetic energy cloud map.

Images