Method for improving simulation efficiency of interlayer body with compression delay under crustal stress change principle

Abstract

The invention relates to a method for improving the simulation efficiency of an interlayer body with compression time delay under a crustal stress change principle, and the method comprises the steps: carrying out the half-thickness grid spatial discretization of the compression time delay interlayer body through a rectangular grid, and obtaining a simulation grid unit of the compression time delay interlayer body; constructing a water balance model according to the compressed time-delay interlayer body simulation grid unit; and combining the water balance model of each node unit to obtain a difference equation set, and simulating compaction or expansion of each node unit in the compression time-delay interlayer body by using the difference equation set. According to the method, the simulation efficiency of the interlayer body with the compression time delay can be remarkably improved, the memory simulation requirement is reduced, and the method has wide application prospects.

Description

Technical Field

[0001] The present invention relates to the technical field of groundwater system simulation, and in particular to a method for improving the simulation efficiency of a sandwich body with compression time lag under the principle of ground stress variation. Background Art

[0002] Land subsidence, a downward movement of the Earth's crust caused by the consolidation and compression of loose underground strata, threatens the sustainable development of the natural environment and the economy and society. Quantitative research on land subsidence is crucial for developing conservation policies and restoring the environment. The continuous decline in groundwater levels is often the primary cause of land subsidence. Currently, many scholars are using monitoring technologies, model tests, artificial intelligence, and other methods and techniques to study the mechanisms and trends of land subsidence. Statistical methods based on monitoring data are widely used to identify and predict land subsidence patterns. However, these methods ignore the deformation mechanisms of land subsidence. Numerical simulation methods based on physical mechanisms can account for the interaction between groundwater and land subsidence, making them an important tool for studying and predicting the evolution of land subsidence.

[0003] Currently, the most widely used MODFLOW software development modules include SUB and CSUB. Both modules have the function of simulating intercalated bodies with compressive time lags by solving the one-dimensional vertical head diffusion equation of the intercalated body. The SUB software package is based on the principle of head change, that is, it is assumed that the change in effective stress at any position in the aquifer is equal to the negative value of the change in the aquifer head. Therefore, when simulating intercalated bodies with compressive time lags, the effective stress acting on the top and bottom of the intercalated body is the same. When solving the diffusion equation, the principle of symmetry can be used to discretize half the thickness of the intercalated body, thereby improving simulation efficiency and reducing memory requirements. This method is called the half-thickness discretization format. The CSUB software package integrated in the latest version of MODFLOW6 integrates all the functions of simulating coarse-grained media in aquifers, intercalated bodies without compressive time lags, intercalated bodies with compressive time lags, and the compaction release and expansion storage of pore water in aquifers. Compared to the SUB module, CSUB's settlement simulations are more realistic when involving phreatic aquifers because the impact of phreatic level fluctuations on geostress can be reflected in CSUB, whereas SUB, based on the principle of varying hydraulic head, cannot, and thus typically simulates significantly larger settlements. However, based on the principle of varying geostress, the effective stress at any location in the aquifer is equal to the difference between the total geostress at that location and the pore water pressure. Since total geostress increases with depth, the effective stress acting on the top and bottom of a compressive time-lag interlayer differs, generally requiring a full-thickness discretization scheme. Consequently, for the same vertical discretization spacing of the interlayer, the computational workload and memory requirements of CSUB when simulating a compressive time-lag interlayer are approximately twice those of SUB. This significantly impacts model runtime and, when the simulation scale is large, places higher demands on the computing device's memory capacity. Summary of the Invention

[0004] In order to solve the problems existing in the above-mentioned prior art, the purpose of the present invention is to provide a method for improving the simulation efficiency of interlayer bodies with compressive time lag under the principle of ground stress change, which can significantly improve the simulation efficiency of interlayer bodies with compressive time lag and reduce the simulation memory requirements, and has broad application prospects.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A method for improving the simulation efficiency of a sandwich body with compression time lag based on the principle of geostress variation, comprising:

[0007] The half-thickness grid space discretization of the compressed time-delay sandwich body is performed using rectangular grids to obtain the grid elements for simulating the compressed time-delay sandwich body.

[0008] Constructing a water balance model based on the compressed time-lag sandwich body simulation grid unit;

[0009] The water balance models of the node units are combined to obtain a differential equation group, and the differential equation group is used to simulate the compaction or expansion of each node unit in the compression time-lag sandwich body.

[0010] Optionally, constructing the water balance model includes:

[0011] A one-dimensional hydraulic head diffusion model is established to describe the compaction water release and expansion water storage of the compressed time-delay interlayer:

[0012]

[0013] The one-dimensional head diffusion model is numerically discretized to obtain the water balance model at the node element i=1:

[0014]

[0015] Among them, K′ v is the vertical permeability coefficient of the interlayer with compression lag, γ w is the specific gravity of water, is the hydraulic head of groundwater grid cell j at the end of the mth calculation period, is the hydraulic head on the first node element at the end of the mth calculation period, Δz is the vertical discrete distance between the centers of two node elements on the equivalent sandwich body, is the water storage rate to be selected at the first node element in the mth calculation period, is the effective stress acting on the first node element at the end of the mth calculation period, is the pre-consolidation stress on the first node element at the end of the m-1th calculation period, It is the elastic water storage rate on the first node unit within the m time period.

[0016] Optionally, the effective stress calculation model acting on the first node unit at the end of the m-th calculation period includes:

[0017]

[0018] Among them, is the total in-situ stress acting on the first node unit at the end of the m-th calculation period, z1 is the central elevation of the first node unit, is the water head on the first node unit at the end of the m-th calculation period.

[0019] Optionally, the system of difference equations includes:

[0020] The first difference model: Substitute the effective stress calculation model acting on the first node unit at the end of the m-th calculation period into the water balance model at node unit i = 1 to obtain the first difference model:

[0021]

[0022] Among them, K′ v is the vertical permeability coefficient of the compressible lag interlayer body, γ w is the unit weight of water, is the water head of the groundwater grid unit j at the end of the m-th calculation period, is the water head on the first node unit at the end of the m-th calculation period, is the water head on the second node unit at the end of the m-th calculation period, Δz is the vertical discrete distance between the centers of two node units on the equivalent interlayer body, is the optional water storage rate on the first node unit in the m-th calculation period, is the effective stress acting on the first node unit at the end of the m-th calculation period, is the pre-consolidation stress on the first node unit at the end of the (m - 1)-th calculation period, is the elastic water storage rate on the first node unit within the m time period, is the total in-situ stress on the first node unit at the end of the m-th calculation period;

[0023] The second difference model: Substitute the effective stress calculation model acting on the first <i<NN node units at the end of the m-th calculation period into the water balance model at node units 1 <i<NN to obtain the second difference model:

[0024]

[0025] Among them, is the water head on the i-1th node element at the end of the mth calculation period, is the water head on the i+1th node element at the end of the mth calculation period, z i is the center position of the i-th node unit, S skei is the elastic water storage rate of the i-th node element in the time period.

[0026] Optionally, the group of difference equations further includes:

[0027] Third differential model: Substitute the effective stress calculation model acting on the i=NNth node element at the end of the mth calculation period into the water balance model at the node element i=NN to obtain the third differential model:

[0028]

[0029] in, is the water head on the NN-1th node element at the end of the mth calculation period, is the water storage rate to be selected on the NNth node unit in the mth calculation period, z NN is the center position of the NNth node unit, S skeNN is the elastic water storage rate of the NNth node element in the time period.

[0030] Optionally, simulating the compaction or expansion of each node element in the compressed time-delay sandwich body by using the differential equation group includes:

[0031] All discrete grid units are combined according to the differential equation group to obtain a matrix model;

[0032] Calculate the mean of the effective stress at the i-th and 2NN-1-i-th node units, and if the distances between the i-th and 2NN-1-i-th node units and the center of the sandwich body are equal, then determine the center of the compression lag sandwich body;

[0033] When in the half-thickness discrete format, according to the mean and the center position, if the head at the i-th node unit is equal to the mean of the heads of the i-th and 2NN-1-i-th node units in the full-thickness discrete format, it meets the equalization expectation, and the matrix model is further used to simulate the compaction or expansion of each node unit in the compressed time-lag sandwich body.

[0034] Optionally, the matrix model includes:

[0035] [A] m [h] m =[r] m

[0036] Among them, [A] m for, [h]m for, [r] m is a one-dimensional known vector of NN elements.

[0037] Optionally, calculating the mean of the effective stresses at the i-th and 2NN-1-i-th node elements includes:

[0038]

[0039] in, is the mean of the effective stresses at the i-th and 2NN-1-i-th node elements, is the effective stress at the 2NN-1-i node element, is the total ground stress acting on the first node element at the end of the mth calculation period, is the water head on the first node element at the end of the mth calculation period, is the water head on the 2NN-1-ith node element at the end of the mth calculation period, z 2NN-1-i is the center position of the 2NN-1-ith node unit.

[0040] Optionally, determining the center position of the compressed time-delay sandwich body includes:

[0041]

[0042] in, is the effective stress at the 2NN-1-i node element, z B It is the vertical position at the center of the sandwich body.

[0043] The beneficial effects of the present invention are:

[0044] In the present invention, the number of discrete elements in the full-thickness discretization scheme for a compressed time-delay sandwich is 2NN-1 node elements. However, using the half-thickness discretization scheme, this number is reduced to NN node elements, directly halving the number of nodes required for calculation. For example, the full-thickness discretization scheme has 1001 discrete elements, while the half-thickness discretization scheme has only 501. This significantly reduces the amount of data required for calculation, thereby significantly improving computational efficiency.

[0045] Due to the reduction in the number of discrete nodes, the size of the differential equation system formed by the simultaneous node units is also reduced accordingly. When solving the differential equation system, the amount of computation and solution time are greatly reduced as the size of the equation system is reduced. For example, the number of variables in the matrix equation solved by the half-thickness discretization method is only half that of the full-thickness discretization method. This significantly reduces the number of iterations and computing resources required during the solution process, thereby improving the overall efficiency of the simulation calculation.

[0046] During computer simulations, each discrete node element in the present invention needs to store relevant data, such as head value, effective stress, and water storage rate. The half-thickness discretization method reduces the number of discrete nodes by half, and accordingly, the memory space required to store this data is also halved. Reducing memory requirements is particularly important for large-scale simulations, as it allows simulations to proceed smoothly with limited computing resources, avoiding computational failures or frequent memory swapping due to insufficient memory, thereby further improving the feasibility and efficiency of simulations.

[0047] In summary, the present invention proposes a half-thickness discretization method for delayed interlayer bodies based on the characteristic that the ground stress is linearly distributed in the vertical direction of the interlayer body. Compared with the existing technology, it has higher computational efficiency and can also achieve good accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0049] Figure 1 This is a flow chart of a method for improving the simulation efficiency of a sandwich body with compression time lag based on the principle of ground stress variation according to an embodiment of the present invention;

[0050] Figure 2 A one-dimensional spatial discrete graph of a sandwich body with a compressed time delay according to an embodiment of the present invention;

[0051] Figure 3 A diagram showing a model configuration for an embodiment of the present invention;

[0052] Figure 4 A comparison chart of the simulated settlement amounts of MODFLOW-CSUB in an embodiment of the present invention and the present invention;

[0053] Figure 5 This is a comparison diagram of the water heads at the nodes calculated by MODFLOW-CSUB and the present invention at the 99.15th day of the embodiment of the present invention;

[0054] Figure 6 This is a comparison diagram of the water heads at the nodes calculated by MODFLOW-CSUB and the present invention at the 501.28th day of the embodiment of the present invention;

[0055] Figure 7 This is a comparison diagram of the water heads at the nodes calculated by MODFLOW-CSUB and the present invention at the 1000th day in an embodiment of the present invention. DETAILED DESCRIPTION

[0056] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. 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.

[0057] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0058] like Figure 1 As shown, this embodiment discloses a method for improving the simulation efficiency of a compressive time-lag interlayer body under the principle of ground stress change, comprising: using a rectangular grid to perform half-thickness gridding spatial discretization on the compressive time-lag interlayer body to obtain grid units for simulating the compressive time-lag interlayer body; constructing a water balance model based on the grid units for simulating the compressive time-lag interlayer body; combining the water balance models of each node unit to obtain a group of differential equations, and using the group of differential equations to simulate the compaction or expansion of each node unit in the compressive time-lag interlayer body.

[0059] Specifically, this embodiment discloses a method for improving the simulation efficiency of a sandwich body with compression time lag based on the principle of geostress variation, including:

[0060] S1. Based on the linear distribution of ground stress in the vertical direction of the interlayer, the compressive time-lag interlayer is spatially discretized using half-thickness grids to obtain the grid units for simulating the compressive time-lag interlayer.

[0061] S2. Construct the water balance equation for discrete grid cells at different locations and modify the effective stress calculation formula.

[0062] S3. The equations of each node unit are combined to form a differential equation group, and the equations are solved to simulate the compaction / expansion process of each node unit of the sandwich body.

[0063] Specifically, step S1 is to discretize the half-thickness of the sandwich body with compression time lag into one-dimensional NN grid units using a rectangular grid, such as Figure 2 .

[0064] Furthermore, constructing the water balance model includes: establishing a one-dimensional head diffusion model to describe the compaction water release and expansion water storage in the compression time-lag interlayer simulation grid unit; numerically discretizing the one-dimensional head diffusion model to obtain the water balance model.

[0065] Specifically, when simulating land subsidence, the water flow in the interlayer with compression lag can be considered vertical. Based on the principle of ground stress change, the compaction water release and expansion water storage of the interlayer can be expressed by the following one-dimensional head diffusion equation:

[0066]

[0067] Where: K′ v is the vertical permeability coefficient of the interlayer with compression lag (L / T); h is the hydraulic head in the interlayer (L); z is the vertical spatial coordinate (L); S s is the water storage rate (-); t is the time (T); γ w is the specific gravity of water (M / L2 / T2); σ' is the effective stress (M / L / T2).

[0068] According to formula (1), the numerical discretization of the partial differential equation is performed. For any m-th calculation period, the water balance equation at node element i=1 is:

[0069]

[0070] Where: Δz is the vertical discrete spacing between the centers of two node elements on the equivalent sandwich body (L); Δt is the time step (T); is the water storage rate to be selected at the first node element in the mth calculation period (1 / L); is the hydraulic head (L) of groundwater grid cell j (assuming the equivalent interlayer body is located on the groundwater grid cell) at the end of the mth calculation period; is the water head on the first node element at the end of the mth calculation period; is the water head on the second node element at the end of the mth calculation period; is the water storage rate to be selected on the first node unit in the m period (1 / L); is the effective stress (in water column height) acting on the first node element at the end of the mth calculation period (L); is the pre-consolidation stress (in water column height) on the first node element at the end of the m-1th calculation period (L); is the elastic water storage rate (1 / L) of the first node element in the m-time period; It is the effective stress (in water column height) acting on the first node element at the end of the m-1th calculation period (L).

[0071] Furthermore, the system of difference equations includes: The first difference model: Substitute the effective stress calculation model acting on the first node unit at the end of the m-th calculation period into the water balance model at node unit i = 1 to obtain the first difference model. The second difference model: Substitute the effective stress calculation model acting on the node units where 1 < i < NN at the end of the m-th calculation period into the water balance model at the node units where 1 < i < NN to obtain the second difference model. The third difference model: Substitute the effective stress calculation model acting on the node unit i = NN at the end of the m-th calculation period into the water balance model at node unit i = NN to obtain the third difference model.

[0072] Specifically, the effective stress acting on the first node unit at the end of the m-th calculation period in formula (2) is calculated as:

[0073]

[0074] Where: is the total in-situ stress acting on the first node unit at the end of the m-th calculation period (in terms of water column height) (L); z1 is the central elevation of the first node unit (L).

[0075] Combining formulas (2) and (3) and arranging gives:

[0076]

[0077] Similarly, establish the difference equation for the node units where 1 < i < NN, and we have:

[0078]

[0079] Arranging gives:

[0080]

[0081] Finally, establish the difference equation for node unit NN. Note that the thickness of node unit NN is only Δz / 2, and we have:

[0082]

[0083] Arranging gives:

[0084]

[0085] S3. Combine the equations of each node unit to form a system of difference equations, and solve to realize the simulation of the compaction / expansion process of each node unit of the interlayer body.

[0086] Furthermore, the use of a group of differential equations to simulate the compaction or expansion of each node unit in the compression time-lag sandwich body includes: combining all discrete grid units according to the group of differential equations to obtain a matrix model; calculating the mean of the effective stresses at the i-th and 2NN-1-i-th node units, and at the same time, the distances between the i-th and 2NN-1-i-th node unit nodes and the center position of the sandwich body are equal, and then determining the center position of the compression time-lag sandwich body; when in a half-thickness discrete format, according to the mean and the center position, determining that the water head at the i-th node unit is approximately equal to the mean of the water heads of the i-th and 2NN-1-i-th node units in the full-thickness discrete format, which meets the equalization expectation, and further using the matrix model to simulate the compaction or expansion of each node unit in the compression time-lag sandwich body.

[0087] By referring to formulas (4), (6) and (8) and combining all discrete grid elements, we can obtain the matrix equation:

[0088] [A] m [h] m =[r] m (9)

[0089] Where: [A] m is a tridiagonal symmetric matrix of NN×NN; [h] m is a one-dimensional head vector with NN elements; [r] m is a one-dimensional known vector of NN elements.

[0090] The elements in the matrix equation are:

[0091]

[0092] in, is the element at the i-th row and j-th column of the tridiagonal symmetric matrix, is the element at the 1st row and 1st column of the tridiagonal symmetric matrix, is the element at the i-th row and i-th column of the tridiagonal symmetric matrix, is the element at the NNth row and NNth column of the tridiagonal symmetric matrix, is the one-dimensional known vector of the first element, r i m is the one-dimensional known vector of the i-th element, is the one-dimensional known vector of the NNth element.

[0093] The calculation formula of water storage rate in the present invention is:

[0094]

[0095] in, The modified effective stress at node element i is calculated as:

[0096]

[0097] Where: z B is the vertical position of the center of the sandwich body (L), is the total ground stress at the center of the interlayer (measured in water column height) (M / L2 / T2).

[0098] For the same sandwich with compression lag, assuming that the unit discretization spacing in the half-thickness discretization format is equal to that in the full-thickness discretization format, the half-thickness discretization format has NN node elements and the full-thickness discretization format has 2NN-1 node elements. For the i-th node in the half-thickness discretization format, we expect its water storage rate to be the average of the water storage rates of the i-th and 2NN-1-i node elements in the full-thickness discretization format. To this end, we calculate the average of the effective stresses at the i-th and 2NN-1-i node elements. have:

[0099]

[0100] Since the total ground stress is linearly distributed in the vertical direction of the sandwich body, and the distances between the i-th and 2NN-1-i-th node unit nodes and the center of the sandwich body are equal, the reference center position of the sandwich body is:

[0101]

[0102] Substituting formula (20) into formula (19), we can see that:

[0103]

[0104] In the half-thickness discretization format, the water head at the i-th node element is approximately equal to the mean of the water heads at the i-th and 2NN-1-i-th node elements in the full-thickness discretization format. Therefore, comparing formula (19) and formula (21), we have In line with our average expectations.

[0105] The compaction / water storage process of each node element of the sandwich body can be simulated by solving the matrix equation (Formula (9)). Since the number of variables in the equation group is only half of the full-thickness discretization method, the simulation efficiency is higher and the memory requirement is less.

[0106] This example simulates the drainage process of a thick interlayer due to the gradual decrease of the aquifer head. It is divided into 1 stress period and the total simulation time is 1000 days. It is divided into 100 calculation periods of different lengths using a time step multiplier of 1.05. The model settings are as follows Figure 3As shown, the model grid consists of one layer, one row, and three columns. The top elevation of the model is 0 meters, the bottom elevation is -1000 meters, and the row and column widths are both 1 meter. A delayed interlayer is located in the middle grid, and the hydraulic head values of the fixed-head grid cells on both sides are both 0 meters. The hydraulic conductivity of the middle grid is set to 1.0E+6 to ensure that the hydraulic head value is also 0 meters. Detailed model parameters are shown in Table 1. For comparison, the full-thickness discretization method of the MODFLOW6-CSUB package was also used to simulate the aquifer in this example.

[0107] Table 1 Model parameters

[0108]

[0109]

[0110] Application effect achieved:

[0111] This embodiment generally does not meet such harsh conditions under actual working conditions, and can be used as an extreme scenario to verify the reliability of the half-thickness discretization method. Figure 4 This is a comparison chart of the settlement calculation results of the present invention and MODFLOW-CSUB. It can be seen that the simulated settlement of the present invention is almost the same as that of MODFLOW-CSUB, and both have high accuracy. Figure 5 、 6 Figures 7 and 8 show a comparison of the full-thickness discretization method, where the hydraulic head values at all symmetric nodes are summed and averaged, with the half-thickness discretization method. The three figures correspond to time points at 99.15 days, 501.28 days, and 1000 days, respectively. It can be seen that the hydraulic head distribution is essentially consistent. Table 2 shows the groundwater system equilibrium table. In this case, the interlayer undergoes both elastic and inelastic compaction. Influenced by the water storage rate, the inelastic compaction is significantly greater than the elastic compaction. The final groundwater equilibrium results are also similar. The relative deviation in interlayer water release between the two methods is only 0.1874%. These results demonstrate that the half-thickness discretization method achieves similar computational accuracy to the full-thickness discretization method in scenarios with extreme thicknesses and a large number of discrete elements, improving computational efficiency while also being mechanistically sound.

[0112] Table 2 Groundwater system balance table

[0113]

[0114] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. A method for improving the simulation efficiency of sandwich structures with compression time lag based on the principle of geostress variation, characterized in that: Comprising: Performing semi - thickness grid - based spatial discretization on the compressed time - lag interlayer body using a rectangular grid to obtain simulation grid units of the compressed time - lag interlayer body; Constructing a water balance model based on the simulation grid units of the compressed time - lag interlayer body; Combining the water balance models of each node unit to obtain a system of difference equations, and using the system of difference equations to simulate the compaction or expansion of each node unit in the compressed time - lag interlayer body.

2. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 1, characterized in that: Constructing the water balance model includes: Establishing a one - dimensional head diffusion model for describing the water release during compaction and water storage during expansion of the compressed time - lag interlayer body: Performing numerical discretization on the one - dimensional head diffusion model to obtain the water balance model at node unit \(i = 1\): Among them, K′ v is the vertical permeability coefficient of the interlayer with compression lag, γ w is the specific gravity of water, is the hydraulic head of groundwater grid cell j at the end of the mth calculation period, is the hydraulic head on the first node element at the end of the mth calculation period, Δz is the vertical discrete distance between the centers of two node elements on the equivalent sandwich body, is the water storage rate to be selected at the first node element in the mth calculation period, is the effective stress acting on the first node element at the end of the mth calculation period, is the pre-consolidation stress on the first node element at the end of the m-1th calculation period, is the elastic water storage rate of the first node element in the m-time period.

3. The method for improving the simulation efficiency of sandwich structures with compression time lag based on the principle of geostress variation according to claim 2, characterized in that: The effective stress calculation model acting on the first node unit at the end of the \(m\) - th calculation period includes: in, is the total ground stress acting on the first node element at the end of the mth calculation period, z1 is the center elevation of the first node element, is the water head on the first node element at the end of the mth calculation period.

4. The method for improving the simulation efficiency of sandwich structures with compression time lag based on the principle of geostress variation according to claim 1, characterized in that: The system of difference equations includes: The first difference model: Substituting the effective stress calculation model acting on the first node unit at the end of the \(m\) - th calculation period into the water balance model at node unit \(i = 1\) to obtain the first difference model: Among them, K′ v is the vertical permeability coefficient of the interlayer with compression lag, γ w is the specific gravity of water, is the hydraulic head of groundwater grid cell j at the end of the mth calculation period, is the water head on the first node element at the end of the mth calculation period, is the hydraulic head on the second node element at the end of the mth calculation period, Δz is the vertical discrete distance between the centers of the two node elements on the equivalent sandwich body, is the water storage rate to be selected at the first node element in the mth calculation period, is the effective stress acting on the first node element at the end of the mth calculation period, is the pre-consolidation stress on the first node element at the end of the m-1th calculation period, is the elastic water storage rate of the first node element in the m-time period, is the total ground stress on the first node element at the end of the mth calculation period; The second difference model: Substituting the effective stress calculation model acting on the node units where \(1\lt i\lt NN\) at the end of the \(m\) - th calculation period into the water balance model at node units where \(1\lt i\lt NN\) to obtain the second difference model: in, is the water head on the i-1th node element at the end of the mth calculation period, is the water head on the i+1th node element at the end of the mth calculation period, z i is the center position of the i-th node unit, S skei is the elastic water storage rate of the i-th node element in the time period.

5. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 4, characterized in that: The system of difference equations further includes: The third difference model: Substituting the effective stress calculation model acting on the node unit \(i = NN\) at the end of the \(m\) - th calculation period into the water balance model at node unit \(i = NN\) to obtain the third difference model: in, is the water head on the NN-1th node element at the end of the mth calculation period, is the water storage rate to be selected on the NNth node unit in the mth calculation period, z NN is the center position of the NNth node unit, S skeNN is the elastic water storage rate of the NNth node element in the time period.

6. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 1, characterized in that: Using the system of difference equations to simulate the compaction or expansion of each node unit in the compressed time - lag interlayer body includes: 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 7. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 6, characterized in that: 联立所有离散网格单元,得到矩阵模型; [A] m [h] m =[r] m Among them, [A] m for, [h] m for, [r] m is a one-dimensional known vector of NN elements.

8. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 6, characterized in that: 联立所有离散网格单元,得到矩阵模型; in, is the mean of the effective stresses at the i-th and 2NN-1-i-th node elements, is the effective stress at the 2NN-1-i node element, is the total ground stress acting on the first node element at the end of the mth calculation period, is the water head on the first node element at the end of the mth calculation period, is the water head on the 2NN-1-ith node element at the end of the mth calculation period, z 2NN-1-i is the center position of the 2NN-1-ith node unit.

9. The method for improving the simulation efficiency of sandwich bodies with compression time lag based on the principle of geostress variation according to claim 6, characterized in that: 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; 联立所有离散网格单元,得到矩阵模型; in, is the effective stress at the 2NN-1-i node element, z B It is the vertical position at the center of the sandwich body.

Images