Soil shrinkage cracking analysis method and device based on near-field dynamics

Abstract

The invention discloses a soil shrinkage cracking analysis method and device based on near-field dynamics. The analysis method comprises the following steps: establishing near-field dynamic models of a soil body, wherein the near-field dynamic models comprise a heat conduction model, a moisture diffusion model and a motion model; setting initial parameters and boundary conditions of the model; solving the heat conduction model and the moisture diffusion model through an explicit algorithm, solving the motion model through an implicit algorithm, and obtaining a displacement field variable value under the current time step; and based on the variable value, updating a scalar value function to obtain local damage and strain energy density, and obtaining a stress field variable value under the current time step through a near-field dynamics differential operator. According to the method, complex behaviors of the soil body under the heat-water-force coupling effect are considered, the explicit-implicit hybrid algorithm is adopted, and the dry shrinkage cracking process of the soil body can be accurately simulated.

Description

Technical Field

[0001] The present invention belongs to the technical field of geotechnical engineering, and particularly relates to a peridynamics method for analyzing soil dry shrinkage cracking. Background Technique

[0002] Near-surface soil is prone to shrinkage cracking due to severe water loss, and this process will significantly reduce the hydraulic and mechanical properties of the soil. As an important factor inducing natural disasters such as land subsidence, dam break, debris flow, landslide, and chemical substance leakage, in-depth and accurate analysis of soil dry shrinkage cracking has important practical significance in reducing engineering risks and environmental impacts.

[0003] Current analysis methods mainly adopt the continuous medium mechanics method, but it is generally based on the local theory and is difficult to accurately simulate the crack initiation and propagation process. Peridynamics (PD), as a new non-local numerical method, replaces the differential equation with an integral equation, successfully avoiding the problem of the derivative singularity at the crack tip, showing unique advantages in the simulation of soil dry shrinkage cracking, and is expected to replace the traditional continuous medium mechanics analysis method.

[0004] However, it still faces some key challenges, including: it needs to construct a coupling model of multi-physical field interaction, and the model construction is difficult; there are significant differences in time steps in the solution processes of the motion equation and the diffusion equation. If the explicit algorithm is directly used, the calculation cost increases significantly, and if the implicit algorithm is completely used, there are technical bottlenecks in matrix storage and assembly. Therefore, it is urgent to solve the above technical problems to obtain an analysis method that can accurately simulate the soil dry shrinkage cracking process through peridynamics. Summary of the Invention

[0005] To solve the problems of the existing technology, the purpose of the present invention is to propose an analysis method and an analysis device that can accurately simulate the soil dry shrinkage cracking process through peridynamics.

[0006] The technical solution of the present invention is as follows:

[0007] A peridynamics-based method for analyzing soil dry shrinkage cracking, which includes:

[0008] S1 Establish a three-dimensional solid model of the soil to be analyzed, perform grid discretization processing on the three-dimensional solid model through the orthogonal uniform discretization method to obtain its discrete model, set the material parameters of the discrete model according to the characteristics of the soil to be analyzed, and construct a multi-field coupling peridynamics model of the temperature field, moisture field, stress field, and displacement field of the discrete model. The peridynamics model includes a peridynamics heat conduction model, a peridynamics moisture diffusion model, and a peridynamics motion model;

[0009] S2 Set the initial parameters and boundary conditions of the peridynamic model according to the external environmental parameters of the soil to be analyzed. Among them, the boundary conditions of the peridynamic moisture diffusion model include the evaporation flux obtained according to the reference crop evapotranspiration formula, and obtain the initialized peridynamic model;

[0010] S3 Introduce a scalar value function and solve the initialized peridynamic model through an explicit-implicit hybrid algorithm. The solution process includes first solving the peridynamic heat conduction model and the peridynamic moisture diffusion model through the explicit algorithm to obtain the temperature field variable value and the moisture field variable value at the current time step. Based on the temperature field variable value and the moisture field variable value, then solve the peridynamic motion model through the implicit algorithm to obtain the displacement field variable value at the current time step;

[0011] S4 Update the scalar value function based on the displacement field variable value at the current time step, obtain the local damage and strain energy density generated by it, and obtain the stress field variable value at the current time step through the peridynamic differential operator;

[0012] S5 Execute S3-S4 in sequence for each time step, loop until all time steps are completed, and output the analysis results.

[0013] Preferably, the analysis results include one or more of the contour maps of the moisture field, temperature field, displacement field, stress field, local damage, and strain energy density.

[0014] Preferably, the peridynamic heat conduction model is:

[0015]

[0016] Preferably, the peridynamic moisture diffusion model is:

[0017]

[0018] Preferably, the peridynamic motion model is:

[0019]

[0020] Among them, (ρc) eff is the effective heat capacity, (ρc) eff =(1 - φ)ρ s c s + φρ ω c ω , where ρ is the density, c is the specific heat capacity, φ is the porosity, and the subscripts s and ω represent the solid phase and the liquid phase respectively; K eff represents the effective thermal conductivity, Among them, K s is the thermal conductivity of the soil, Kω is the thermal conductivity of water; T and T′ are the temperatures of any two material points x and x′ respectively; ω and ω′ are the volumetric water contents of the material points x and x′ respectively; g(ξ) is the peridynamic function, and G(ξ) is the matrix composed of peridynamic functions; t is a certain moment; H x represents the peridynamic domain of the material point x; <·> represents the peridynamic dot product; S is the volumetric heat source term; D T is the water diffusion coefficient affected by temperature, where the subscripts T and To represent any temperature and the reference temperature respectively, η T is the dynamic viscosity of water at any temperature, η To is the dynamic viscosity of water at the reference temperature, ξ is the relative position vector in the reference configuration, ξ = x′ - x; V x′ is the volume of the material point x′; D T and D′ T are the water diffusion coefficients of the material points x and x′ respectively; R is the volumetric source-sink term; ü is the acceleration; c(ξ, δ) is the improved micro-elastic modulus function, δ is the radius of the peridynamic domain; s is the bond elongation; α and β are the linear water shrinkage coefficient and thermal expansion coefficient respectively; and represent the changes in the average water content and temperature relative to the initial moment respectively; ξ + η is the relative position vector in the current configuration, ξ + η = y′ - y, where η is the relative displacement vector, and y and y′ represent the deformed positions of x and x′; n is the unit vector; b is the body force density vector;

[0021] Among them, the improved micro-elastic modulus function at non-boundaries is:

[0022]

[0023] The improved micro-elastic modulus function at boundaries is:

[0024]

[0025] Among them, E is the elastic modulus; υ is the Poisson's ratio, υ = 1 / 3 in plane stress and υ = 1 / 4 in three-dimensional and plane strain.

[0026] Preferably, the boundary conditions of the peridynamic water diffusion model are set as:

[0027]

[0028] Among them, Δ bc is the boundary layer thickness, n is the unit normal vector, q *is the evaporation flux converted from the reference evapotranspiration ET0, and ET0 is calculated as follows:

[0029]

[0030] Wherein, is the gradient of the saturation pressure curve; R n is the net radiation of the earth's surface; T a is the ambient temperature; u2 is the wind speed at 2 meters above the ground surface; e s is the saturation vapor pressure; e a is the actual vapor pressure.

[0031] More preferably, wherein:

[0032] The saturation vapor pressure is:

[0033] The actual vapor pressure is:

[0034] The gradient of the saturation pressure curve is:

[0035] The net radiation of the earth's surface is:

[0036]

[0037] Wherein, χ is the albedo of the earth; T R is the relative sunshine duration; RH is the relative humidity; S0 is the short-wave radiation at the top of the atmosphere.

[0038] Preferably, the boundary conditions set in S2 further include the boundary conditions of the peridynamic heat conduction model, as follows:

[0039] Or

[0040] Wherein, S is the volume heat source term; T bc is the temperature at the side material point; h c is the convective heat transfer coefficient; T a is the known ambient temperature; is the boundary emissivity, represents the Stefan-Boltzmann constant.

[0041] Preferably, in the solution of the peridynamic heat conduction model and the peridynamic moisture diffusion model by the explicit algorithm, the solution model of the peridynamic heat conduction model is as follows:

[0042]

[0043] Preferably, the solution model of the peridynamic moisture diffusion model is as follows:

[0044]

[0045] where Δt is the time step of the calculation, the superscripts n and n + 1 represent the nth and (n + 1)th time steps, and the subscripts j and i represent the jth and ith material points respectively; μ ij is a scalar-valued function; N i is the total number of material points within the peridynamic range of material point i; V j is the volume of material point j, which is V j =(Δx) 3 in three dimensions, and V j =h(Δx) 2 in two dimensions, and V j =A(Δx) in one dimension, where Δx is the distance between material points.

[0046] Preferably, in the solution of the peridynamic motion model by the implicit algorithm, the solution model of the peridynamic motion model is as follows:

[0047]

[0048] KU + F ext = 0;

[0049] where K is the global stiffness matrix, which is composed of the global stiffness matrices; U is the displacement column vector of all material points; F ext is the external load column vector of all material points; G is the known coefficient matrix, λ is the Lagrange multiplier, U * is the known displacement boundary condition, and T represents the transpose;

[0050] The calculation of the global stiffness matrix is as follows:

[0051]

[0052] where k is the global stiffness matrix; l, m, and n are the direction cosines in the global coordinate system, with l = (x j - x i ) / |ξ|, m = (y j - y i ) / |ξ|, and n = (z j - z i ) / |ξ|.

[0053] According to the analysis method described in claim 1, wherein the scalar-valued function is set as follows:

[0054]

[0055] where s is the elongation rate of the bond, and s = (|ξ + η| - |ξ|) / |ξ|; s c is the critical elongation rate;

[0056] The local damage is obtained through the following calculation model:

[0057]

[0058] The strain energy density is obtained through the following calculation model:

[0059]

[0060] Preferably, the stress field variable value includes one or more of the following:

[0061] The deformation gradient is calculated as:

[0062] The strain tensor is calculated as:

[0063] The stress tensor is calculated as: σ = C:(ε - ε T - ε ω );

[0064] where F is the peridynamic deformation gradient; g(ξ) is the peridynamic function; denotes the dyadic product; ε T and ε ω respectively represent the strain tensors caused by temperature and moisture content changes; Ι is the identity matrix; C is the fourth-order elastic tensor.

[0065] The present invention further provides a peridynamics-based soil dry shrinkage and cracking analysis device, which includes a storage medium storing an executable program for implementing the above analysis method.

[0066] The present invention has the following beneficial effects:

[0067] (1) The present invention adopts the peridynamics non-local method, considers the thermo-hydro-mechanical coupling process of the soil, and incorporates complex environmental conditions into the model by introducing the reference crop evapotranspiration formula, and can more accurately simulate the soil dry shrinkage and cracking phenomenon.

[0068] (2) The present invention adopts an explicit-implicit hybrid algorithm, in which the diffusion equation is solved explicitly and the motion equation is solved implicitly, which greatly improves the calculation efficiency and significantly reduces the calculation cost, and is suitable for simulating soil dry shrinkage and cracking. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 is a flowchart of the soil dry shrinkage and cracking analysis method in the specific implementation manner;

[0070] Figure 2 It is a schematic diagram of the soil body model in the specific implementation manner;

[0071] Figure 3 It is a schematic diagram of the evolution of the water content of the soil body at different stages in the embodiment;

[0072] Figure 4 It is a schematic diagram of the evolution of the temperature of the soil body at different stages in the embodiment;

[0073] Figure 5 It is a schematic diagram of the evolution of the displacement of the soil body in three directions at different stages in the embodiment. Among them, Figure 5 (a) is the displacement in the X direction, Figure 5 (b) is the displacement in the Y direction, Figure 5 (c) is the displacement in the Z direction;

[0074] Figure 6 It is a schematic diagram of the evolution of the damage of the soil body at different stages in the embodiment;

[0075] Figure 7 It is a schematic diagram of the evolution of the normal stress of the soil body at different stages in the embodiment;

[0076] Figure 8 It is a schematic diagram of the evolution of the strain energy density of the soil body at different stages in the embodiment. Specific implementation manner

[0077] The present invention will be described in detail below in conjunction with the embodiments and the drawings. However, it should be understood that the embodiments and the drawings are only used for an exemplary description of the present invention, and cannot constitute any limitation to the protection scope of the present invention. All reasonable transformations and combinations within the scope of the inventive concept of the present invention fall within the protection scope of the present invention.

[0078] Referring to the attached Figure 1 , in some specific implementation manners, the soil body dry shrinkage cracking analysis method of the present invention includes:

[0079] S1 Establish a three-dimensional solid model of the soil body to be analyzed, perform grid discretization processing on the three-dimensional solid model by the orthogonal uniform discretization method to obtain its discrete model, set the material parameters of the discrete model according to the characteristics of the soil body to be analyzed, and construct a multi-field coupling peridynamic model of the temperature field, moisture field, stress field and displacement field of the discrete model. The peridynamic model includes a peridynamic heat conduction model, a peridynamic moisture diffusion model and a peridynamic motion model.

[0080] In a specific embodiment, referring to the attached Figure 2, the three-dimensional solid model is a cuboid with a length of 100 mm, a width of 100 mm, and a height of 20 mm. Each grid point formed in its discrete model represents a material point that makes up the continuous medium, and the distance between material points is Δx = Δ bc = 2 mm, and the near-field radius is δ = 3Δx, where Δ bc represents the boundary layer thickness.

[0081] In a specific embodiment, the set material parameters include: Young's modulus E = 1 MPa, Poisson's ratio υ = 0.25, porosity φ = 0.4, liquid-phase density ρ ω = 1000 kg / m 3 , solid-phase density ρ s = 2680 kg / m 3 , liquid-phase thermal conductivity K ω = 0.598 W / (m·K), solid-phase thermal conductivity K s = 0.3 W / (m·K), liquid-phase specific heat capacity c ω = 4190 J / (kg·K), solid-phase specific heat capacity c s = 900 J / (kg·K), moisture diffusion coefficient D To = 2×10 -9 m 2 / s, linear moisture shrinkage coefficient α = 0.47, linear thermal expansion coefficient β = 1×10 -5 ℃ -1 .

[0082] In some specific embodiments, the peridynamic model includes:

[0083] Peridynamic heat conduction model:

[0084]

[0085] Peridynamic moisture diffusion model:

[0086]

[0087] and peridynamic motion model:

[0088]

[0089] In the above models, (ρc) eff is the effective heat capacity, (ρc) eff = (1 - φ)ρ s c s + φρ ω c ω , where ρ is the density, c is the specific heat capacity, φ is the porosity, and the subscripts s and ω respectively represent the solid phase and the liquid phase; K eff represents the effective thermal conductivity, where K s is the thermal conductivity of soil, and K ω is the thermal conductivity of water; T and T′ are the temperatures of any two material points x and x′ respectively; ω and ω′ are the volumetric water contents of material points x and x′ respectively; g(ξ) is the peridynamic function, and G(ξ) is a matrix composed of peridynamic functions; t is a certain moment; H x represents the peridynamic domain of material point x; <·> represents the peridynamic dot product; S is the volumetric heat source term; D T is the water diffusion coefficient affected by temperature, where the subscripts T and To represent any temperature and the reference temperature respectively, and η T is the dynamic viscosity of water at any temperature, η To is the dynamic viscosity of water at the reference temperature, ξ is the relative position vector in the reference configuration, ξ = x′ - x; V x′ is the volume of material point x′; D T and D′ T are the water diffusion coefficients of material points x and x′ respectively; R is the volumetric source-sink term; ü is the acceleration; c(ξ, δ) is the improved micro-elastic modulus function, δ is the radius of the peridynamic domain; s is the bond elongation; α and β are the linear water shrinkage coefficient and thermal expansion coefficient respectively; and represent the changes in the average water content and temperature relative to the initial moment respectively; ξ + η is the relative position vector in the current configuration, ξ + η = y′ - y, where η is the relative displacement vector, and y and y′ represent the deformed positions of x and x′; n is the unit vector; b is the body force density vector;

[0090] where the improved micro-elastic modulus function is expressed as:

[0091]

[0092] In the formula, E is the elastic modulus; υ is the Poisson's ratio, υ = 1 / 3 in plane stress and υ = 1 / 4 in three-dimensional and plane strain; the selected kernel function is (1 - (|ξ| / δ) 2 ) 2 , which reflects the attenuation of the long-range force with space and is more in line with the distribution of non-local forces; by equating the peridynamic strain energy density with the strain energy of continuum mechanics, the specific expression of the micro-elastic modulus function under different problem conditions is obtained; only the initial relative position |ξ| in the equation is the unknown, and only the initial positions of each material point need to be considered;

[0093] Preferably, the peridynamic domain at the boundary is incomplete. To correct the error, the micro-elastic modulus function at the boundary is set as:

[0094]

[0095] S2 sets the initial parameters and boundary conditions of the peridynamic model according to the external environmental parameters of the soil to be analyzed, wherein the moisture boundary condition includes the evaporation flux obtained according to the reference crop evapotranspiration formula, and an initialized peridynamic model is obtained.

[0096] In a specific embodiment, the external environmental parameters include: environmental temperature T a = 20 °C, relative humidity RH = 70%, short-wave radiation S0 at the top of the atmosphere = 37.7 MJ / (m 2 ·day), wind speed u2 at 2 m above the ground surface = 1 m / s, earth albedo χ = 0.23, relative sunshine duration T R = 0.7.

[0097] In a specific embodiment, the reference crop evapotranspiration formula is calculated as follows:

[0098]

[0099]

[0100] wherein, ET0 represents the reference evapotranspiration, and its unit is the same as that of the evaporation flux and can be converted to each other; is the gradient of the saturation pressure curve; R n is the net radiation at the earth's surface; χ is the earth albedo; T R is the relative sunshine duration; RH is the relative humidity; S0 is the short-wave radiation at the top of the atmosphere; T a is the environmental temperature; u2 is the wind speed at 2 m above the ground surface; e s is the saturation vapor pressure; e a is the actual vapor pressure.

[0101] According to this embodiment, the reference evapotranspiration is converted into the evaporation flux q * = 4×10 -8 m / s.

[0102] In a specific embodiment, the initial parameters include: soil temperature 15 °C, soil moisture content 40%.

[0103] In some specific embodiments, the boundary conditions of the peridynamic moisture diffusion model in the peridynamic model are set as:

[0104]

[0105] The boundary conditions of the peridynamic heat conduction model are set as any one of the following types:

[0106] Dirichlet boundary condition: T = T * (x, t);

[0107] Neumann boundary condition:

[0108] Robin boundary condition:

[0109] Preferably, it is set as the Robin boundary condition;

[0110] wherein, T bc is the temperature at the edge material point; S is the volume heat source term; T * is the known temperature; q * is the known heat flux vector; Δ bc is the boundary layer thickness; n is the unit normal vector; h c is the convective heat transfer coefficient; T a is the known ambient temperature; is the boundary emissivity, represents the Stefan-Boltzmann constant; the boundary condition of the moisture diffusion equation is similar to that of the heat conduction equation;

[0111] Considering the influence of heat convection and heat radiation on the soil mass in reality, it is set as the Robin boundary condition. Other boundary conditions include: moisture evaporation, heat convection, and solar radiation only act on the surface, that is, the evaporation flux is q * = 4×10 -8 m / s, the convective heat transfer coefficient is h c = 3.7W·(m 2 ·K) -1 , the boundary divergence rate is The Stefan-Boltzmann constant is The fluxes of other surfaces are zero, that is, q * = S = R = 0. In addition, the bottom surface of the soil block cannot move horizontally or vertically, that is, u x = u y = u z = 0, and the soil block cannot move normally around, that is, u x = 0 or u y = 0.

[0112] S3 introduces a scalar-valued function and solves the initialized peridynamic model through an explicit-implicit hybrid algorithm. The solution process includes first solving the peridynamic heat conduction model and the peridynamic moisture diffusion model through the explicit algorithm to obtain the temperature field variable value and the moisture field variable value at the current time step. Based on the temperature field variable value and the moisture field variable value, then solve the peridynamic motion model through the implicit algorithm to obtain the displacement field variable value at the current time step.

[0113] In some specific embodiments, the solution model of the peridynamic heat conduction model is as follows:

[0114]

[0115] The solution model of the peridynamic moisture diffusion model is as follows:

[0116]

[0117] In the above solution process of the implicit algorithm, the forces caused by the changes in temperature and moisture content can be regarded as known external loads. Considering the slow changes in temperature and moisture, the soil deformation is regarded as a quasi-static process, and the acceleration term can be ignored; then, the motion equation is expressed under the interaction of two material points as:

[0118]

[0119] In the formula, the internal force vector f int and the external force vector f ext are respectively expressed as:

[0120]

[0121] On this basis, under the small deformation assumption, the mechanical behavior of the interaction between two material points can be expressed in matrix form as:

[0122]

[0123]

[0124] In the formula, μ is a scalar-valued function; k is the global stiffness matrix; u is the displacement column vector; l, m, and n are the direction cosines in the global coordinate system, with l = (x j - x i ) / |ξ|, m = (y j - y i ) / |ξ|, and n = (z j - z i ) / |ξ|;

[0125] Assembling it into the matrix form of all material points, it is expressed as:

[0126] KU + F ext = 0;

[0127] where K is the overall stiffness matrix, which is composed of the global stiffness matrix; U is the displacement column vector of all material points; F ext is the external load column vector of all material points;

[0128] Then the solution model of the corresponding peridynamic motion model is as follows:

[0129]

[0130] where, in the formula, Δt is the computational time step; the superscripts n and n + 1 represent time steps; the subscripts i and j represent the i-th and j-th material points respectively; μ ij is a scalar-valued function; N i is the total number of material points within the near field of material point i; V j is the volume of material point j, in three dimensions V j = (Δx) 3 in two dimensions V j = h(Δx) 2 in one dimension V j = A(Δx), where Δx is the distance between material points, h is the thickness of the plate, A is the cross-sectional area; G is a known coefficient matrix, λ is the Lagrange multiplier, U * is the known displacement boundary condition, and T represents the transpose.

[0131] S4 updates the scalar-valued function based on the displacement field variable values at the current time step, obtains the local damage and strain energy density generated by it, and obtains the stress field variable values at the current time step through the peridynamic differential operator;

[0132] Specifically, using the displacement field results, the damage failure and strain energy density are calculated, including the following formulas:

[0133] The scalar-valued function is:

[0134] The local damage is:

[0135] The strain energy density is:

[0136] where s is the elongation rate of the bond, s = (|ξ + η| - |ξ|) / |ξ|; s c is the critical elongation rate

[0137] To calculate the stress field, the strain field and deformation gradient results are required. The following formula is used:

[0138] The deformation gradient is:

[0139] The strain tensor is as follows:

[0140] The stress tensor is: σ = C : (ε - ε T - ε ω ));

[0141] where F is the peridynamic deformation gradient; g(ξ) is the peridynamic function; denotes the dyadic product; ε T and ε ω respectively represent the strain tensors caused by temperature and moisture content changes; Ι is the identity matrix; C is the fourth-order elastic tensor.

[0142] In each time step of S5, S3 - S4 are sequentially executed, and the loop continues until all time steps are completed, and the contour maps of the moisture field, temperature field, displacement field, stress field, local damage, and strain energy density are output.

[0143] In a specific embodiment, the present invention can obtain the attached Figures 3 - 8 analysis results, which show the changes in soil moisture content, temperature, displacement, damage evolution, normal stress, and strain energy density. It can be seen that the analysis method of the present invention can accurately simulate the initiation and development of cracks, which not only provides strong support for in-depth understanding of the mechanical behavior of soil, but also provides a scientific basis for solving soil cracking problems in practical engineering.

[0144] Embodiments of the present invention can be implemented in the form of a method, system, or computer program product. Therefore, the present invention can be implemented in a completely hardware manner, a completely software manner, or a combination of software and hardware. In addition, the present invention can also be stored in a computer-usable storage medium including but not limited to disk memory, CD-ROM, optical memory, etc. in the form of a computer program product for computer execution.

[0145] The present invention describes its implementation method, device (system), and computer program product by combining the steps in the flowcharts or block diagrams. It should be understood that each process or block in the flowchart or block diagram, as well as the combination of multiple processes and blocks, can be implemented by computer program instructions. These computer program instructions can be transmitted to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices, so as to construct a system capable of performing specific functions according to the steps defined in the flowchart or block diagram by executing these instructions.

[0146] These computer program instructions can also be stored in a computer-readable storage medium and can direct a computer or other programmable device to operate in a predetermined manner, ultimately forming a device capable of performing the functions described in the present invention. By loading these instructions, a computer or other device can execute in accordance with the steps to achieve the various functions specified in the flowchart or block diagram.

[0147] The above embodiments are only used to illustrate the technical solutions of the present application and are not intended to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features; and these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for analyzing soil dry shrinkage cracking based on peridynamics, comprising: S1 Establish a three-dimensional solid model of the soil to be analyzed, perform grid discretization processing on the three-dimensional solid model through orthogonal uniform discretization to obtain its discrete model, set the material parameters of the discrete model according to the characteristics of the soil to be analyzed, and construct a multi-field coupling peridynamics model of the temperature field, moisture field, stress field, and displacement field of the discrete model. The peridynamics model includes a peridynamics heat conduction model, a peridynamics moisture diffusion model, and a peridynamics motion model; S2 Set the initial parameters and boundary conditions of the peridynamics model according to the external environment parameters of the soil to be analyzed. Among them, the boundary conditions of the peridynamics moisture diffusion model include the evaporation flux obtained according to the reference crop evapotranspiration formula to obtain an initialized peridynamics model; S3 Introduce a scalar value function, and solve the initialized peridynamics model through an explicit-implicit hybrid algorithm. The solution process includes first solving the peridynamics heat conduction model and the peridynamics moisture diffusion model through the explicit algorithm to obtain the temperature field variable value and moisture field variable value at the current time step. Based on the temperature field variable value and moisture field variable value, then solve the peridynamics motion model through the implicit algorithm to obtain the displacement field variable value at the current time step; S4 Based on the displacement field variable value at the current time step, update the scalar value function to obtain the local damage and strain energy density generated by it, and obtain the stress field variable value at the current time step through the peridynamics differential operator; S5 Execute S3 - S4 in sequence for each time step, loop until all time steps are completed, and output the analysis results.

2. The analysis method according to claim 1, characterized in that, The analysis results include one or more of the following: contour maps of the moisture field, temperature field, displacement field, stress field, local damage, and strain energy density.

3. The analysis method according to claim 1, wherein Among them, The peridynamics heat conduction model is: The peridynamics moisture diffusion model is: And / or, the peridynamics motion model is: where (ρc) eff is the effective heat capacity, and (ρc) eff =(1 - φ)ρ s c s + φρ ω c ω , where ρ is the density, c is the specific heat capacity, φ is the porosity, and the subscripts s and ω represent the solid phase and the liquid phase respectively; K eff represents the effective thermal conductivity, where, K s is the thermal conductivity of the soil, and K ω is the thermal conductivity of water; T and T′ are the temperatures of any two material points x and x′ respectively; ω and ω′ are the volumetric water contents of the material points x and x′ respectively; g(ξ) is the peridynamic function, and G(ξ) is a matrix composed of peridynamic functions; t is a certain moment; represents the peridynamic domain of the material point x; <·> represents the peridynamic dot product; S is the volumetric heat source term; D T is the water diffusion coefficient affected by temperature, where the subscripts T and To represent any temperature and the reference temperature respectively, and η T is the dynamic viscosity of water at any temperature, η To is the dynamic viscosity of water at the reference temperature, ξ is the relative position vector in the reference configuration, ξ = x′ - x; V x′ is the volume of the material point x′; D T and D′ T are the water diffusion coefficients of the material points x and x′ respectively; R is the volumetric source-sink term; ü is the acceleration; c(ξ, δ) is the improved micro-elastic modulus function, and δ is the radius of the peridynamic domain; s is the bond elongation; α and β are the linear water shrinkage coefficient and thermal expansion coefficient respectively; and represent the changes in the average water content and temperature respectively relative to the initial moment; ξ + η is the relative position vector in the current configuration, ξ + η = y′ - y, where η is the relative displacement vector, and y and y′ represent the deformed positions of x and x′; n is the unit vector; b is the body force density vector; Among them, the improved micro-elastic modulus function is: The improved micro-elastic modulus function at the boundary is: Among them, E is the elastic modulus; υ is the Poisson's ratio, υ = 1 / 3 in plane stress, and υ = 1 / 4 in three-dimensional and plane strain.

4. The analysis method according to claim 1, wherein The boundary conditions of the peridynamics moisture diffusion model are set as: where, Δ bc is the boundary layer thickness, n is the unit normal vector, and q * is the evaporation flux converted from the reference evapotranspiration ET0, and ET0 is calculated as follows: wherein, is the gradient of the saturation pressure curve; R n is the net radiation at the Earth's surface; T a is the ambient temperature; u2 is the wind speed at 2 meters above the ground surface; e s is the saturation vapor pressure; e a is the actual vapor pressure.

5. The analysis method according to claim 1, characterized in that The boundary conditions set in S2 also include the boundary conditions of the peridynamics heat conduction model, as follows: or Among them, S is the volumetric heat source term; T bc is the temperature at the edge material point; h c is the convective heat transfer coefficient; T a is the known ambient temperature; is the boundary emissivity, represents the Stefan-Boltzmann constant.

6. The analysis method according to claim 1, characterized in that, In the process of solving the peridynamics heat conduction model and the peridynamics moisture diffusion model through the explicit algorithm, the solution model of the peridynamics heat conduction model is as follows: The solution model of the peridynamics moisture diffusion model is as follows: where Δt is the calculated time step, the superscripts n and n + 1 represent the nth and (n + 1)th time steps, and the subscripts i and j represent the ith and jth material points respectively; μ ij is a scalar-valued function; N i is the total number of material points within the near field of material point i; V j is the volume of material point j, which is V j = (Δx) 3 in three dimensions, and V j = h(Δx) 2 in two dimensions, and V j = A(Δx) in one dimension, where Δx is the distance between material points, h is the thickness of the plate, and A is the cross-sectional area.

7. The analysis method according to claim 6, characterized in that In the process of then solving the peridynamics motion model through the implicit algorithm, the solution model of the peridynamics motion model is as follows: KU+F ext = 0; Among them, K is the overall stiffness matrix, which is composed of the global stiffness matrix; U is the displacement column vector of all material points; F ext is the external load column vector of all material points; G is the known coefficient matrix, λ is the Lagrange multiplier, and U * is the known displacement boundary condition, and T represents the transpose; The global stiffness matrix is calculated as follows: Among them, k is the global stiffness matrix; l, m, and n are the direction cosines in the global coordinate system, with l = (x j - x i ) / |ξ|, m = (y j - y i ) / |ξ|, and n = (z j - z i ) / |ξ|.

8. The analysis method according to claim 1, wherein The scalar value function is set as follows: where s is the elongation rate of the key, s = (|ξ + η| - |ξ|) / |ξ|; s c is the critical elongation rate; The local damage is obtained through the following calculation model: The strain energy density is obtained through the following calculation model:

9. The analysis method according to claim 8, wherein The stress field variable values include one or more of the following: The deformation gradient is calculated as: The strain tensor is calculated as: The stress tensor is calculated as: σ = C:(ε - ε T - ε ω ); where \(F\) is the peridynamic deformation gradient; \(g(\xi)\) is the peridynamic function; denotes the dyadic product; \(\varepsilon\) T and \(\varepsilon\) ω denote the strain tensors caused by the changes in temperature and moisture content, respectively; \(I\) is the identity matrix; \(C\) is the fourth-order elastic tensor.

10. An analysis device for soil dry shrinkage cracking based on peridynamics, which includes a storage medium storing an executable program for implementing the analysis method described in any one of claims 1-9.

Images