Simulation calculation method for ground gas leakage diffusion

Abstract

The invention relates to the field of urban gas pipeline integrity management, and discloses a ground gas leakage diffusion simulation calculation method, which comprises the following steps of: initializing environmental parameters and leakage physical conditions, constructing a rotation alignment wind field calculation grid, and dividing dense and sparse time sequences by adopting a dynamic time step strategy; a leakage flow dynamic sequence is generated in combination with linear interpolation, diffusion coefficients are calculated based on atmospheric stability level sections, and an intermediate coefficient matrix is pre-calculated to optimize the calculation efficiency of a Gaussian puff model; an RK4 numerical solution is introduced to process an unsteady state diffusion process, and fusion of an analytical model and a numerical method is realized; and performing standard condition conversion, anomaly correction and structured storage on the concentration data to generate a multi-dimensional visualization result. According to the method, minute-level high-resolution grid leakage diffusion simulation is realized, leaked gas volume concentration distribution and diffusion trend are accurately output, and reliable technical support is provided for smart city operation, city gas safety management, real-time early warning and emergency decision making.

Description

Technical Field

[0001] The invention relates to the technical field of urban gas pipeline integrity management, and in particular to a ground gas leakage diffusion simulation calculation method. Background Art

[0002] With the acceleration of global urbanization, gas, as the core carrier of urban energy supply, has become a key issue in the field of public safety. The frequent occurrence of gas leakage accidents not only threatens the safety of people's lives and property, but may also cause secondary environmental disasters, which have a profound impact on urban functions and social stability. The diffusion behavior after a gas leak is constrained by multiple dynamic factors, including the strength of the leakage source, meteorological conditions (wind speed, wind direction, atmospheric stability), topography and building interference, among which the dynamic changes in wind speed and wind direction will significantly change the diffusion path. The turbulent effect and obstacle shielding in the complex urban environment further aggravate the nonlinear characteristics of the concentration distribution. Therefore, accurately predicting the three-dimensional spatiotemporal evolution of gas diffusion is the core prerequisite for formulating emergency response strategies, and it is also the technical bottleneck for optimizing urban gas safety management.

[0003] In the field of gas leak diffusion simulation, existing technologies generally face the challenge of balancing computational efficiency and simulation accuracy. Although numerical simulation methods based on computational fluid dynamics (CFD) can provide high-precision results, their large grid size and strict stability conditions (such as Courant number restrictions) lead to a surge in demand for computing resources. For example, for a typical urban leak scenario (3 km×3 km×30 m computational domain), the number of nodes exceeds 2 million when a 5-meter resolution grid is used. Combined with the hour-level simulation time and the second-level time step, it takes several hours to several days of computing time, which is difficult to meet the timeliness requirements of real-time response to accidents. In addition, the CFD method's high reliance on the user's professional experience also limits its popularization and application.

[0004] Although the simplified method based on the Gaussian puff model can achieve fast calculations, its assumptions (such as constant wind speed, uniform flow field, steady-state leakage) deviate significantly from the actual scenario, especially when the wind direction changes dynamically, the leakage flow fluctuates, or complex terrain interferes, the prediction error increases significantly. For example, the existing Gaussian model does not take into account the strong time dependence of the leakage source, and cannot effectively characterize the shielding effect of buildings on the diffusion path, which limits its applicability in urban environments. At the same time, although empirical models that rely on experimental data or statistical regression can provide fast predictions in specific scenarios, their generalization ability is insufficient and it is difficult to adapt to the combined effects of variable environmental parameters and leakage conditions.

[0005] In general, existing technologies have significant shortcomings in terms of high computing efficiency, adaptability to complex environments, and dynamic process simulation capabilities: CFD methods are difficult to achieve minute-level responses due to computing resource bottlenecks; Gaussian models lead to distortion in complex scene predictions due to idealized assumptions; and empirical models cannot cover the coupling effects of multiple factors due to data limitations. How to build a gas diffusion simulation method that takes into account both efficiency and accuracy, supports real-time adjustment of dynamic parameters, and adapts to complex urban environments has become a technical problem that needs to be overcome in this field. Summary of the invention

[0006] In view of the shortcomings of the existing technology, the present invention provides a ground gas leakage diffusion simulation calculation method, which solves the technical problems of low calculation efficiency, insufficient dynamic process simulation capability and difficulty in balancing accuracy and real-time performance of the existing gas leakage diffusion simulation method in complex urban environments.

[0007] To achieve the above objectives, the present invention is implemented through the following technical solutions: a ground gas leakage diffusion simulation calculation method, comprising the following steps: Initializing environmental parameters, leakage source parameters, and calculation domain parameters, wherein the environmental parameters include wind speed, wind direction, and atmospheric stability level; Obtain flow and concentration data of underground leaks over time; Generate a computational grid, and adjust the grid direction based on the wind direction in the environmental parameter so that the grid coordinates are aligned with the wind field direction; Generate time series based on preset simulation duration and dynamic time step strategy; Interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points; dynamically calculating a diffusion coefficient based on the atmospheric stability level and the time series; Initialize the concentration matrix and intermediate coefficient matrix storing the gas concentration distribution; The gas concentration distribution is calculated based on an improved Gaussian puff model, wherein the improved Gaussian puff model introduces dynamic leakage parameters, wind speed-driven puff center displacement, and dynamic diffusion coefficient, and combines a numerical iteration method to update the concentration value; Perform standard conversion on the concentration calculation results and output the visualized data.

[0008] Preferably, the leakage source parameters include the initial flow rate, initial concentration and ground coordinates of the leakage point; the calculation domain parameters include the simulation area range, initial time step and total simulation time.

[0009] Preferably, the step of obtaining flow and concentration data of underground leakage points varying with time comprises: Extract discrete time series data of flow rate and concentration at the leakage point from the simulation results of underground pipeline leakage; The discrete data are checked for rationality, including flow non-negativity check and concentration physical range verification.

[0010] Preferably, the step of generating a calculation grid and adjusting the grid direction based on the wind direction in the environmental parameter comprises: Generate a uniform three-dimensional grid according to the preset grid size; The grid coordinates are rotated according to the wind direction angle. The transformation formula is: ; in, is the wind direction angle; is the grid coordinate before rotation; is the rotated grid coordinate.

[0011] Preferably, the step of generating a time series according to a preset simulation duration and a dynamic time step strategy comprises: Divide the total simulation time into at least two consecutive periods, where: The first time step is used in the first period , used to capture rapid concentration changes at the beginning of a leak; The second time step is used for the second period ,in , used to reduce redundancy in later calculations; The time step of subsequent periods is gradually increased according to the preset rules and meets the ,in The time period number.

[0012] Preferably, the step of interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points comprises: Obtain discrete time series data of flow and concentration from underground leakage simulation results; According to the flow rate and concentration values ​​at adjacent time points, the following formula is used to calculate Traffic flow at the moment and concentration : ; ; in, and Represents two adjacent known underground data time points; Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage concentration at Indicates at a point in time Underground leakage concentration at .

[0013] Preferably, the step of dynamically calculating the diffusion coefficient according to the atmospheric stability level and the time series comprises: A parameterized formula is selected according to the atmospheric stability level, and when the time t is less than a preset threshold, the diffusion coefficient is calculated according to the first formula; When the time t is greater than or equal to the preset threshold, the diffusion coefficient is calculated according to the second formula; The coefficients and exponents in the second formula are smaller than those in the first formula.

[0014] Preferably, the step of calculating the gas concentration distribution based on the improved Gaussian puff model includes: The improved Gaussian puff model is used to calculate the concentration distribution under unit flow rate. The formula is: ; in, express Time in space coordinates Gas mass concentration at express Dynamic leakage flow at each moment; , , It represents the diffusion coefficient in the horizontal, longitudinal and vertical directions; Indicates wind speed; Represents the rotated spatial grid coordinates; represents the simulation time; Unit flow concentration coefficient Precomputed as intermediate coefficient matrix , and based on real-time traffic Update actual concentration value: ; in, Indicates time Time grid point The corresponding unit flow concentration coefficient.

[0015] Preferably, the step of calculating the gas concentration distribution based on the improved Gaussian puff model further includes: In areas where the concentration gradient exceeds the preset threshold, the RK4 numerical method is used to correct the concentration value, including: ; in, is the diffusion coefficient tensor; is the Laplace operator of concentration; is the convection term caused by wind speed; Update the concentration according to the RK4 iterative formula: ; in, is the dynamic time step; , , , is the fourth slope estimate.

[0016] Preferably, the step of converting the concentration calculation result to standard conditions is to convert the mass concentration into the standard volume concentration.

[0017] The present invention provides a method for simulating and calculating ground gas leakage and diffusion, which has the following beneficial effects: 1. The present invention significantly improves the prediction accuracy of gas leakage diffusion by dynamically adjusting the diffusion coefficient and comprehensively considering multiple environmental factors, and can more realistically reflect the concentration distribution and diffusion range in complex urban environments. This effect helps to timely identify high-risk areas, reduce safety hazards caused by prediction deviations, and provide more reliable protection for the safety of residents' lives and property.

[0018] 2. Compared with the time-consuming traditional computational fluid dynamics methods, this invention adopts a fast iteration algorithm and a dynamic time step strategy, which greatly shortens the simulation calculation time, enabling relevant departments to quickly obtain diffusion prediction results after a gas leak accident occurs. This rapid response capability helps to formulate timely and effective evacuation and disposal plans, and minimize the economic losses and social impact of the accident.

[0019] 3. Aiming at the characteristics of changeable wind direction and complex terrain in urban environment, the present invention optimizes the adaptability of the model and can accurately simulate the diffusion behavior of gas in densely built areas or non-uniform terrain. This feature enables urban gas safety management personnel to more efficiently assess leakage risks, optimize pipeline inspections and emergency plans, and improve overall management efficiency.

[0020] 4. By converting the simulation results into standard volume concentration and supporting three-dimensional visualization output, the present invention provides emergency commanders with intuitive diffusion range and concentration distribution diagrams. Compared with the abstract data of traditional methods, the output form of this method simplifies the decision-making process, allowing non-professionals to quickly understand and apply the prediction results, and improving the operability of emergency response.

[0021] 5. The present invention is not only applicable to a single leakage scenario, but can also adapt to gas leakage accidents of different scales and conditions by adjusting parameters, providing a scientific basis for the planning, operation and maintenance of urban gas systems. Its high efficiency and accuracy make it valuable for promotion, and it helps to promote the overall improvement of energy security level in the process of urbanization. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 It is a schematic diagram of the method flow of the present invention; Figure 2 It is a schematic diagram of the leakage center point of the present invention; Figure 3 This is a diagram of ground diffusion results of an embodiment of the present invention. DETAILED DESCRIPTION

[0023] The following will be combined with the drawings in the specification of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0024] Please refer to the attached Figure 1 -Attached Figure 2 The present invention provides a simulation calculation method for ground gas leakage diffusion, which realizes efficient and accurate simulation of gas leakage diffusion by dynamically coupling environmental parameters and leakage source data, combining grid adaptive transformation and improved Gaussian smoke model.

[0025] like Figure 1 As shown, the ground gas leakage diffusion simulation calculation method may include the following steps: S1, initializing environmental parameters, leakage source parameters and calculation domain parameters; S2, obtain the flow and concentration data of underground leakage points changing with time; S3, generating a computational grid, and adjusting the grid direction based on the wind direction in the environmental parameters so that the grid coordinates are aligned with the wind field direction; S4, generating a time series according to a preset simulation duration and a dynamic time step strategy; S5, interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at continuous time points; S6. Dynamic calculation of diffusion coefficient based on atmospheric stability level and time series; S7, initializing a concentration matrix and an intermediate coefficient matrix storing gas concentration distribution; S8. Calculate gas concentration distribution based on improved Gaussian puff model; S9. Perform standard conversion on the concentration calculation results and output visual data.

[0026] The following is a detailed description of each step in the method of the present invention, which comprehensively describes the specific implementation principle, technical details and process of each step.

[0027] Regarding step S1, in this embodiment, initializing the environmental parameters, leakage source parameters, and calculation domain parameters is to provide complete input conditions of the physical scene for subsequent simulation calculations, ensuring the consistency between the model and the real environment.

[0028] The configuration of environmental parameters includes: Wind speed :Acquire real-time wind speed data in the area where the leak is located based on meteorological monitoring equipment, or set it to a fixed value based on historical meteorological data, in meters per second (m / s), with a range covering typical wind conditions for gas diffusion, such as 0 to 30 m / s. The wind speed parameter directly affects the calculation of the downwind displacement of the center of the smoke puff in subsequent steps.

[0029] wind direction : Based on the north direction, the wind direction angle is defined in radians, for example, the north direction is 0 radians and the east direction is π / 2 radians. This parameter is used to rotate the computational grid to ensure that the grid coordinate system is aligned with the wind field direction, thereby avoiding the diffusion direction deviation caused by the fixed coordinates of the traditional Gaussian model.

[0030] Atmospheric stability level: Determined according to the Pasquill classification method, the specific basis includes temperature gradient, solar radiation intensity and wind speed data, and is divided into A (extremely unstable) to F (extremely stable). This parameter is used to dynamically select the diffusion coefficient calculation formula to reflect the impact of different atmospheric stratification states on turbulent diffusion.

[0031] Preferably, wind speed and direction parameters can be dynamically updated through a real-time data interface to support real-time simulation in emergency scenarios.

[0032] The configuration of the leakage source parameters includes: Initial flow : The gas mass flow rate at the initial moment output by the underground pipeline leakage simulation software, in kilograms per second (kg / s), is determined by the pipeline pressure, leakage aperture and gas physical parameters. The initial flow rate is used as an input parameter of the Gaussian puff model to determine the initial diffusion intensity of the leak.

[0033] Initial concentration : It indicates the initial volume concentration of the gas at the leakage point, in parts per million (ppm). Its value needs to be set according to the explosion limit range of the gas components (such as methane and propane). For example, natural gas is usually set to a volume concentration of 5% to 15%.

[0034] Leakage point coordinates: define the three-dimensional coordinates with the ground leakage point as the center ,in Indicates that the leak source is located on the ground. This coordinate is used as the starting point for diffusion calculation and to locate the initial position of the center of the smoke puff.

[0035] The configuration of the computational domain parameters includes: Simulation area range: Set the three-dimensional space size according to the maximum impact range of gas diffusion, for example, the length, width and height are , , , unit: m, to ensure coverage of dangerous areas where gas may spread.

[0036] Dynamic time step strategy: Set the time step in segments according to the time-varying characteristics of the diffusion process to balance the calculation accuracy and efficiency. Preferably, a smaller step size (e.g. 0.1 second) is used in the early stage of leakage (0-60 seconds) to capture the rapid diffusion process; the step size is gradually increased (e.g. 1 second) in the middle stage (60-600 seconds); and a larger step size (e.g. 10 seconds) is used in the late stage (600-3600 seconds) to reduce redundant calculations in the low dynamic change stage.

[0037] Preferably, the computational domain parameters support custom region ranges, such as dynamically adjusting the simulation region according to geographic information system (GIS) data to meet the needs of different urban environments.

[0038] Through the above parameter initialization operation, this embodiment provides accurate physical scene input for subsequent steps and ensures the reliability and adaptability of model calculation.

[0039] For step S2, in this embodiment, the operation of obtaining the flow and concentration data of the underground leakage point that varies with time is to provide highly reliable dynamic input parameters for the subsequent diffusion simulation, ensuring that the model can accurately reflect the time-varying characteristics of the leakage process.

[0040] The flow rate and concentration data of the underground leakage point are generated by pipeline leakage simulation software, which simulates the transient leakage process after pipeline rupture based on fluid mechanics equations (such as Navier-Stokes equations) and outputs discrete time series data.

[0041] The data consists of a time series (satisfy ,in is the total simulation time), and the corresponding flow value (Unit: kg / s) and concentration value Unit: ppm). Preferably, the time point interval is set according to the dynamic change rate of the leakage, for example, dense sampling (interval of 1 second) is used in the early stage of leakage, and sparse sampling (interval of 10 seconds) is used in the later stage.

[0042] Before inputting data into the diffusion model, the following checks are performed to ensure that the data conforms to physical laws: Traffic non-negativity check: Detect traffic sequence Whether there is a negative value. If detected , it is set to zero and triggers an alarm signal to avoid model calculation failure due to negative flow.

[0043] Concentration range calibration: verify the concentration value Whether it is within the flammable range of the gas (for example, natural gas is 5% to 15% volume concentration). If the concentration exceeds the preset safety threshold, the abnormal data is marked and manual intervention is prompted to ensure that the input parameters are consistent with the actual leakage scenario.

[0044] The verified data is stored in a structured format (such as CSV or JSON), including three columns of data: timestamp, flow rate, and concentration.

[0045] Through the above operations, this embodiment ensures the physical rationality and timing continuity of the input data, and provides high-precision dynamic driving parameters for subsequent diffusion simulation.

[0046] For step S3, in this embodiment, the operation of generating a computational grid and adjusting the grid direction based on the wind direction in the environmental parameters is to construct a computational coordinate system consistent with the wind field direction, thereby optimizing the numerical calculation accuracy of the diffusion model and avoiding the diffusion path error caused by the wind direction deviation of the traditional fixed grid.

[0047] The generation of the computational grid is achieved through the following logic: Based on the calculation domain parameters defined in step S1 (such as the simulation area range ), generate a uniformly distributed three-dimensional structured grid. Preferably, the size of the grid in the horizontal direction (x, y) and the vertical direction (z) is set according to the simulation accuracy requirements, for example, the horizontal grid size can be set to 5 meters × 5 meters, and the vertical grid size can be set to 5 meters. The grid node coordinates are initially defined as a set of discrete points in a Cartesian coordinate system, covering the entire simulation area.

[0048] like Figure 2 Shown, as an example: Taking the underground pipeline leakage point as the center, the calculation domain is defined as 3km long ( =3000m、width 3km( =3000m), height 30m ( =30m), covering the main impact area of ​​gas diffusion in the urban environment. The height of 30m takes into account the ground diffusion height. A uniform three-dimensional grid is created in the calculation domain. Assuming that the grid unit size is set to 5m×5m×5m, the entire space is gridded, and the number of grids is calculated as follows: Number of grids in X direction: 3000m / 5m = 600; Number of grids in Y direction: 3000m / 5m = 600; Number of grids in Z direction: 30m / 5m = 6; The total number of grids is 600×600×6=2,160,000 (about 2.16 million). The grid size affects the calculation efficiency. The finer the space is divided, the more grids are formed, the higher the calculation accuracy, and the slower the calculation speed. This part directly affects the calculation speed and accuracy of the entire model.

[0049] Dynamic adjustment of the grid orientation is achieved through rotation transformation: According to the wind direction angle parameter obtained in step S1 , the horizontal coordinates of the grid nodes are rotated so that the transformed grid The axis is aligned with the wind direction. The specific transformation formula is: ; in: is the grid node coordinate in the original Cartesian coordinate system; is the coordinate in the new coordinate system after rotation.

[0050] Through this transformation, the wind field direction and the grid The axis direction is consistent, thus simplifying the calculation logic of the downwind displacement of the smoke puff in the subsequent diffusion model.

[0051] After rotation, the leakage diffusion direction is consistent with the wind direction to avoid concentration distribution deviation caused by wind direction deviation. For example, if θ=π / 4 (northeast wind, θ=0 is set to due north), the grid will be redistributed along the 45° direction. Ensure that the grid does not exceed , , range; at the same time, ensure that the coordinates of the grid point and the leakage point , , The total number of grids is about 2.16 million, which is sufficient to support the fast iterative calculation of the Gaussian puff model and solve the computational burden of tens of millions of grids in the CFD method.

[0052] For step S4, in this embodiment, the operation of generating a time series according to the preset simulation duration and dynamic time step strategy is to adjust the time step in stages to improve the calculation efficiency while ensuring the simulation accuracy, so as to adapt to the dynamic characteristics differences at different stages in the gas leakage and diffusion process.

[0053] The construction of the dynamic time step strategy is implemented through the following logic: Based on the total simulation time defined in step S1 , the simulation process is divided into multiple stages, each stage using a different time step.

[0054] Preferably, the stage division is based on the physical characteristics of leakage diffusion, including an initial rapid diffusion stage, a mid-term stable diffusion stage, and a late slow decay stage. The time step of each stage is dynamically adjusted according to the diffusion rate, with a smaller step size in the initial stage to capture rapid changes and a larger step size in the late stage to reduce redundant calculations.

[0055] In this embodiment, the method for generating a time series includes the following steps: 1. Stage division and step size setting: Initial stage : Use a smaller time step , which is used to accurately describe the dramatic changes in flow and concentration caused by the sudden pressure drop at the initial stage of leakage.

[0056] Mid-term : Use a medium time step , balancing computational accuracy and efficiency.

[0057] Late Stage : Use a longer time step , adapting to the smooth change process after the diffusion rate slows down.

[0058] in, , is the preset time threshold, satisfying .

[0059] 2. Time point generation and merging: In each stage, a sequence of equally spaced time points is generated according to the corresponding step length, and the formula is: ; Finally, the time points of all stages are merged into a global time series in ascending order. ,in , .

[0060] 3. Time point deduplication and verification: Perform deduplication on the merged time series to ensure , and check whether the time points cover the entire simulation duration. If there are gaps (such as due to step size switching ), additional time points are inserted to ensure continuity.

[0061] As an example: Simulation total duration setting: The simulation duration is initialized to T = 3600 seconds (1 hour) to construct a time series covering the complete leakage diffusion process.

[0062] Phase division and step size setting: 1. Initial stage ( Second): Step Length: Second Time point generation: by formula Generate a discrete time series consisting of Seconds, total points indivual.

[0063] In the early stage of gas leakage, the concentration gradient changes dramatically due to the sudden drop in pressure. A small step size can accurately capture the transient characteristics of the diffusion front (such as the peak concentration position).

[0064] 2. Mid-term stage ( Second): Time point generation: by formula Generate a sequence containing Seconds, total points indivual.

[0065] The diffusion enters a stable stage, the concentration changes slowly, and the medium step size reduces redundant calculations while ensuring accuracy.

[0066] 3. Late stage ( Second): Step Length: Second Time point generation: by formula Generate a sequence containing Seconds, total points indivual.

[0067] After the smoke puff has fully diffused, the concentration gradient is significantly reduced. A large step size can accelerate the simulation of the dilution process and improve computational efficiency.

[0068] Time series merging and verification: Merge operation: Merge the three-stage time points into a global sequence in ascending order, with a total number of points N=600+540+300=1440.

[0069] Deduplication and gap processing: Check whether the intervals between adjacent time points are strictly increasing. If there are gaps (such as the phase switching points t=60 seconds and t=600 seconds), ensure continuity by inserting intermediate time points.

[0070] The small step size in the initial stage of dynamic step size ensures accuracy, and the large step size in the middle and late stages reduces the calculation burden, so that the total simulation time is controlled within seconds to minutes to meet real-time emergency needs.

[0071] For step S5, in this embodiment, the operation of interpolating underground leakage flow and concentration data based on time series is to convert the leakage parameters of discrete time points into continuous time series, provide dynamic input for the subsequent Gaussian puff model, and ensure that the diffusion simulation can reflect the real physical process of leakage parameters changing over time.

[0072] The extraction and preprocessing of discrete data points are based on the original leakage data obtained in step S2, including the time point series , corresponding flow value And concentration value . Time series must meet strict time order ( ,in is the total simulation time), and the data verification module is used to ensure that the flow rate is non-negative and the concentration is within the flammable range. Preferably, if there is an uneven distribution of time points or missing data, a pre-interpolation completion strategy is adopted, such as linear extrapolation of adjacent data or insertion of virtual time points.

[0073] The dynamic determination of the interpolation interval is achieved through a binary search algorithm. , quickly locate adjacent known time points in a time series and ,satisfy The algorithm compares the target time The search range is gradually narrowed down until the midpoint value of the time series is found. Preferably, for large-scale time series, an index optimization strategy is adopted to accelerate the search process through a hash table or a skip table.

[0074] Linear interpolation calculations are performed based on parameter values ​​at adjacent time points. and concentration , and their interpolation formulas are: ; ; in: and Represents two adjacent known underground data time points; Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage concentration at Indicates at a point in time Underground leakage concentration at .

[0075] The linear interpolation method can effectively describe the typical characteristics of the linear change of leakage parameters over time while ensuring the computational efficiency, such as the gradual decrease in flow rate caused by pipeline pressure decay.

[0076] Linear interpolation assumes that flow and concentration change linearly in a short period of time. If the interval between underground data points is too large (e.g., >30 s), the error can be reduced by increasing the number of sampling points or using a higher-order interpolation (e.g., quadratic interpolation).

[0077] For step S6, in this embodiment, the diffusion coefficient of the Gaussian puff model is dynamically calculated based on the atmospheric stability level and time series, which combines meteorological conditions with time evolution characteristics, and accurately describes the three-dimensional spatial broadening characteristics of gas diffusion through a piecewise parameterized formula, thereby reflecting the influence of different atmospheric turbulence intensities on the concentration distribution morphology.

[0078] The mapping and parameter selection of atmospheric stability levels are based on the environmental parameters defined in step S1. According to the Pasquier classification method, atmospheric stability is divided into multiple levels (e.g., level A to level F), and each level corresponds to a set of predefined diffusion coefficient calculation parameters. Preferably, the parameters include coefficients , , With index , , , whose value is calibrated according to experimental data or historical observation results. For example: Extremely unstable atmosphere (such as Class A): The downwind diffusion coefficient can be preferably set to ,index , to reflect the rapid lateral diffusion caused by strong turbulence; Neutral stable atmosphere (such as Class D): The vertical diffusion coefficient can be preferably configured as ,index , characterizing the diffusion characteristics under the confined mixing layer height.

[0079] The segmented calculation logic of the diffusion coefficient is dynamically switched by presetting the time threshold. Preferably, the time threshold can be set The diffusion process is divided into two stages: initial rapid diffusion and later stable diffusion: 1. Initial stage : Using higher coefficients and exponents, for example, for Class A stability, the downwind diffusion coefficient is calculated as: ; 2. Late Stage : Switch to a lower coefficient and exponent. For example, the adjusted downwind diffusion coefficient is: ; Parameter switching is achieved through table lookup, and different stability levels correspond to independent segmented parameter tables.

[0080] In this embodiment, the segmented parameter table corresponding to the atmospheric stability level is as follows:

[0081] The general formula for the three-dimensional diffusion coefficient is expressed as: ; ; ; in, , , , is the diffusion intensity coefficient in each direction, , , is the diffusion rate exponent, and its value is dynamically selected according to the stability level. The power function form can simultaneously describe the magnitude (controlled by the coefficient) and rate (controlled by the exponent) of the diffusion range growing over time.

[0082] For step S7, in this embodiment, the operation of initializing the concentration result matrix and the intermediate coefficient matrix is ​​to provide an efficient data storage structure and pre-calculated parameters for the Gaussian puff model, reduce the real-time computing overhead through the space-for-time strategy, and ensure the integrity and traceability of the concentration data during the diffusion simulation process.

[0083] The concentration result matrix is ​​constructed and initialized based on the computational grid generated in step S3 and the time series defined in step S4. ) is a four-dimensional tensor structure with dimensions ,in is the number of time steps (1440 in this example), , , The number of nodes in the downwind, acrosswind and vertical directions of the calculation grid are respectively (the number of nodes in the example grid is ). Each element Indicates that at time step , grid coordinates after rotation Gas volume concentration at the location (unit: kg / m 3 ). Preferably, the initial values ​​of the matrices are all set to zero, indicating that the gas has not yet diffused into the calculation domain when the simulation starts.

[0084] The pre-calculation logic of the intermediate coefficient matrix is ​​based on the diffusion coefficients dynamically calculated in step S6 ( ) and the grid coordinates rotated in step S3. The intermediate coefficient matrix (denoted as ) is consistent with the concentration matrix dimension, and its elements Stores precomputed time- and space-dependent terms for the Gaussian puff model.

[0085] For step S8, in this embodiment, the operation of calculating the gas concentration distribution based on the improved Gaussian puff model and the numerical solution method combines the computational efficiency of the analytical model with the high-precision characteristics of the numerical method, and updates the grid point concentration value through multi-stage iterations to accurately simulate the spatiotemporal evolution process of gas diffusion.

[0086] The analytical calculation of the improved Gaussian puff model is based on the intermediate coefficient matrix pre-calculated in step S7 and the leakage flow interpolated in step S5 The concentration calculation formula is: ; in: express Time in space coordinates Gas mass concentration at express Dynamic leakage flow at each moment; , , represents the diffusion coefficient in the horizontal, longitudinal and vertical directions, which is the diffusion coefficient dynamically calculated in step S6; represents the wind speed, which is defined in step S1; Represents the rotated spatial grid coordinates; Represents the simulation time, which is the time point in the time series of step S4.

[0087] First calculate the concentration coefficient per unit flow , which represents the concentration contribution of each kilogram of flow in space, and is stored in the intermediate coefficient matrix of step S7 in advance. , obtained from step S5 ,and In Multiply to get the actual concentration: ; in, is the grid point index. Through the above, the real-time calculation overhead can be significantly reduced, so that the concentration update is simplified to the matrix element multiplication operation.

[0088] The numerical solution of the dynamic diffusion process is implemented by the fourth-order Runge-Kutta (RK4) method, which is used to deal with unsteady diffusion or complex boundary conditions. The convection-diffusion equation is expressed as: ; in, is the diffusion coefficient tensor; is the Laplace operator of concentration; is the convection term caused by wind speed.

[0089] The implementation steps of the RK4 method are as follows: 1. Slope estimation: Calculate the current time step The slope : ; in is the right-hand side of the convection-diffusion equation.

[0090] 2. Calculate the intermediate slopes one by one ; 3. Concentration update: Update the concentration value at the next time step using the weighted average slope: ; Iteratively update the concentration at each grid point There are four estimates, each with a time step of The fourth slope estimate within represents the rate of change of concentration Their role is to improve the concentration change trend from the current moment by trying it multiple times. Until the next moment Finally, RK4 uses weighted average Updates the concentration value.

[0091] For step S9, in this embodiment, the operations of correcting the concentration results, converting units and storing data are to convert the physical concentration data output by the simulation into volume concentration under standard working conditions, and provide intuitive and parseable data support for emergency decision-making through structured storage and visualization interface.

[0092] The correction process of the concentration data is based on the concentration matrix calculated in step S8. , including the following operations: 1. Unit conversion: Convert mass concentration (unit: kg / m 3 ) is converted to volume concentration (unit: ppm).

[0093] 2. Boundary condition correction: at the ground boundary The volume concentration is forced not to exceed the lower flammable limit (such as 5% volume concentration for methane) to avoid numerical overflow due to model errors.

[0094] 2. Outlier filtering: For grid points whose volume concentration exceeds the preset safety threshold (such as the upper flammable limit of 15%), median filtering or Gaussian smoothing is used to eliminate isolated noise points.

[0095] Data storage and structured output are achieved through the following methods: Standardized storage format: The corrected volume concentration matrix and the original mass concentration matrix are stored as HDF5 or NetCDF format files, containing the following data layers: Concentration data layer: stores volume concentration values ​​by time step and spatial coordinates; Metadata layer: records simulation parameters (wind speed, stability level, leakage flow sequence); Grid description layer: stores the calculation grid coordinates and rotation matrix parameters defined in step S3.

[0096] Lightweight caching: To meet real-time visualization requirements, the concentration data of key sections is extracted and stored in JSON or CSV format, supporting fast reading and rendering.

[0097] The visual output interface includes the following functions: 3D dynamic rendering: Based on OpenGL or WebGL engine, the volume concentration data is mapped to the computational grid to generate spatiotemporal evolution animation. Preferably, isosurface drawing (such as 1% LEL, 5% LEL isosurface) and thermal map overlay are supported.

[0098] 2D cross-section analysis: Provides interactive tools to select any spatial cross-section (such as horizontal plane, vertical plane), generate concentration distribution contour map or pseudo-color map, and mark the peak concentration position and diffusion range.

[0099] Timeline control: allows users to drag the timeline to view the concentration distribution at a specific moment, or export the concentration change curve within a specified time interval.

[0100] In general, the present invention forms a full-process technical solution covering data preprocessing, dynamic simulation and result analysis by initializing environmental parameters, constructing a computing grid for rotating and aligned wind fields, dynamic time step strategy, leakage data interpolation, segmented diffusion coefficient calculation, pre-calculated matrix optimization, hybrid simulation calculation combining improved Gaussian smoke puff model and RK4 numerical solution, as well as standard condition conversion and visualization output. It can efficiently simulate the three-dimensional spatiotemporal diffusion process of gas leakage, and significantly improve the calculation efficiency through matrix pre-calculation, dynamic parameter switching and parallel optimization strategy. Finally, it outputs the volume concentration distribution and interactive visualization results under standard conditions, providing accurate data support for real-time monitoring, risk assessment and emergency decision-making of gas leakage accidents.

[0101] Example: This embodiment provides an implementation example of a gas leakage diffusion simulation method, and verifies the feasibility and efficiency of the method through actual leakage scenarios and computing performance tests, as follows: Implementation scenarios and parameter configuration: A buried gas pipeline leakage accident was used as the simulation object. The pipeline operating pressure was 1 MPa, the leakage aperture diameter was 100 mm, the pipeline buried depth was 5 m, and the leakage direction was vertically upward. The simulation environment parameters included: northeast wind direction (wind speed 3.5 m / s, wind direction angle 45°), atmospheric stability level D, and the calculation domain covered a spatial range of 3 km×3 km×30 m. The leakage simulation was divided into two stages, above ground and underground, each lasting 300 seconds, and the total simulation time was 3600 seconds.

[0102] The implementation steps are as follows: 1. Environment initialization and grid construction: Rotate the computational domain according to the wind direction angle of 45° to align the grid along the wind direction with the wind field direction and eliminate coordinate deviation; Generate a computational grid with a spatial resolution of 5 m × 5 m × 5 m (with a total of approximately 2.16 million nodes), covering the leak point coordinates (0,0,0) to the far boundary; In the dynamic time step strategy, a 0.1 second step size is used in the initial stage (0-60 seconds) to capture rapid diffusion, and a 10 second step size is switched to accelerate the calculation in the later stage (600-3600 seconds). The total number of time steps is about 1440.

[0103] 2. Leakage data and diffusion parameter calculation: Subsurface leakage model output flow series , the initial flow rate is about 0.5 kg / s, which decays linearly with time (e.g., the flow rate is 0.45 kg / s at 300 seconds); Based on D-level stability, the diffusion coefficient is dynamically calculated. In the initial stage ( seconds) The downwind diffusion coefficient is (For example, when t=300 seconds, ), the vertical coefficient is significantly smaller than the horizontal coefficient, reflecting the diffusion inhibition characteristics of the stable atmospheric stratification.

[0104] 3. Pre-calculation and simulation execution: Initialize the four-dimensional concentration matrix (time × space) and the intermediate coefficient matrix to store the concentration contribution parameters under unit flow rate; The improved Gaussian puff model is used to calculate the concentration distribution, combined with the pre-calculated matrix to achieve real-time simulation (such as processing about 100,000 grid nodes per second), and supports switching to the RK4 method to handle dynamic diffusion processes; The calculated results are converted to standard conditions, the mass concentration is converted into volume concentration (unit: ppm), and the ground boundary concentration is limited to not exceed the lower flammable limit (such as 5% volume concentration).

[0105] 4. Visualization and output: Generate a visualization of the ground diffusion results such as Figure 3 shown.

[0106] Computational performance verification: In the hardware environment of 12th generation Intel i7-12700H processor (20 cores) and 32GB memory, the operating efficiency under different grid densities is tested: 10-meter grid (about 270,000 nodes): The full process runs in about 1.7 minutes, suitable for quick preliminary evaluation; 5-meter grid (about 2.16 million nodes): running time is about 2.1 minutes, and the accuracy meets the needs of real-time emergency response; Finer grid (e.g. 3 meters, about 33.75 million nodes): The running time is increased to about 4 minutes, which is suitable for high-precision post-mortems.

[0107] This example realizes minute-level simulation of leakage diffusion at a 5-meter grid resolution through dynamic parameter adjustment (time step, diffusion coefficient) and pre-calculation optimization, accurately depicts the diffusion path of gas from the leakage point to the northeast and the concentration attenuation trend, and supports multi-level grid density to adapt to different emergency scenarios, verifying the method's ability to balance real-time performance, accuracy and resource consumption.

[0108] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for simulating and calculating ground gas leakage and diffusion, characterized in that: The following steps are involved: Initializing environmental parameters, leakage source parameters, and calculation domain parameters, wherein the environmental parameters include wind speed, wind direction, and atmospheric stability level; Obtain flow and concentration data of underground leaks over time; Generate a computational grid, and adjust the grid direction based on the wind direction in the environmental parameter so that the grid coordinates are aligned with the wind field direction; Generate time series based on preset simulation duration and dynamic time step strategy; Interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points; dynamically calculating a diffusion coefficient based on the atmospheric stability level and the time series; Initialize the concentration matrix and intermediate coefficient matrix storing the gas concentration distribution; The gas concentration distribution is calculated based on an improved Gaussian puff model, wherein the improved Gaussian puff model introduces dynamic leakage parameters, wind speed-driven puff center displacement, and dynamic diffusion coefficient, and combines a numerical iteration method to update the concentration value; Perform standard conversion on the concentration calculation results and output the visualized data.

2. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The leakage source parameters include the initial flow rate, initial concentration and ground coordinates of the leakage point; the calculation domain parameters include the simulation area range, initial time step and total simulation time.

3. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The step of obtaining the flow rate and concentration data of the underground leakage point varying with time comprises: Extract discrete time series data of flow rate and concentration at the leakage point from the simulation results of underground pipeline leakage; The discrete data are checked for rationality, including flow non-negativity check and concentration physical range verification.

4. The method for simulating and calculating the ground gas leakage diffusion according to claim 1, characterized in that: The step of generating a computational grid and adjusting the grid direction based on the wind direction in the environmental parameter comprises: Generate a uniform three-dimensional grid according to the preset grid size; The grid coordinates are rotated according to the wind direction angle. The transformation formula is: ； in, is the wind direction angle; is the grid coordinate before rotation; is the rotated grid coordinate.

5. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The step of generating a time series according to a preset simulation duration and a dynamic time step strategy comprises: Divide the total simulation time into at least two consecutive periods, where: The first time step is used in the first period , used to capture rapid concentration changes at the beginning of a leak; The second time step is used in the second period ,in , used to reduce redundancy in later calculations; The time step of subsequent periods is gradually increased according to the preset rules and meets the ,in The time period number.

6. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The step of interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points comprises: Obtain discrete time series data of flow and concentration from underground leakage simulation results; According to the flow rate and concentration values ​​at adjacent time points, the following formula is used to calculate Traffic flow at the moment and concentration : ； ； in, and Represents two adjacent known underground data time points; Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage flow at Indicates at a point in time Underground leakage concentration at Indicates at a point in time Underground leakage concentration at .

7. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The step of dynamically calculating the diffusion coefficient according to the atmospheric stability level and the time series comprises: A parameterized formula is selected according to the atmospheric stability level, and when the time t is less than a preset threshold, the diffusion coefficient is calculated according to the first formula; When the time t is greater than or equal to the preset threshold, the diffusion coefficient is calculated according to the second formula; The coefficients and exponents in the second formula are smaller than those in the first formula.

8. The ground gas leakage diffusion simulation calculation method according to claim 1 is characterized in that: The step of calculating the gas concentration distribution based on the improved Gaussian puff model includes: ； in, express Time in space coordinates Gas mass concentration at express Dynamic leakage flow at each moment; , , It represents the diffusion coefficient in the horizontal, longitudinal and vertical directions; Indicates wind speed; Represents the rotated spatial grid coordinates; represents the simulation time; Unit flow concentration coefficient Precomputed as intermediate coefficient matrix , and based on real-time traffic Update actual concentration value: ； in, Indicates time Time grid point The corresponding unit flow concentration coefficient.

9. The method for simulating and calculating the ground gas leakage diffusion according to claim 8, characterized in that: The step of calculating the gas concentration distribution based on the improved Gaussian puff model also includes: In areas where the concentration gradient exceeds the preset threshold, the RK4 numerical method is used to correct the concentration value, including: ； in, is the diffusion coefficient tensor; is the Laplace operator of concentration; is the convection term caused by wind speed; Update the concentration according to the RK4 iterative formula: ； in, is the dynamic time step; , , , is the fourth slope estimate.

10. The ground gas leakage diffusion simulation calculation method according to claim 1, characterized in that: The step of converting the concentration calculation result to standard conditions is to convert the mass concentration into the standard volume concentration.

Images