Method for calculating rheological numerical value of deep covering layer panel rock-fill dam under consideration of seepage action

Abstract

The invention discloses a deep covering layer panel rock-fill dam rheological numerical calculation method considering the seepage effect. The method comprises the steps that firstly, contact deformation between a panel and a dam body is simulated; then simulating the unstable seepage process of the rock-fill dam and the foundation, and solving the seepage problem with a free surface; carrying out seepage-stress coupling calculation on the dam body material, and establishing a seepage-stress coupling control equation of the dam body material; carrying out seepage-stress coupling calculation on the overburden foundation, establishing a seepage-stress coupling control equation of the foundation, and finally obtaining a relationship among stress, strain and pore pressure to realize the seepage-stress coupling calculation; numerical simulation is conducted, and a concrete faced rockfill dam rheological numerical simulation result considering the seepage effect is obtained. The problems that in the prior art, the stress deformation rule of the concrete faced rockfill dam under the seepage-rheology coupling effect is not clear yet, and meanwhile, a seepage-rheology coupling calculation method is not perfect enough are solved.

Description

Technical Field

[0001] The present invention belongs to the technical field of rheological numerical calculation of face rockfill dams on deep overburden, and particularly relates to a rheological numerical calculation method for face rockfill dams on deep overburden considering seepage action. Background Art

[0002] A primary concern when constructing a concrete face rockfill dam on a deep overburden layer is dam deformation. This deformation primarily consists of transient and time-dependent deformation. Transient deformation is primarily caused by the dam's own weight and water loads during the impoundment period, while time-dependent deformation is primarily rheological deformation resulting from normal dam operation. Deformation measurements from displacement monitoring equipment installed on existing, functioning concrete face rockfill dams, as well as field test results, indicate that rheological deformation plays a key role in the deformation of concrete face rockfill dams. It is also the primary factor contributing to long-term dam deformation and failure of anti-seepage systems. Because the overburden foundation is compressible and generally distributed over a large area, concrete face rockfill dams constructed on overburden layers experience saturated seepage in the overburden foundation and adjacent dam bodies during normal operation. During this period, the rheological properties of the concrete face rockfill dam, influenced by pore water pressure, are particularly pronounced.

[0003] Current research on seepage-stress coupling of face rockfill dams on deep overburden foundations focuses mainly on the impact of the seepage field on the instantaneous deformation of the dam, and has not considered the impact of the seepage field in long-term deformation analysis. The stress-deformation law of face rockfill dams under the action of seepage-rheological coupling is still unclear, and the calculation method of seepage-rheological coupling is not yet perfect. Summary of the Invention

[0004] The purpose of the present invention is to provide a numerical calculation method for the rheology of a deep overburden face rockfill dam under the action of seepage, which solves the problems in the prior art that the stress and deformation laws of face rockfill dams under the action of seepage-rheology coupling are still unclear and the calculation method of seepage-rheology coupling is not yet perfect.

[0005] The technical solution adopted by the present invention is a numerical calculation method for the rheological properties of a deep overburden face rockfill dam under seepage, which is specifically implemented in the following steps:

[0006] Step 1: Simulate the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material to obtain the total strain at any point in the rockfill; simulate the contact deformation between the face plate and the dam body.

[0007] Step 2: Simulate the unstable seepage process of the rockfill dam and foundation, solve the seepage problem with free surface, overcome the instability of the overflow point of the potential seepage surface in the free surface iteration step during the seepage process, and accurately locate the seepage overflow point and free surface point;

[0008] Step 3: Perform seepage-stress coupling calculation on the dam material and establish the seepage-stress coupling control equation of the dam material; perform seepage-stress coupling calculation on the overburden foundation and establish the seepage-stress coupling control equation of the foundation. Finally, obtain the relationship between stress, strain and pore pressure to achieve seepage-stress coupling calculation;

[0009] Step 4: Perform numerical simulation to obtain the rheological numerical simulation results of the face rockfill dam considering the seepage effect.

[0010] The present invention is also characterized in that:

[0011] Step 1 is implemented as follows:

[0012] Step 1.1, calculate the plastic strain matrix;

[0013] Step 1.2, calculate the total strain at any point in the rockfill;

[0014] Step 1.3: Use thin layer contact surface elements to simulate the deformation behavior of the contact surface.

[0015] Step 1.1 is implemented as follows:

[0016] The Drucker-Prager plasticity model is used to describe the instantaneous deformation of rockfill and overburden gravel. The yield surface function of the Drucker-Prager model is:

[0017]

[0018] in,

[0019] q=σ1-σ3 (3)

[0020]

[0021] Where: F is the yield function; b is the generalized shear stress parameter; p is the average value of the three normal stresses; is the internal friction angle of the material; c is the cohesion of the material; q is the generalized shear stress; a is the ratio of the triaxial tensile and compressive yield stresses; r is the third invariant of the deviatoric stress; σ1 is the first principal stress; σ3 is the third principal stress;

[0022] Where c is calculated differently depending on the stress state. When hardening is defined by the uniaxial compressive yield stress:

[0023]

[0024] When hardening is defined by the uniaxial tension yield stress:

[0025]

[0026] When defining the hardening law, the cohesion is expressed as cohesion, and the cohesion of formula (5) and formula (6) are equal; is the material friction angle; d is the material parameter, 0.778≤d≤1.0; c, σ c , σ t is the isotropic hardening parameter;

[0027] For the non-associated flow law, the direction of the plastic strain is perpendicular to the plastic potential G:

[0028]

[0029] Where: is the plastic strain increment; C is a constant defined according to the hardening parameter; G is the plastic potential function; σ is the principal stress.

[0030] Where G = b-ptanψ; The method is as follows:

[0031]

[0032] Where: ψ is the expansion angle. The range of the expansion angle is determined by the flow law: ψ>71.50, that is, tanψ>3;

[0033] The total strain component consists of two parts: elastic strain component and plastic strain component, which can be expressed as follows:

[0034] dε=(dε ij ) e +(dε ij ) p (9)

[0035] Where: dε is the total strain component; (dε ij ) e is the elastic strain component; (dε ij ) p is the plastic strain component;

[0036] The elastic strain increment is calculated using the following formula:

[0037]

[0038] Where dσ ij , dσ kk is the stress component; μ is the Lame constant; ν is the Poisson's ratio; E is the elastic modulus; δ ij is the Kronecker function; D e is the elastic stress-strain matrix; dσ is the total stress component;

[0039] Solution of the elastic stress-strain matrix:

[0040] The stress and strain components are:

[0041] σ T ={σ x σ y σ z τ xy τ yz τ xz}ε T ={ε x ε y ε z ε xy ε yz ε xz} (11)

[0042] Where: σ is the principal stress; τ is the shear stress; ε is the principal strain;

[0043] Therefore, the elastic stress-strain matrix is expressed as:

[0044]

[0045] Where:

[0046]

[0047] The plastic strain increment is calculated using the following formula:

[0048] Assuming that the plastic strain increment is proportional to the plastic potential Q, and its deformation direction is orthogonal to the plastic potential surface, we have:

[0049]

[0050] Where dλ is a non-negative proportional constant, called the plastic multiplier; Q is the plastic potential;

[0051] Solution of the plastic stress-strain matrix:

[0052] Differentiating the function of the DP yield criterion yields:

[0053] a T dσ-Adλ=0 (17)

[0054] Among them, a T and A are defined as follows:

[0055]

[0056] The A value in the formula is determined by uniaxial test. After calculation, A=H, that is, stress is equal to the slope on the plastic deformation curve, and we get:

[0057]

[0058] Where: H represents the slope of the plastic deformation curve; d D =D e a.

[0059] Therefore, the plastic strain matrix is:

[0060]

[0061] Step 1.2 is implemented as follows:

[0062] The time-hardening rheological model is used to describe the rheological deformation of rockfill materials and overburden gravel materials. The time-hardening rheological model is as follows:

[0063]

[0064] Where: ε cr is the equivalent rheological strain; A, n, m are rheological parameters obtained by experimental fitting; σ cr The definition of σ comes from the hardening type of the meridian DP model, corresponding to different types of σ cr It is related to q,p, The functional form of is different; t is the rheological time;

[0065] The rheological model coupled with DP plasticity adopts the hyperbolic rheological flow potential function. The direction of the hyperbolic flow potential function to maintain the rheological deformation is determined by the following formula:

[0066]

[0067] Where: G cr is the hyperbolic flow potential function; q is the generalized shear stress; p is the average normal stress; θ is the internal friction angle of the material;

[0068] Combined with step 1.1, the total strain at any point in the rockfill is obtained as:

[0069] ε=(ε ij ) e +(ε ij ) p +ε cr (26)

[0070] Where: ε is the total strain; (ε ij ) e is the elastic strain; (ε ij ) p is the plastic strain.

[0071] Step 1.3 is implemented as follows:

[0072] The thin layer contact surface element is used to simulate the deformation behavior of the contact surface. The constitutive matrix of the thin layer element is:

[0073]

[0074] in,

[0075]

[0076]

[0077] Where: dσ i is the normal stress component, i=x,y,z; dε i is the positive strain component, i=x,y,z; dτ i is the shear stress component, i=x,y,z; dγ ij are shear strain components, i, j = x, y, z; D1, D2, D3, k zz 、k yz 、k zx Both represent the elastic coefficients of the elastic matrix in the anisotropic elastic constitutive model; E and μ represent the elastic modulus and Poisson's ratio, respectively; k zz0 is the initial normal stiffness of the contact surface, which can be taken as D1; v is the normal relative deformation; t is the thickness of the contact surface; k zx0 and k yz0 is the initial tangential stiffness, which can be taken as D3; c is the cohesion; f is the friction coefficient; α is the friction shear displacement angle; r is the failure ratio; τ is the contact surface shear stress, and γ is the contact surface shear strain.

[0078] Step 2 is implemented as follows:

[0079] Step 2.1, saturated seepage motion equation:

[0080] The Darcy law is used to solve the saturated seepage field. The saturated seepage motion equation based on Darcy law is expressed as:

[0081]

[0082] Where: q represents the seepage flow rate; v represents the seepage velocity; k represents the permeability coefficient; Represents the hydraulic gradient.

[0083] Step 2.2, unsaturated seepage motion equation:

[0084] The unsaturated seepage field is solved based on Darcy's law. The unsaturated seepage motion equation based on Darcy's law is expressed as:

[0085]

[0086] Step 2.3, basic differential equation of unsaturated seepage:

[0087] Take a sample with a thickness of d in unsaturated soil z Differential unit, water passes through the differential unit in the x and y directions, and the unsaturated continuity equation is obtained according to the principle of mass conservation:

[0088]

[0089] Where: θ w is the volumetric water content; v wx 、v wy represent the components of the seepage velocity in the x and y directions respectively;

[0090] The basic differential equation of unsaturated seepage is obtained by combining equations (35):

[0091]

[0092] Where: k wx 、k wy is the component of the unsaturated permeability coefficient in the x and y directions; H is the total hydraulic head;

[0093] Step 2.4, saturated-unsaturated seepage control equation:

[0094] The basic seepage equation is:

[0095]

[0096] The unsaturated seepage stress state is related to two stress variables, namely the normal stress (σ-u a ) and (u a -u w ) matrix suction, and the volumetric water content will also change as a result, expressed as:

[0097]

[0098] Where: θ w is the volumetric moisture content; is the volume change coefficient of water related to the change of normal stress; σ is the total stress; u a is the pore gas pressure; is the volume change coefficient related to the change of matrix suction; u w is the pore water pressure.

[0099] There is a relationship:

[0100]

[0101] The combined equation (37) yields:

[0102]

[0103] When the gas in the soil is continuous and there is no external load, Substituting into the above formula we get:

[0104]

[0105] Where: It represents the absolute value of the slope of the soil-water characteristic curve;

[0106] Also because and Substituting into the above formula we get:

[0107]

[0108] The above equation is the governing equation for saturated-unsaturated seepage, where ρ w is the density of water; g is the acceleration due to gravity;

[0109] Step 2.5, solution conditions:

[0110] When calculating saturated-unsaturated seepage, the corresponding initial conditions and boundary conditions are:

[0111]

[0112] Where: H = H1 is the boundary head; Γ1 is the known boundary of the head distribution; is the direction cosine of the normal line outside the boundary surface; q n is the normal flow rate per unit time boundary; Γ2 is the known flow boundary; z is the node coordinate of the seepage surface; Γ3 is the seepage surface boundary; H0 is the initial condition;

[0113] Step 2.6: Use the double-parameter penalty function cutoff negative pressure method to solve the seepage problem with free surface, and solve the problems of calculation accuracy and convergence in the seepage calculation. According to the variational principle, we can get:

[0114]

[0115] Where: I(h) is the variational functional; Ω is the spatial domain of the calculation area, including the geometric range of the entire dam body and foundation; k ij is the permeability tensor; h ,i 、h ,j is the head relative to the spatial coordinate x i 、x j The partial derivative of S s is the water storage rate; h is the total water head; Γ2∪Γ3 is the set of boundary conditions; q is the flow rate per unit area;

[0116] Using the improved cut-off negative pressure method, Substitute into the above formula and assume that z does not change with time, that is So the above formula is equivalent to:

[0117]

[0118] Where: I(p) is the energy function related to the pore water pressure field p; p is the pore water pressure; γ is the specific gravity of water; z is the position head;

[0119] Using the Galerkin method and isoparametric 8-node hexahedral elements, The governing equation of the unit is obtained:

[0120]

[0121] in

[0122]

[0123] Where: [K] e is the element stiffness matrix; {p} e is the element node pressure vector; [S] e The unit storage matrix; [G] e is the unit boundary condition matrix; {F} e is the element load vector; is the corresponding element of nodes a and b in the element stiffness matrix; γ is the specific gravity of water; N a,i 、N b,j is the unit shape function, a, b = 1 to 8, i, j = 1 to 3; k ij is the component of the permeability tensor; dΩ is the integral of the small volume within the unit; The elements of the matrix are stored for the unit; Ω c is the spatial area of the current finite element; S s is the water storage rate; N a 、N b is the shape function; is the coupling matrix element at the unit level; Γ3 is the seepage surface boundary; μ is the penalty function parameter; cosθ is the cosine of the angle between the boundary normal and the vertical direction; dΓ is the integral element on the boundary; is the equivalent nodal force of unit node a; H ε (p) is the penalty function; z ,j is the position coordinate z relative to the spatial coordinate x j The partial derivative of ; q is the flow rate per unit area per unit time;

[0124] Integrate the unit governing equations to obtain the overall finite element governing equations:

[0125]

[0126] Where: [K] is the stiffness matrix; {p} is the node pressure vector; [S] is the storage matrix; [G] is the boundary condition matrix; {F} is the load vector;

[0127] An implicit finite difference scheme is used for time, namely:

[0128]

[0129] Substituting into the above formula we get:

[0130]

[0131] Where: is the partial derivative of pore water pressure with respect to time; Δt is the time step; p t+Δt is the pore water pressure at time t+Δt; p t is the known pore water pressure value at the current time step t; {p} t+Δt is the pressure vector at time t+Δt; {p} t is the pressure vector at time t;

[0132] The incremental iteration method is used for iterative calculation. The iterative format suitable for calculation is:

[0133]

[0134] in

[0135]

[0136] Where: [A] is the system matrix; is the pressure increment vector; is the residual vector; is the equivalent nodal force increment caused by the pressure field change; {p} t is the converged pressure field vector at the end of the previous time step t; is the node force deviation caused by the coupling of the seepage field and the stress field in the current iteration step; Ω c is the spatial region of the current finite element; is the penalty function; N a,j is the shape function N a For the space coordinate x i The partial derivative of k ij are the components of the permeability tensor; is the potential energy gradient of the seepage; dΩ is the integral of the small volume in the unit; Γ2 is the known flow boundary; q is the flow rate per unit area per unit time; N a is the value of the shape function of node a on the boundary; dΓ is the integral element on the boundary;

[0137] For each calculation period, in order to make the iterative solution have sufficient accuracy, two convergence criteria are adopted, namely:

[0138] ① Unbalanced external force ΔB k / ΔB 0 ≤ε;

[0139] ②Pressure increment Δp k / Δp 0 ≤ε

[0140] Where ε is the convergence accuracy. By determining the above boundary conditions and convergence conditions, the seepage calculation results of the rockfill and foundation are obtained.

[0141] Step 3 is implemented as follows:

[0142] The fluid continuity governing equation of the Richard coupling model is:

[0143]

[0144] Where: ρ is the liquid density; C m is the water content; g is the acceleration due to gravity; S e is saturation; S is water storage coefficient; p is pressure; k s is the saturated permeability; η is the fluid viscosity; k r is the relative permeability; D is the position head; Q m is the liquid source and sink term;

[0145] The solid mechanics steady-state equilibrium control equation based on the Biot principle in the model is:

[0146]

[0147] Where: S is the water storage model; S ext is the saturation coefficient; F v is the pore pressure; S e is saturation; p is absolute pressure; p ref Reference stress level;

[0148] The continuity equation for the fluid in the poroelastic storage model is:

[0149]

[0150] Where: is the porous storage term; ρ is the liquid density; u is the Darcy velocity; Q m is the liquid source and sink term;

[0151] Poroelasticity relates fluid flow and particle deformation in porous media, resulting in the constitutive relationship between stress, strain, and pore pressure:

[0152]

[0153] Where: σ is the Cauchy stress tensor; C is the elastic matrix; ε is the strain tensor; α B is the Biot consolidation coefficient; p f is the fluid pore pressure; I is the unit matrix; M is the modulus; ζ is the change in fluid content in the porous medium; ε vol is the volume strain;

[0154] The convergence criterion is that the deformation and head of the two iterations meet the specified tolerances:

[0155]

[0156] Where: is the total water head; i is the number of iterations; ε is the stress deformation; ε1 and ε2 are the specified allowable errors, and ε1=ε2=10 -3 ;

[0157] In the process of seepage-stress coupling, when the numerical simulation at a certain moment is completed, the stress tensor is calculated as an additional isotropic term at the next moment and is regarded as an initial stress at the next moment. The entire coupling process is iteratively solved in each time period until the settlement is completed.

[0158] Step 4 is implemented as follows:

[0159] The rheological model of the rockfill material and overburden gravel material, the seepage model of the rockfill dam and foundation, and the seepage-deformation coupling analysis method of this project were calculated on the ABAQUS software platform, and a rheological calculation model of the overburden face rockfill dam considering the seepage effect was established: first, the three-dimensional geometric model of the face rockfill dam and deep overburden foundation was created using ABAQUS software, and the Drucker-Prager plasticity model of the rockfill material and overburden gravel material was set. The time-hardening rheological model of the rockfill material and overburden gravel material was implemented through user subroutines. Secondly, ABAQUS software was used to simulate the unsteady seepage process of the rockfill dam and foundation. Finally, according to the proposed seepage-rheology coupling analysis method, a coupled analysis step was created in ABAQUS software to realize the seepage-rheology coupling process of the dam material and the overburden gravel foundation material. The rheological deformation of the face rockfill dam under seepage was simulated through the established seepage-rheology coupling analysis method, and the final stress-deformation was obtained.

[0160] The beneficial effect of this invention is that it provides a numerical calculation method for the rheological behavior of deep-overburden face rockfill dams. Previous studies have not yet fully understood the stress and deformation patterns of face rockfill dams under the influence of seepage-rheological coupling, and the calculation method for seepage-rheological coupling is also incomplete. This invention considers the influence of seepage on the rheological deformation of face rockfill dams and proposes a method for indirect seepage-rheological coupling. This method can effectively solve the problems of slow calculation speed and high iteration error. Based on the numerical simulation of the rheological deformation of deep-overburden face rockfill dams, the seepage and deformation characteristics of face rockfill dams can be more accurately obtained. BRIEF DESCRIPTION OF THE DRAWINGS

[0161] Figure 1 This is a flow chart of a numerical calculation method for rheological properties of a deep overburden face rockfill dam considering seepage effects according to the present invention;

[0162] Figure 2 This is a calculation flow chart of the cutoff negative pressure method of the dual-parameter penalty function of the present invention;

[0163] Figure 3 It is a comparison chart of the calculation results and measured results of the dam settlement of the present invention. DETAILED DESCRIPTION

[0164] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0165] This invention considers the rheological numerical calculation method of a deep overburden face rockfill dam under seepage. First, the Drucker-Prager plasticity model and the time-hardening rheological model are used to describe the instantaneous and rheological deformations of the rockfill material and the overburden gravel material, respectively, to obtain the total strain at any point within the rockfill. The Desai contact surface thin layer element is then used to simulate the contact deformation between the two materials. Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. A cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. Simultaneously, an incremental iteration method is employed to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, enabling accurate positioning of the seepage overflow point and the free surface point. Based on the porous media seepage theory, the Richard coupling equation is used to perform coupled seepage-stress calculations on the dam material, establishing the seepage-stress coupling governing equation for the dam material. The poroelastic storage equation is used to perform coupled seepage-stress calculations on the overburden foundation, establishing the seepage-stress coupling governing equation for the foundation. Finally, the relationship between stress, strain, and pore pressure is obtained, realizing the coupled calculation of seepage-stress. This invention reveals the rheological mechanism of deep overburden face rockfill dams under seepage, providing a theoretical basis and technical support for the construction of face rockfill dams.

[0166] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0167] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0168] Step 1 is implemented as follows:

[0169] Step 1.1, calculate the plastic strain matrix;

[0170] Step 1.1 is implemented as follows:

[0171] The Drucker-Prager plasticity model is used to describe the instantaneous deformation of rockfill and overburden gravel. The yield surface function of the Drucker-Prager model is:

[0172]

[0173] in,

[0174] q=σ1-σ3 (3)

[0175]

[0176] Where: F is the yield function; b is the generalized shear stress parameter; p is the average value of the three normal stresses; is the internal friction angle of the material; c is the cohesion of the material; q is the generalized shear stress; a is the ratio of the triaxial tensile and compressive yield stresses; r is the third invariant of the deviatoric stress; σ1 is the first principal stress; σ3 is the third principal stress;

[0177] Where c is calculated differently depending on the stress state. When hardening is defined by the uniaxial compressive yield stress:

[0178]

[0179] When hardening is defined by the uniaxial tension yield stress:

[0180]

[0181] When defining the hardening law, the cohesion is expressed as cohesion, and the cohesion of formula (5) and formula (6) are equal; is the material friction angle; d is the material parameter, 0.778≤d≤1.0; c, σ c , σ t is the isotropic hardening parameter, which is determined by the equivalent plastic strain.

[0182] For the non-associated flow law, the direction of the plastic strain is perpendicular to the plastic potential G:

[0183]

[0184] Where: is the plastic strain increment; C is a constant defined according to the hardening parameter; G is the plastic potential function; σ is the principal stress.

[0185] Where G = b-ptanψ; The method is as follows:

[0186]

[0187] Where: ψ is the expansion angle. The range of the expansion angle is determined by the flow law: ψ>71.50, that is, tanψ>3;

[0188] The total strain component consists of two parts: elastic strain component and plastic strain component, which can be expressed as follows:

[0189] dε=(dε ij ) e +(dε ij ) p (9)

[0190] Where: dε is the total strain component; (dε ij ) e is the elastic strain component; (dε ij ) p is the plastic strain component;

[0191] The elastic strain increment is calculated using the following formula:

[0192]

[0193] Where dσ ij , dσ kk is the stress component; μ is the Lame constant; ν is the Poisson's ratio; E is the elastic modulus; δ ij is the Kronecker function; D e is the elastic stress-strain matrix; dσ is the total stress component;

[0194] Solution of the elastic stress-strain matrix:

[0195] The stress and strain components are:

[0196] σ T ={σ x σ y σ z τ xy τyz τ xz}ε T ={ε x ε y ε z ε xy ε yz ε xz} (11)

[0197] Where: σ is the principal stress; τ is the shear stress; ε is the principal strain;

[0198] Therefore, the elastic stress-strain matrix is expressed as:

[0199]

[0200] Where:

[0201]

[0202] The plastic strain increment is calculated using the following formula:

[0203] Assuming that the plastic strain increment is proportional to the plastic potential Q, and its deformation direction is orthogonal to the plastic potential surface, we have:

[0204]

[0205] Where dλ is a non-negative proportional constant, called the plastic multiplier; Q is the plastic potential;

[0206] Solution of the plastic stress-strain matrix:

[0207] Differentiating the function of the DP yield criterion yields:

[0208] a T dσ-Adλ=0 (17)

[0209] Among them, a T and A are defined as follows:

[0210]

[0211] The A value in the formula is determined by uniaxial test. After calculation, A=H, that is, stress is equal to the slope on the plastic deformation curve, and we get:

[0212]

[0213] Where: H represents the slope of the plastic deformation curve; d D =D e a.

[0214] Therefore, the plastic strain matrix is:

[0215]

[0216] Step 1.2, calculate the total strain at any point in the rockfill;

[0217] Step 1.2 is implemented as follows:

[0218] The time-hardening rheological model is used to describe the rheological deformation of rockfill materials and overburden gravel materials. The time-hardening rheological model is as follows:

[0219]

[0220] Where: ε cr is the equivalent rheological strain; A, n, m are rheological parameters obtained by experimental fitting; σ cr The definition of σ comes from the hardening type of the meridian DP model, corresponding to different types of σ cr It is related to q,p, The functional form is different; t is the rheological time.

[0221] The time hardening model is a law obtained from experiments assuming that the stress remains constant, and is obtained by controlling the independent variable σ cr The law of strain rate is obtained. In the process of time change, stress changes accordingly. Assuming that the stress in each time period remains unchanged, the time duration is divided into many time intervals. At the same time, it is assumed that the stress in each interval is a constant value. Finally, the strain rate is obtained by segmentation. Without considering the influence of the previous rheology on the subsequent rheology, assuming that the rheological strain transitions in different time periods, ε cr Mutations can occur while t remains constant. As t increases, ε cr It continues to accumulate and increase, while the strain rate is only related to the current stress state and the current time.

[0222] The rheological model coupled with DP plasticity adopts the hyperbolic rheological flow potential function. The direction of the hyperbolic flow potential function to maintain the rheological deformation is determined by the following formula:

[0223]

[0224] Where: G cr is the hyperbolic flow potential function; q is the generalized shear stress; p is the average normal stress; θ is the internal friction angle of the material;

[0225] Combined with step 1.1, the total strain at any point in the rockfill is obtained as:

[0226] ε=(ε ij ) e +(ε ij ) p +ε cr (26)

[0227] Where: ε is the total strain; (ε ij ) e is the elastic strain; (ε ij ) p is the plastic strain.

[0228] Step 1.3: Use thin layer contact surface elements to simulate the deformation behavior of the contact surface.

[0229] Step 1.3 is implemented as follows:

[0230] The thin layer contact surface element is used to simulate the deformation behavior of the contact surface. The constitutive matrix of the thin layer element is:

[0231]

[0232] in,

[0233]

[0234] Where: dσ i is the normal stress component, i=x,y,z; dε i is the positive strain component, i=x,y,z; dτ i is the shear stress component, i=x,y,z; dγ ij are shear strain components, i, j = x, y, z; D1, D2, D3, k zz 、k yz 、k zx Both represent the elastic coefficients of the elastic matrix in the anisotropic elastic constitutive model; E and μ represent the elastic modulus and Poisson's ratio, respectively; k zz0 is the initial normal stiffness of the contact surface, which can be taken as D1; v is the normal relative deformation; t is the thickness of the contact surface; k zx0 and k yz0 is the initial tangential stiffness, which can be taken as D3; c is the cohesion; f is the friction coefficient; α is the friction shear displacement angle; r is the failure ratio; τ is the contact surface shear stress, and γ is the contact surface shear strain.

[0235] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0236] like Figure 2 As shown, step 2 is implemented according to the following steps:

[0237] Step 2.1, saturated seepage motion equation:

[0238] The Darcy law is used to solve the saturated seepage field. The saturated seepage motion equation based on Darcy law is expressed as:

[0239]

[0240] Where: q represents the seepage flow rate; v represents the seepage velocity; k represents the permeability coefficient; Represents the hydraulic gradient.

[0241] Step 2.2, unsaturated seepage motion equation:

[0242] In unsaturated seepage, water spontaneously flows from high to low water potential, which is consistent with saturated seepage. Assuming that the air space is incompressible during seepage, the air space is considered a solid phase consistent with the soil, and water flows only along its own space. In this case, the seepage calculation formula for saturated soil can be applied to the seepage of unsaturated soil. Therefore, the unsaturated seepage field is solved based on Darcy's law. The unsaturated seepage motion equation based on Darcy's law is expressed as:

[0243]

[0244] Although the above seepage formula is similar to the saturated seepage formula (34), the parameters in the formula have different meanings. In unsaturated seepage, the total head is roughly equal to the position head plus the matrix potential head, where the position head is related to the height relative to the reference surface, and the matrix potential head is related to the saturation of the soil. In addition, the effect of saturation on the permeability coefficient of unsaturated soil cannot be ignored.

[0245] Step 2.3, basic differential equation of unsaturated seepage:

[0246] Take a sample with a thickness of d in unsaturated soil z Differential unit, water passes through the differential unit in the x and y directions, ignoring the deformation of the soil skeleton. Water is incompressible, so according to the principle of conservation of mass, the unsaturated continuity equation can be obtained:

[0247]

[0248] Where: θ w is the volumetric water content; v wx 、v wy represent the components of the seepage velocity in the x and y directions respectively.

[0249] The basic differential equation of unsaturated seepage is obtained by combining equations (35):

[0250]

[0251] Where: k wx 、kwy are the components of the unsaturated permeability coefficient in the x and y directions; H is the total hydraulic head.

[0252] Step 2.4, saturated-unsaturated seepage control equation:

[0253] The basic seepage equation is:

[0254]

[0255] The unsaturated seepage stress state is related to two stress variables, namely the normal stress (σ-u a ) and (u a -u w ) matrix suction, and the volumetric water content will also change as a result, expressed as:

[0256]

[0257] Where: θ w is the volumetric moisture content; is the volume change coefficient of water related to the change of normal stress; σ is the total stress; u a is the pore gas pressure; is the volume change coefficient related to the change of matrix suction; u w is the pore water pressure.

[0258] because can be treated as constants, so there is a relationship:

[0259]

[0260] The combined equation (37) yields:

[0261]

[0262] When the gas in the soil is continuous and there is no external load, Substituting into the above formula we get:

[0263]

[0264] Where: It represents the absolute value of the slope of the soil-water characteristic curve;

[0265] Also because and Substituting into the above formula we get:

[0266]

[0267] The above equation is the governing equation for saturated-unsaturated seepage. w is the density of water; g is the acceleration due to gravity.

[0268] Step 2.5, solution conditions:

[0269] When calculating saturated-unsaturated seepage, the corresponding initial conditions and boundary conditions are:

[0270]

[0271] Where: H = H1 is the boundary head; Γ1 is the known boundary of the head distribution; is the direction cosine of the normal line outside the boundary surface; q n is the normal flow rate per unit time boundary; Γ2 is the known flow boundary; z is the node coordinate of the seepage surface; Γ3 is the seepage surface boundary; H0 is the initial condition;

[0272] Step 2.6: Use the double-parameter penalty function cutoff negative pressure method to solve the seepage problem with free surface, and solve the problems of calculation accuracy and convergence in the seepage calculation. According to the variational principle, we can get:

[0273]

[0274] Where: I(h) is the variational functional; Ω is the spatial domain of the calculation area, including the geometric range of the entire dam body and foundation; k ij is the permeability tensor; h ,i 、h ,j is the head relative to the spatial coordinate x i 、x j The partial derivative of S s is the water storage rate; h is the total water head; Γ2∪Γ3 is the set of boundary conditions; q is the flow rate per unit area.

[0275] Using the improved cut-off negative pressure method, Substitute into the above formula and assume that z does not change with time, that is So the above formula is equivalent to:

[0276]

[0277] Where I(p) is the energy function related to the pore water pressure field p; p is the pore water pressure; γ is the specific gravity of water; and z is the positional head.

[0278] Using the Galerkin method and isoparametric 8-node hexahedral elements, The unit governing equation is obtained:

[0279]

[0280] in

[0281]

[0282] Where: [K] e is the element stiffness matrix; {p} e is the element node pressure vector; [S] e The unit storage matrix; [G] e is the unit boundary condition matrix; {F} e is the element load vector; is the corresponding element of nodes a and b in the element stiffness matrix; γ is the specific gravity of water; N a,i 、N b,j is the unit shape function, a, b = 1 to 8, i, j = 1 to 3; k ij is the component of the permeability tensor; dΩ is the integral of the small volume within the unit; The elements of the matrix are stored for the unit; Ω c is the spatial area of the current finite element; S s is the water storage rate; N a 、N b is the shape function; is the coupling matrix element at the unit level; Γ3 is the seepage surface boundary; μ is the penalty function parameter; cosθ is the cosine of the angle between the boundary normal and the vertical direction; dΓ is the integral element on the boundary; is the equivalent nodal force of unit node a; H ε (p) is the penalty function; z ,j is the position coordinate z relative to the spatial coordinate x j The partial derivative of ; q is the flow rate per unit area per unit time;

[0283] Integrate the unit governing equations to obtain the overall finite element governing equations:

[0284]

[0285] Where: [K] is the stiffness matrix; {p} is the node pressure vector; [S] is the storage matrix; [G] is the boundary condition matrix; {F} is the load vector;

[0286] An implicit finite difference scheme is used for time, namely:

[0287]

[0288] Substituting into the above formula we get:

[0289]

[0290] Where: is the partial derivative of pore water pressure with respect to time; Δt is the time step; p t+Δt is the pore water pressure at time t+Δt; p t is the known pore water pressure value at the current time step t; {p}t+Δt is the pressure vector at time t+Δt; {p} t is the pressure vector at time t;

[0291] The incremental iteration method is used for iterative calculation. The iterative format suitable for calculation is:

[0292]

[0293] in

[0294]

[0295] Where: [A] is the system matrix; is the pressure increment vector; is the residual vector; is the equivalent nodal force increment caused by the pressure field change; {p} t is the converged pressure field vector at the end of the previous time step t; is the node force deviation caused by the coupling of the seepage field and the stress field in the current iteration step; Ω c is the spatial region of the current finite element; is the penalty function; N a,j is the shape function N a For the space coordinate x i The partial derivative of k ij are the components of the permeability tensor; is the potential energy gradient of the seepage; dΩ is the integral of the small volume in the unit; Γ2 is the known flow boundary; q is the flow rate per unit area per unit time; N a is the value of the shape function of node a on the boundary; dΓ is the integral element on the boundary.

[0296] For each calculation period, in order to make the iterative solution have sufficient accuracy, two convergence criteria are adopted, namely:

[0297] ① Unbalanced external force ΔB k / ΔB 0 ≤ε;

[0298] ②Pressure increment Δp k / Δp 0 ≤ε

[0299] Where ε is the convergence accuracy, which is usually 10 -2 .

[0300] By determining the above boundary conditions and convergence conditions, the seepage calculation results of the rockfill and foundation are obtained.

[0301] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0302] Step 3 is implemented as follows:

[0303] The fluid continuity governing equation of the Richard coupling model is:

[0304]

[0305] Where: ρ is the liquid density; C m is the water content; g is the acceleration due to gravity; S e is saturation; S is water storage coefficient; p is pressure; k s is the saturated permeability; η is the fluid viscosity; k r is the relative permeability; D is the position head; Q m is the liquid source and sink term;

[0306] The solid mechanics steady-state equilibrium control equation based on the Biot principle in the model is:

[0307]

[0308] Where: S is the water storage model; S ext is the saturation coefficient; F v is the pore pressure; S e is saturation; p is absolute pressure; p ref Reference pressure level.

[0309] The continuity equation for the fluid in the poroelastic storage model is:

[0310]

[0311] Where: is the porous storage term; ρ is the liquid density; u is the Darcy velocity; Q m is the liquid source and sink term.

[0312] Poroelasticity relates fluid flow and particle deformation in porous media, resulting in the constitutive relationship between stress, strain, and pore pressure:

[0313]

[0314] Where: σ is the Cauchy stress tensor; C is the elastic matrix; ε is the strain tensor; α B is the Biot consolidation coefficient; p f is the fluid pore pressure; I is the unit matrix; M is the modulus; ζ is the change in fluid content in the porous medium; ε vol is the volume strain.

[0315] The convergence criterion is that the deformation and head of the two iterations meet the specified tolerances:

[0316]

[0317] Where: is the total water head; i is the number of iterations; ε is the stress deformation; ε1 and ε2 are the specified allowable errors, and ε1=ε2=10 -3 ;

[0318] In the process of seepage-stress coupling, when the numerical simulation at a certain moment is completed, the stress tensor is calculated as an additional isotropic term at the next moment and is regarded as an initial stress at the next moment. The entire coupling process is iteratively solved in each time period until the settlement is completed.

[0319] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0320] Step 4 is implemented as follows:

[0321] The rheological model of the rockfill material and overburden gravel material, the seepage model of the rockfill dam and foundation, and the seepage-deformation coupling analysis method of this project were calculated on the ABAQUS software platform, and a rheological calculation model of the overburden face rockfill dam considering the seepage effect was established: first, the three-dimensional geometric model of the face rockfill dam and deep overburden foundation was created using ABAQUS software, and the Drucker-Prager plasticity model of the rockfill material and overburden gravel material was set. The time-hardening rheological model of the rockfill material and overburden gravel material was implemented through user subroutines. Secondly, ABAQUS software was used to simulate the unsteady seepage process of the rockfill dam and foundation. Finally, according to the proposed seepage-rheology coupling analysis method, a coupled analysis step was created in ABAQUS software to realize the seepage-rheology coupling process of the dam material and the overburden gravel foundation material. The rheological deformation of the face rockfill dam under seepage was simulated through the established seepage-rheology coupling analysis method, and the final stress-deformation was obtained.

[0322] Example 1

[0323] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, and is specifically implemented according to the following steps:

[0324] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0325] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0326] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0327] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0328] Example 2

[0329] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0330] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0331] Step 1 is implemented as follows:

[0332] Step 1.1, calculate the plastic strain matrix;

[0333] Step 1.1 is implemented as follows:

[0334] The Drucker-Prager plasticity model is used to describe the instantaneous deformation of rockfill and overburden gravel. The yield surface function of the Drucker-Prager model is:

[0335]

[0336] in,

[0337] q=σ1-σ3 (3)

[0338]

[0339] Where: F is the yield function; b is the generalized shear stress parameter; p is the average value of the three normal stresses; is the internal friction angle of the material; c is the cohesion of the material; q is the generalized shear stress; a is the ratio of the triaxial tensile and compressive yield stresses; r is the third invariant of the deviatoric stress; σ1 is the first principal stress; σ3 is the third principal stress;

[0340] Where c is calculated differently depending on the stress state. When hardening is defined by the uniaxial compressive yield stress:

[0341]

[0342] When hardening is defined by the uniaxial tension yield stress:

[0343]

[0344] When defining the hardening law, the cohesion is expressed as cohesion, and the cohesion of formula (5) and formula (6) are equal; is the material friction angle; d is the material parameter, 0.778≤d≤1.0; c, σ c , σ t is the isotropic hardening parameter, which is determined by the equivalent plastic strain.

[0345] For the non-associated flow law, the direction of the plastic strain is perpendicular to the plastic potential G:

[0346]

[0347] Where: is the plastic strain increment; C is a constant defined according to the hardening parameter; G is the plastic potential function; σ is the principal stress.

[0348] Where G = b-ptanψ; The method is as follows:

[0349]

[0350] Where: ψ is the expansion angle. The range of the expansion angle is determined by the flow law: ψ>71.50, that is, tanψ>3;

[0351] The total strain component consists of two parts: elastic strain component and plastic strain component, which can be expressed as follows:

[0352] dε=(dε ij ) e +(dε ij ) p (9)

[0353] Where: dε is the total strain component; (dε ij ) e is the elastic strain component; (dε ij ) p is the plastic strain component;

[0354] The elastic strain increment is calculated using the following formula:

[0355]

[0356] Where dσ ij , dσ kk is the stress component; μ is the Lame constant; ν is the Poisson's ratio; E is the elastic modulus; δ ij is the Kronecker function; D e is the elastic stress-strain matrix; dσ is the total stress component;

[0357] Solution of the elastic stress-strain matrix:

[0358] The stress and strain components are:

[0359] σ T ={σ x σ y σ z τ xy τ yz τ xz}ε T ={ε x ε y ε z ε xy ε yz ε xz} (11)

[0360] Where: σ is the principal stress; τ is the shear stress; ε is the principal strain;

[0361] Therefore, the elastic stress-strain matrix is expressed as:

[0362]

[0363] Where:

[0364]

[0365] The plastic strain increment is calculated using the following formula:

[0366] Assuming that the plastic strain increment is proportional to the plastic potential Q, and its deformation direction is orthogonal to the plastic potential surface, we have:

[0367]

[0368] Where dλ is a non-negative proportional constant, called the plastic multiplier; Q is the plastic potential;

[0369] Solution of the plastic stress-strain matrix:

[0370] Differentiating the function of the DP yield criterion yields:

[0371] a T dσ-Adλ=0 (17)

[0372] Among them, a T and A are defined as follows:

[0373]

[0374] The A value in the formula is determined by uniaxial test. After calculation, A=H, that is, stress is equal to the slope on the plastic deformation curve, and we get:

[0375]

[0376] Where: H represents the slope of the plastic deformation curve; d D =D e a.

[0377] Therefore, the plastic strain matrix is:

[0378]

[0379] Step 1.2, calculate the total strain at any point in the rockfill;

[0380] Step 1.2 is implemented as follows:

[0381] The time-hardening rheological model is used to describe the rheological deformation of rockfill materials and overburden gravel materials. The time-hardening rheological model is as follows:

[0382]

[0383] Where: ε cris the equivalent rheological strain; A, n, m are rheological parameters obtained by experimental fitting; σ cr The definition of σ comes from the hardening type of the meridian DP model, corresponding to different types of σ cr It is related to q,p, The functional form is different; t is the rheological time.

[0384] The time hardening model is a law obtained from experiments assuming that the stress remains constant, and is obtained by controlling the independent variable σ cr The law of strain rate is obtained. In the process of time change, stress changes accordingly. Assuming that the stress in each time period remains unchanged, the time duration is divided into many time intervals. At the same time, it is assumed that the stress in each interval is a constant value. Finally, the strain rate is obtained by segmentation. Without considering the influence of the previous rheology on the subsequent rheology, assuming that the rheological strain transitions in different time periods, ε cr Mutations can occur while t remains constant. As t increases, ε cr It continues to accumulate and increase, while the strain rate is only related to the current stress state and the current time.

[0385] The rheological model coupled with DP plasticity adopts the hyperbolic rheological flow potential function. The direction of the hyperbolic flow potential function to maintain the rheological deformation is determined by the following formula:

[0386]

[0387] Where: G cr is the hyperbolic flow potential function; q is the generalized shear stress; p is the average normal stress; θ is the internal friction angle of the material;

[0388] Combined with step 1.1, the total strain at any point in the rockfill is obtained as:

[0389] ε=(ε ij ) e +(ε ij ) p +ε cr (26)

[0390] Where: ε is the total strain; (ε ij ) e is the elastic strain; (ε ij ) p is the plastic strain.

[0391] Step 1.3: Use thin layer contact surface elements to simulate the deformation behavior of the contact surface.

[0392] Step 1.3 is implemented as follows:

[0393] The thin layer contact surface element is used to simulate the deformation behavior of the contact surface. The constitutive matrix of the thin layer element is:

[0394]

[0395] in,

[0396]

[0397] Where: dσ i is the normal stress component, i=x,y,z; dε i is the positive strain component, i=x,y,z; dτ i is the shear stress component, i=x,y,z; dγ ij are shear strain components, i, j = x, y, z; D1, D2, D3, k zz 、k yz 、k zx Both represent the elastic coefficients of the elastic matrix in the anisotropic elastic constitutive model; E and μ represent the elastic modulus and Poisson's ratio, respectively; k zz0 is the initial normal stiffness of the contact surface, which can be taken as D1; v is the normal relative deformation; t is the thickness of the contact surface; k zx0 and k yz0 is the initial tangential stiffness, which can be taken as D3; c is the cohesion; f is the friction coefficient; α is the friction shear displacement angle; r is the failure ratio; τ is the contact surface shear stress, and γ is the contact surface shear strain.

[0398] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0399] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0400] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0401] Example 3

[0402] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0403] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0404] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0405] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0406] Step 3 is implemented as follows:

[0407] The fluid continuity governing equation of the Richard coupling model is:

[0408]

[0409] Where: ρ is the liquid density; C m is the water content; g is the acceleration due to gravity; S e is saturation; S is water storage coefficient; p is pressure; k s is the saturated permeability; η is the fluid viscosity; k r is the relative permeability; D is the position head; Q m is the liquid source and sink term;

[0410] The solid mechanics steady-state equilibrium control equation based on the Biot principle in the model is:

[0411]

[0412] Where: S is the water storage model; S ext is the saturation coefficient; F vis the pore pressure; S e is saturation; p is absolute pressure; p ref Reference pressure level.

[0413] The continuity equation for the fluid in the poroelastic storage model is:

[0414]

[0415] Where: is the porous storage term; ρ is the liquid density; u is the Darcy velocity; Q m is the liquid source and sink term.

[0416] Poroelasticity relates fluid flow and particle deformation in porous media, resulting in the constitutive relationship between stress, strain, and pore pressure:

[0417]

[0418] Where: σ is the Cauchy stress tensor; C is the elastic matrix; ε is the strain tensor; α B is the Biot consolidation coefficient; p f is the fluid pore pressure; I is the unit matrix; M is the modulus; ζ is the change in fluid content in the porous medium; ε vol is the volume strain.

[0419] The convergence criterion is that the deformation and head of the two iterations meet the specified tolerances:

[0420]

[0421] Where: is the total water head; i is the number of iterations; ε is the stress deformation; ε1 and ε2 are the specified allowable errors, and ε1=ε2=10 -3 ;

[0422] In the process of seepage-stress coupling, when the numerical simulation at a certain moment is completed, the stress tensor is calculated as an additional isotropic term at the next moment and is regarded as an initial stress at the next moment. The entire coupling process is iteratively solved in each time period until the settlement is completed.

[0423] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0424] Example 4

[0425] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0426] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0427] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0428] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0429] Step 3 is implemented as follows:

[0430] The fluid continuity governing equation of the Richard coupling model is:

[0431]

[0432] Where: ρ is the liquid density; C m is the water content; g is the acceleration due to gravity; S e is saturation; S is water storage coefficient; p is pressure; k s is the saturated permeability; η is the fluid viscosity; k r is the relative permeability; D is the position head; Q m is the liquid source and sink term;

[0433] The solid mechanics steady-state equilibrium control equation based on the Biot principle in the model is:

[0434]

[0435] Where: S is the water storage model; S ext is the saturation coefficient; F v is the pore pressure; Se is saturation; p is absolute pressure; p ref Reference pressure level.

[0436] The continuity equation for the fluid in the poroelastic storage model is:

[0437]

[0438] Where: is the porous storage term; ρ is the liquid density; u is the Darcy velocity; Q m is the liquid source and sink term.

[0439] Poroelasticity relates fluid flow and particle deformation in porous media, resulting in the constitutive relationship between stress, strain, and pore pressure:

[0440]

[0441] Where: σ is the Cauchy stress tensor; C is the elastic matrix; ε is the strain tensor; α B is the Biot consolidation coefficient; p f is the fluid pore pressure; I is the unit matrix; M is the modulus; ζ is the change in fluid content in the porous medium; ε vol is the volume strain.

[0442] The convergence criterion is that the deformation and head of the two iterations meet the specified tolerances:

[0443]

[0444] Where: is the total water head; i is the number of iterations; ε is the stress deformation; ε1 and ε2 are the specified allowable errors, and ε1=ε2=10 -3 ;

[0445] In the process of seepage-stress coupling, when the numerical simulation at a certain moment is completed, the stress tensor is calculated as an additional isotropic term at the next moment and is regarded as an initial stress at the next moment. The entire coupling process is iteratively solved in each time period until the settlement is completed.

[0446] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0447] Step 4 is implemented as follows:

[0448] The rheological model of the rockfill material and overburden gravel material, the seepage model of the rockfill dam and foundation, and the seepage-deformation coupling analysis method of this project were calculated on the ABAQUS software platform, and a rheological calculation model of the overburden face rockfill dam considering the seepage effect was established: first, the three-dimensional geometric model of the face rockfill dam and deep overburden foundation was created using ABAQUS software, and the Drucker-Prager plasticity model of the rockfill material and overburden gravel material was set. The time-hardening rheological model of the rockfill material and overburden gravel material was implemented through user subroutines. Secondly, ABAQUS software was used to simulate the unsteady seepage process of the rockfill dam and foundation. Finally, according to the proposed seepage-rheology coupling analysis method, a coupled analysis step was created in ABAQUS software to realize the seepage-rheology coupling process of the dam material and the overburden gravel foundation material. The rheological deformation of the face rockfill dam under seepage was simulated through the established seepage-rheology coupling analysis method, and the final stress-deformation was obtained.

[0449] Example 5

[0450] The present invention considers the rheological numerical calculation method of deep overburden face rockfill dam under seepage action, the flow chart is as follows Figure 1 As shown, please follow the steps below:

[0451] Step 1: Use the Drucker-Prager plasticity model and the time-hardening rheological model to describe the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material, respectively, and obtain the total strain at any point in the rockfill. Use the Desai contact surface thin layer element to simulate the contact deformation between the face plate and the dam body.

[0452] Step 1 is implemented as follows:

[0453] Step 1.1, calculate the plastic strain matrix;

[0454] Step 1.2, calculate the total strain at any point in the rockfill;

[0455] Step 1.3: Use thin layer contact surface elements to simulate the deformation behavior of the contact surface.

[0456] Step 2: Based on the basic theory of saturated-unsaturated seepage, the unstable seepage process of the rockfill dam and foundation is simulated. The cutoff negative pressure method with a two-parameter penalty function is used to solve the seepage problem with a free surface. At the same time, the incremental iteration method is used to overcome the instability of the overflow point of the potential seepage surface during the free surface iteration step during the seepage process, and the overflow point and free surface point are accurately located.

[0457] Step 3: Based on the deformation field calculation results obtained in step 1 and the seepage field calculation results obtained in step 2, and based on the porous media seepage theory, the Richard coupling equation is used to perform seepage-stress coupling calculation on the dam material, and the seepage-stress coupling control equation of the dam material is established; the porous elastic storage equation is used to perform seepage-stress coupling calculation on the overburden foundation, and the seepage-stress coupling control equation of the foundation is established. Finally, the relationship between stress, strain and pore pressure is obtained, and the seepage-stress coupling calculation is realized;

[0458] Step 4: The rheological deformation of the rockfill material and the overburden gravel material obtained in step 1 and the seepage field calculation results obtained in step 2 are numerically simulated on the ABAQUS software platform according to the coupling analysis method proposed in step 3 to obtain the rheological numerical simulation results of the face rockfill dam under the action of seepage.

[0459] Step 4 is implemented as follows:

[0460] The rheological model of the rockfill material and overburden gravel material, the seepage model of the rockfill dam and foundation, and the seepage-deformation coupling analysis method of this project were calculated on the ABAQUS software platform, and a rheological calculation model of the overburden face rockfill dam considering the seepage effect was established: first, the three-dimensional geometric model of the face rockfill dam and deep overburden foundation was created using ABAQUS software, and the Drucker-Prager plasticity model of the rockfill material and overburden gravel material was set. The time-hardening rheological model of the rockfill material and overburden gravel material was implemented through user subroutines. Secondly, ABAQUS software was used to simulate the unsteady seepage process of the rockfill dam and foundation. Finally, according to the proposed seepage-rheology coupling analysis method, a coupled analysis step was created in ABAQUS software to realize the seepage-rheology coupling process of the dam material and the overburden gravel foundation material. The rheological deformation of the face rockfill dam under seepage was simulated through the established seepage-rheology coupling analysis method, and the final stress-deformation was obtained.

[0461] Example 6

[0462] In order to verify the rationality and accuracy of the seepage-rheological coupling model of the face rockfill dam proposed in this paper, the Miaojiaba face rockfill dam project was used as a research example to carry out numerical simulation calculations. The stress and deformation calculation results were compared with the actual monitoring data. The comparison results are shown in Figure 2. Figure 3 .

[0463] from Figure 3As can be seen from the data, the dam's rheological deformation primarily occurs before water storage is complete, after which the rate of rheological deformation gradually decreases. Considering seepage-rheological coupling, the numerical calculation results are most consistent with the measured results, with the numerical results for the dam's internal settlement deformation index differing by 0.3 cm from the measured results. This indicates that the seepage effect also has a certain impact on the dam's long-term deformation. In this case, the incremental internal settlement caused by seepage accounts for 3.5% of the total deformation. Therefore, when considering seepage-rheological coupling, the calculated dam deformation has a small error compared to the measured results, indicating that the present invention can relatively accurately predict the seepage-rheological coupling characteristics of face rockfill dams.

[0464] This paper considers the impact of seepage on the rheology of face rockfill dams. The rheology of face rockfill dams is influenced by the rheology of the rockfill material in the dam body, and seepage directly affects the rheological properties of the foundation and part of the dam body below the infiltration line. The proposed indirect seepage-rheology coupling method effectively addresses the issues of slow calculation speed and high iteration error. Through numerical simulation of the rheological deformation of face rockfill dams with deep overburden and high depth, the seepage and deformation characteristics of face rockfill dams can be accurately determined.

Claims

1. A numerical calculation method for the rheological properties of a deep overburden face rockfill dam under seepage, characterized by: Please follow the steps below to implement it: Step 1: Simulate the instantaneous deformation and rheological deformation of the rockfill material and the overburden gravel material to obtain the total strain at any point in the rockfill; simulate the contact deformation between the face plate and the dam body; Step 2: Simulate the unstable seepage process of the rockfill dam and foundation, solve the seepage problem with free surface, and accurately locate the seepage overflow point and free surface point; Step 3: perform seepage-stress coupling calculation on the dam material and establish the seepage-stress coupling control equation of the dam material; Perform seepage-stress coupling calculation on the overburden foundation, establish the seepage-stress coupling control equation of the foundation, and finally obtain the relationship between stress, strain and pore pressure to realize the seepage-stress coupling calculation; Step 4: Perform numerical simulation to obtain the rheological numerical simulation results of the face rockfill dam considering the seepage effect.

2. The numerical calculation method for rheological properties of deep overburden face rockfill dam considering seepage according to claim 1 is characterized in that: The step 1 is specifically implemented according to the following steps: Step 1.1, calculate the plastic strain matrix; Step 1.2, calculate the total strain at any point in the rockfill; Step 1.3: Use thin layer contact surface elements to simulate the deformation behavior of the contact surface.

3. The numerical calculation method for rheological properties of deep overburden face rockfill dam considering seepage according to claim 2 is characterized in that: The step 1.1 is specifically implemented according to the following steps: The Drucker-Prager plasticity model is used to describe the instantaneous deformation of rockfill and overburden gravel. The yield surface function of the Drucker-Prager model is: in, q=σ1-σ3 (3) Where: F is the yield function; b is the generalized shear stress parameter; p is the average value of the three normal stresses; is the internal friction angle of the material; c is the cohesion of the material; q is the generalized shear stress; a is the ratio of the triaxial tensile and compressive yield stresses; r is the third invariant of the deviatoric stress; σ1 is the first principal stress; σ3 is the third principal stress; Where c is calculated differently depending on the stress state. When hardening is defined by the uniaxial compressive yield stress: When hardening is defined by the uniaxial tension yield stress: When defining the hardening law, the cohesion is expressed as cohesion, and the cohesion of formula (5) and formula (6) are equal; is the material friction angle; d is the material parameter, 0.778≤d≤1.0; c, σ c , σ t is the isotropic hardening parameter; For the non-associated flow law, the direction of the plastic strain is perpendicular to the plastic potential G: Where: dε pl is the plastic strain increment; C is a constant defined according to the hardening parameter; G is the plastic potential function; σ is the principal stress; Where G = b-ptanψ; dε pl The method is as follows: Where: ψ is the expansion angle. The range of the expansion angle is determined by the flow law: ψ>71.50, that is, tanψ>3; The total strain component consists of two parts: elastic strain component and plastic strain component, which can be expressed as follows: dε=(dε ij ) e +(dε ij ) p (9) Where: dε is the total strain component; (dε ij ) e is the elastic strain component; (dε ij ) p is the plastic strain component; The elastic strain increment is calculated using the following formula: Where dσ ij , dσ kk is the stress component; μ is the Lame constant; ν is the Poisson's ratio; E is the elastic modulus; δ ij is the Kronecker function; D e is the elastic stress-strain matrix; dσ is the total stress component; Solution of the elastic stress-strain matrix: The stress and strain components are: s T ={σ x s y s z t xy t yz t xz }e T ={e x e y e z e xy e yz e xz } (11) Where: σ is the principal stress; τ is the shear stress; ε is the principal strain; Therefore, the elastic stress-strain matrix is expressed as: Where: The plastic strain increment is calculated using the following formula: Assuming that the plastic strain increment is proportional to the plastic potential Q, and its deformation direction is orthogonal to the plastic potential surface, we have: Where dλ is a non-negative proportional constant, called the plastic multiplier; Q is the plastic potential; Solution of the plastic stress-strain matrix: Differentiating the function of the DP yield criterion yields: a T dσ-Adλ=0 (17) Among them, a T and A are defined as follows: The A value in the formula is determined by uniaxial test. After calculation, A=H, that is, stress is equal to the slope on the plastic deformation curve, and we get: Where: H represents the slope of the plastic deformation curve; d D =D e a; Therefore, the plastic strain matrix is:

4. The numerical calculation method for rheological properties of deep overburden face rockfill dams considering seepage according to claim 3 is characterized by: The step 1.2 is specifically implemented according to the following steps: The time-hardening rheological model is used to describe the rheological deformation of rockfill materials and overburden gravel materials. The time-hardening rheological model is as follows: Where: ε cr is the equivalent rheological strain; A, n, m are rheological parameters obtained by experimental fitting; σ cr The definition of σ comes from the hardening type of the meridian DP model, corresponding to different types of σ cr It is related to q,p, The functional form of is different; t is the rheological time; The rheological model coupled with DP plasticity adopts the hyperbolic rheological flow potential function. The direction of the hyperbolic flow potential function to maintain the rheological deformation is determined by the following formula: Where: G cr is the hyperbolic flow potential function; q is the generalized shear stress; p is the average normal stress; θ is the internal friction angle of the material; Combined with step 1.1, the total strain at any point in the rockfill is obtained as: e=(e ij ) e +(e ij ) p +e cr (26) Where: ε is the total strain; (ε ij ) e is the elastic strain; (ε ij ) p is the plastic strain.

5. The numerical calculation method for rheological properties of deep overburden face rockfill dam considering seepage according to claim 4 is characterized in that: The step 1.3 is specifically implemented according to the following steps: The thin layer contact surface element is used to simulate the deformation behavior of the contact surface. The constitutive matrix of the thin layer element is: in, Where: dσ i is the normal stress component, i=x,y,z; dε i is the positive strain component, i=x,y,z; dτ i is the shear stress component, i=x,y,z; dγ ij are shear strain components, i, j = x, y, z; D1, D2, D3, k zz 、k yz 、k zx Both represent the elastic coefficients of the elastic matrix in the anisotropic elastic constitutive model; E and μ represent the elastic modulus and Poisson's ratio, respectively; k zz0 is the initial normal stiffness of the contact surface, which can be taken as D1; v is the normal relative deformation; t is the thickness of the contact surface; k zx0 and k yz0 is the initial tangential stiffness, which can be taken as D3; c is the cohesion; f is the friction coefficient; α is the friction shear displacement angle; r is the failure ratio; τ is the contact surface shear stress, and γ is the contact surface shear strain.

6. The numerical calculation method for rheological properties of deep overburden face rockfill dam considering seepage according to claim 5 is characterized in that: The step 2 is specifically implemented according to the following steps: Step 2.1, saturated seepage motion equation: Darcy's law is used to solve the saturated seepage field. The saturated seepage motion equation based on Darcy's law is expressed as: Where: q represents the seepage flow rate; v represents the seepage velocity; k represents the permeability coefficient; represents the hydraulic gradient; Step 2.2, unsaturated seepage motion equation: The unsaturated seepage field is solved based on Darcy's law. The unsaturated seepage motion equation based on Darcy's law is expressed as: Step 2.3, basic differential equation of unsaturated seepage: Take a sample with a thickness of d in unsaturated soil z Differential unit, water passes through the differential unit in the x and y directions, and the unsaturated continuity equation is obtained according to the principle of mass conservation: Where: θ w is the volumetric water content; v wx 、v wy represent the components of the seepage velocity in the x and y directions respectively; The basic differential equation of unsaturated seepage is obtained by combining equations (35): Where: k wx 、k wy is the component of the unsaturated permeability coefficient in the x and y directions; H is the total hydraulic head; Step 2.4, saturated-unsaturated seepage control equation: The basic seepage equation is: The unsaturated seepage stress state is related to two stress variables, namely the normal stress (σ-u a ) and (u a -u w ) matrix suction, and the volumetric water content will also change as a result, expressed as: Where: θ w is the volumetric moisture content; is the volume change coefficient of water related to the change of normal stress; σ is the total stress; u a is the pore gas pressure; is the volume change coefficient related to the change of matrix suction; u w is the pore water pressure; There is a relationship: The combined equation (37) yields: When the gas in the soil is continuous and there is no external load, Substituting into the above formula we get: Where: It represents the absolute value of the slope of the soil-water characteristic curve; Also because and Substituting into the above formula we get: The above equation is the governing equation for saturated-unsaturated seepage, where ρ w is the density of water; g is the acceleration due to gravity; Step 2.5, solution conditions: When calculating saturated-unsaturated seepage, the corresponding initial conditions and boundary conditions are: Where: H = H1 is the boundary head; Γ1 is the known boundary of the head distribution; is the direction cosine of the normal line outside the boundary surface; q n is the normal flow rate per unit time boundary; Γ2 is the known flow boundary; z is the node coordinate of the seepage surface; Γ3 is the seepage surface boundary; H0 is the initial condition; Step 2.6: Use the double-parameter penalty function cutoff negative pressure method to solve the seepage problem with free surface, and solve the problems of calculation accuracy and convergence in the seepage calculation. According to the variational principle, we can get: Where: I(h) is the variational functional; Ω is the spatial domain of the calculation area, including the geometric range of the entire dam body and foundation; k ij is the permeability tensor; h ,i 、h ,j is the head relative to the spatial coordinate x i 、x j The partial derivative of S s is the water storage rate; h is the total water head; Γ2∪Γ3 is the set of boundary conditions; q is the flow rate per unit area; Using the improved cut-off negative pressure method, Substitute into the above formula and assume that z does not change with time, that is So the above formula is equivalent to: Where: I(p) is the energy function related to the pore water pressure field p; p is the pore water pressure; γ is the specific gravity of water; z is the position head; Using the Galerkin method and isoparametric 8-node hexahedral elements, The governing equation of the unit is obtained: in Where: [K] e is the element stiffness matrix; {p} e is the element node pressure vector; [S] e The unit storage matrix; [G] e is the unit boundary condition matrix; {F} e is the element load vector; is the corresponding element of nodes a and b in the element stiffness matrix; γ is the specific gravity of water; N a,i 、N b,j is the unit shape function, a, b = 1 to 8, i, j = 1 to 3; k ij is the component of the permeability tensor; dΩ is the integral of the small volume within the unit; The elements of the matrix are stored for the unit; Ω c is the spatial area of the current finite element; S s is the water storage rate; N a 、N b is the shape function; is the coupling matrix element at the unit level; Γ3 is the seepage surface boundary; μ is the penalty function parameter; cosθ is the cosine of the angle between the boundary normal and the vertical direction; dΓ is the integral element on the boundary; is the equivalent nodal force of unit node a; H ε (p) is the penalty function; z ,j is the position coordinate z relative to the spatial coordinate x j The partial derivative of ; q is the flow rate per unit area per unit time; Integrate the unit governing equations to obtain the overall finite element governing equations: Where: [K] is the stiffness matrix; {p} is the node pressure vector; [S] is the storage matrix; [G] is the boundary condition matrix; {F} is the load vector; An implicit finite difference scheme is used for time, namely: Substituting into the above formula we get: Where: is the partial derivative of pore water pressure with respect to time; Δt is the time step; p t+Δt is the pore water pressure at time t+Δt; p t is the known pore water pressure value at the current time step t; {p} t+Δt is the pressure vector at time t+Δt; {p} t is the pressure vector at time t; The incremental iteration method is used for iterative calculation. The iterative format suitable for calculation is: in Where: [A] is the system matrix; is the pressure increment vector; is the residual vector; is the equivalent nodal force increment caused by the pressure field change; {p} t is the converged pressure field vector at the end of the previous time step t; is the node force deviation caused by the coupling of the seepage field and the stress field in the current iteration step; Ω c is the spatial region of the current finite element; is the penalty function; N a,j is the shape function N a For the space coordinate x i The partial derivative of k ij are the components of the permeability tensor; is the potential energy gradient of the seepage; dΩ is the integral of the small volume in the unit; Γ2 is the known flow boundary; q is the flow rate per unit area per unit time; N a is the value of the shape function of node a on the boundary; dΓ is the integral element on the boundary; For each calculation period, in order to make the iterative solution have sufficient accuracy, two convergence criteria are adopted, namely: ① Unbalanced external force ΔB k / ΔB 0 ≤ε; ②Pressure increment Δp k / Δp 0 ≤ε Where ε is the convergence accuracy. By determining the above boundary conditions and convergence conditions, the seepage calculation results of the rockfill and foundation are obtained.

7. The numerical calculation method for rheological properties of deep overburden face rockfill dams considering seepage according to claim 6 is characterized in that: The step 3 is specifically implemented according to the following steps: The fluid continuity governing equation of the Richard coupling model is: Where: ρ is the liquid density; C m is the water content; g is the acceleration due to gravity; S e is saturation; S is water storage coefficient; p is pressure; k s is the saturated permeability; η is the fluid viscosity; k r is the relative permeability; D is the position head; Q m is the liquid source and sink term; The solid mechanics steady-state equilibrium control equation based on the Biot principle in the model is: Where: S is the water storage model; S ext is the saturation coefficient; F v is the pore pressure; S e is saturation; p is absolute pressure; p ref Reference stress level; The continuity equation for the fluid in the poroelastic storage model is: Where: is the porous storage term; ρ is the liquid density; u is the Darcy velocity; Q m is the liquid source and sink term; Poroelasticity relates fluid flow and particle deformation in porous media, resulting in the constitutive relationship between stress, strain, and pore pressure: Where: σ is the Cauchy stress tensor; C is the elastic matrix; ε is the strain tensor; α B is the Biot consolidation coefficient; p f is the fluid pore pressure; I is the unit matrix; M is the modulus; ζ is the change in fluid content in the porous medium; ε vol is the volume strain; The convergence criterion is that the deformation and head of the two iterations meet the specified tolerances: Where: is the total water head; i is the number of iterations; ε is the stress deformation; ε1 and ε2 are the specified allowable errors, and ε1=ε2=10 -3 ; In the process of seepage-stress coupling, when the numerical simulation at a certain moment is completed, the stress tensor is calculated as an additional isotropic term at the next moment and is regarded as an initial stress at the next moment. The entire coupling process is iteratively solved in each time period until the settlement is completed.

8. The numerical calculation method for rheological properties of deep overburden face rockfill dams considering seepage according to claim 7 is characterized in that: The step 4 is specifically implemented according to the following steps: The rheological model of the rockfill material and overburden gravel material, the seepage model of the rockfill dam and foundation, and the seepage-deformation coupling analysis method of this project were calculated on the ABAQUS software platform, and a rheological calculation model of the overburden face rockfill dam considering the seepage effect was established: first, the three-dimensional geometric model of the face rockfill dam and deep overburden foundation was created using ABAQUS software, and the Drucker-Prager plasticity model of the rockfill material and overburden gravel material was set. The time-hardening rheological model of the rockfill material and overburden gravel material was implemented through user subroutines. Secondly, ABAQUS software was used to simulate the unsteady seepage process of the rockfill dam and foundation. Finally, according to the proposed seepage-rheology coupling analysis method, a coupled analysis step was created in ABAQUS software to realize the seepage-rheology coupling process of the dam material and the overburden gravel foundation material. The rheological deformation of the face rockfill dam under seepage was simulated through the established seepage-rheology coupling analysis method, and the final stress-deformation was obtained.

Images