A simulation calculation method for ground gas leakage and diffusion

Abstract

This invention relates to the field of urban gas pipeline integrity management and discloses a method for simulating and calculating ground gas leak diffusion. The method comprises the following steps: initializing environmental parameters and leakage physical conditions, constructing a rotationally aligned wind field computational grid, and employing a dynamic time step strategy to divide dense and sparse time series; generating a dynamic series of leakage flow using linear interpolation, calculating diffusion coefficients based on atmospheric stability levels in segments, and precalculating intermediate coefficient matrices to optimize the computational efficiency of the Gaussian puff model; introducing the RK4 numerical solution to handle the non-steady-state diffusion process, integrating analytical models with numerical methods; and performing standard condition conversion, anomaly correction, and structured storage on concentration data to generate multi-dimensional visualization results. This invention implements minute-level, high-resolution grid-based leakage diffusion simulation, accurately outputting the volume concentration distribution and diffusion trend of leaked gas, and provides reliable technical support for smart city operations, urban gas safety management, real-time early warning, and emergency decision-making.

Description

Technical Field

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

[0002] With the acceleration of global urbanization, the safe operation of 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 also may trigger secondary environmental disasters, with profound impacts on urban functions and social stability. The diffusion behavior of gas leaks is constrained by multiple dynamic factors, including the intensity of the leak source, meteorological conditions (wind speed, wind direction, atmospheric stability), topography, and interference from buildings. Dynamic changes in wind speed and direction can significantly change the diffusion path, and turbulent effects and obstruction in complex urban environments further exacerbate the nonlinear characteristics of concentration distribution. Therefore, accurately predicting the three-dimensional spatiotemporal evolution of gas diffusion is a core prerequisite for formulating emergency response strategies and a 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. While numerical simulation methods based on computational fluid dynamics (CFD) can provide highly accurate results, their large grid sizes and stringent stability requirements (such as Courant number limits) lead to a surge in computational resource requirements. For example, for a typical urban leak scenario (a 3 km × 3 km × 30 m computational domain), using a 5-meter resolution grid would require over 2 million nodes. Combined with hourly simulation times and second-level time steps, this would require computational timescales ranging from several hours to several days, making it difficult to meet the timeliness requirements of real-time incident response. Furthermore, CFD methods' heavy reliance on user expertise limits their widespread application.

[0004] While simplified methods based on the Gaussian puff model enable rapid calculations, their assumptions (such as constant wind speed, uniform flow field, and steady-state leakage) deviate significantly from real-world scenarios. This is particularly true when wind direction changes dynamically, leakage flow fluctuates, or complex terrain interferes, significantly increasing prediction errors. For example, existing Gaussian models fail to account for the strong temporal dependence of leakage sources and cannot effectively characterize the obstruction effects of buildings on diffusion paths, limiting their applicability in urban environments. Furthermore, while empirical models relying on experimental data or statistical regression can provide rapid predictions in specific scenarios, their generalization capabilities are insufficient, making it difficult to adapt to the combined influence of variable environmental parameters and leakage conditions.

[0005] In summary, existing technologies have significant shortcomings in high computational efficiency, adaptability to complex environments, and the ability to simulate dynamic processes. CFD methods struggle to achieve minute-level responses due to computational resource bottlenecks; Gaussian models distort predictions in complex scenarios due to idealized assumptions; and empirical models, due to data limitations, are unable to account for the coupled effects of multiple factors. Developing a gas diffusion simulation method that balances efficiency and accuracy, supports real-time dynamic parameter adjustment, and adapts to complex urban environments has become a pressing technical challenge in this field. Summary of the Invention

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

[0007] To achieve the above objectives, the present invention is implemented through the following technical solutions: A ground gas leakage diffusion simulation calculation method includes the following steps:

[0008] Initializing environmental parameters, leakage source parameters, and calculation domain parameters, wherein the environmental parameters include wind speed, wind direction, and atmospheric stability level;

[0009] Obtain flow and concentration data of underground leak points over time;

[0010] 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;

[0011] Generate time series based on preset simulation duration and dynamic time step strategy;

[0012] interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points;

[0013] dynamically calculating a diffusion coefficient based on the atmospheric stability level and the time series;

[0014] Initialize the concentration matrix and intermediate coefficient matrix storing the gas concentration distribution;

[0015] Calculating gas concentration distribution based on an improved Gaussian puff model that incorporates dynamic leakage parameters, wind-driven puff center displacement, and dynamic diffusion coefficients, and combining a numerical iterative method to update concentration values;

[0016] Perform standard conversion on the concentration calculation results and output the visual data.

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

[0018] Preferably, the step of obtaining flow and concentration data of underground leakage points varying with time includes:

[0019] Extract discrete time series data of flow rate and concentration at the leakage point from the underground pipeline leakage simulation results;

[0020] The discrete data is checked for rationality, including flow non-negativity check and concentration physical range verification.

[0021] Preferably, the step of generating a calculation grid and adjusting the grid direction based on the wind direction in the environmental parameter includes:

[0022] Generate a uniform three-dimensional grid according to the preset grid size;

[0023] The grid coordinates are rotated according to the wind direction angle. The transformation formula is:

[0024] ;

[0025] in, is the wind direction angle; is the grid coordinate before rotation; is the rotated grid coordinate.

[0026] Preferably, the step of generating a time series according to a preset simulation duration and a dynamic time step strategy includes:

[0027] Divide the total simulation time into at least two consecutive periods, where:

[0028] The first time step is used in the first period , used to capture rapid concentration changes at the initial stage of a leak;

[0029] The second time step is used for the second period ,in , used to reduce redundancy in later calculations;

[0030] The time step of subsequent periods increases gradually according to the preset rules and satisfies ,in The time period number.

[0031] Preferably, the step of interpolating the underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points includes:

[0032] Obtain discrete time series data of flow and concentration from underground leakage simulation results;

[0033] According to the flow rate and concentration values at adjacent time points, the following formula is used to calculate Traffic at any time and concentration :

[0034] ;

[0035] ;

[0036] 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 the site.

[0037] Preferably, the step of dynamically calculating the diffusion coefficient according to the atmospheric stability level and the time series includes:

[0038] A parameterized formula is selected based on the atmospheric stability level. When the time t is less than the preset threshold, the diffusion coefficient is calculated using the first formula.

[0039] When the time t is greater than or equal to the preset threshold, the diffusion coefficient is calculated according to the second formula;

[0040] The coefficients and exponents in the second formula are smaller than those in the first formula.

[0041] Preferably, the step of calculating the gas concentration distribution based on the improved Gaussian puff model includes:

[0042] The improved Gaussian puff model is used to calculate the concentration distribution under unit flow rate. The formula is:

[0043] ;

[0044] 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;

[0045] Unit flow concentration coefficient Precomputed as an intermediate coefficient matrix , and based on real-time traffic Update actual concentration value:

[0046] ;

[0047] in, Indicates time Time grid points The corresponding unit flow concentration coefficient.

[0048] Preferably, the step of calculating the gas concentration distribution based on the improved Gaussian puff model further includes:

[0049] In areas where the concentration gradient exceeds a preset threshold, the RK4 numerical method is used to correct the concentration value, including:

[0050] ;

[0051] in, is the diffusion coefficient tensor; is the Laplace operator of concentration; is the convection term caused by wind speed;

[0052] Update the concentration according to the RK4 iterative formula:

[0053] ;

[0054] in, is the dynamic time step; , , , is the fourth-order slope estimate.

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

[0056] The present invention provides a simulation calculation method for ground gas leakage and diffusion, which has the following beneficial effects:

[0057] 1. This invention significantly improves the accuracy of gas leak diffusion prediction by dynamically adjusting the diffusion coefficient and comprehensively considering multiple environmental factors. It can more accurately reflect the concentration distribution and diffusion range in complex urban environments. This helps to promptly identify high-risk areas, reduce safety hazards caused by prediction errors, and provide more reliable protection for residents' lives and property.

[0058] 2. Compared to the time-consuming nature of traditional computational fluid dynamics methods, this invention utilizes a fast iteration algorithm and a dynamic time-stepping strategy, significantly reducing simulation time. This allows relevant departments to quickly obtain diffusion prediction results after a gas leak accident. This rapid response capability facilitates the development of timely and effective evacuation and disposal plans, minimizing the economic losses and social impact of the accident.

[0059] 3. This invention optimizes the model's adaptability to address the variable wind directions and complex terrain of urban environments, accurately simulating the diffusion of gas in densely built-up areas or uneven terrain. This feature enables urban gas safety managers to more efficiently assess leakage risks, optimize pipeline inspections and emergency response plans, and improve overall management efficiency.

[0060] By converting simulation results into standard volume concentrations and supporting 3D visualization, this method provides emergency response personnel with intuitive maps of diffusion range and concentration distribution. Compared to the abstract data generated by traditional methods, this output simplifies the decision-making process, enabling even non-specialists to quickly understand and apply the prediction results, thus improving the operability of emergency response.

[0061] 5. This invention is not only applicable to single leak scenarios but can also adapt to gas leaks of varying 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 highly valuable for promotion and will help improve energy security in the overall urbanization process. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 Schematic diagram of the method flow of the present invention;

[0063] Figure 2 This is a schematic diagram of the leakage center point of the present invention;

[0064] Figure 3 This is a diagram of the ground diffusion results of an embodiment of the present invention. DETAILED DESCRIPTION

[0065] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the drawings in the present specification. 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 making creative efforts are within the scope of protection of the present invention.

[0066] Please see the attached Figure 1 -Attached Figure 2The 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 puff model.

[0067] like Figure 1 As shown, the ground gas leakage diffusion simulation calculation method may include the following steps:

[0068] S1, initialize environmental parameters, leakage source parameters and calculation domain parameters;

[0069] S2. Obtain flow and concentration data of underground leakage points that change over time;

[0070] S3. Generate a computational grid and adjust the grid direction based on the wind direction in the environmental parameters so that the grid coordinates are aligned with the wind field direction;

[0071] S4, generating a time series according to a preset simulation duration and a dynamic time step strategy;

[0072] S5. interpolating underground leakage flow and concentration data based on the time series to generate dynamic leakage parameters at consecutive time points;

[0073] S6. Dynamically calculate the diffusion coefficient based on the atmospheric stability level and time series;

[0074] S7, initializing the concentration matrix and intermediate coefficient matrix storing the gas concentration distribution;

[0075] S8. Calculate gas concentration distribution based on the improved Gaussian puff model;

[0076] S9. Perform standard condition conversion on the concentration calculation results and output visual data.

[0077] The following is a detailed description of each step in the method of the present invention, which comprehensively explains the specific implementation principles, technical details and processes of each step.

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

[0079] The configuration of environmental parameters includes:

[0080] wind speed Obtain real-time wind speed data in the leak area from meteorological monitoring equipment, or set a fixed value based on historical meteorological data. The value range is in meters per second (m / s), covering typical wind conditions for gas diffusion, for example, 0 to 30 m / s. This wind speed parameter directly impacts the calculation of the downwind displacement of the puff center in subsequent steps.

[0081] wind direction : Defines wind direction angles in radians, relative to true north. For example, true north is 0 radians, and due east is π / 2 radians. This parameter rotates the computational grid to ensure alignment of the grid coordinate system with the wind direction, thus avoiding the diffusion direction deviation caused by the fixed coordinates of the traditional Gaussian model.

[0082] Atmospheric stability class: Determined according to the Pasquill classification system, it ranges from A (very unstable) to F (very stable) based on temperature gradient, solar radiation intensity, and wind speed data. This parameter is used to dynamically select the diffusion coefficient calculation formula to reflect the impact of different atmospheric stratification states on turbulent diffusion.

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

[0084] The configuration of the leakage source parameters includes:

[0085] Initial flow Initial gas mass flow rate, expressed in kilograms per second (kg / s), is output by underground pipeline leak simulation software. This value is determined by pipeline pressure, leak aperture, and gas physical properties. This initial flow rate serves as an input parameter for the Gaussian puff model, determining the initial diffusion intensity of the leak.

[0086] Initial concentration : Indicates the initial volume concentration of the gas at the leak 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%.

[0087] Leak 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 calculations and to locate the initial position of the center of the smoke puff.

[0088] The configuration of the computational domain parameters includes:

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

[0090] Dynamic time step strategy: The time step is set in segments based on the time-varying characteristics of the diffusion process to balance computational accuracy and efficiency. Preferably, a smaller step size (e.g., 0.1 second) is used in the early stages of a leak (0-60 seconds) to capture rapid diffusion processes; a gradually larger step size (e.g., 1 second) is used in the middle stages (60-600 seconds); and a larger step size (e.g., 10 seconds) is used in the late stages (600-3600 seconds) to reduce redundant computations during periods of low dynamics.

[0091] Preferably, the calculation domain parameters support custom area ranges, such as dynamically adjusting the simulation area according to geographic information system (GIS) data to adapt to the needs of different urban environments.

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

[0093] Regarding step S2, in this embodiment, the operation of obtaining the flow and concentration data of the underground leakage point that changes 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.

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

[0095] The data consists of a time series of points (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.

[0096] Before inputting data into the diffusion model, the following checks are performed to ensure that the data conforms to physical laws:

[0097] Traffic non-negativity check: detecting traffic sequence Is there a negative value? , it is set to zero and triggers an alarm signal to avoid model calculation failure due to negative flow.

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

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

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

[0101] For step S3, in this embodiment, the operation of generating a calculation grid and adjusting the grid direction based on the wind direction in the environmental parameters is to construct a calculation 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 offset of the traditional fixed grid.

[0102] The generation of the computational grid is achieved through the following logic:

[0103] Based on the calculation domain parameters defined in step S1 (such as the simulation area range ), generating a uniformly distributed 3D structured grid. Preferably, the grid dimensions in the horizontal (x, y) and vertical (z) directions are set based on the simulation accuracy requirements. For example, the horizontal grid size can be set to 5 meters by 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 that covers the entire simulation area.

[0104] like Figure 2 Shown, as an example:

[0105] With the underground pipeline leakage point as the center, the calculation domain is defined as 3km long ( =3000m), width 3km ( =3000m)、30m high( =30m), covering the main impact areas of gas diffusion in urban environments. The 30m height accounts for the ground diffusion height. A uniform 3D grid is created within the computational domain. Assuming a grid cell size of 5m × 5m × 5m, the entire space is meshed. The number of cells is calculated as follows:

[0106] Number of grids in the X direction: 3000m / 5m = 600;

[0107] Number of grids in the Y direction: 3000m / 5m = 600;

[0108] Number of grids in the Z direction: 30m / 5m = 6;

[0109] The total number of grid cells is 600 × 600 × 6 = 2,160,000 (approximately 2.16 million). Grid size affects computational efficiency. Finer spatial divisions create more grid cells, resulting in higher accuracy but slower computational speed. This directly impacts both the speed and accuracy of the entire model.

[0110] Dynamic adjustment of the grid orientation is achieved through rotation transformation:

[0111] 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:

[0112] ;

[0113] in:

[0114] is the grid node coordinate in the original Cartesian coordinate system;

[0115] is the coordinate in the new coordinate system after rotation.

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

[0117] 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 points and the leakage points are within the ... , , The total number of grids is approximately 2.16 million, which is sufficient to support the rapid iterative calculation of the Gaussian puff model while solving the computational burden of tens of millions of grids in CFD methods.

[0118] Regarding 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 improve the calculation efficiency while ensuring the simulation accuracy by adjusting the time step in stages, so as to adapt to the dynamic characteristics differences at different stages in the gas leakage and diffusion process.

[0119] The construction of the dynamic time step strategy is achieved through the following logic:

[0120] Based on the total simulation time defined in step S1 , the simulation process is divided into multiple stages, each stage uses a different time step.

[0121] Preferably, the phases are divided according to the physical properties of leakage diffusion, including an initial rapid diffusion phase, a mid-term stable diffusion phase, and a late slow decay phase. The time step size of each phase is dynamically adjusted based on the diffusion rate, with a smaller step size in the early phase to capture rapid changes and a larger step size in the later phase to reduce redundant calculations.

[0122] In this embodiment, the time series generation method includes the following steps:

[0123] 1. Stage division and step size setting:

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

[0125] mid-term stage : Use a medium time step , balancing computational accuracy and efficiency.

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

[0127] in, 、 is the preset time threshold, which satisfies .

[0128] 2. Time point generation and merging:

[0129] In each stage, a sequence of equally spaced time points is generated according to the corresponding step size. The formula is:

[0130] ;

[0131] Finally, the time points of all stages are merged into a global time series in ascending order ,in , .

[0132] 3. Time point deduplication and verification:

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

[0134] As an example:

[0135] Simulation total duration setting:

[0136] The simulation duration is initialized to T = 3600 seconds (1 hour) to construct a time series covering the complete leakage diffusion process.

[0137] Stage division and step size setting:

[0138] 1. Initial stage ( Second):

[0139] Step length: Second

[0140] Time point generation: according to the formula Generate a discrete time series consisting of Seconds, total points indivual.

[0141] 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).

[0142] 2. Mid-term stage ( Second):

[0143] Time point generation: according to the formula Generate a sequence containing Seconds, total points indivual.

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

[0145] 3. Late stage ( Second):

[0146] Step length: Second

[0147] Time point generation: according to the formula Generate a sequence containing Seconds, total points indivual.

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

[0149] Time series merging and verification:

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

[0151] Deduplication and gap handling: Check that the intervals between adjacent time points are strictly increasing. If there are gaps (e.g., between phase switching points t=60 seconds and t=600 seconds), insert intermediate time points to ensure continuity.

[0152] Small step sizes in the initial stage of dynamic step size ensure accuracy, while large step sizes in the middle and late stages reduce the computational burden, so that the total simulation time is controlled within seconds to minutes, meeting real-time emergency needs.

[0153] 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, providing dynamic input for the subsequent Gaussian puff model, and ensuring that the diffusion simulation can reflect the real physical process of the leakage parameters changing over time.

[0154] 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 values . Time series must meet strict time order ( ,in = total simulation time), and a data validation 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 strategy is used, such as linear extrapolation of adjacent data or insertion of virtual time points.

[0155] 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 and the midpoint value of the time series, gradually narrowing the search range until the 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.

[0156] Linear interpolation calculations are performed based on the parameter values at adjacent time points. and concentration , and the interpolation formulas are:

[0157] ;

[0158] ;

[0159] in:

[0160] and Represents two adjacent known underground data time points;

[0161] Indicates at a point in time Underground leakage flow at

[0162] Indicates at a point in time Underground leakage flow at

[0163] Indicates at a point in time Underground leakage concentration at

[0164] Indicates at a point in time Underground leakage concentration at the site.

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

[0166] 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., >30s), the error can be reduced by increasing the number of sampling points or using a higher-order interpolation (e.g., quadratic interpolation).

[0167] Regarding 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.

[0168] 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 based on experimental data or historical observation results. For example:

[0169] 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;

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

[0171] 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:

[0172] 1. Initial stage : Using higher coefficients and exponents, for example, for Class A stability, the downwind diffusion coefficient is calculated as:

[0173] ;

[0174] 2. Late Stage : Switch to a lower coefficient and exponent. For example, the adjusted downwind diffusion coefficient is:

[0175] ;

[0176] Parameter switching is achieved through table lookup, and different stability levels correspond to independent segmented parameter tables.

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

[0178]

[0179] The general formula for the three-dimensional diffusion coefficient is expressed as:

[0180] ;

[0181] ;

[0182] ;

[0183] in, , , , is the diffusion intensity coefficient in each direction, , , is the diffusion rate exponent, whose value is dynamically selected based on the stability level. The power function form can simultaneously describe the magnitude (controlled by the coefficient) and the rate (controlled by the exponent) of the diffusion range over time.

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

[0185] 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, crosswind 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.

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

[0187] Regarding 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 accurately simulates the spatiotemporal evolution process of gas diffusion by iteratively updating the grid point concentration values in multiple stages.

[0188] 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:

[0189] ;

[0190] in:

[0191] express Time in space coordinates Gas mass concentration at

[0192] express Dynamic leakage flow at each moment;

[0193] 、 、 represents the diffusion coefficients in the horizontal, longitudinal, and vertical directions, which are the diffusion coefficients dynamically calculated in step S6;

[0194] represents the wind speed, which is defined in step S1;

[0195] Represents the rotated spatial grid coordinates;

[0196] Represents the simulation time, which is the time point in the time series of step S4.

[0197] 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:

[0198] ;

[0199] in, is the grid point index. Through the above, the real-time computation overhead can be significantly reduced, and the concentration update is simplified to a matrix element multiplication operation.

[0200] 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:

[0201] ;

[0202] in, is the diffusion coefficient tensor; is the Laplace operator of concentration; is the convection term caused by wind speed.

[0203] The implementation steps of the RK4 method are as follows:

[0204] 1. Slope estimation:

[0205] Calculate the current time step The slope :

[0206] ;

[0207] in is the right-hand side term of the convection-diffusion equation.

[0208] 2. Calculate the intermediate slopes sequentially ;

[0209] 3. Concentration update:

[0210] Update the concentration value at the next time step using the weighted average slope:

[0211] ;

[0212] Iteratively update the concentration at each grid point There are four estimates, each time step is 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 multiple times Until the next moment Finally, RK4 uses weighted average Updates the concentration value.

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

[0214] The correction process of the concentration data is based on the concentration matrix calculated in step S8. , including the following operations:

[0215] 1. Unit conversion: Convert mass concentration (unit: kg / m 3 ) is converted to volume concentration (unit: ppm).

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

[0217] 2. Outlier filtering: For grid points where the 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.

[0218] Data storage and structured output are achieved through the following methods:

[0219] 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:

[0220] Concentration data layer: stores volume concentration values by time step and spatial coordinates;

[0221] Metadata layer: records simulation parameters (wind speed, stability level, leakage flow sequence);

[0222] Grid description layer: stores the calculated grid coordinates and rotation matrix parameters defined in step S3.

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

[0224] The visual output interface includes the following functions:

[0225] 3D dynamic rendering: Based on OpenGL or WebGL engines, volume concentration data is mapped to a computational grid to generate spatiotemporal evolution animations. Optimally, isosurface rendering (e.g., 1% LEL, 5% LEL isosurfaces) and thermal map overlays are supported.

[0226] 2D Cross-Section Analysis: Provides interactive tools to select any spatial cross-section (such as horizontal or vertical planes) to generate concentration distribution contour maps or pseudo-color maps, and annotate the peak concentration location and diffusion range.

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

[0228] 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 computational grid for rotating aligned wind fields, dynamic time step strategy, leakage data interpolation, segmented diffusion coefficient calculation, pre-calculation matrix optimization, hybrid simulation calculation combining an improved Gaussian puff model with the 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 computational efficiency through matrix pre-calculation, dynamic parameter switching and parallel optimization strategy. It finally outputs the volume concentration distribution and interactive visualization results under standard working conditions, providing accurate data support for real-time monitoring, risk assessment and emergency decision-making of gas leakage accidents.

[0229] Example:

[0230] This embodiment provides an implementation example of a gas leakage and diffusion simulation method. The feasibility and efficiency of the method are verified through actual leakage scenarios and computational performance tests, as follows:

[0231] Implementation scenarios and parameter configuration:

[0232] A buried gas pipeline leak accident was simulated. The pipeline operating pressure was 1 MPa, the leak aperture diameter was 100 mm, the pipeline was buried at a depth of 5 m, and the leak direction was vertically upward. Simulation environmental parameters included northeasterly wind direction (3.5 m / s wind speed, 45° wind angle), atmospheric stability level D, and a computational domain covering a spatial area of 3 km × 3 km × 30 m. The leak simulation was divided into two phases, aboveground and underground, each lasting 300 seconds, for a total simulation time of 3600 seconds.

[0233] The implementation steps are as follows:

[0234] 1. Environment initialization and grid construction:

[0235] Rotate the computational domain according to the wind direction angle of 45° to align the grid with the wind direction and eliminate coordinate deviation.

[0236] 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;

[0237] 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 approximately 1440.

[0238] 2. Leakage data and diffusion parameter calculation:

[0239] Underground 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);

[0240] Based on the 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.

[0241] 3. Pre-calculation and simulation execution:

[0242] Initialize the four-dimensional concentration matrix (time × space) and the intermediate coefficient matrix to store the concentration contribution parameters under unit flow;

[0243] An improved Gaussian puff model is used to calculate concentration distribution, combined with pre-calculated matrices to achieve real-time simulation (e.g., processing approximately 100,000 grid nodes per second), and supports switching to the RK4 method to handle dynamic diffusion processes;

[0244] Perform standard condition conversion on the calculation results, convert the mass concentration into volume concentration (unit: ppm), and limit the ground boundary concentration to not exceed the lower flammable limit (such as 5% volume concentration).

[0245] 4. Visualization and output:

[0246] Generate a visualization of the ground diffusion results as shown below Figure 3 shown.

[0247] Computational performance verification:

[0248] In the hardware environment of the 12th generation Intel i7-12700H processor (20 cores) and 32GB of memory, the operating efficiency under different grid densities was tested:

[0249] 10-meter grid (approximately 270,000 nodes): The entire process runs in about 1.7 minutes, suitable for quick preliminary assessments;

[0250] 5-meter grid (approximately 2.16 million nodes): Run time is approximately 2.1 minutes, and the accuracy meets real-time emergency response requirements;

[0251] Finer grids (e.g., 3 meters, approximately 33.75 million nodes): Run time increases to approximately 4 minutes, making it suitable for high-precision post-mortems.

[0252] This example achieves minute-level simulation of leak diffusion at a 5-meter grid resolution through dynamic parameter adjustment (time step, diffusion coefficient) and pre-calculation optimization. It accurately depicts the diffusion path and concentration attenuation trend of the gas from the leak point to the northeast, 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.

[0253] While 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 these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A method for simulating the diffusion of ground gas leakage, 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 leak points 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; Calculating gas concentration distribution based on an improved Gaussian puff model that incorporates dynamic leakage parameters, wind-driven puff center displacement, and dynamic diffusion coefficients, and combining a numerical iterative method to update concentration values; Perform standard condition conversion on the concentration calculation results and output visual data; 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 underground pipeline leakage simulation results; Performing rationality checks on the discrete data, including flow non-negativity checks and concentration physical range verification; The step of generating a time series according to a preset simulation duration and a dynamic time step strategy includes: Divide the total simulation time into at least two consecutive periods, where: The first time period uses the first time step Δt1 to capture the rapid concentration changes at the initial stage of leakage; The second time period uses the second time step Δt2, where Δt2>Δt1, to reduce redundancy in later calculations; The time step of the subsequent period is gradually increased according to the preset rules, and satisfies Δt n+1 ≥Δt n , where n is the time period number; The step of dynamically calculating the diffusion coefficient according to the atmospheric stability level and the time series comprises: A parameterized formula is selected based on the atmospheric stability level. When the time t is less than the preset threshold, the diffusion coefficient is calculated using 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.

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 generating a calculation grid and adjusting the grid direction based on the wind direction in the environmental parameter includes: 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: Among them, θ is the wind direction angle; (x, y) is the grid coordinate before rotation; (x ′ ,y ′ ) are the rotated grid coordinates.

4. 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 includes: Obtain discrete time series data of flow and concentration from underground leakage simulation results; According to the flow and concentration values at adjacent time points, the flow q(t) and concentration c(t) at time t are calculated using the following formula: Among them, t i and t i+1 represents two adjacent known underground data time points; q(t i ) represents the time point t i Underground leakage flow at q(t i+1 ) represents the time point t i+1 Underground leakage flow at the location; c(t i ) represents the time point t i Underground leakage concentration at i+1 ) represents the time point t i+1 Underground leakage concentration at the site.

5. 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: Where c(t,x,y,z) represents the gas mass concentration at the spatial coordinate (x,y,z) at time t; q(t) represents the dynamic leakage flow at time t; σ x , σ y , σ z represents the diffusion coefficient in the horizontal, longitudinal and vertical directions; u represents the wind speed; x, y, z represent the spatial grid coordinates after rotation; t represents the simulation time; The unit flow concentration coefficient c / q(t) is pre-calculated as the intermediate coefficient matrix M, and the actual concentration value is updated according to the real-time flow q(t): c(t,x,y,z)=q(t)·M[t,x,y,z] Where M[t,x,y,z] represents the unit flow concentration coefficient corresponding to the grid point (x,y,z) at time t.

6. The ground gas leakage diffusion simulation calculation method according to claim 5 is characterized in that: The step of calculating the gas concentration distribution based on the improved Gaussian puff model further includes: In areas where the concentration gradient exceeds a preset threshold, the RK4 numerical method is used to correct the concentration value, including: Where D 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: Where Δt is the dynamic time step; k1, k2, k3, k4 are the fourth slope estimates.

7. The ground gas leakage diffusion simulation calculation method according to claim 1 is 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