Method for improving simulation efficiency of body with compression time delay sandwich under ground stress change principle

Abstract

The present application relates to a kind of methods for improving the simulation efficiency of compression time lag sandwich body under the principle of ground stress variation, comprising: using rectangular grid to carry out half-thickness grid space dispersion to compression time lag sandwich body, obtain compression time lag sandwich body simulation grid unit;According to the water balance model of compression time lag sandwich body simulation grid unit, obtain difference equation set by combining the water balance model of each node unit, and the compaction or expansion of each node unit in compression time lag sandwich body is simulated using difference equation set.The present application can significantly improve the simulation efficiency of compression time lag sandwich body, and reduce the simulation memory requirement, with wide application prospect.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of groundwater system simulation, in particular to a method for improving simulation efficiency of a compression time lag interlayer under the principle of ground stress change. BACKGROUND

[0002] Land subsidence is a downward movement caused by the consolidation compression of underground loose strata, which threatens the sustainable development of the natural environment and the economy and society. The quantitative study of land subsidence is a key issue for formulating protection policies and restoring the environment. The continuous decline of the groundwater level is usually the main cause of land subsidence. At present, many scholars use monitoring technology, model test, artificial intelligence and other methods and technologies to study the subsidence mechanism and trend prediction. The method of using statistical methods based on monitoring data to identify the development law of land subsidence and to make prediction analysis is widely used, but this method ignores the deformation mechanism of land subsidence. The numerical simulation method based on the physical mechanism can consider the interaction process of groundwater and land subsidence, and has become an important means to study and predict the evolution law of land subsidence.

[0003] At present, the most widely used MODFLOW software has two modules, SUB and CSUB, which can simulate compression time lag interlayers by solving one-dimensional vertical water head diffusion equations. The SUB software package is based on the principle of water head change, that is, it is assumed that the change of effective stress at any position in the aquifer is equal to the negative value of the change of water head in the aquifer. Therefore, when simulating compression time lag interlayers, the effective stresses acting on the top and bottom of the interlayer are the same. When solving the diffusion equation, the symmetry principle can be used to discretize half the thickness of the interlayer, thereby improving the simulation efficiency and reducing the memory requirement. This method is called half-thickness discretization format. The latest version of MODFLOW6 integrates the CSUB software package, which integrates all functions of simulating coarse-grained media in aquifers, non-compression time lag interlayers, compression time lag interlayers and the compaction release and swelling storage of pore water in aquifers. Compared with the SUB module, the subsidence simulation of CSUB involving phreatic aquifers is more realistic because the influence of phreatic water level fluctuations on ground stress can be reflected in CSUB, while the SUB based on the water head change principle cannot. Generally, the simulated subsidence will be significantly larger. However, based on the principle of ground stress change, the effective stress at any position in the aquifer is equal to the difference between the total ground stress at that position and the pore water pressure, and the total ground stress increases with the increase of the burial depth. Therefore, for compression time lag interlayers, the effective stresses acting on the top and bottom of the interlayer are different. Generally, the full-thickness discretization format is used, which means that the calculation workload and memory requirement of CSUB in simulating compression time lag interlayers are about twice that of SUB under the same vertical discretization interval of the interlayer. This significantly affects the model running time and requires higher memory capacity of the computing device when simulating large-scale models. SUMMARY

[0004] In order to solve the problems existing in the prior art, the purpose of the present application is to provide a method for improving the simulation efficiency of a compression time lag sandwich body under the principle of ground stress change, which can significantly improve the simulation efficiency of the compression time lag sandwich body and reduce the simulation memory requirement, and has a wide application prospect.

[0005] To achieve the above purpose, the present application provides the following scheme:

[0006] The method for improving the simulation efficiency of a compression time lag sandwich body under the principle of ground stress change comprises the following steps:

[0007] The compression time lag sandwich body is semi-thickness meshed and spatially dispersed by using a rectangular grid to obtain a compression time lag sandwich body simulation grid unit;

[0008] According to the compression time lag sandwich body simulation grid unit, a water balance model is constructed;

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

[0010] Optionally, the construction of the water balance model comprises:

[0011] A one-dimensional water head diffusion model is established to describe the compaction water release and swelling water storage of the compression time lag sandwich body:

[0012]

[0013] The one-dimensional water head diffusion model is numerically dispersed to obtain a water balance model at node unit i=1:

[0014]

[0015] wherein, K' v is the vertical permeability coefficient of the compression time lag sandwich body, γ w is the unit weight of water, is the water head of the groundwater grid unit j at the end of the mth calculation period, is the water head on the first node unit at the end of the mth calculation period, Δz is the vertical discrete interval between the centers of two node units on the equivalent sandwich body, is the selected water storage rate on the first node unit at the mth calculation period, is the effective stress acting on the first node unit at the end of the mth calculation period, is the pre-consolidation stress on the first node unit at the end of the m-1th calculation period, is the elastic water storage rate of the first node unit in the mth period.

[0016] Optionally, the effective stress calculation model acting on the first node unit at the end of the mth calculation period comprises:

[0017]

[0018] wherein, is the total ground stress acting on the first node unit at the end of the mth calculation period, and z1 is the central elevation of the first node unit, is the water head of the first node unit at the end of the mth calculation period.

[0019] Optionally, the difference equation set comprises:

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

[0021]

[0022] wherein, K' v is the vertical permeability coefficient of the compression time 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 mth calculation period, is the water head of the first node unit at the end of the mth calculation period, is the water head of the second node unit at the end of the mth calculation period, and Δz is the vertical discrete spacing between the centers of the two node units on the equivalent interlayer body, is the selected water storage rate of the first node unit in the mth period, is the effective stress acting on the first node unit at the end of the mth calculation period, is the pre-consolidation stress of the first node unit at the end of the m-1th calculation period, is the elastic water storage rate of the first node unit in the mth period, is the total ground stress of the first node unit at the end of the mth calculation period;

[0023] The second difference model: the effective stress calculation model acting on the first <i<NN node unit at the end of the mth calculation period is substituted into the water balance model at the node unit 1<i<NN to obtain the second difference model:

[0024]

[0025] wherein, a water head on the i-1th node unit at the end of the mth calculation period, a water head on the i+1th node unit at the end of the mth calculation period, i a center position of the i th node unit, skei a elastic water storage rate on the i th node unit in the period.

[0026] Optionally, the difference equation set further comprises:

[0027] a third difference model: substituting a effective stress calculation model acting on the i=NNth node unit at the end of the mth calculation period into a water quantity balance model at the i=NNth node unit, to obtain the third difference model:

[0028]

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

[0030] Optionally, simulating compaction or expansion of each node unit in the compression time lag interlayer body by using the difference equation set comprises:

[0031] According to the difference equation set, all discrete grid units are solved to obtain a matrix model;

[0032] When the mean value of effective stresses at the i th and the 2NN-1-i th node units is calculated, and the distance between the i th and the 2NN-1-i th node units and the center position of the interlayer body is equal, the center position of the compression time lag interlayer body is determined;

[0033] When in a half-thickness discrete format, according to the mean value and the center position, if the water head at the i th node unit is equal to the mean value of water heads at the i th and the 2NN-1-i th node units in a full-thickness discrete format, it meets the homogenization expectation, and further simulating compaction or expansion of each node unit in the compression time lag interlayer body by using the matrix model.

[0034] Optionally, the matrix model comprises:

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

[0036] wherein, [A] m is, [h]m is the i-th node element, m is a one-dimensional known vector with NN elements.

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

[0038]

[0039] wherein, is the mean of the effective stress at the i-th and the 2NN-1-i-th node elements, is the effective stress at the 2NN-1-i-th node element, is the total ground stress acting on the 1st node element at the end of the m-th calculation period, is the water head on the 1st node element at the end of the m-th calculation period, is the water head on the 2NN-1-i-th node element at the end of the m-th calculation period, 2NN-1-i is the central position of the 2NN-1-i-th node element.

[0040] Optionally, determining the central position of the compression time lag sandwich body comprises:

[0041]

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

[0043] The present application has the following beneficial effects:

[0044] Under the full-thickness discrete format, the number of discrete elements of the compression time lag sandwich body is 2NN-1 node elements, while under the half-thickness discrete format, the number of discrete elements is reduced to NN node elements, directly reducing the number of nodes that need to be calculated by half. For example, the full-thickness discrete format has 1001 discrete elements, while the half-thickness discrete format has only 501 discrete elements, greatly reducing the amount of data that needs to be processed in the calculation process, thereby significantly improving the calculation efficiency.

[0045] Due to the reduction in the number of discrete nodes, the size of the difference equation system formed by each node element is also correspondingly reduced. When solving the difference equation system, the amount of calculation and the solving time will be greatly reduced with the reduction of the size of the equation system. For example, the number of variables of the matrix equation solved by the half-thickness discrete method is only half of that of the full-thickness discrete method, which significantly reduces the number of iterations and computing resources required in the solving process, thereby improving the overall efficiency of the simulation calculation.

[0046] In the computer simulation process, each discrete node unit needs to store relevant data such as water head value, effective stress, water storage rate and the like. The half-thickness discrete method reduces the number of discrete nodes by half, and accordingly the memory space required for storing these data is also reduced by half. For large-scale simulation problems, the reduction of memory requirement is particularly important, which can make the simulation smoothly proceed under limited computing resources, avoid the calculation failure or frequent memory exchange operation caused by insufficient memory, and further improve the feasibility and efficiency of the simulation.

[0047] In summary, the present application proposes a half-thickness discretization method of the delayed interlayer based on the linear distribution of the ground stress in the vertical direction of the interlayer, which has higher calculation efficiency compared with the prior art, and can also have good accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0048] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments will be briefly introduced as follows. Obviously, the drawings in the following description only constitute some embodiments of the present application, and other drawings can also be obtained by those skilled in the art without creative labor under the premise of the drawings.

[0049] Figure 1 A flow chart of a method for improving simulation efficiency of a compressed time-lag interlayer under a ground stress change principle according to an embodiment of the present application;

[0050] Figure 2 A one-dimensional space discretization diagram of a compressed time-lag interlayer according to an embodiment of the present application;

[0051] Figure 3 A model setting diagram according to an embodiment of the present application;

[0052] Figure 4 A comparison diagram of MODFLOW-CSUB and the simulated settlement amount according to an embodiment of the present application;

[0053] Figure 5 A comparison diagram of MODFLOW-CSUB and the calculated node water head at 99.15d according to an embodiment of the present application;

[0054] Figure 6 A comparison diagram of MODFLOW-CSUB and the calculated node water head at 501.28d according to an embodiment of the present application;

[0055] Figure 7 A comparison diagram of MODFLOW-CSUB and the calculated node water head at 1000d according to an embodiment of the present application. DETAILED DESCRIPTION

[0056] The technical solutions in the embodiments of the present application will be apparently and completely described below with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all the other embodiments obtained by those skilled in the art without creative work belong to the protection scope of the present application.

[0057] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with the drawings and specific embodiments.

[0058] As shown in the drawings, the present embodiment discloses a method for improving the simulation efficiency of a compression time-lag interlayer body under the principle of ground stress change, comprising: using a rectangular grid to perform half-thickness grid spatial discretization on the compression time-lag interlayer body to obtain a compression time-lag interlayer body simulation grid unit; constructing a water balance model according to the compression time-lag interlayer body simulation grid unit; combining the water balance models of each node unit to obtain a difference equation set, and using the difference equation set to simulate the compaction or swelling of each node unit in the compression time-lag interlayer body. Figure 1 Specifically, the present embodiment discloses a method for improving the simulation efficiency of a compression time-lag interlayer body under the principle of ground stress change, comprising:

[0059] S1, based on the characteristic that the ground stress is linearly distributed in the vertical direction of the interlayer body, performing half-thickness grid spatial discretization on the compression time-lag interlayer body to obtain a compression time-lag interlayer body simulation grid unit.

[0060] S2, constructing a water balance equation for the discrete grid units at different positions, and modifying the effective stress calculation formula.

[0061] S3, combining the equations of each node unit to form a difference equation set, and solving to realize the compaction / swelling process simulation of each node unit of the interlayer body.

[0062] Specifically, step S1 is specifically that the compression time-lag interlayer body is half-thickness discretized into one-dimensional NN grid units by using a rectangular grid, as shown in the drawings.

[0063] . Figure 2

[0064] Further, constructing the water balance model comprises: establishing a one-dimensional water head diffusion model for describing the compaction water release and swelling water storage in the compression time-lag interlayer body simulation grid unit; and performing numerical discretization on the one-dimensional water head diffusion model to obtain the water balance model.

[0065] ​Specifically, in land subsidence simulation, the water flow within the compression-delayed interlayer can be considered vertical. Based on the principle of geostress variation, the compaction-induced water release and expansion-induced water storage within the interlayer can be expressed by the following one-dimensional hydraulic head diffusion equation:

[0066]

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

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

[0069]

[0070] Where: Δz is the vertical discrete distance (L) between the centers of two node elements on the equivalent sandwich structure; Δt is the time step (T); Let L be the candidate water storage rate (1 / L) at the first node unit in the m-th calculation period; Let L be the water head (L) of groundwater grid cell j at the end of the m-th calculation period (assuming the equivalent interlayer is located on the groundwater grid cell); The head of the first node unit at the end of the m-th calculation period; The water head at the second node unit at the end of the m-th calculation period; The water storage rate (1 / L) to be selected on the first node unit within time period m; The effective stress (in terms of water column height) acting on the first node element at the end of the m-th calculation period is (L). The preconsolidation stress (in terms of water column height) on the first node element at the end of the (m-1)th calculation period is (L). Let be the elastic water storage ratio (1 / L) of the first node unit within time period m; The effective stress (in terms of water column height) acting on the first node element at the end of the (m-1)th calculation period is (L).

[0071] Further, the difference equation set comprises: a first difference model: the effective stress calculation model acting on the first node unit at the end of the mth calculation period is substituted into the water balance model at the node unit i=1 to obtain the first difference model. A second difference model: the effective stress calculation model acting on the first to (i-1)th node units at the end of the mth calculation period is substituted into the water balance model at the node unit 1

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

[0073]

[0074] wherein: is the total ground stress (in terms of water column height) (L) acting on the first node unit at the end of the mth calculation period; z1 is the central elevation (L) of the first node unit.

[0075] By combining formula (2) and (3), we have:

[0076]

[0077] Similarly, the difference equation of the first to (i-1)th node units is established, and we have:

[0078]

[0079] By combining formula (2) and (3), we have:

[0080]

[0081] Finally, the difference equation of the node unit NN is established, and it is noted that the thickness of the node unit NN is only Δz / 2, and we have:

[0082]

[0083] By combining formula (2) and (3), we have:

[0084]

[0085] S3, the equations of each node unit are combined to form a difference equation set, and the compaction / expansion process simulation of each node unit of the interlayer body is realized by solving.

[0086] Further, simulating compaction or expansion of each node unit in the compression time-delay sandwich body by using the difference equation set comprises: obtaining a matrix model according to the difference equation set in combination with all discrete grid units; calculating the average of effective stresses at the i th and the 2 N N-1-i th node units, and simultaneously the distances of the i th and the 2 N N-1-i th node units to the center position of the sandwich body are equal, so that the center position of the compression time-delay sandwich body is determined; when being in a half-thickness discrete format, according to the average and the center position, it is determined that the water head at the i th node unit is equal to the average of the water heads of the i th and the 2 N N-1-i th node units in a full-thickness discrete format, so that the homogenization expectation is met, and further, the matrix model is used to simulate the compaction or expansion of each node unit in the compression time-delay sandwich body.

[0087] According to the formula (4), (6) and (8) in combination with all discrete grid units, the matrix equation can be obtained:

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

[0089] Wherein: [A] m is a three-diagonal symmetric matrix of NN*NN; [h] m is a one-dimensional water head vector to be solved of NN elements; [r] m is a one-dimensional known vector of NN elements.

[0090] Each element in the matrix equation is:

[0091]

[0092] Wherein, is the element at the i th row and the j th column position of the three-diagonal symmetric matrix, is the element at the 1 th row and the 1 th column position of the three-diagonal symmetric matrix, is the element at the i th row and the i th column position of the three-diagonal symmetric matrix, is the element at the N N th row and the N N th column position of the three-diagonal symmetric matrix, is a one-dimensional known vector of the 1 st element, r i m is a one-dimensional known vector of the i th element, is a one-dimensional known vector of the N N th element.

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

[0094]

[0095] Wherein, The calculation formula of the corrected effective stress at the node unit i is:

[0096]

[0097] where: z B is the vertical position (L) of the center of the interlayer, is the total earth stress (in terms of water column height) at the center of the interlayer (M / L2 / T2).

[0098] For the same interlayer with compressive time lag, assuming the half-thickness discretization has the same element discretization interval as the full-thickness discretization, the half-thickness discretization has NN nodes, while the full-thickness discretization has 2NN-1 nodes. For the ith node in the half-thickness discretization, we expect that its water storage rate should be the average of the water storage rates of the ith and the 2NN-1-i th nodes in the full-thickness discretization. To this end, we calculate the average of the effective stresses at the ith and the 2NN-1-i th nodes has:

[0099]

[0100] Since the total earth stress is linearly distributed vertically in the interlayer, and the ith and the 2NN-1-i th nodes are equidistant from the center of the interlayer, we have

[0101]

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

[0103]

[0104] For the half-thickness discretization, the water head at the ith node is approximately equal to the average of the water heads at the ith and the 2NN-1-i th nodes in the full-thickness discretization. Therefore, comparing equation (19) and equation (21), we have which is consistent with our expectation of homogenization.

[0105] The compaction / water storage process of the interlayer nodes can be simulated by solving the matrix equation (equation (9)). Since the number of variables in the equation set is only half of that in the full-thickness discretization method, the simulation is more efficient and requires less memory.

[0106] This example simulates the process of thick interlayer drainage due to the gradual decrease of aquifer water head. It is divided into 1 stress period, with a total simulation time of 1000 days, and divided into 100 calculation periods of different lengths using a time step multiplier equal to 1.05. The model settings are as follows: Figure 3As shown, the model grid is composed of 1 layer, 1 row and 3 columns, the model top elevation is 0 m, the bottom elevation is -1000 m, and the row width and column width are both 1 m. The delayed interlayer body is located in the middle grid, the water head values of the two side constant head grid units are both 0 m, the permeability coefficient of the middle grid is set to 1.0E+6 to ensure that the water head value is also 0 m, and the detailed model parameters are shown in Table 1. As a comparison, the full-thickness discretization method of the MODFLOW6-CSUB package is also used to simulate the aquifer of the embodiment.

[0107] Table 1 Model parameters

[0108]

[0109]

[0110] Achieved application effect:

[0111] The actual working condition of the embodiment generally does not reach such harsh conditions, and can be used as an extreme scenario to verify the reliability of the half-thickness discretization method. Figure 4 For the comparison chart of the simulation settlement amount of the embodiment and the settlement amount calculated by MODFLOW-CSUB, it can be seen that the simulation settlement amount of the embodiment is almost the same as that of MODFLOW-CSUB, both of which have high precision. Figure 5 、 6 、7 is a comparison chart of the half-thickness discretization method and the full-thickness discretization method after adding and averaging the water head values of all symmetric nodes, and the three charts correspond to three time nodes of 99.15d, 501.28d and 1000d respectively. It can be seen that the water head distribution is basically consistent. Table 2 is a groundwater system balance table. At this time, the interlayer body has both elastic compaction amount and non-elastic compaction amount, and the non-elastic compaction amount is obviously larger than the elastic compaction amount due to the size of the water storage rate. The final groundwater balance result also has little difference, and the relative deviation of the interlayer body water release amount under the two methods is only 0.1874%. The above results show that the half-thickness discretization method has similar calculation precision to the full-thickness discretization method in the extreme thickness and large number of discretization unit scenario, and is reasonable in mechanism while improving the calculation efficiency.

[0112] Table 2 Groundwater system balance table

[0113]

[0114] The above-described embodiments are only descriptions of the preferred modes of the present application and do not limit the scope of the present application. Without departing from the design spirit of the present application, various modifications and improvements to the technical solutions of the present application made by those skilled in the art shall fall within the protection scope determined by the claims of the present application.

Claims

1. A method for improving the simulation efficiency of a body with a compressive time lag sandwich under the principle of ground stress change, characterized in that, The method comprises the steps of: Semi-thickness grid space discretization is performed on the compressed time-lag sandwich body by using a rectangular grid to obtain a simulation grid unit of the compressed time-lag sandwich body; A water balance model is constructed according to the simulation grid unit of the compressed time-lag sandwich body; The water balance model is constructed by: A one-dimensional water head diffusion model is established to describe the release of water and the storage of water in the compressed time-lag sandwich body: ； The one-dimensional water head diffusion model is numerically discretized to obtain a water balance model at node unit i=1: ； wherein, Kz is the vertical hydraulic conductivity of the compressed time-lag interlayer, ρw is the water density, hj(m) is the water head in the grid cell j at the end of the mth calculation period, h1(m) is the water head in the first node cell at the end of the mth calculation period, L is the vertical distance between the centers of two node cells in the equivalent interlayer, S1(m) is the selected water storage rate in the first node cell at the mth calculation period, σ1(m) is the effective stress in the first node cell at the end of the mth calculation period, σ1(m-1) is the pre-consolidation stress in the first node cell at the end of the (m-1)th calculation period, S1(m) is the elastic water storage rate in the first node cell during the mth period. The water balance models of all node units are combined to obtain a difference equation set, and the difference equation set is used to simulate the compaction or expansion of each node unit in the compressed time-lag sandwich body; The difference equation set comprises: A first difference model is obtained by substituting an effective stress calculation model acting on the first node unit at the end of the mth calculation period into the water balance model at node unit i=1: ； wherein, Kz is the vertical hydraulic conductivity of the compressed interlayer, ρw is the water density, hj(m) is the head in the jth grid cell of the groundwater at the end of the mth calculation period, h1(m) is the head in the first node cell at the end of the mth calculation period, h2(m) is the head in the second node cell at the end of the mth calculation period, d is the vertical separation between the centers of the two node cells of the equivalent interlayer, S1(m) is the selected water storage rate in the first node cell for the mth calculation period, σ1(m) is the effective stress acting in the first node cell at the end of the mth calculation period, σ1(m-1) is the preconsolidation stress in the first node cell at the end of the (m-1)th calculation period, S1(m) is the elastic water storage rate in the first node cell for the mth period, σ1(m) is the total geostatic stress in the first node cell at the end of the mth calculation period. A second difference model is obtained by substituting effective stress calculation models acting on the first to (NN-1)th node units at the end of the mth calculation period into the water balance models at node units 1 to (NN-1): ； wherein, is the head at the end of the mth computation period at the i th node element, is the head at the end of the mth computation period at the i th node element, is the head at the end of the mth computation period at the i th node element, is the head at the end of the mth computation period at the i th node element, is the center position of the i th node element, is the elastic water storage rate at the i th node element during the period.

2. The method for improving the simulation efficiency of a compressed delaminated body under the principle of ground stress variation according to claim 1, characterized in that, The effective stress calculation model acting on the first node unit at the end of the mth calculation period comprises: ； wherein, is the total ground stress acting on the first node element at the end of the mth calculation period, is the central elevation of the first node element, is the water head on the first node element at the end of the mth calculation period.

3. The method of claim 1, wherein the method is used to improve the simulation efficiency of a compressed delaminated body with time lag under the principle of ground stress variation. The difference equation set further comprises: A third difference model is obtained by substituting an effective stress calculation model acting on the i=NNth node unit at the end of the mth calculation period into the water balance model at node unit i=NN: ； wherein, is the water head on the mth node unit at the end of the mth calculation period, is the water head on the mth node unit at the end of the mth calculation period, is the water head on the mth node unit at the end of the mth calculation period, is the water head on the mth node unit at the end of the mth calculation period, is the center position of the mth node unit, is the elastic water storage rate on the mth node unit within the period.

4. The method for improving the simulation efficiency of a compressed delaminated body under the principle of ground stress variation according to claim 1, characterized in that, The simulation of the compaction or expansion of each node unit in the compressed time-lag sandwich body by using the difference equation set comprises: All discrete grid units are combined according to the difference equation set to obtain a matrix model; The average value of the effective stresses at the ith and (2NN-1-i)th node units is calculated, and the distance between the ith and (2NN-1-i)th node units and the center position of the sandwich body is equal, so that the center position of the compressed time-lag sandwich body is determined; When in the semi-thickness discretization format, according to the average value and the center position, if the water head at the ith node unit is equal to the average value of the water heads at the ith and (2NN-1-i)th node units in the full-thickness discretization format, the homogenization expectation is met, and the matrix model is further used to simulate the compaction or expansion of each node unit in the compressed time-lag sandwich body.

5. The method of claim 4, wherein the method further comprises, The matrix model comprises: ； wherein is, is, is a one-dimensional known vector of NN elements.

6. The method for improving the simulation efficiency of a compressed delaminated body under the principle of ground stress variation according to claim 4, characterized in that, The calculation of the average value of the effective stresses at the ith and (2NN-1-i)th node units comprises: ； wherein, is the average of the effective stress at the 1st and 2nd NN -1 -i node elements, is the effective stress at the 2nd NN -1 -i node element, is the total geostatic stress acting on the 1st node element at the end of the mth calculation period, is the water head on the 1st node element at the end of the mth calculation period, is the water head on the 2nd NN -1 -i node element at the end of the mth calculation period, is the center position of the 2nd NN -1 -i node element.

7. The method of claim 4, wherein the method further comprises, The determination of the center position of the compressed time-lag sandwich body comprises: ； wherein, is the effective stress at the 2NN-1-i node element, is the vertical position at the center of the interlayer.

Images