Self-adaptive grid dynamic encryption method for underground water pollution migration simulation

Abstract

The invention discloses a self-adaptive grid dynamic encryption method for underground water pollution migration simulation, which belongs to the technical field of underground water pollution migration simulation, and comprises the following steps of: forming 0.05-meter super-resolution grids at the periphery of a dispersed pollution source through a local grid dynamic encryption technology triggered by a pollution source intensity gradient threshold; the pollutant frontal surface tracking error is greatly reduced, and the high-precision pollutant migration simulation capability is realized. Moreover, based on an intelligent grid reconstruction algorithm of a flux threshold, the number of calculation units is reduced by 62%, the time consumed by single simulation is shortened to 3.3 hours, and the overall mass conservation error is kept to be smaller than 1% while the calculation efficiency is improved. Meanwhile, for a heterogeneous aquifer, through a pollution flux and stratum permeability coefficient dual-drive encryption strategy, the solute transport prediction precision of the model at the boundary of a clay lenticular body is greatly improved, and the problems that in the prior art, the simulation precision of the periphery of a pollution source is insufficient, and the calculation efficiency is low due to global encryption are effectively solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of groundwater pollution migration simulation, and in particular relates to an adaptive grid dynamic encryption method for groundwater pollution migration simulation. Background Art

[0002] The Pearl River Delta region exhibits a fragmented distribution pattern characterized by "one village, multiple enterprises, distributed along roads," resulting in a diffuse pollution emission pattern characterized by multiple locations, low intensity, and wide coverage. Unlike centralized pollution sources, these dispersed sources exhibit significant spatiotemporal heterogeneity: 1) Spatially, 63% of these point sources are located in shallow aquifers 5-15 meters below the surface, influenced by hidden pathways such as abandoned seepage wells and damaged pipelines, forming a pollution network that spans administrative boundaries. 2) In terms of emission intensity, pollutant fluxes fluctuate by two to three orders of magnitude due to the intermittent production of small and medium-sized enterprises in rural towns (annual operating rates fluctuate by ±35%). 3) In terms of migration paths, the interweaving of high-permeability sand layers and clay lenses in the Pearl River Delta's sedimentary layers creates a complex three-dimensional heterogeneous seepage field. This complex pollution pattern of "dispersed point sources, hidden pathways, and heterogeneous media" results in pollutant capture rates of less than 35% using traditional monitoring networks, and even lower than 28% in rural industrial parks. There is an urgent need to develop high-resolution numerical simulation techniques tailored to the characteristics of dispersed pollution sources.

[0003] Current groundwater pollution models generally use static structured grids. However, their application in rural industrial parks in the Pearl River Delta faces three major bottlenecks:

[0004] 1. In terms of spatial resolution contradiction, a uniform grid across the entire region must balance computational efficiency and accuracy, with a typical grid spacing of 1-5 meters. However, the gradient variation area of ​​the actual pollution plume (e.g., the 2-meter range around the pollution source) requires a resolution of 0.1 meters, resulting in a local mass conservation error of 12-18%.

[0005] 2. In terms of computing resource waste, in order to capture the migration details of dispersed pollution sources, the fully encrypted grid will cause the number of computing units to surge 30 times, and a single simulation will take more than 72 hours, far exceeding the time requirements of actual engineering applications.

[0006] 3. Existing technologies lack a dynamic response mechanism. Faced with sudden changes in pollution flux caused by production fluctuations (such as a fluctuation of more than 200% in the concentration of benzene series within 24 hours), the fixed grid model has a delayed response of more than 48 hours, resulting in a risk assessment deviation of more than 40%, which seriously restricts the effectiveness of pollution warning and emergency response. Summary of the Invention

[0007] The purpose of the present invention is to provide an adaptive grid dynamic encryption method for groundwater pollution migration simulation, which dynamically adjusts the grid density to Implementing local encryption to 0.05 meters can greatly improve the accuracy of front tracking; triggering grid reconstruction based on pollutant flux thresholds can significantly reduce the overall computing units.

[0008] The technical solution adopted in the present invention is as follows:

[0009] An adaptive grid dynamic refinement method for groundwater pollution migration simulation includes the following steps:

[0010] (1) Based on the real-time monitoring data obtained by the IoT perception layer, the DBSCAN clustering algorithm is used to identify discrete pollution source clusters, and three levels of encryption zones are established with the cluster centroid as the origin: the pollution source core zone, the transition zone, and the peripheral zone;

[0011] (2) By calculating the pollutant concentration gradient value of each grid unit in real time, when the pollutant concentration gradient value is detected to be greater than the threshold When , the local grid encryption operation is triggered, and the encryption depth is controlled by an iterative formula, constraining the adjacent grid size change rate to ≤ 50% to avoid numerical oscillation;

[0012] (3) A dual-driving factor weight model of pollution flux and formation permeability coefficient was established. The variance analysis of historical data was performed through orthogonal test method to determine the weight ratio of pollution flux and permeability coefficient. The dual-driving factor weight model dynamically adjusted the grid density based on the encryption priority index P calculated in real time. The pH, electrical conductivity EC and heavy metal concentration data were used to establish a grid update cycle control mechanism of no more than 6 hours. When the encryption priority index P exceeded the preset threshold, the full grid reconstruction was triggered. The Delaunay triangulation was used to generate an unstructured grid. Local encryption was implemented at the unit boundary where the permeability coefficient difference was greater than 100 times. The Laplace smoothing algorithm was used to optimize the grid quality so that the Jacobian matrix condition number was less than 10 to ensure numerical stability.

[0013] Furthermore, in step (1), the three-level encryption zone is established with the cluster centroid as the origin, specifically: the radius of the pollution source core area is 20 meters, and the initial grid size is 0.05 meters; the radius of the transition area is 50 meters, and the initial grid size is 0.2 meters; the radius of the peripheral area is 200 meters, and the initial grid size is 1 meter; the background area uses a basic grid of 5 meters × 5 meters × 2 meters for full coverage.

[0014] Furthermore, in step (2), the maximum encryption depth is limited to level 5, corresponding to a minimum grid size of 0.05 meters.

[0015] Furthermore, in step (2), the specific process of calculating the pollutant concentration gradient is as follows:

[0016] Based on the spatial discretization principle of the finite volume method, the formula for calculating the pollutant concentration gradient is as follows:

[0017]

[0018] wherein, represents the magnitude of the pollutant concentration gradient at the three-dimensional grid point (i,j,k), reflecting the degree of spatial variation of the pollutant concentration at the point, C represents the pollutant concentration per unit volume at a point (x,y,z) in space, respectively represent the rate of change of the pollutant concentration in the x and y horizontal directions, represents the rate of change of the pollutant concentration in the vertical z direction;

[0019] In the Cartesian coordinate system, the calculation domain is discretized, and the concentration value C of each grid point (i,j,k) is obtained by interpolation of the values of its six adjacent cells. The partial derivatives in each direction are calculated using the central difference method:

[0020]

[0021] wherein, C i+1,j,k , C i-1,j,k represent the concentration values at the two grid points adjacent to the current grid point in the x direction, Δx represents the grid spacing in the x direction, C i,j+1,k , C i,j-1,k represent the concentration values at the two grid points adjacent to the current grid point in the y direction, Δy represents the grid spacing in the y direction, C i,j,k+1 , C i,j,k-1 represent the concentration values at the two grid points adjacent to the current grid point in the z direction, Δz represents the grid spacing in the z direction.

[0022] Further, in the step (2), the specific process of calculating the encryption depth is:

[0023] The encryption depth has a logarithmic relationship with the amplitude of the concentration gradient threshold value:

[0024]

[0025] wherein, L enc represents the strength of the encryption operation, which is a dimensionless integer; is the critical gradient value triggering the encryption operation; the relationship between the encryption depth L and the grid size is defined as:

[0026] h L = h0 / 2 L

[0027] wherein, h L represents the grid cell size after L times of encryption operation, which is a direct measure of spatial discretization, h0 represents the reference grid size before encryption, and L represents the iteration number of grid encryption;

[0028] According to the numerical stability requirements, the encrypted grid must satisfy the local Péclet number < 2, that is:

[0029]

[0030] Where v is the pore water velocity and D is the diffusion coefficient. The required infill depth can be obtained by combining the above equations:

[0031]

[0032] Establish an exponential relationship between the gradient ratio and the encryption depth, introduce a rounding function to ensure that the encryption depth is an integer, and set the maximum encryption depth L max =5, to prevent computing resource exhaustion caused by excessive encryption.

[0033] Furthermore, in step (3), the specific formula of the dual-driving factor weight model is:

[0034]

[0035] Among them, P represents the grid encryption priority index, ranging from 0 to 1, and Q represents the real-time pollution flux. max represents the maximum pollution flux, K represents the formation permeability coefficient, K max Indicates the maximum permeability coefficient of the formation;

[0036] According to Darcy's law and Fick's law,

[0037]

[0038] Where n is the porosity, A is the cross-sectional area, v is the Darcy velocity, D is the diffusion tensor, and C is the pollutant concentration. The formation database is normalized by the range.

[0039]

[0040] Among them, K norm represents the normalized formation permeability coefficient, K represents the formation permeability coefficient, K min Indicates the minimum permeability coefficient of the formation, K max Indicates the maximum permeability coefficient of the formation;

[0041] The pollution flux weight, permeability coefficient weight, update cycle, and encryption threshold are orthogonalized, and the simulation accuracy is used as the response variable. The optimal combination is obtained through range analysis:

[0042]

[0043] Finally, take w Q =0.6, w K =0.4, where R QIndicates the extreme value of the pollution flux weight factor, R K Indicates the range of the permeability coefficient weight factor, w Q Indicates the final distribution value of pollution flux weight, w K Indicates the final assigned value of the permeability coefficient weight.

[0044] Furthermore, in step (3), the grid real-time update cycle control mechanism is specifically as follows:

[0045] The grid reconstruction period Δt is determined by the fluctuation rate of the monitoring data:

[0046]

[0047] Wherein, ΔC represents the change in pollutant concentration, C0 represents the initial concentration reference value, ΔK represents the change in permeability coefficient, and K0 represents the initial permeability coefficient;

[0048] The key parameter change rate thresholds are defined as the concentration change rate threshold α < 15% / h, that is, the pollutant concentration change rate does not exceed 15% per hour, and the permeability coefficient change rate β < 10% / h, that is, the permeability coefficient change rate does not exceed 10% per hour. The maximum change rate is taken according to the most unfavorable principle:

[0049]

[0050] Among them, γ represents the maximum change rate, which is used to comprehensively evaluate the most unfavorable changes in concentration and permeability coefficient. The larger value of the relative change rate of the two is taken. When γ exceeds the threshold α = 15% / h or β = 10% / h, an early warning is triggered and the parameters are adjusted. C t 、C t-1 Represent the concentration values ​​at the current time t and the previous time t-1, respectively. Reflects the intensity of concentration fluctuation, K t , K t-1 Represent the medium penetration capacity at the current moment t and the previous moment t-1, quantify the degree of penetrant mutation;

[0051] The simultaneous stability condition γ·Δt≤α yields: Δt≤α / γ. To prevent excessive calculation frequency due to a small Δt, an upper limit of 6 hours is set.

[0052] Furthermore, in step (3), the specific process of optimizing the grid quality is as follows:

[0053] The Jacobian condition number is used to quantify the degree of mesh cell distortion:

[0054]

[0055] Where K(J) represents the Jacobian matrix condition number, which is used to measure the degree of mesh unit distortion. When K(J)>10, the mesh needs to be repaired. J is the Jacobian matrix, which represents the mapping quality of the mesh unit from the parameter space to the physical space. max Represents the maximum singular value, reflecting the maximum tensile deformation of the grid unit, σ min Represents the minimum singular value, which is used to detect the risk of grid degradation;

[0056] Define the unit mapping function: map the physical space unit (x, y, z) to the reference unit (ξ, η, ζ), and the Jacobian matrix is:

[0057]

[0058] Compute the singular value decomposition (SVD),

[0059] J=U∑V T

[0060] ∑=diag(σ1,σ2,σ3)

[0061] Where U is the left singular vector matrix, V is the right singular vector matrix, ∑ is the singular value matrix, σ i represents the i-th singular value.

[0062] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0063] 1. In the present invention, the pollution source intensity gradient The triggered local grid dynamic encryption technology forms a 0.05-meter super-resolution grid around scattered pollution sources (traditional methods are ≥1 meter), reducing the pollutant front tracking error from 18.7% to ≤4.8%, and has the ability to simulate pollutant migration with high precision.

[0064] 2. In the present invention, the intelligent grid reconstruction algorithm based on flux threshold can reduce the number of computing units by 62%, and the time required for a single simulation is compressed to 3.3 hours (72 hours for the traditional method). While improving the computing efficiency, the overall mass conservation error is kept below 1%.

[0065] 3. In this invention, for the typical heterogeneous aquifer in the Pearl River Delta (with permeability coefficients varying by three orders of magnitude), a dual-drive encryption strategy is used to improve the solute transport prediction accuracy of the model at the clay lens boundary by 76%. DETAILED DESCRIPTION

[0066] In order to make the objects, technical solutions and advantages of the present application clearer, the technical solutions in the embodiments of the present application are described clearly and completely below in combination with embodiments. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Meanwhile, the specific embodiments described herein are only used to explain the present application, and are not used to limit the present application.

[0067] An adaptive grid dynamic encryption method for groundwater pollution migration simulation, comprising the following steps:

[0068] (1) Firstly, based on the real-time monitoring data obtained by the perception layer of the Internet of Things, the DBSCAN clustering algorithm (neighborhood radius ε=50 meters, minimum sample number MinPts=3) is used to identify the discrete pollution source cluster, and a three-level encryption area is established with the cluster center as the origin: core area (radius 20 meters, initial grid size 0.05 meters), transition area (radius 50 meters, 0.2 meters), and peripheral area (radius 200 meters, 1 meter). The background area is covered by the basic grid of 5 meters x 5 meters x 2 meters (X x Y x Z).

[0069] (2) Then, the pollution concentration gradient value of each grid unit is calculated in real time, and when it is detected that the gradient exceeds the threshold , the local grid encryption operation is triggered, and the encryption depth is controlled by the iterative formula, and the maximum encryption depth is limited to 5 levels, corresponding to the minimum grid size of 0.05 meters, and the size change rate of adjacent grids is limited to ≤50% to avoid numerical oscillation.

[0070] Specifically, the pollution concentration gradient is the core criterion for triggering grid encryption, which is calculated based on the spatial discretization principle of the finite volume method:

[0071]

[0072] Wherein, represents the pollution concentration gradient at the three-dimensional grid point (i, j, k), reflecting the spatial variation of the pollution concentration at the point. The larger the calculation result is, the faster the concentration changes in different directions. C represents the pollution concentration per unit volume at a point (x, y, z) in space. represents the change rate of the pollution concentration in the x and y horizontal directions, respectively. represents the change rate of the pollution concentration in the vertical (z) direction.

[0073] In the Cartesian coordinate system, the calculation domain is discretely structured, and the concentration value C of each grid point (i, j, k) is obtained by interpolating the values of its adjacent 6 units, and the central difference method is used to calculate the partial derivative in each direction,

[0074]

[0075] C i+1,j,k 、C i-1,j,k Indicates the concentration values ​​at the two grid points adjacent to the current grid point in the x direction. Δx represents the grid spacing (step size) along the x direction. C i,j+1,k 、C i,j-1,k Indicates the concentration value at the two grid points adjacent to the current grid point in the y direction. Δy represents the grid spacing (step length) along the y direction. C i,j,k+1 、C i,j,k-1 Indicates the concentration values ​​at the two grid points adjacent to the current grid point in the z direction. Δz represents the grid spacing (step size) along the z direction.

[0076] Specifically, the density depth is logarithmically related to the magnitude of the concentration gradient exceeding the threshold:

[0077]

[0078] Among them, L enc Indicates the intensity or spatial resolution level of the encryption operation. It is usually a dimensionless integer. The larger the value, the higher the precision of the calculation required for the area. It is the critical gradient value that triggers the encryption operation, which is set by the specific problem.

[0079] Define the relationship between encryption depth L and grid size:

[0080] h L =h0 / 2 L

[0081] Where h L represents the size of the mesh element after L refinement operations and is a direct measure of spatial discretization. h0 represents the baseline mesh size before refinement. L represents the number of mesh refinement iterations.

[0082] According to the numerical stability requirements, the refined mesh must satisfy the local Péclet number < 2:

[0083]

[0084] Where v is the pore water velocity and D is the diffusion coefficient. Combining the above two equations, the required infill depth can be obtained:

[0085]

[0086] Establish an exponential relationship between the gradient ratio and the encryption depth, introduce a rounding function to ensure that the encryption depth is an integer, and set the maximum encryption depth L max =5, to prevent computing resource exhaustion caused by excessive encryption.

[0087] (3) In response to the dual complexity of the heterogeneous aquifer and dynamic pollution sources in the Pearl River Delta, a spatiotemporal multi-dimensional grid optimization system was constructed. First, a dual-driving factor weight model of pollution flux and formation permeability was established, in which the pollution flux weight accounted for 60% and the permeability weight accounted for 40%. This ratio was determined by variance analysis of 100 sets of historical data using the orthogonal test method. The model dynamically adjusts the grid density based on the real-time calculated encryption priority index P. When P>0.8, a full grid reconstruction is triggered. The pH, electrical conductivity (EC), and heavy metal concentration data are updated every 6 hours to drive parameter iteration. For complex geological structures such as clay lenses, Delaunay triangulation is used to generate unstructured grids. Local encryption is implemented at the unit boundaries where the permeability coefficient difference is greater than 100 times. The Laplace smoothing algorithm is used to optimize the grid quality, making the Jacobian matrix condition number <10 to ensure numerical stability.

[0088] Specifically, the dual-driving factor priority index formula is:

[0089] The encryption priority index P quantitatively characterizes the grid:

[0090]

[0091] Among them, P represents the grid encryption priority index, ranging from 0 to 1, and Q represents the real-time pollution flux. max represents the maximum pollution flux, K represents the formation permeability coefficient, K max Indicates the maximum permeability of the formation.

[0092] According to Darcy's law and Fick's law,

[0093]

[0094] Where n is the porosity, A is the cross-sectional area of ​​the flow path, v is the Darcy velocity, D is the diffusion tensor, and C is the pollutant concentration.

[0095] The Pearl River Delta stratigraphic database was normalized.

[0096]

[0097] Among them, K norm represents the normalized formation permeability coefficient, K represents the formation permeability coefficient, K min Indicates the minimum permeability coefficient of the formation, K max Indicates the maximum permeability of the formation.

[0098] The four factors and three levels are orthogonalized (the four factors are: pollution flux weight, permeability coefficient weight, update cycle, and encryption threshold), and the simulation accuracy is used as the response variable. The optimal combination is obtained through range analysis:

[0099]

[0100] Finally, take w Q =0.6, w K =0.4. Where R Q Indicates the extreme value of the pollution flux weight factor, R K Indicates the range of the permeability coefficient weight factor, w Q Indicates the final distribution value of pollution flux weight, w K Indicates the final assigned value of the permeability coefficient weight.

[0101] (2) Grid quality evaluation indicators

[0102] The Jacobian condition number is used to quantify the degree of mesh cell distortion:

[0103]

[0104] Where K(j) represents the Jacobian matrix condition number, which is used to measure the degree of mesh unit distortion. When K(J)>10, the mesh needs to be repaired. J is the Jacobian matrix, which represents the mapping quality of the mesh unit from parameter space to physical space. σ max Represents the maximum singular value, reflecting the maximum tensile deformation of the grid unit. min Represents the minimum singular value, which is used to detect the risk of grid degradation (if σ min →0 will cause the calculation to diverge).

[0105] Define the unit mapping function: map the physical space unit (x, y, z) to the reference unit (ξ, η, ζ), and the Jacobian matrix is:

[0106]

[0107] Compute the singular value decomposition (SVD),

[0108] J=U∑C T

[0109] ∑=diag(σ1,σ2,σ3)

[0110] Where U is the left singular vector matrix, V is the right singular vector matrix, ∑ is the singular value matrix, σ i represents the i-th singular value.

[0111] (3) Real-time grid update cycle formula

[0112] The grid reconstruction period Δt is determined by the fluctuation rate of monitoring data.

[0113]

[0114] Where ΔC represents the change in pollutant concentration, C0 represents the initial concentration reference value, ΔK represents the change in permeability coefficient, and K0 represents the initial permeability coefficient.

[0115] The key parameter change rate thresholds are defined as the concentration change rate threshold α < 15% / h (the pollutant concentration change rate does not exceed 15% per hour), the permeability coefficient change rate β < 10% / h (the permeability coefficient change rate does not exceed 10% per hour), and the maximum change rate is taken according to the most unfavorable principle:

[0116]

[0117] Where γ represents the maximum rate of change, which is used to comprehensively evaluate the most unfavorable changes in concentration and permeability coefficient, taking the larger value of the relative change rate of the two. When γ exceeds the threshold α = 15% / h or β = 10% / h, an early warning is triggered and parameters are adjusted. t 、C t-1 Represent the concentration values ​​at the current time t and the previous time t-1 respectively. Reflects the intensity of concentration fluctuations. A larger value indicates a more intense pollution diffusion or reaction. t , K t-1 Represent the medium penetration capacity at the current moment t and the previous moment t-1, Quantifies the degree of permeability mutation. Large values ​​may be caused by medium blockage or crack expansion.

[0118] The simultaneous stability condition γ·Δt≤α yields: Δt≤α / γ. To prevent excessive calculation frequency due to a small Δt, an upper limit of 6 hours is set.

[0119] The above are the embodiments of the present invention. The foregoing are the preferred embodiments of the present invention. If the preferred implementation methods in each preferred embodiment are not obviously self-contradictory or based on a preferred implementation method, each preferred implementation method can be arbitrarily superimposed and used in combination. The embodiments and the specific parameters in the embodiments are only for the purpose of clearly describing the verification process of the invention, and are not intended to limit the scope of patent protection of the present invention. The scope of patent protection of the present invention is still subject to its claims. Any equivalent structural changes made by using the contents of the description of the present invention should also be included in the scope of protection of the present invention.

Claims

1. An adaptive grid dynamic encryption method for groundwater pollution migration simulation, characterized in that: The following steps are involved: (1) Based on the real-time monitoring data obtained by the IoT perception layer, the DBSCAN clustering algorithm is used to identify discrete pollution source clusters, and three levels of encryption zones are established with the cluster centroid as the origin: the pollution source core zone, the transition zone, and the peripheral zone; (2) By calculating the pollutant concentration gradient value of each grid unit in real time, when the pollutant concentration gradient value is detected to be greater than the threshold When , the local grid encryption operation is triggered, and the encryption depth is controlled by an iterative formula, constraining the adjacent grid size change rate to ≤ 50% to avoid numerical oscillation; (3) A dual-driving factor weight model of pollution flux and formation permeability coefficient was established. The variance analysis of historical data was performed through orthogonal test method to determine the weight ratio of pollution flux and permeability coefficient. The dual-driving factor weight model dynamically adjusted the grid density based on the encryption priority index P calculated in real time. The pH, electrical conductivity EC and heavy metal concentration data were used to establish a grid update cycle control mechanism of no more than 6 hours. When the encryption priority index P exceeded the preset threshold, the full grid reconstruction was triggered. The Delaunay triangulation was used to generate an unstructured grid. Local encryption was implemented at the unit boundary where the permeability coefficient difference was greater than 100 times. The Laplace smoothing algorithm was used to optimize the grid quality so that the Jacobian matrix condition number was less than 10 to ensure numerical stability.

2. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (1), a three-level encryption zone is established with the cluster centroid as the origin, specifically: the radius of the pollution source core area is 20 meters, and the initial grid size is 0.05 meters; the radius of the transition area is 50 meters, and the initial grid size is 0.2 meters; the radius of the peripheral area is 200 meters, and the initial grid size is 1 meter; the background area uses a basic grid of 5 meters × 5 meters × 2 meters for full coverage.

3. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (2), the maximum encryption depth is limited to 5 levels, corresponding to a minimum grid size of 0.05 meters.

4. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (2), the specific process of calculating the pollutant concentration gradient is as follows: Based on the spatial discretization principle of the finite volume method, the formula for calculating the pollutant concentration gradient is as follows: in, It represents the pollutant concentration gradient at the three-dimensional grid point (i, j, k), reflecting the intensity of the spatial variation of the pollutant concentration at that point. C represents the pollutant concentration per unit volume at a point (x, y, z) in space. represent the rate of change of pollutant concentration in the x and y horizontal directions, respectively. It represents the rate of change of pollutant concentration in the vertical z direction; In the Cartesian coordinate system, the computational domain is structured and discretized. The concentration value C of each grid point (i, j, k) is obtained by interpolating the values ​​of its six adjacent cells. The partial derivatives in each direction are calculated using the central difference method: Among them, C i+1,j,k 、C i-1,j,k represents the concentration value at two grid points adjacent to the current grid point in the x direction, Δx represents the grid spacing along the x direction, C i,j+1,k 、C i,j-1,k represents the concentration value at two grid points adjacent to the current grid point in the y direction, Δy represents the grid spacing along the y direction, C i,j,k+1 、C i,j,k-1 represents the concentration values ​​at the two grid points adjacent to the current grid point in the z direction, and Δz represents the grid spacing along the z direction.

5. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 or 3, characterized in that: In step (2), the specific process of calculating the encryption depth is as follows: The density depth is logarithmically related to the magnitude of the concentration gradient exceeding the threshold: Among them, L enc Indicates the strength of the encryption operation, which is a dimensionless integer; Is the critical gradient value that triggers the encryption operation; defines the relationship between the encryption depth L and the grid size: h L =h0 / 2 L Among them, h L It represents the size of the grid unit after L encryption operations, which is a direct measure of spatial discretization. h0 represents the base grid size before encryption, and L represents the number of iterations of grid encryption. According to the numerical stability requirements, the encrypted grid must satisfy the local Péclet number < 2, that is: Where v is the pore water velocity and D is the diffusion coefficient. The required infill depth can be obtained by combining the above equations: Establish an exponential relationship between the gradient ratio and the encryption depth, introduce a rounding function to ensure that the encryption depth is an integer, and set the maximum encryption depth L max =5, to prevent computing resource exhaustion caused by excessive encryption.

6. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (3), the specific formula of the dual-driving factor weight model is: Among them, P represents the grid encryption priority index, ranging from 0 to 1, and Q represents the real-time pollution flux. max represents the maximum pollution flux, K represents the formation permeability coefficient, K max Indicates the maximum permeability coefficient of the formation; According to Darcy's law and Fick's law, Among them, n is porosity, A is flow cross-sectional area, v is pore water velocity, D is diffusion coefficient, and C is pollutant concentration. The formation database is normalized by range, Among them, K norm represents the normalized formation permeability coefficient, K represents the formation permeability coefficient, K min Indicates the minimum permeability coefficient of the formation, K max Indicates the maximum permeability coefficient of the formation; The pollution flux weight, permeability coefficient weight, update cycle, and encryption threshold are orthogonalized, and the simulation accuracy is used as the response variable. The optimal combination is obtained through range analysis: R Q ＝0.32， Finally, take w Q =0.6, w K =0.4, where R Q Indicates the range of the pollution flux weight factor, R K Indicates the range of the permeability coefficient weight factor, w Q Indicates the final distribution value of pollution flux weight, w K Indicates the final assigned value of the permeability coefficient weight.

7. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (3), the grid real-time update cycle control mechanism is specifically as follows: The grid reconstruction period Δt is determined by the fluctuation rate of the monitoring data: Wherein, ΔC represents the change in pollutant concentration, C0 represents the initial concentration reference value, ΔK represents the change in permeability coefficient, and K0 represents the initial permeability coefficient; The key parameter change rate thresholds are defined as the concentration change rate threshold α < 15% / h, that is, the pollutant concentration change rate does not exceed 15% per hour, and the permeability coefficient change rate β < 10% / h, that is, the permeability coefficient change rate does not exceed 10% per hour. The maximum change rate is taken according to the most unfavorable principle: Among them, γ represents the maximum change rate, which is used to comprehensively evaluate the most unfavorable changes in concentration and permeability coefficient. The larger value of the relative change rate of the two is taken. When γ exceeds the threshold α = 15% / h or β = 10% / h, an early warning is triggered and the parameters are adjusted. C t 、C t-1 Represent the concentration values ​​at the current time t and the previous time t-1 respectively, Reflects the intensity of concentration fluctuation, K t , K t-1 Represent the medium penetration capacity at the current moment t and the previous moment t-1, quantify the degree of penetrant mutation; The simultaneous stability condition γ·Δt≤α yields: Δt≤α / γ. To prevent Δt from being too small, which leads to an excessively high calculation frequency, an upper limit of 6 hours is set.

8. The adaptive grid dynamic encryption method for groundwater pollution migration simulation according to claim 1 is characterized in that: In step (3), the specific process of optimizing the grid quality is as follows: The Jacobian condition number is used to quantify the degree of mesh cell distortion: Where K(J) represents the Jacobian matrix condition number, which is used to measure the degree of mesh unit distortion. When K(J)>10, the mesh needs to be repaired. J is the Jacobian matrix, which represents the mapping quality of the mesh unit from the parameter space to the physical space. max Represents the maximum singular value, reflecting the maximum tensile deformation of the grid unit, σ min Represents the minimum singular value, which is used to detect the risk of grid degradation; Define the unit mapping function: map the physical space unit (x, y, z) to the reference unit (ξ, η, ζ), and the Jacobian matrix is: Compute the singular value decomposition (SVD), J=U∑V T ∑=diag(σ1,σ2,σ3) Where U is the left singular vector matrix, V is the right singular vector matrix, ∑ is the singular value matrix, σ i represents the i-th singular value.

Images