Direct current GIL temperature rise distribution rapid calculation method combining transient state and steady state

Abstract

The invention relates to the technical field of numerical simulation of power equipment, and discloses a transient and steady state combined direct current GIL temperature rise distribution rapid calculation method, which comprises the following steps: firstly, constructing a multi-physics field coupling model and obtaining an initial Joule heat source density; then, an artificial viscosity model controlled by a flow field characteristic monitoring criterion is introduced, self-adaptive transient pre-calculation based on viscosity residual collaborative relaxation is executed, and a flow field skeleton is established while numerical oscillation is inhibited; then space mapping transfer based on the momentum energy ratio is executed, space reconstruction and momentum compensation are conducted on the transient velocity field through the constructed correction operator, and an equivalent steady-state solving initial value is generated; and finally, starting a steady-state solver by using the initial value, and executing closed-loop iteration including conductor resistivity updating and closed cavity air pressure correction in a real physical environment until convergence, thereby effectively solving the problem that high Rayleigh number flow field calculation is difficult to converge and consumes long time, and considering both calculation stability and efficiency.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for power equipment, specifically a rapid calculation method for the temperature rise distribution of DC GIL combined with metastable state. Background Technology

[0002] Gas-insulated direct current (GIL) transmission lines are high-capacity power transmission equipment. Their internal temperature rise distribution affects the insulation strength of the insulating gas, the electrical performance of the insulators, and the sealing performance of the metal casing. Obtaining the steady-state temperature rise distribution of a DC GIL under rated current-carrying conditions is of great significance for equipment structural optimization and safe operation assessment.

[0003] Currently, multiphysics numerical simulation is the primary method for analyzing the temperature rise characteristics of gas-insulated liquids (GILs). GIL devices are large in size, and the insulating gas filling them has a low dynamic viscosity. Driven by Joule heating from the high-voltage conductor, natural convection in the gas region is at a high Rayleigh number. Existing numerical methods face the challenge of balancing computational stability and efficiency when dealing with such high Rayleigh number natural convection coupled heat transfer problems.

[0004] If a steady-state solver is used directly, the buoyancy driving force is greater than the fluid viscous drag under high Rayleigh number conditions, resulting in strong flow field nonlinearity. Starting the calculation with zero initial velocity can easily cause numerical oscillations, leading to solver non-convergence. To address the convergence problem, existing techniques often use transient solvers to simulate the process of the system from cold start to thermal equilibrium. Due to the large heat capacity of the metallic conductor and the insulating shell, the system's thermal time constant is long, while the time step of transient calculations is less constrained by flow field stability. This difference in time scales means that transient calculations require a large number of time step iterations, resulting in excessive computational time, consuming computational resources, and failing to meet the needs of rapid solutions in engineering analysis.

[0005] Furthermore, existing computational models simplify the physical field coupling mechanism. Conventional methods often neglect the characteristic that conductor resistivity increases with temperature, employing a heat source model with a fixed power density; or they ignore the characteristic that the GIL is a constant-volume closed container, and its internal air pressure changes with the average temperature, using only a constant pressure boundary. During the operation of a DC GIL, the electric field, thermal field, flow field, and air pressure field interact with each other. These simplifications lead to deviations between the calculated temperature rise distribution and the actual operating conditions, affecting the accuracy of the simulation results. Summary of the Invention

[0006] To address the challenges of convergence difficulties in steady-state calculations for multiphysics coupled calculations of DC gas-insulated transmission lines (GILs) due to high Rayleigh number natural convection effects, and the excessively long time consumption of purely transient calculations, which fail to balance computational stability and time efficiency, this invention provides a rapid calculation method for the temperature rise distribution of DC GILs that combines metastable and steady-state approaches. This method introduces an artificial viscosity model for transient pre-calculation to establish the flow field topology, and then uses a spatial mapping operator to transfer the results to the steady-state solver, achieving rapid and stable solutions under complex flow fields.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A rapid calculation method for the temperature rise distribution of DC GIL combined with metastable state includes the following steps:

[0009] Step S1: Construct a three-dimensional multiphysics coupling model of the DC gas-insulated transmission line and perform mesh discretization processing, set the basic material parameters, solve the steady-state current conservation equation based on the ambient temperature conditions, and obtain the initial Joule heat source density.

[0010] Step S2: Load the initial Joule heat source density obtained in step S1 into the transient solver, introduce the artificial viscosity model controlled by the flow field characteristic monitoring criterion, perform adaptive transient pre-calculation based on viscosity residual collaborative relaxation, establish the flow field skeleton while suppressing numerical oscillations, and output transient terminal field data including velocity field, temperature field and pressure field.

[0011] Step S3: Read the transient end field data output in step S2, perform spatial mapping transfer based on momentum-energy ratio, and use the correction operator constructed by viscosity ratio and flow index to perform spatial reconstruction and momentum compensation on the transient velocity field in the transient end field data to generate physically equivalent steady-state solution initial values.

[0012] Step S4: Start the steady-state solver using the initial steady-state solution value generated in step S3. Perform calculations in a real physical environment with the artificial viscosity model removed, and execute closed-loop iterations including conductor resistivity updates and closed cavity gas property pressure corrections until the calculation converges. Output the final corrected steady-state temperature field and flow field distribution.

[0013] Preferably, in step S1, the construction process of the three-dimensional multiphysics coupling model includes: for the fluid-structure interaction interface, generating a multi-layer prismatic mesh in the region of the fluid domain near the solid wall, analyzing the velocity gradient and temperature gradient near the wall to meet the requirements for solving the dimensionless wall distance; and using an unstructured tetrahedral mesh in the mainstream region of the insulating gas. The governing equations of the three-dimensional multiphysics coupling model include: fluid mass conservation equation and fluid momentum conservation equation for the insulating gas region, and multiphysics energy conservation equation for the entire domain; the fluid momentum conservation equation includes pressure gradient force, viscous stress term, and buoyancy source term.

[0014] Preferably, in step S1, the specific process of obtaining the initial Joule heat source density is as follows: assuming that the system as a whole is under a uniform initial ambient temperature, the potential distribution inside the high-voltage conductor is calculated using the steady-state current conservation equation; based on the calculated potential distribution, the initial Joule heat source density at the initial moment is calculated by the relationship between the electric field strength and the conductor conductivity at the initial ambient temperature.

[0015] Preferably, in step S2, the artificial viscosity model is constructed using an effective dynamic viscosity model. The effective dynamic viscosity model defines the effective dynamic viscosity at each time step as the product of the true physical viscosity of the insulating gas and an artificial gain term. The artificial gain term consists of a unit value and a dimensionless artificial viscosity relaxation factor. The initial value of the dimensionless artificial viscosity relaxation factor is set to a preset high value, so that the fluid at the initial moment exhibits the target viscosity characteristics and the Reynolds number is in a laminar state.

[0016] Preferably, in step S2, the specific process of performing the adaptive transient pre-calculation based on viscosity residual collaborative relaxation includes: calculating the convective heat transfer intensity on the surface of the high-voltage conductor in real time using the average Nusselt number calculation formula; calculating the filtered Nusselt number time-varying rate characteristic value using the convective heat transfer intensity and the Nusselt number derivative smoothing characteristic value formula; substituting the Nusselt number time-varying rate characteristic value into the artificial viscosity relaxation factor evolution equation, and updating the dimensionless artificial viscosity relaxation factor at each time step using the flow field stability gating function. The flow field stability gating function dynamically outputs a first state value to forcibly lock the current artificial viscosity, or outputs a second state value to initiate the viscosity decay process, based on the relationship between the Nusselt number time-varying rate characteristic value and a preset flow field stability criterion threshold. When the Nusselt number time-varying rate characteristic value is less than the preset Nusselt number change rate tolerance, the transient pre-calculation is determined to meet the termination condition. The flow field stability criterion threshold is a preset value that characterizes the tolerance for the rate of change of the Nusselt number; the Nusselt number change rate tolerance is a preset tolerance value used to determine whether the flow field evolution tends to a quasi-steady state.

[0017] Preferably, in step S3, the specific process of constructing the correction operator includes: using the ratio between the effective dynamic viscosity and the true physical viscosity at the end of the transient calculation as the viscosity ratio, and combining it with the flow index, determining the theoretical amplification factor of the velocity field using the theoretical correction operator calculation formula; comparing the theoretical amplification factor with a preset numerical stability upper limit threshold of the correction operator using the applied correction operator limiting formula, and taking the smaller value to obtain the final applied correction operator applied to the velocity field; defining the applied correction operator as the correction operator. Wherein, the numerical stability upper limit threshold of the correction operator is a preset value determined based on the Courant number constraint relationship between the mesh size and the time step.

[0018] Preferably, in step S3, the specific process of generating the initial value for steady-state solution based on the spatial mapping transfer of momentum-energy ratio includes: substituting the transient velocity field vector in the transient terminal field data output in step S2 and the application correction operator determined in step S3 into the steady-state initial velocity field reconstruction formula for product operation to obtain the initial velocity field vector passed to the steady-state solver; combining the initial velocity field vector with the temperature field and pressure field directly transferred in the transient terminal field data to form the initial value for steady-state solution.

[0019] Preferably, in step S4, the calculation under the real physical environment of removing the artificial viscosity model specifically means: forcibly setting the dimensionless artificial viscosity relaxation factor to zero, so that the fluid dynamic viscosity is strictly equal to the real physical viscosity at the current temperature, and solving the steady-state Navier-Stokes equations, continuity equations, and energy equations.

[0020] Preferably, in step S4, the specific process of performing the conductor resistivity update in the closed-loop iteration includes: extracting the volume average temperature of the high-voltage conductor during the steady-state iteration; and using the volume average temperature and the conductor resistivity temperature change correction formula to update the conductor resistivity and Joule heat source power density.

[0021] Preferably, in step S4, the specific process of performing the closed-loop iteration to correct the gas property pressure of the closed cavity until the calculation converges includes: using the ratio of the average temperature of the insulating gas domain in the current iteration step to the ambient temperature at the initial inflation, updating the background pressure of the closed cavity using the constant-volume pressure correction formula; the solver using the updated background pressure as the operating pressure to automatically update the gas density field of the entire domain; based on the updated resistivity and Joule heat source power density and the updated background pressure in this step, calculating the relative error using a dual relative error convergence criterion; when the maximum relative error of both is less than a preset convergence tolerance threshold, the calculation is determined to be converged. The convergence tolerance threshold is a preset tolerance value used to determine whether the electro-thermal-gas multiphysics coupling iteration has reached an equilibrium state.

[0022] This invention provides a rapid calculation method for the temperature rise distribution of DC GILs based on metastable state conditions. It offers the following advantages:

[0023] 1. This invention introduces an artificial viscosity model controlled by flow field characteristic monitoring criteria in the transient pre-calculation stage. By increasing the effective dynamic viscosity of the fluid in the early stage of calculation, the Reynolds number is controlled within the laminar range, which suppresses the non-physical numerical oscillations caused by excessive buoyancy driving force under natural convection at high Rayleigh numbers. This mechanism ensures that the numerical calculation can start stably and quickly establishes a flow field convection topology that conforms to physical laws while suppressing oscillations, thus solving the problem of easy divergence when directly performing steady-state calculations at high Rayleigh numbers.

[0024] 2. This invention proposes a spatial mapping transfer method based on momentum-energy ratio. It utilizes a constructed correction operator to spatially reconstruct and compensate for momentum in the transient calculated velocity field, eliminating the difference in momentum transfer efficiency caused by artificial viscosity. This method transforms flow field data under high-viscosity artificial environments into physically equivalent steady-state initial values, enabling the steady-state solver to start close to the true solution. This effectively shortens the computation time required for subsequent iteration convergence, achieving a balance between computational efficiency and stability.

[0025] 3. This invention constructs a closed-loop iterative mechanism that includes conductor resistivity updates and closed-cavity gas property pressure corrections. During the steady-state solution phase, the conductor resistivity and Joule heat source intensity are dynamically adjusted based on real-time temperature field feedback, and the background gas pressure and gas density field are corrected based on the constant-volume closed-loop characteristics. This strongly coupled multi-physics solution strategy eliminates the errors of single-physics calculations, ensuring that the final output temperature and flow field distributions accurately reflect the thermo-hydraulic characteristics of the DC GIL under operating conditions, thus improving the accuracy of the calculation results. Attached Figure Description

[0026] Figure 1 This is a flowchart illustrating the overall process of a rapid calculation method for DC GIL temperature rise distribution based on metastable state conditions, according to an embodiment of the present invention.

[0027] Figure 2 This is a schematic diagram of the three-dimensional geometric model and mesh generation of the DC GIL according to an embodiment of the present invention;

[0028] Figure 3 This is a schematic diagram of the adaptive transient pre-calculation process based on viscosity residual collaborative relaxation in an embodiment of the present invention;

[0029] Figure 4 This is a schematic diagram of the spatial mapping state transfer process based on momentum-energy ratio according to an embodiment of the present invention;

[0030] Figure 5This is a schematic diagram of the steady-state solution and closed-loop correction process for multiple physical parameters (electricity, heat, and gas) according to an embodiment of the present invention. Detailed Implementation

[0031] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] See attached document Figure 1 This invention provides a rapid calculation method for the temperature rise distribution of DC GIL (Gas Inlet Liquidity) based on metastable state conditions, comprising the following steps:

[0033] S1. Construct a three-dimensional multiphysics coupling model of a DC gas-insulated transmission line and perform mesh discretization. Set the basic material parameters and solve the steady-state current conservation equation based on the ambient temperature conditions to obtain the initial Joule heat source density.

[0034] S2 loads the initial Joule heat source density into the transient solver, introduces an artificial viscosity model controlled by the flow field characteristic monitoring criterion, performs adaptive transient pre-calculation based on viscosity residual collaborative relaxation, establishes the flow field skeleton while suppressing numerical oscillations, and outputs transient terminal field data including velocity field, temperature field and pressure field.

[0035] S3 reads transient end field data, performs spatial mapping transfer based on momentum-energy ratio, and uses the correction operator constructed by viscosity ratio and flow index to perform spatial reconstruction and momentum compensation of transient velocity field, generating physically equivalent steady-state solution initial values.

[0036] S4. Start the steady-state solver using the initial values ​​of the steady-state solution, perform calculations in a real physical environment where artificial viscosity is removed, and execute closed-loop iterations including conductor resistivity updates and pressure corrections for gas properties in the closed cavity until the calculation converges, outputting the final corrected steady-state temperature field and flow field distribution.

[0037] The technical details of each of the above steps will be explained in detail below, combining mathematical models and specific implementation logic.

[0038] See attached document Figure 2 In step S1, which involves constructing and initializing a multiphysics coupling model, the specific execution process can be further divided into the following sub-steps:

[0039] S101. Construct a three-dimensional geometric model of the DC gas-insulated transmission line and execute a mesh discretization strategy. Based on the actual structural parameters of the DC gas-insulated transmission line to be analyzed, a three-dimensional solid model is established, including the high-voltage conductor, insulating gas, basin insulators, and metal shell. The high-voltage conductor is typically set to aluminum alloy, the insulating gas to sulfur hexafluoride (SF6) or a mixture thereof, and the metal shell to ground potential. For the fluid-structure interaction interface, especially the outer surface of the high-voltage conductor and the inner surface of the metal shell, a boundary layer mesh refinement strategy is adopted. Specifically, a multi-layer prism layer is generated in the fluid domain near the solid wall to analyze the near-wall velocity and temperature gradients, satisfying the requirement for solving dimensionless wall distances; an unstructured tetrahedral mesh is used in the mainstream region of the insulating gas to ensure that the spatial discretization accuracy of the computational domain meets the requirements for solving high Rayleigh number natural convection. The specific operations of geometric modeling and basic mesh generation are conventional techniques well-known to those skilled in the art and will not be elaborated here.

[0040] S102, Establish the multiphysics control equations for fluid and solid.

[0041] Given the density changes of the gas and the natural convection effect driven by thermal buoyancy during the temperature rise process inside a DC GIL, the insulating gas is considered as a compressible Newtonian fluid. A set of partial differential equations obeying the laws of conservation of mass, momentum, and energy are established in both the fluid and solid domains, respectively, as the physical governing equations for numerical calculations.

[0042] First, for the insulating gas region, the compressible Navier-Stokes equations are adopted as the governing equations for the flow field. The fluid mass conservation equation is as follows:

[0043] ;

[0044] In the formula, This indicates the density of the insulating gas, expressed in kilograms per cubic meter (kg / m³). ); Time is expressed in seconds. ); This represents the fluid velocity vector, with units of meters per second (m / s). ); This represents the Nabla operator (differential operator), here This represents divergence calculation. The fluid momentum conservation equation is:

[0045] ;

[0046] In the formula, This represents the fluid pressure field variable, with units of Pascals (Pa). ); Expresses the dynamic viscosity of a fluid, measured in Pascals per second (Pa). ); Represents a second-order unit tensor; superscript Represents the matrix transpose symbol; This represents the gravitational acceleration vector, with units of meters per second squared (m²). The terms on the right side of the equation correspond to the pressure gradient force, the viscous stress term (including the bulk viscosity term), and the buoyancy source term, respectively.

[0047] Secondly, a multiphysics energy conservation equation is established for the entire domain (including the fluid and solid domains) to describe the evolution of the temperature field:

[0048] ;

[0049] In the formula, For general material density, take the following in the gas domain: In the solid domain, take the density of the corresponding solid material. ; This indicates the specific heat capacity at constant pressure, expressed in joules per kilogram of Kelvin (kJ / kg). ); Represents the temperature field variable, with units of Kelvin (K). ); This represents the general thermal conductivity, measured in watts per meter Kelvin (W / m). ); Joule heat source power density, expressed in watts per cubic meter (W / m³). In the solid domain, due to the velocity vector... If the result is zero, the above equation automatically degenerates into the solid heat conduction equation.

[0050] S103, perform initial current field solution and steady-state estimation of the starting heat source. At the initial moment before the transient calculation starts ( Assuming the system as a whole is in a uniform initial environment temperature Next, we temporarily disregard the effect of uneven temperature distribution on the material's conductivity and use the steady-state current conservation equation to calculate the current and potential distribution inside the high-voltage conductor. The steady-state current conservation equation is:

[0051] ;

[0052] In the formula, This represents the conductor conductivity at the initial ambient temperature, expressed in Siemens units per meter (S / m). ); Electric potential is expressed in volts (V). ).

[0053] Based on the calculated potential distribution, the Joule heat source density at the initial moment is calculated using the relationship between electric field strength and conductivity. :

[0054] ;

[0055] In the formula, That is, the initial Joule heat source density; This represents the magnitude of the vector. The initial heat source density will serve as the initial heat load input for the transient solver in subsequent step S2, driving the flow field to generate natural convection due to the temperature difference.

[0056] See attached document Figure 3 After completing the multiphysics model construction and initial heat source calculation, the process proceeds to step S2, which is the adaptive transient pre-calculation based on viscosity residual cooperative relaxation. This step constructs an artificial physical environment that dynamically evolves with the flow field development state, thereby suppressing numerical oscillations during the initiation of natural convection at high Rayleigh numbers while rapidly establishing the convection structure of the flow field. The specific implementation process includes the following sub-steps:

[0057] S201, construct a dynamic artificial viscosity model controlled by the flow field state.

[0058] In the transient calculation of the momentum equation, the dynamic viscosity term of the insulating gas is replaced with the effective dynamic viscosity. This effective dynamic viscosity is obtained by multiplying the gas's true physical viscosity by a time-varying artificial gain term. The fluid viscosity at each time step is calculated using the effective dynamic viscosity model, which is as follows:

[0059] ;

[0060] In the formula: This represents the effective dynamic viscosity actually substituted into the Navier-Stokes equations during transient calculations, expressed in Pascal-seconds (Pa). ); Indicates based on the current temperature field variables The true physical viscosity of the insulating gas, obtained by referring to tables or calculating using the Sutherland formula, is expressed in Pascals per second (Pa). ); This represents the dimensionless artificial viscosity relaxation factor, with its initial value. The value is set to be greater than 100, so that the fluid at the initial moment exhibits the target high viscosity characteristics, thereby keeping the Reynolds number at an extremely low level and ensuring that the calculation can be stably started from a low Reynolds number laminar flow state.

[0061] S202 performs flow field characteristic monitoring and smoothing based on Nusselt number derivatives.

[0062] To achieve closed-loop control of artificial viscosity, the average Nusselt number on the high-pressure conductor surface is selected as a global scalar index characterizing the intensity of natural convection evolution in the flow field. The convective heat transfer intensity on the high-pressure conductor surface is calculated in real time using the average Nusselt number calculation formula, which is:

[0063] ;

[0064] In the formula: The average Nusselt number (dimensionless) represents the surface of a high-voltage conductor. This represents the outer surface area of ​​a high-voltage conductor, expressed in square meters (m²). ); This represents the local convective heat transfer coefficient of a conductor surface, measured in watts per square kelvin (W / m²). ); The characteristic length is represented by its diameter for cylindrical conductors, and the unit is meters (m). The thermal conductivity of an insulating gas is expressed in watts per meter Kelvin (W / m). ); This represents the area integral operation performed on the outer surface of the high-voltage conductor.

[0065] Due to the discrete errors in numerical calculations, directly... Differentiation introduces high-frequency noise, causing instability in the control system. Therefore, it is necessary to calculate the characteristic value of the rate of change after filtering.

[0066] The flow field change rate is calculated using the Nusselt number derivative smoothing eigenvalue formula, which is:

[0067] ;

[0068] In the formula: This represents the characteristic value of the rate of change of the Nusselt number over time after smoothing, in reciprocal seconds (Rs). ); This represents the direct derivative of the mean Nusselt number with respect to time. This represents a numerical low-pass filter smoothing operator. In practice, it employs an exponentially weighted moving average (EWMA) algorithm or a sliding window average algorithm to filter out high-frequency numerical noise and extract the trend characteristics of flow field evolution.

[0069] S203, Execution of viscosity residual collaborative feedback evolution mechanism and gating function design.

[0070] To maintain high viscosity during periods of drastic flow field adjustment to suppress numerical oscillations, and to rapidly release viscosity during periods of flow field stabilization to approximate the true physical state, a relaxation factor is established. With flow field eigenvalues The negative feedback evolution relationship between them.

[0071] The relaxation factor at each time step is updated using an artificial viscosity relaxation factor evolution equation, which is as follows:

[0072] ;

[0073] In the formula: This represents the rate of change of the relaxation factor over time; This represents the viscosity decay rate constant, in reciprocating seconds (π / 2). ), used to control the base rate of viscosity release; This represents the flow field stability gate function, used to dynamically start and stop the viscosity decay process based on the flow field state. The definition of the flow field stability gate function is as follows:

[0074] ;

[0075] In the formula: This represents the preset flow field stability criterion threshold, in counts of seconds (p / s). ).when When this occurs, it indicates that a violent topological reconstruction is taking place inside the flow field, and the function output is 0, making... That is, to forcibly lock the current artificial viscosity value and prevent it from decreasing; when When the flow field evolution is stable, the function output is 1. At this time, the relaxation factor decays exponentially, reducing the effective dynamic viscosity. Gradually approaching the true physical viscosity .

[0076] S204, Determine the adaptive termination condition for transient calculations. At the end of each time step of the transient solution, check whether the system state meets the termination criterion. To minimize computation time while ensuring the correctness of the flow field structure, the transient pre-calculation can be terminated when the following quasi-steady-state condition of flow field evolution is met: the flow field evolution tends to a quasi-steady state:

[0077] (in (This is the preset tolerance for the rate of change of the Nusselt number). At this point, if the artificial viscosity relaxation factor... It has not yet decayed to zero (i.e.) Then, corrections are made based on subsequent mapping steps;

[0078] like It has decayed to the preset trace level. (For example, 0.01), then the condition is naturally met.

[0079] After meeting the above conditions, extract the flow field data at the current moment and output transient terminal field data including velocity field, temperature field and pressure field. , , ), which serves as the input for the subsequent step S3.

[0080] See attached document Figure 4 After completing the adaptive transient pre-calculation and obtaining the transient terminal field data, proceed to step S3. Because artificial viscosity is introduced in the transient calculation to ensure convergence stability at high Rayleigh numbers, the flow field at its end, while possessing the correct convective topology, deviates from the flow velocity amplitude in the real physical environment. To enable the subsequent steady-state solver to converge smoothly and quickly, a space mapping state transfer based on the momentum-energy ratio is performed, reconstructing the energy levels of the velocity field through a physical mapping operator. The specific implementation process includes the following sub-steps:

[0081] S301, construct a theoretical correction operator based on the flow index and viscosity ratio.

[0082] By utilizing the ratio between the effective dynamic viscosity and the true physical viscosity at the end of the transient calculation, and combining it with the current flow characteristics, a dimensionless theoretical correction operator is constructed to quantify the difference in momentum transport efficiency caused by viscosity differences.

[0083] The theoretical magnification factor of the velocity field is determined using the theoretical correction operator calculation formula, which is as follows:

[0084] ;

[0085] In the formula: This represents the calculated theoretical momentum correction operator, which is a dimensionless scalar. Indicates the end time of transient calculation The effective dynamic viscosity used in the last iteration, in Pascals per second (Pa). ); Represents the temperature field distribution at the end of the transient. The calculated physical viscosity of the real gas is expressed in Pascals per second (Pa). ); The flow index represents the nonlinear relationship between fluid momentum dissipation and viscosity. For laminar flow, which is common in DC GIL, the value ranges from 0.9 to 1.1. Represented by Grashof numbers The correction factor is used to compensate for the difference in nonlinear response of buoyancy driving force under different viscosities, and can be set to 1.0 in simplified calculations.

[0086] S302 performs operator safety limiting to prevent numerical divergence.

[0087] In numerical calculations, if If the numerical value is too large, directly doubling the velocity field will lead to excessively high local flow velocities, thereby violating the Courant-Friedrich-Levy (CFL) condition of the computational grid and causing subsequent steady-state calculations to diverge at startup. Therefore, a numerical stability limiting mechanism is introduced.

[0088] The final applied correction coefficient is determined using the applied correction operator limiting formula, which is as follows:

[0089] ;

[0090] In the formula: This represents the applied correction operator actually applied to the velocity field; This represents the preset upper limit threshold for the numerical stability of the correction operator. This threshold is determined based on the CFL constraint relationship between the mesh size and the time step, and the preferred range is 5.0 to 10.0. This indicates that the minimum value function operation is performed to ensure that the final correction ratio does not exceed the safety threshold. .

[0091] S303 performs the reconstruction of the initial velocity field and the transfer of physical fields for the steady-state solver.

[0092] Using the calculated application correction operator The transient velocity field is reconstructed globally to compensate for momentum deficit while maintaining the continuity of the scalar field, generating physically equivalent steady-state initial values.

[0093] The velocity vector field is updated using the steady-state initial velocity field reconstruction formula, which is:

[0094] ;

[0095] In the formula: This represents the initial velocity field vector passed to the steady-state solver, in meters per second. ; This represents the velocity field vector at the end of the transient calculation, in meters per second. .

[0096] Simultaneously, scalar field data is directly transmitted using the following steady-state initial scalar field transmission formula:

[0097] ;

[0098] ;

[0099] In the formula: and Represent the initial temperature field of the steady-state solver (units). ) and initial pressure field (unit) ); and These represent the temperature field and pressure field at the end of the transient calculation, respectively.

[0100] Through the above steps, the state mapping from a high-viscosity artificial field to a low-viscosity real field is completed in physical space, enabling the steady-state solver to start on a flow field that has been fully developed and whose momentum level is close to that of reality.

[0101] See attached document Figure 5 After obtaining the initial field for steady-state calculation using state-space mapping, the process proceeds to step S4. This step iteratively resolves the strong nonlinear coupling relationships between the electric field, flow field, temperature field, and pressure field in the DC GIL device, particularly the coupling characteristics of conductor resistivity changing with temperature and gas pressure within the sealed cavity changing with the average gas temperature. The specific implementation process includes the following sub-steps:

[0102] S401 performs a steady-state solution for the real physical field after removing artificial viscosity.

[0103] The output in step S3 , and The initial values ​​are substituted into the steady-state solver. At this stage, the artificial viscosity model introduced in step S2 is completely removed; that is, the relaxation factor is included in the solver configuration. This is forced to be 0, ensuring that the hydrodynamic viscosity is exactly equal to the actual physical viscosity at the current temperature. The solver solves the steady-state Navier-Stokes equations, continuity equations, and energy equations under real physical parameters to obtain the intermediate physical field distribution under the current heat source and pressure boundary conditions.

[0104] S402 performs dynamic correction of conductor resistivity based on temperature feedback.

[0105] The resistivity of high-voltage conductors increases with increasing temperature, leading to an increase in the power density of the Joule heat source. During steady-state iteration, the volume average temperature of the high-voltage conductor is extracted, and the material's conductivity and heat source intensity are updated accordingly.

[0106] The resistivity and Joule heat source of the conductor are updated using the conductor resistivity temperature change correction formula. The conductor resistivity temperature change correction formula is as follows:

[0107] ;

[0108] ;

[0109] In the formula: Indicates the current average temperature of the conductor. The resistivity of the sample is expressed in ohmmeters. ; Indicates reference temperature (The reference resistivity is taken at 20 degrees Celsius); The temperature coefficient of resistance of a conductor material is expressed in inverted Kelvin (TK). ); This represents the volume average temperature of the conductor in the current calculation step, in Kelvin (K). ); Indicates the first The Joule heat source power density corrected by the next coupling iteration step, in watts per cubic meter. ; This indicates the rated load current of the DC GIL, expressed in amperes (A). This represents the cross-sectional area of ​​a conductor, expressed in square meters (m²). ).

[0110] S403 performs closed-loop correction of background pressure and physical properties based on constant volume closed-loop characteristics.

[0111] DC GIL is a constant-volume closed system, and the insulating gas inside follows the law of conservation of mass. As the operating temperature increases, the background absolute pressure in the closed gas chamber increases, which in turn changes the operating density of the gas and affects the Rayleigh number of natural convection.

[0112] The background pressure is updated using the constant volume pressure correction formula, which is:

[0113] ;

[0114] In the formula: Indicates the first The background absolute pressure after the next iteration is expressed in Pascals (Pa). This indicates the initial inflation pressure (absolute pressure), and the unit is Pascal (Pa). Indicates the first The average temperature of the insulating gas domain obtained from the next iteration step is in Kelvin (K). ); This indicates the ambient temperature at the time of initial inflation, expressed in Kelvin (K). ).

[0115] In updating the background absolute pressure Then, the solver uses the updated background pressure as the operating pressure, automatically updates the global gas density field based on the compressible fluid state equation, and feeds the updated pressure field and density field back into the fluid control equation.

[0116] S404 executes the dual convergence criterion of heat source pressure and the iterative termination strategy.

[0117] To ensure the accuracy of multiphysics coupling calculations, a dual convergence index of heat source intensity and background pressure is set.

[0118] The double relative error convergence criterion is used to determine whether the iteration terminates. The double relative error convergence criterion is as follows:

[0119] ;

[0120] In the formula: This represents the maximum relative error of the current iteration step; and These represent the Joule heat source power densities of the current iteration step and the previous iteration step, respectively. and These represent the background absolute pressures of the current iteration step and the previous iteration step, respectively. This represents the function operation to find the maximum value. This represents the preset convergence tolerance threshold, with a preferred range of [value missing]. to .

[0121] like Then the corrected heat source and pressure will be used as the new boundary conditions, and the process will return to step S401 to continue the steady-state solution.

[0122] like If the calculation is successful, the iteration is terminated, and the final corrected steady-state temperature field, flow field, and pressure distribution are output.

[0123] To verify the effectiveness of the proposed DC GIL multiphysics adaptive coupling calculation method and ensure the feasibility of the technical solution, the following details the hardware and software environment configuration required for implementing this invention and the recommended value range of key empirical parameters of the algorithm.

[0124] Hardware and software implementation environment configuration:

[0125] Given that the three-dimensional DC GIL multiphysics coupling model involves multi-domain coupling of fluid, solid, electric, and thermal fields, and the number of meshes typically reaches millions, the computational process involves iterative solutions of large-scale sparse matrices, requiring the support of a high-performance computing (HPC) hardware platform. The hardware platform preferably employs a multi-core parallel processing workstation or server cluster, with a central processing unit (CPU) configured with a server-grade processor with a clock speed of 3.0 GHz or higher and at least 16 physical cores to support domain decomposition parallel computation; random access memory (RAM) capacity is configured to be at least 64 GB to meet the storage requirements of large-scale mesh data and intermediate iterative matrices. Regarding the software environment, the numerical computation platform is based on secondary development of general multiphysics simulation software using the finite volume method (FVM) or finite element method (FEM) architecture. The core algorithm of this invention, including an artificial viscosity model, gating functions, and state mapping operators, is written through user-defined function (UDF) interfaces or script interfaces (such as C, C++, or Python) provided by the simulation software and embedded into the iterative loop of the solver to achieve dynamic control of fluid viscosity, boundary conditions, and initial fields. In the solver's underlying settings, the fluid flow equations employ the SIMPLE or PISO pressure-velocity coupled algorithm, and the preferred spatial discretization scheme is the second-order upwind scheme to ensure computational accuracy.

[0126] Recommended values ​​for key empirical parameters of the algorithm:

[0127] To ensure the stable convergence of the adaptive transient pre-calculation based on viscosity residual cooperative relaxation (step S2) and the spatial mapping based on momentum-energy ratio (step S3), based on numerical experiments of DC GIL models for different voltage levels and pipe diameters, the following recommended ranges and physical meanings of the parameters are given:

[0128] Regarding the parameter settings for the artificial viscosity model, the initial value of the dimensionless artificial viscosity relaxation factor is... The set range is to The preferred value is 1000. The purpose of setting this value is to reduce the initial equivalent Reynolds number to below 1, forcing the flow field into a low Reynolds number laminar state, thereby suppressing non-physical oscillations in the early stages of the calculation. Viscosity decay rate constant. The set range is to The preferred value is This parameter determines the release rate of artificial viscosity. If the value is too large, the sudden drop in viscosity may cause secondary oscillations; if the value is too small, it will significantly reduce the time efficiency of transient calculations.

[0129] Regarding the parameter settings of the flow field stability gating function and the threshold of the flow field stability criterion. The set range is to The preferred value is This threshold characterizes the tolerance to the rate of change of the Nusselt number; the smaller the value, the more stringent the gating mechanism and the smoother the viscosity decay process. For artificial viscosity at minute levels... The value is set to 0.01, meaning that when the artificial viscosity component drops to 1% of the true viscosity, its impact on the flow field structure is considered negligible. For the quasi-steady-state criterion... The value is .

[0130] Regarding the parameter settings of the state mapping operator, the flow index It depends on the Rayleigh number range within the DC GIL. For most operating conditions involving laminar natural convection, a value of [value to be specified] is recommended. If local turbulence effects may occur in large-size UHV pipelines, the threshold can be adjusted to 0.6 to 0.8. Adjust the upper limit threshold for operator numerical stability. The set range is 5.0 to 10.0, with a preferred value of 8.0. This limit is used to prevent excessive multiplication factors due to numerical noise in the low-speed stagnation region, which could lead to local velocity overflow and divergence. It serves as the convergence tolerance threshold for steady-state calculations. , set as This is to ensure that the electro-thermal-gas multi-physics field coupling iteration achieves a balanced state with engineering precision.

Claims

1. A rapid calculation method for the temperature rise distribution of DC GIL combined with metastable state, characterized in that, Includes the following steps: S1. Construct a three-dimensional multiphysics coupling model of a DC gas-insulated transmission line and perform mesh discretization. Set the basic material parameters and solve the steady-state current conservation equation based on the ambient temperature conditions to obtain the initial Joule heat source density. S2. Load the initial Joule heat source density obtained in step S1 into the transient solver, introduce the artificial viscosity model controlled by the flow field feature monitoring criterion, perform adaptive transient pre-calculation based on viscosity residual collaborative relaxation, establish the flow field skeleton while suppressing numerical oscillations, and output transient terminal field data including velocity field, temperature field and pressure field. S3. Read the transient end field data output in step S2, perform spatial mapping transfer based on momentum-energy ratio, and use the correction operator constructed by viscosity ratio and flow index to perform spatial reconstruction and momentum compensation on the transient velocity field in the transient end field data to generate physically equivalent steady-state solution initial values. S4. Start the steady-state solver using the initial steady-state solution value generated in step S3. Perform calculations in a real physical environment with the artificial viscosity model removed. Execute a closed-loop iteration that includes conductor resistivity updates and pressure corrections for gas properties in the closed cavity until the calculation converges. Output the final corrected steady-state temperature field and flow field distribution.

2. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 1, characterized in that, In step S1, the construction process of the three-dimensional multiphysics coupling model includes: For fluid-structure interaction interfaces, multi-layer prismatic meshes are generated in the fluid domain near the solid wall to analyze the velocity and temperature gradients near the wall, satisfying the requirement for solving dimensionless wall distances; unstructured tetrahedral meshes are used in the mainstream region of the insulating gas. The governing equations of the three-dimensional multiphysics coupling model include: fluid mass conservation equations and fluid momentum conservation equations for the insulating gas region, and multiphysics energy conservation equations for the entire domain; the fluid momentum conservation equations include pressure gradient force, viscous stress term and buoyancy source term.

3. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 1, characterized in that, In step S1, the specific process of obtaining the initial Joule heat source density is as follows: Assuming the system is under a uniform initial ambient temperature, the potential distribution inside the high-voltage conductor is calculated using the steady-state current conservation equation. Based on the calculated potential distribution, the initial Joule heat source density at the initial moment is calculated by the relationship between the electric field strength and the conductor conductivity at the initial ambient temperature.

4. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 1, characterized in that, In step S2, the artificial viscosity model is specifically constructed using an effective dynamic viscosity model; The effective dynamic viscosity model defines the effective dynamic viscosity at each time step as the product of the true physical viscosity of the insulating gas and the artificial gain term. The artificial gain term consists of a unit value and a dimensionless artificial viscosity relaxation factor. The initial value of the dimensionless artificial viscosity relaxation factor is set to a preset high value, so that the fluid at the initial moment exhibits the target viscosity characteristics and the Reynolds number is in a laminar state.

5. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 4, characterized in that, In step S2, the specific process of performing the adaptive transient pre-calculation based on viscosity residual collaborative relaxation includes: The convective heat transfer intensity on the surface of the high-voltage conductor is calculated in real time using the average Nusselt number calculation formula; Using the convective heat transfer intensity and the Nusselt number derivative smoothing eigenvalue formula, the filtered Nusselt number time change rate eigenvalue is calculated. The Nusselt number time change rate characteristic value is substituted into the artificial viscosity relaxation factor evolution equation, and the dimensionless artificial viscosity relaxation factor at each time step is updated in combination with the flow field stability gating function. The flow field stability gating function dynamically outputs a first state value to forcibly lock the current artificial viscosity, or outputs a second state value to start the viscosity decay process, based on the relationship between the Nusselt number time change rate characteristic value and the preset flow field stability criterion threshold. When the characteristic value of the time rate of change of the Nusselt number is less than the preset tolerance of the rate of change of the Nusselt number, it is determined that the transient pre-calculation meets the termination condition. The flow field stability criterion threshold is a preset value that characterizes the tolerance for the rate of change of the Nusselt number; the Nusselt number change rate tolerance is a preset tolerance value used to determine whether the flow field evolution tends to a quasi-steady state.

6. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 4, characterized in that, In step S3, the specific process of constructing the correction operator includes: Using the ratio between the effective dynamic viscosity and the true physical viscosity at the end of the transient calculation as the viscosity ratio, and combining it with the flow index, the theoretical amplification factor of the velocity field is determined by the theoretical correction operator calculation formula. The theoretical amplification factor is compared with the preset upper limit threshold of numerical stability of the correction operator by using the applied correction operator limiting formula, and the smaller value is taken to obtain the final applied correction operator applied to the velocity field; The application correction operator is defined as the correction operator; The upper limit threshold for the numerical stability of the modified operator is a preset value determined based on the Courant number constraint relationship between the grid size and the time step.

7. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 6, characterized in that, In step S3, the specific process of generating the initial values ​​for the steady-state solution based on the spatial mapping transfer of the momentum-energy ratio includes: The transient velocity field vector in the transient end field data output in step S2 and the application correction operator determined in step S3 are substituted into the steady-state initial velocity field reconstruction formula for product operation to obtain the initial velocity field vector passed to the steady-state solver. The initial velocity field vector is combined with the temperature and pressure fields directly transmitted from the transient end field data to form the initial values ​​for the steady-state solution.

8. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 4, characterized in that, In step S4, the calculation under the real physical environment where the artificial viscosity model is removed specifically refers to: The dimensionless artificial viscosity relaxation factor is forcibly set to zero, so that the hydrodynamic viscosity is strictly equal to the actual physical viscosity at the current temperature. The steady-state Navier-Stokes equations, continuity equations, and energy equations are then solved.

9. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 1, characterized in that, In step S4, the specific process of performing the conductor resistivity update in the closed-loop iteration includes: During the steady-state iteration process, the volume average temperature of the high-voltage conductor is extracted; The resistivity of the conductor and the power density of the Joule heat source are updated using the volume average temperature and the conductor resistivity temperature change correction formula.

10. The method for rapid calculation of DC GIL temperature rise distribution based on metastable state as described in claim 9, characterized in that, In step S4, the specific process of performing the closed-loop gas property pressure correction in the closed-loop iteration until the calculation converges includes: The background pressure of the closed cavity is updated using the ratio of the average temperature of the insulating gas domain in the current iteration step to the ambient temperature at the initial inflation. The solver uses the updated background pressure as the operating pressure to automatically update the gas density field across the entire domain. Based on the resistivity and Joule heat source power density updated in step S9 and the background pressure updated in this step, the relative error is calculated using a dual relative error convergence criterion. When the maximum relative error of both is less than the preset convergence tolerance threshold, the calculation is determined to be converged. The convergence tolerance threshold is a preset tolerance value used to determine whether the electro-thermal-gas multiphysics coupling iteration has reached an equilibrium state.

Images