阅读说明：本技术 一种荷载和地下水开采引发地面沉降三维变参数全耦合模拟计算方法 (Load and ground settlement caused by underground water exploitation three-dimensional variable parameter fully-coupled simulation calculation method ) 是由 骆祖江 付延玲 王鑫 李兆 张也 成磊 李子义 于 2021-05-28 设计创作，主要内容包括：本发明提出一种荷载和地下水开采引发地面沉降和地裂缝三维变参数全耦合模拟计算方法,所解决的是在土体上覆荷载和地下水开采双重作用下的地面沉降量和地裂缝的计算问题。本发明根据地下水开采的实际背景,建立地质概念模型。基于Biot固结理论,实现了渗流场与应力场全耦合,构建了用于预测地面沉降量以及地裂缝发育的三维全耦合数值模型。本计算方法考虑岩土体压缩过程中的机理,并耦合弹性、黏弹性和黏塑性的本构模型,结合岩土体参数的变化,将渗流场与应力场紧密联系,可保证能同时模拟地面荷载和地下水动态下降情况下的地面沉降量和地裂缝延伸量,为研究地下水开采过程中的渗流场演化和地面沉降演化规律以及流固耦合作用机理提供了依据。(The invention provides a method for simulating and calculating ground settlement and ground cracks caused by load and underground water exploitation in a three-dimensional variable parameter fully-coupled mode, and solves the problem of calculating ground settlement and ground cracks under the double effects of earth covering load and underground water exploitation. The invention builds a geological conceptual model according to the actual background of underground water exploitation. Based on a Biot consolidation theory, the full coupling of a seepage field and a stress field is realized, and a three-dimensional full coupling numerical model for predicting the ground settlement and the ground crack development is constructed. The calculation method considers the mechanism in the compression process of the rock-soil body, couples the constitutive models of elasticity, viscoelasticity and viscoplasticity, combines the change of the parameters of the rock-soil body, closely connects the seepage field with the stress field, can ensure that the ground settlement and the ground crack extension under the conditions of ground load and groundwater dynamic decline can be simulated simultaneously, and provides a basis for researching the seepage field evolution and the ground settlement evolution rule and the fluid-solid coupling mechanism in the groundwater exploitation process.)

1. A three-dimensional variable parameter fully-coupled simulation calculation method for ground settlement and ground cracks caused by load and underground water exploitation is characterized by comprising the following steps: the method comprises the following steps:

step 1: determining the range, elevation, vertical stratum thickness and lithology conditions of the research area according to the range, hydrogeological conditions and engineering geological conditions of the research area, determining the layering characteristics of each stratum, and establishing a three-dimensional hydrogeological and engineering geological conceptual model of the research area;

step 2: carrying out finite element mesh subdivision on hydrogeology and engineering geology conceptual models, dividing the conceptual models into structural unit meshes, establishing finite element mesh models, carrying out parameter partition on the finite element mesh models, deriving node and unit information of effective meshes of all partitions, respectively selecting various boundary units and nodes in the finite element mesh models, respectively deriving the boundary units and the node information, and manufacturing different boundary types;

and step 3: assigning different attribute values to the different types of boundaries in the step 2, writing the boundary parameters of the model into the derived unit and node information of each boundary, and generating boundary conditions;

and 4, step 4: writing the space coordinate positions and the initial water levels of all the nodes into the node information by using all the effective node information derived in the step 2 to generate an initial water level condition;

and 5: screening out unit and node information of all pumping wells and observation wells in underground water exploitation by using all effective unit and node information of the model derived in the step 2, generalizing the pumping wells into a drainage substructure, writing the flow of the pumping wells into the node information of the pumping wells, and finally writing the water level value of the observation wells into the node information of the observation wells to generate pumping well conditions and observation well conditions;

step 6: screening out rock-soil body node information for bearing loads by using the effective node information derived in the step 2, writing the additional loads of all nodes into the node information, and generating rock-soil body load conditions, wherein the ground load of each node is determined by the overlying load and the additional stress coefficient;

and 7: inputting permeability coefficient, porosity, water supply degree/water storage rate data with nonlinear change in each parameter partition into the parameter partition by using all effective nodes and parameter partition information derived in the step 2 to generate hydrogeological conditions, and writing the elastic modulus and Poisson ratio with nonlinear change in each parameter partition into the parameter partition to generate engineering geological conditions;

and 8: importing the files in the steps 3 to 7 into a groundwater system (GWS) based on a Biot consolidation theory to complete the establishment of a numerical model, carrying out time dispersion by using a Galerkin weighted allowance method, and carrying out finite element solution on the established numerical model by using a pretreatment conjugate gradient method (PCG);

and step 9: and (3) comparing the underground water level change obtained by calculation in the step (8) with the water level of the actual pumping test and the water level of the hydrological pilot hole: if the difference value between the calculated value and the observed value is more than 10% of the water level change value, returning to the step 8 to adjust the parameters until the result is in accordance, and storing the model; and finally, predicting the water level and the ground settlement of the seepage field by using the adjusted model.

2. The method for simulating and calculating the three-dimensional variable parameters of the ground settlement and the ground cracks caused by the load and the underground water exploitation according to claim 1, which is characterized in that: the well treatment method in the step 4 is a drainage substructure method, nodes of the drainage substructure are divided into 4 layers and are numbered again, the number of nodes of an inner layer is 9-12, the number of nodes of a secondary inner layer is 5-8, the number of nodes of an outer layer is 1-4, the number of nodes of an outermost layer is 13-16, 9-12 are drainage hole nodes of a known water head, 1-8 nodes are unknown water head nodes of the secondary inner layer and the outer layer, and the nodes form an internal conduction matrix based on well flow; the outermost 13-16 nodes are an external conductive matrix for connecting the internal and external structures, and an overall conductive matrix [ K ] composed of the internal conductive matrix and the external conductive matrix^s]And corresponding flow arrays q^sFinally, the total conduction matrix of the drainage substructure and the corresponding flow array are:

in the formula: [ K ]_ii]A conduction matrix among unknown water head nodes in the drainage substructure; [ K ]_ib]、[K_bi]Respectively are conduction matrixes between nodes of an internal unknown water head and outlet nodes; [ K ]_bb]A conduction matrix between unknown water head nodes at the outlet of the drainage substructure; { h_iIs an internal known node of the drainage substructureA waterhead array; { h_bThe drainage substructure is a node water head array at an outlet of the drainage substructure; { q ] q_iThe flow rate of the known node water head to the internal unknown node water head is obtained; { q ] q_bAnd the flow rate of the known node water head to the unknown node water head of the outlet is obtained.

3. The method for simulating and calculating the three-dimensional variable parameters of the ground settlement and the ground cracks caused by the load and the underground water exploitation according to claim 1, which is characterized in that: in the step 6, the overburden load of the rock and soil mass is accurately described, an additional stress coefficient is introduced, and the final overburden load calculation formula of the rock and soil mass is as follows:

P＝P_d·a_sr

in the formula: p_dBuilding load for the ground; a is_srFor additional stress coefficients, the calculation is as follows:

wherein, l is the long side of the rectangular load, b is the short side of the rectangular load, x, y and z are the coordinate values of the calculation point, a_srIs an additional stress factor.

4. The method for calculating the three-dimensional variable parameter fully-coupled simulation of the ground settlement and the ground fracture caused by the load and the underground water exploitation as claimed in claim 3, wherein the parameters of the porosity, the permeability coefficient, the elastic modulus and the Poisson's ratio generated in the step 7 are characterized by nonlinear changes, and the parameter change equation is as follows:

in the formula: n is porosity; n is₀Is the initial porosity; Δ P is the pore water pressure change; k_sThe volume elastic compression modulus of the solid particles of the porous medium skeleton; epsilon_vIs a bulk strain; k is the permeability coefficient; k₀Is the initial permeability coefficient; e_tIs the tangent modulus of elasticity; p_aThe atmospheric pressure is adopted, k is the initial elastic modulus when the confining pressure is 100KPa, and m is the slope of the curve of the elastic modulus and the consolidation pressure and is between 1.0 and 0.2; c is the cohesive force of the soil body;the internal friction angle of the soil body; r_fAs a destruction ratio; g₁For confining pressure of 100KP_aAn initial poisson's ratio of (a); D. f is a soil test constant; epsilon₁Is the axial strain; sigma₁Is the maximum principal stress; sigma₃Is the minimum principal stress.

5. The method for simulating and calculating the three-dimensional variable parameters of the ground settlement and the ground cracks caused by the load and the underground water exploitation according to claim 1, which is characterized in that: in step 8, the GWS system comprising the Biot consolidation theory is improved in the following way when the combined action of overburden load and pore water pressure reduction of rock and soil mass is considered:

in the formula: g is the shear modulus,e is elastic modulus, and v is Poisson's ratio; w is a_x,w_y,w_zIs the displacement in x, y, z direction; k is a radical of_x,k_y,k_zRespectively permeability coefficient in three directions, gamma_wIs the severity of the water; gamma is the rock soil gravity;

initial conditions of the ground stress, displacement and pore water pressure of the numerical model were as follows:

in the formula: sigma_xAnd σ_zThe initial horizontal and vertical stress of the soil body; gamma is the rock soil gravity; z is the calculated point depth; k₁Is the static side pressure coefficient; w is a₀(x, y, z) is the initial displacement within the region of interest; u. of₀(x, y, z) is the initial pore water pressure within the study area;

the pore water pressure boundary conditions for the numerical model are as follows:

in the formula: u (x, y, z) is the pore water pressure at each point; u. of_sIs the head boundary gamma₁Known pore water pressure above;

the boundary conditions for the flow of the numerical model are as follows:

in the formula: k is the permeability coefficient; h is a water head; n is₁Is a normal vector;is a boundary F₂A known flow per unit area;

the free-surface boundary conditions of the numerical model are as follows:

u＝γ_wz₁

in the formula: u is free surface pore water pressure; gamma ray_wIs the severity of the water; z is a radical of₁Is the elevation of the free surface; q is the flow per unit area through the free surface boundary; mu is the soil body water supply degree; theta is an intersection angle of the normal direction outside the free surface and the vertical line;

the displacement boundary conditions for the numerical model are as follows:

in the formula (I), the compound is shown in the specification,is a displacement boundary gamma₄Known displacements in the upper three directions.

6. The method for simulating and calculating the three-dimensional variable parameters of the ground settlement and the ground cracks caused by the load and the underground water exploitation according to claim 1, which is characterized in that: step 2 also comprises the step of checking the quality of the rock-soil body grid after the grid subdivision of the research area is completed, if the quality is qualified, step 3 is started, and if the quality is not qualified, the rock-soil body grid needs to be subdivided again until the quality is qualified.

Technical Field

The invention relates to the field of geological environment protection and restoration, in particular to a three-dimensional variable parameter fully-coupled simulation calculation method technology for ground settlement and ground cracks caused by load and underground water exploitation.

Background

The ground settlement is a phenomenon that the ground elevation is slowly reduced within a certain range due to the compression of an underground rock-soil body framework under the influence of natural conditions or human activities, and the ground settlement seriously influences the safety of buildings such as high-rise buildings, highways, high-speed rails and the like. Ground fractures are also a derivative of ground subsidence. Among the many factors that influence, the load increase of underground water over-exploitation and engineering is the most dominant factor that causes ground subsidence. The existing prediction method for ground settlement mainly comprises a regression prediction model based on the statistics of underground water exploitation amount and ground settlement amount and a partial fluid-solid coupling model based on the coupling of underground water and a rock-soil body framework. However, the two models have certain advantages in researching the influence of underground water mining on ground settlement, and when the overlying load of the rock-soil mass exists, the precision is difficult to guarantee. The traditional fluid-solid coupling model mainly considers the superposition effect of a seepage field and a stress field, and does not fully consider the three-dimensional full-coupling model of rock-soil body parameter change and constitutive model change in the processes of load and underground water migration. Therefore, in order to more accurately predict the evolution of the seepage field and the change of the ground settlement amount in the underground water mining process, a fully coupled model considering the variation of geotechnical parameters and the variation of a constitutive model is used for predicting the ground settlement amount.

Disclosure of Invention

The invention provides a seepage field and stress field full-coupling simulation method considering dual effects of rock-soil body overburden load and underground water mining, aiming at solving the problem of ground settlement caused by underground water level reduction and rock-soil body overburden load increase in the process of continuous underground water mining.

In order to realize the current, the invention adopts the following technical scheme: a three-dimensional variable parameter fully-coupled simulation calculation method for ground settlement and ground cracks caused by load and underground water exploitation comprises the following steps:

step 1: determining the range, elevation, vertical stratum thickness and lithology conditions of the research area according to the range, hydrogeological conditions and engineering geological conditions of the research area, determining the layering characteristics of each stratum, and establishing a three-dimensional hydrogeological and engineering geological conceptual model of the research area;

step 2: carrying out finite element mesh subdivision on hydrogeology and engineering geology conceptual models, subdividing the conceptual models into structural unit meshes, establishing finite element mesh models, carrying out parameter partition on the finite element mesh models, deriving node and unit information of effective meshes of all partitions, respectively selecting various boundary units and nodes in the finite element mesh models, respectively deriving the boundary units and the node information, and manufacturing different boundary types;

and step 3: assigning different attribute values to the different types of boundaries in the step 2, writing the boundary parameters of the model into the derived unit and node information of each boundary, and generating boundary conditions;

and 4, step 4: writing the space coordinate positions and the initial water levels of all the nodes into the node information by using all the effective node information derived in the step 2 to generate an initial water level condition;

and 5: screening out unit and node information of all pumping wells and observation wells in underground water exploitation by using all effective unit and node information of the model derived in the step 2, generalizing the pumping wells into a drainage substructure, writing the flow of the pumping wells into the node information of the pumping wells, and finally writing the water level value of the observation wells into the node information of the observation wells to generate pumping well conditions and observation well conditions;

step 6: screening out rock-soil body node information for bearing loads by using the effective node information derived in the step 2, writing the additional loads of all nodes into the node information, and generating rock-soil body load conditions, wherein the ground load of each node is determined by the overlying load and the additional stress coefficient;

and 7: inputting permeability coefficient, porosity, water supply degree/water storage rate data with nonlinear change in each parameter partition into the parameter partition by using all effective nodes and parameter partition information derived in the step 2 to generate hydrogeological conditions, writing the elastic modulus and Poisson ratio with nonlinear change in each parameter partition into the parameter partition, and generating engineering geological conditions;

and 8: importing the files in the steps 3 to 7 into a groundwater system (GWS) based on a Biot consolidation theory to complete the establishment of a numerical model, carrying out time dispersion by using a Galerkin weighted allowance method, and carrying out finite element solution on the established numerical model by using a pretreatment conjugate gradient method (PCG);

and step 9: and (3) comparing the underground water level change obtained by calculation in the step (8) with the water level of the actual pumping test and the hydrological long observation hole water level: if the difference value between the calculated value and the observed value is more than 10% of the water level change value, returning to the step 8 to adjust the parameters until the result is in accordance, and storing the model; and finally, predicting the water level and the ground settlement of the seepage field by using the adjusted model.

Further, it is characterized in that: the well processing method in the step 4 is a drainage substructure method, nodes of the drainage substructure are divided into 4 layers and are numbered again, the number of nodes on the inner layer is 9-12, the number of nodes on the secondary inner layer is 5-8, the number of nodes on the outer layer is 1-4, the number of nodes on the outermost layer is 13-16, 9-12 are drainage hole nodes of a known water head, 1-8 nodes are unknown water head nodes on the secondary inner layer and the outer layer, and the nodes form an internal conduction matrix based on well flow; the outermost 13-16 nodes are an external conductive matrix for connecting the internal and external structures, and an overall conductive matrix [ K ] composed of the internal conductive matrix and the external conductive matrix^s]And corresponding flow arrays q^sFinally, the total conduction matrix of the drainage substructure and the corresponding flow array are:

in the formula: [ K ]_ii]A conduction matrix among unknown water head nodes in the drainage substructure; [ K ]_ib]、[K_bi]Respectively are conduction matrixes between nodes of an internal unknown water head and outlet nodes; [ K ]_bb]A conduction matrix between unknown water head nodes at the outlet of the drainage substructure; { h_iThe drainage substructure is an internal known node water head array; { h_bThe drainage substructure is a node water head array at an outlet of the drainage substructure; { q ] q_iThe flow rate of the known node water head to the internal unknown node water head is obtained; { q ] q_bAnd the flow rate of the known node water head to the unknown node water head of the outlet is obtained.

Further, it is characterized in that: in the step 6, the overburden load of the rock and soil mass is accurately described, an additional stress coefficient is introduced, and the final overburden load calculation formula of the rock and soil mass is as follows:

P＝P_d·a_sr

in the formula: p_dBuilding load for the ground; a is_srFor additional stress coefficients, the calculation is as follows:

wherein, l is the long side of the rectangular load, b is the short side of the rectangular load, x, y and z are the coordinate values of the calculation point, a_srIs an additional stress factor.

Further, the porosity, permeability coefficient, elastic modulus and poisson's ratio parameters generated in the step 7 are characterized by nonlinear change, and the parameter change equation is as follows:

in the formula: n is porosity; n is₀Is the initial porosity; Δ P is the pore water pressure change; k_sThe bulk elastic compression modulus of the solid particles of the porous medium skeleton; epsilon_vIs a bulk strain; k is the permeability coefficient; k₀Is the initial permeability coefficient; e_tIs the tangent modulus of elasticity; p_aThe atmospheric pressure is adopted, k is the initial elastic modulus when the confining pressure is 100KPa, and m is the slope of the curve of the elastic modulus and the consolidation pressure and is between 1.0 and 0.2; c is the cohesive force of the soil body;the inner friction angle of the soil body; r_fAs a destruction ratio; g₁For confining pressure of 100KP_aAn initial poisson's ratio of (a); D. f is a soil test constant; epsilon₁Is the axial strain; sigma₁Is the maximum principal stress; sigma₃Is the minimum principal stress.

Further: in step 8, the GWS system comprising the Biot consolidation theory is improved in the following way when the combined action of overburden load and pore water pressure reduction of rock and soil mass is considered:

in the formula: g is the shear modulus,e is elastic modulus, and v is Poisson's ratio; w is a_x,w_y,w_zIs the displacement in x, y, z direction; k is a radical of_x,k_y,k_zRespectively permeability coefficient in three directions, gamma_wIs the severity of the water; gamma is the rock soil gravity;

initial conditions of the ground stress, displacement and pore water pressure of the numerical model were as follows:

in the formula: sigma_xAnd σ_zThe initial horizontal and vertical stress of the soil body; gamma is the rock soil gravity; z is the calculated point depth; k₁Is the static side pressure coefficient; w is a₀(x, y, z) is the initial displacement within the region of interest; u. of₀(x, y, z) is the initial pore water pressure within the study area;

the pore water pressure boundary conditions for the numerical model are as follows:

u(x,y,z)|_Γ1＝u_s

in the formula: u (x, y, z) is the pore water pressure at each point; u. of_sIs the head boundary gamma₁Known pore water pressure above;

the boundary conditions for the flow of the numerical model are as follows:

in the formula: k is the permeability coefficient; h is a water head; n is₁Is a normal vector; q. q.s_LIs a boundary F₂A known flow per unit area;

the free-surface boundary conditions of the numerical model are as follows:

u＝γ_wz₁

in the formula: u is free surface pore water pressure; gamma ray_wIs the severity of the water; z is a radical of₁Is the elevation of the free surface; q is the flow per unit area passing through the self-contained surface boundary; mu is the soil body water supply degree; theta is an intersection angle of the normal direction outside the free surface and the vertical line;

the displacement boundary conditions for the numerical model are as follows:

in the formula (I), the compound is shown in the specification,is a displacement boundary gamma₄Known displacements in the upper three directions.

Further: step 2 also comprises the steps of checking the quality of the rock-soil body grids after the grids of the research area are divided, starting step 3 if the grids are qualified, and dividing again if the grids are unqualified until the grids are qualified.

Has the advantages that: compared with the prior art, the technical scheme of the invention has the following technical effects:

1. the calculation method of the invention establishes a ground settlement three-dimensional fluid-solid coupling numerical model by utilizing a Biot consolidation theory and a groundwater seepage theory. The calculation method fully considers the influence of the overburden load of the rock-soil body and the underground water mining coupling effect on the skeleton property of the rock-soil body, and the deformation rule of the rock-soil body in the ground settlement process is different from the single superposition of the seepage field and the stress field in the fixed parameter, so that the complete coupling of the stress field and the seepage field is realized from the angle of the change of the property of the rock-soil body caused by the underground water mining. The invention provides a basis for researching seepage field evolution and ground settlement evolution rules and fluid-solid coupling action mechanisms in the underground water exploitation process.

2. The calculation method of the invention describes the property and the seepage field rule of the rock-soil mass finely, effectively ensures the mass conservation of the rock-soil mass during overburden load and underground water mining, has higher quality of the subdivision grid, and improves the precision of the simulation calculation result.

Drawings

FIG. 1 is a schematic flow chart of a load and ground water mining induced ground settlement fully-coupled simulation calculation method according to an embodiment of the invention;

FIG. 2 is a schematic diagram of a drainage substructure of a well in accordance with an embodiment of the present invention.

Detailed Description

A model diagram of ground settlement caused by load and underground water exploitation is shown in figure 2, and in the process of underground water exploitation, the pore water pressure of a diving aquifer is gradually reduced, and the pressure borne by a rock-soil framework is increased. In addition, the earth surface is influenced by the overlying load, the rock-soil framework can be further compressed, and ground settlement is formed under the double action. Because the underground water and the ground load start synchronously, the interaction between the underground water and the rock-soil body is not simply superposed, but is a coupling effect, namely, the rock-soil body is compressed to reduce the porosity of the aquifer, the circulation of the underground water is influenced, and the deformation rule of the rock-soil body is further influenced after the underground water path-supplementing discharge condition is changed. In addition, when the overlying loads are different, the deformation characteristics of rock-soil bodies are greatly different. Thus, in different situations, the interaction between different physical fields creates different fluid-solid coupling problems. The invention is further illustrated below with reference to the figures and examples.

Example 1: the embodiment provides a load and ground settlement caused by underground water exploitation and a ground crack three-dimensional variable parameter fully-coupled simulation calculation method, and a model diagram of the load and ground settlement caused by the underground water exploitation is shown in figure 2 aiming at the influence of water flow and load on a rock-soil body structure in the underground exploitation. The overall flow of the calculation method of the seepage field-stress field full coupling is shown in fig. 1 and is realized by the following steps.

Step 1: selecting a certain underground water mining area, and establishing a hydrogeology and engineering geology conceptual model of the research area. And carrying out finite element mesh subdivision on the aquifer, subdividing the aquifers into a plurality of aquifer structure unit nets, and carrying out parameter partition on the subdivided effective aquifer structure unit nets. And constructing a finite element numerical model and deriving node and unit information of each partition.

Step 2: and (3) writing the boundary conditions, the pumping well conditions, the observation well conditions, the initial conditions and the rock-soil body load conditions into the node and unit information derived in the step (1). And then writing the hydrogeological parameters and the engineering geological parameters into each parameter partition. The rock-soil mass parameters are shown in table 1.

TABLE 1 summary of the main parameters for the calculation of ground subsidence

And step 3: and (3) introducing all boundary conditions, initial conditions and hydrogeological engineering geological conditions in the step (2) into a numerical model, setting a stress period and a time step, and carrying out finite element solution on the model by adopting a pretreatment conjugate gradient method (PCG) to obtain the change values of the water level and the ground settlement.

And 4, step 4: comparing the underground water level change obtained by calculation in the step 3 with the water level of the actual pumping test and the hydrological long observation hole water level: if the difference between the calculated value and the observed value is large, returning to the step (2) to adjust the parameters until the result is within the error range, and storing the model; if the calculated value and the observed value are within the error range, the model is saved.

And 5: and (4) predicting the evolution situation of the seepage field and the evolution situation of ground settlement caused under the combined action of the overlying load of the rock and soil mass and the underground water mining by using the model in the step (4).

According to the calculation method, through researching the change of hydrogeology and engineering geological conditions under the conditions of underground water exploitation and overlying load, a model for seepage field-stress field three-dimensional fully-coupled numerical calculation is established, the ground settlement in the underground water exploitation process is predicted, and a referable basis is provided for effectively controlling the ground settlement.

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and scope of the present invention are intended to be covered thereby.