Cross-fault site large deformation reaction simulation method based on particle finite element method

Abstract

The invention discloses a cross-fault site large deformation reaction simulation method based on a particle finite element method, and the method comprises the following steps: S1, building a cross-fault site two-dimensional model through Python, and submitting an ABAQUS solver to carry out the first small deformation analysis calculation; s2, establishing a grid distortion early warning mechanism, accessing an ODB file generated by the ABAQUS solver, and exporting a deformed node coordinate; s3, identifying the deformed model boundary by using an Alpha Shape algorithm, and outputting boundary node coordinates to a list file; s4, creating a deformed two-dimensional model based on the obtained boundary node coordinate list file; s5, on the basis of a boundary recognition result, grid division is carried out again in the new computational domain; s6, grid information is extracted, calculation information is mapped to a new grid from an old grid through a variable mapping technology, fault dislocation displacement continues to be applied to the deformed two-dimensional model, and small deformation analysis calculation is conducted again; and S7, the step S2 to the step S6 are repeated until the preset fault dislocation amount is reached.

Description

Technical Field

[0001] The present invention relates to the field of numerical simulation and calculation in geotechnical engineering, and particularly relates to a method for simulating large deformation response of a cross-fault site based on the particle finite element method. Background Art

[0002] The simulation of seismic response of cross-fault sites is the primary problem often faced in the seismic safety analysis and evaluation of cross-fault projects. Reasonably considering the large deformation of the overlying soil layer caused by fault dislocation is the key. From the perspective of research methods, model tests can reproduce the process of disaster occurrence, but are not convenient for parametric research, while numerical simulation can effectively control the parameter values and flexibly simulate different working conditions. Among them, the finite element method is widely used. However, the traditional finite element method is proposed based on the small deformation assumption, and problems such as mesh distortion, decreased calculation accuracy, and even non-convergence are likely to occur when the analysis domain undergoes large deformation. Although the discrete element method is also applicable to solving problems with complex contacts and large deformations, the number of granular materials is often too large, resulting in very low computational analysis efficiency.

[0003] To solve the large deformation problem, various numerical methods have been proposed in recent decades. For example, the Arbitrary Lagrangian-Eulerian framework (ALE) is the most commonly used method for dealing with large deformations in geotechnical engineering. Three ALE-based methods widely used in geotechnical engineering are the small strain interpolation technique (RITSS) proposed by Hu and Randolph, the efficient ALE method (EALE) developed by Nazem, and the built-in coupled Euler-Lagrangian (CEL) method in ABAQUS. Another commonly used method is the particle-based method, in which a set of particles is used to replace the mesh used in the ALE-based method to avoid mesh distortion. Various particle-based large deformation analysis methods have been proposed and applied to geotechnical engineering, such as Smoothed Particle Hydrodynamics SPH, Material Point Method (MPM), and Particle Finite Element Method (PFEM). Among them, the Particle Finite Element Method uses particles to represent the structure, and the calculation process is performed on the mesh, combining the solid mathematical foundation of the finite element method and the advantages of the meshless method in dealing with large deformations and complex boundary problems, and is particularly suitable for simulating large deformation problems in geotechnical engineering. Although the Particle Finite Element Method has begun to show potential in the large deformation analysis of geotechnical engineering, its application to the study of the problem of rupture of overlying soil layers caused by fault dislocation in cross-fault sites is still blank. This type of problem has distinct particularities: the fault dislocation effect will induce complex non-linear mechanical responses in the overlying soil layer, involving significant shear deformation, soil strain softening, formation of localized shear bands, and dynamic expansion of the rupture surface. Currently, the application of the Particle Finite Element Method in geotechnical engineering mainly focuses on the simulation of plastic flow in continuous media, and the simulation of the soil rupture process under fault dislocation conditions has not been reported. Therefore, carrying out research on the PFEM numerical simulation technology suitable for large deformation rupture analysis in cross-fault sites has important theoretical significance and engineering value, and can provide a more accurate analysis method for the seismic design and safety assessment of major cross-fault projects. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for simulating large deformation responses in cross-fault sites based on the Particle Finite Element Method to solve the problems raised in the above background technology.

[0005] To achieve the above purpose, the present invention provides a method for simulating large deformation responses in cross-fault sites based on the Particle Finite Element Method, including the following steps:

[0006] S1. Use Python to establish a two-dimensional model of the cross-fault site and submit it to the ABAQUS solver. By applying fault dislocation displacements on the two-dimensional model to simulate the fracture dislocation process, perform the first small deformation analysis calculation;

[0007] S2. Establish a mesh distortion warning mechanism, access the ODB file generated by the ABAQUS solver, and export the deformed node coordinates;

[0008] S3. Use the Alpha Shape algorithm to identify the boundary of the deformed model and output the coordinates of the boundary nodes to a list file.

[0009] S4. Based on the obtained list file of the boundary node coordinates, create the deformed two-dimensional model.

[0010] S5. Based on the boundary recognition result, re-perform mesh generation in the new computational domain.

[0011] S6. Extract the mesh information, map the calculation information from the old mesh to the new mesh through variable mapping technology, and continue to apply the fault dislocation displacement on the deformed two-dimensional model, and perform the small deformation analysis calculation again.

[0012] S7. Repeat steps S2 to S6 until the preset fault dislocation amount is reached.

[0013] In a preferred embodiment, in step S1, simulate the fracture dislocation process and submit it to the ABAQUS solver for the first small deformation analysis calculation, including:

[0014] Before performing the fracture dislocation analysis, first complete the in-situ stress balance calculation: apply the gravity load in the vertically downward direction, create a static analysis step for stress balance calculation, and define the initial stress field that conforms to the in-situ stress distribution on site; after completing the in-situ stress balance calculation, create a new static analysis step to apply the fault dislocation displacement, with the initial step size taken as 0.1% of the model height, and use the penetration of a certain amplitude of equivalent plastic strain as the standard for soil rupture, and dynamically adjust the displacement increment step size according to the plastic zone expansion rate. The mesh element type adopts the 6-node modified quadratic plane strain triangular element, complete the settings in the ABAQUS input file, and submit it to the ABAQUS solver for the first small deformation analysis calculation.

[0015] Among them, the method of applying the fault dislocation displacement is: apply horizontal constraints and no vertical constraints on the boundary of the passive side of the fault, constrain the horizontal and vertical degrees of freedom at the bottom, apply horizontal displacement and no vertical constraints on the boundary of the active side of the fault, and apply horizontal and vertical displacements in both directions at the bottom simultaneously, and simulate different bedrock fracture dips by controlling the ratio of the vertical displacement amount to the horizontal displacement amount.

[0016] In a preferred embodiment, in step S2, establish a mesh distortion warning mechanism, including: adopt a dual-trigger mechanism for mesh distortion judgment. When the set increment step threshold is calculated or it is judged that the mesh is about to be distorted, use a script program written in the Python language to access the ODB document generated by the ABAQUS solver.

[0017] In a preferred embodiment, in step S2, access the ODB file generated by the ABAQUS solver and export the deformed node coordinates, including: S21. Open the ODB file; S22. Obtain the instances in the model; S23. Obtain the analysis step and the last frame; S24. Extract the original node coordinates and deformation data; S25. Close the ODB file; S26. Calculate the deformed node coordinates, where the deformed node coordinates = the original node coordinates + the deformation data; S27. Close the ODB file and write the deformed node coordinates to a list file.

[0018] In a preferred embodiment, in step S3, use the Alpha Shape algorithm to perform boundary recognition based on the exported deformed node coordinates, visualize the boundary points and output a list of boundary point coordinates, including: S31. Read the deformed node coordinate file exported in step S2 and change the format to the form of ID, X, Y; S32. Identify the boundary and output the boundary points, including reading the point set, selecting the α value, extracting the boundary points, visualizing the results, plotting the original points and plotting the boundary; S33. Output the boundary point coordinate list file and configure the file path.

[0019] In a preferred embodiment, in step S3, the process of selecting the α value is as follows: First, preliminarily determine the reasonable range of the α value, and then use the bisection method to iteratively optimize within this range. Calculate the boundary recognition results corresponding to the intermediate values in turn, and through the comparison error analysis with the boundary before deformation, gradually narrow the range of the α value until the accuracy requirement is met, and control the boundary recognition error within 5%.

[0020] In a preferred embodiment, in step S4, based on the obtained list file of the deformed boundary point coordinates, use the Python language to write a script program to create a deformed two-dimensional model, including: S41. Read the node coordinate file and change the format to the form of X, Y; S42. Create a sketch, draw points and connect lines in the sketch to create a new two-dimensional model.

[0021] In a preferred embodiment, in step S5, based on the boundary recognition result, re-perform mesh division in the new computational domain, including: Based on the boundary recognition result, re-perform mesh division through the built-in mesh generation method in ABAQUS. Use global seeding when re-performing mesh division, re-set the global seeding size, and the element type is a 6-node modified quadratic plane strain triangular element.

[0022] In a preferred embodiment, in step S6, grid information is extracted, calculation information is mapped from the old grid to the new grid through variable mapping technology, and fault dislocation displacement is continuously applied to the deformed two-dimensional model, and small deformation analysis and calculation are performed again, including: checking the maximum incremental number of steps in the previous analysis and calculation, calling Map Solution, modifying the input file to specify the variable mapping command, and specifying the new job and the old job in the ABAQUS command box to perform small deformation analysis and calculation again.

[0023] In a preferred embodiment, the two-dimensional model includes the thickness of the overlying soil layer, the fracture position, the dip angle of the bedrock fracture, and the fault dislocation amount. In order to reduce the influence of the boundary on the deformation of the soil body center, the width of the two-dimensional model is set to 3-4 times the height.

[0024] Compared with the prior art, the beneficial effects of the present invention are:

[0025] The present invention regards the nodes in the finite element grid as free particles, updates the positions of the particles according to the displacement increment of the analysis step, and solves the problem of grid distortion in the traditional finite element method when simulating large deformations by periodically calculating domain boundary recognition and grid re - meshing, can effectively realize the numerical simulation of large deformations of rock and soil masses, and realizes this numerical simulation technology on the finite element software ABAQUS platform through secondary development. Compared with the traditional small - deformation finite element method, the present invention can simulate a 50% increase in the fault dislocation amount, successfully breaking through the bottleneck of grid distortion and calculation non - convergence caused by large deformations. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 is the method flow chart of the present invention;

[0027] Figure 2 is the schematic diagram of the two - dimensional model of the cross - fault site of the present invention;

[0028] Figure 3 is the boundary recognition result when different values of α are taken in the present invention;

[0029] Figure 4A is the maximum vertical dislocation that can be carried out by the traditional finite element method (the deformation scaling factor is 3);

[0030] Figure 4B is the maximum vertical dislocation that can be carried out by the method of the present invention (the deformation scaling factor is 3);

[0031] Figure 5 is the boundary recognition result when the vertical dislocation in the present invention is 50 cm;

[0032] Figure 6 is the grid re - meshing when the vertical dislocation in the present invention is 50 cm;

[0033] Figure 7A The stress nephogram before mapping when the vertical dislocation of the present invention is 50 cm;

[0034] Figure 7B The stress nephogram after mapping when the vertical dislocation of the present invention is 50 cm;

[0035] Figure 8 The boundary recognition result when the vertical dislocation of the present invention is 100 cm;

[0036] Figure 9 The mesh re - division when the vertical dislocation of the present invention is 100 cm;

[0037] Figure 10A The stress nephogram before mapping when the vertical dislocation of the present invention is 100 cm;

[0038] Figure 10B The stress nephogram after mapping when the vertical dislocation of the present invention is 100 cm;

[0039] Figure 11 The stress nephogram of the overlying soil when the vertical dislocation of the present invention is 150 cm;

[0040] Figure 12 The mesh diagram after the last small - deformation calculation of the present invention. Detailed implementation manners

[0041] Next, the technical solutions in the embodiments of the present invention will be described clearly and completely. For the embodiments in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0042] Embodiment 1

[0043] As Figures 1 to 3 shown, the large - deformation response simulation method for a cross - fault site based on the particle finite element method of the present invention includes the following steps:

[0044] Step S1: Use Python to establish a two - dimensional model of the cross - fault site and submit it to the ABAQUS solver. Considering that the stiffness of the bedrock is much greater than that of the overlying soil layer, simulate the bedrock as a rigid body and input the fault dislocation displacement at the top of the bedrock to simulate the fault dislocation process.

[0045] Among them, the two - dimensional model includes the thickness of the overlying soil layer, the fracture position, the fracture dip angle of the bedrock, and the fault dislocation amount. And in order to reduce the influence of the boundary on the central deformation of the soil, set the width of the two - dimensional model to 3 - 4 times its height.

[0046] Before performing the fault dislocation analysis, the in-situ stress balance calculation is first completed, which specifically includes: applying a gravity load in the vertically downward direction, defining the initial stress field that conforms to the in-situ stress distribution using the *INITIAL CONDITIONS keyword, and creating a static analysis step containing the "Geostatic" option for stress balance calculation. After completing the in-situ stress balance calculation, a new static analysis step is created to apply the fault dislocation displacement. The initial step size is taken as 0.1% of the model height, and a certain amplitude of equivalent plastic strain penetration is used as the criterion for soil rupture. The displacement increment step size is dynamically adjusted according to the plastic zone expansion rate. The mesh is divided using 6-node modified quadratic plane strain triangular elements, which are set in the ABAQUS input file, and then the ABAQUS solver is submitted for the first small deformation analysis calculation. The method of applying the fault dislocation displacement is as Figure 2 shown. Horizontal constraints are applied to the side boundary of the passive disk of the fault and no vertical constraints are applied (allowing vertical settlement). The horizontal and vertical degrees of freedom are constrained at the bottom (fixed boundary). Horizontal displacement is applied to the side boundary of the active disk of the fault and no vertical constraints are applied (allowing vertical settlement). Displacements in both the horizontal and vertical directions are applied simultaneously at the bottom. Different bedrock fracture dip angles are simulated by controlling the ratio of the vertical displacement to the horizontal displacement.

[0047] Since there are differences in the materials of geotechnical bodies under large deformation conditions and small deformation conditions, for large deformation conditions, the material parameters based on the traditional small deformation process should be corrected and optimized to better conform to actual engineering. The model material properties of the present invention are realized through the user-defined material development interface UMAT of ABAQUS. Specifically, the Mohr-Coulomb model is used to simulate the nonlinear characteristics of geotechnical bodies, and the strain-softening phenomenon of geotechnical bodies under large deformation conditions is considered by introducing the strain-dependent friction angle and dilation angle. In the following equations and ψ mob represent the friction angle and dilation angle that vary with the plastic octahedral shear strain of the geotechnical body respectively, and is the corrected shear strain after the softening of the geotechnical body.

[0048]

[0049] Taking a certain amplitude of equivalent plastic strain penetration as the criterion for soil rupture, for Mohr-Coulomb materials, the shear plastic strain increment is related to the dilation angle, and the critical effective plastic strain can be approximated as:

[0050]

[0051] Step S2: Establish a mesh distortion warning mechanism, access the ODB file generated by the ABAQUS solver, and export the deformed node coordinates as the input data for subsequent boundary recognition.

[0052] Specifically, a grid distortion warning mechanism is established, including: when the set increment step is calculated or it is judged that the grid is about to be distorted, a script program is written in Python language to access the ODB document generated by the ABAQUS solver (as shown in Appendix 1).

[0053] The present invention establishes a dual-trigger mechanism for grid distortion judgment (increment step threshold + grid distortion criterion):

[0054] (1) Technical framework design

[0055]

[0056]

[0057] Conflict handling: When both triggering conditions are met, the grid distortion criterion is preferentially responded to;

[0058] Delay protection: Dynamically adjust the fault dislocation increment to prevent frequent triggering;

[0059] (2) Increment step threshold trigger:

[0060] Set *Restart,write,frequency=N in the ABAQUS input file to automatically save the restart file every N increment steps. The cumulative displacement generated by each small deformation calculation (N increment steps) does not exceed half of the minimum side length of the element. After each small deformation analysis, the grid still has relatively high quality, thus ensuring the calculation accuracy.

[0061] (3) Grid distortion criterion trigger:

[0062]

[0063] Export the deformed node coordinates and write them into a list file as the input data for subsequent boundary recognition. Specifically, it includes: S21. Open the ODB file; S22. Obtain the instances in the model; S23. Obtain the analysis step and the last frame; S24. Extract the original node coordinates and deformation data; S25. Close the ODB file; S26. Calculate the deformed node coordinates, where the deformed node coordinates = the original node coordinates + the deformation data; S27. Close the ODB file and write the deformed node coordinates into the list file.

[0064] Step S3: Use the Alpha Shape algorithm (as shown in Appendix 2) to identify the boundary of the deformed model, and output the boundary node coordinates to a list file for secondary modeling in ABAQUS. Specifically, it includes: S31: Read the deformed node coordinate file exported in Step S2 and change the format to the form of ID, X, Y; S32: Identify the boundary and output the boundary points, including reading the point set, selecting the α value, extracting the boundary points, visualizing the results, plotting the original points, and plotting the boundary; S33: Output the boundary point coordinate list file and configure the file path.

[0065] During the boundary recognition process, the selection of the α value directly affects the accuracy of boundary recognition: when the α value is too small, the boundary contour will show serrated distortion, while when the α value is too large, the true boundary points in the sparse area will be wrongly removed. To obtain the optimal α value, the present invention develops an adaptive α value selection algorithm based on the dichotomy method. The specific implementation process is as follows: First, initially determine a reasonable range of the α value (for example, 2.0 - 3.0), and then use the dichotomy method to iteratively optimize within this range. Calculate the boundary recognition results corresponding to the intermediate values α = 2.5, 2.25, 2.75, etc. in turn. Through the comparative error analysis with the boundary before deformation, gradually narrow the range of the α value until the accuracy requirement is met. Usually, iterating 3 - 5 times can control the boundary recognition error within 5%.

[0066] Step S4: Based on the obtained deformed boundary point coordinate list file, use the Python language to write a script program to create a deformed two-dimensional model. Specifically, it includes: S41: Read the node coordinate file and change the format to the form of X, Y; S42: Create a sketch, draw points and connect lines in the sketch, and create a new two-dimensional model.

[0067] Step S5: Based on the boundary recognition results, re - perform mesh division in the new calculation domain. It should be noted that for boundary recognition, it is required that the particle set used for recognition is evenly distributed. The density of grid nodes (i.e., particles) has a great impact on the quality of boundary recognition. Therefore, global seeding should be adopted during mesh division to reduce the non - uniformity of node distribution. Fault dislocation causes great differences in the strain conditions of different areas of the overlying soil layer. In the rupture area, the soil body deforms violently with a high strain gradient, and a fine mesh is required to accurately capture the deformation and stress concentration; in the area with small strain, the soil body deforms weakly, and a sparse mesh can meet the calculation accuracy requirements to improve the overall analysis efficiency. When re - performing mesh division, adopt global seeding, re - set the global seeding size, and the mesh type is a 6 - node modified quadratic plane strain triangular element.

[0068] Step S6: Extract grid information, map the calculation information from the old grid to the new grid through variable mapping technology, and continue to apply the fault dislocation displacement on the deformed two-dimensional model. Apply a small deformation again and perform small deformation analysis and calculation. During the solution process of large deformation problems, multiple variable mappings are required, which are implemented through the Map Solution technology built in ABAQUS. Since grid nodes may be shared by multiple elements, it is necessary to average the results extrapolated by different elements for the same node. After obtaining the calculation results of all nodes, interpolate them to the Gaussian points of the new grid through the shape function as the initial value for the next calculation. This technology cannot be directly completed through ABAQUS CAE and needs to be completed by modifying the input file after establishing the analysis model (as shown in Appendix 3). Specifically, it includes: checking the maximum number of increment steps in the previous analysis and calculation, calling Map Solution, modifying the input file to specify the variable mapping command, and specifying a new job and an old job in the ABAQUS command box to perform a small deformation analysis and calculation again.

[0069] Step S7: Repeat Step S2 to Step S6 until the preset fault dislocation amount is reached.

[0070] Embodiment 2

[0071] As Figure 1 、 Figure 4A 、 4B to Figure 12 shown, this embodiment mainly studies the process of the deformation and rupture propagation of the overlying soil layer caused by the fault dislocation of the bedrock. Since the rupture propagation of a single layer of soil has a good reference effect for understanding the deformation and rupture of multi-layer soil, the actual stratigraphic distribution is not considered temporarily in the simulation. Homogeneous single silty clay is used, and the finite element calculation geotechnical parameters are shown in Table 1.

[0072] Table 1 Parameters of the overlying soil layer in the model calculation

[0073]

[0074] A two-dimensional model of a cross-fault site is established, with a width of 240 m, a thickness of the overlying soil layer of 50 m, and a dip angle of the fractured bedrock of 45°. When the traditional small deformation finite element method is used, when the analysis step time reaches 0.9274 and the vertical dislocation amount reaches 1 m, the calculation shows non-convergence, as Figure 4A shown.

[0075] When the method of the present invention is used to export the node coordinates after the model deformation when the vertical dislocation amounts reach 50 cm and 100 cm respectively, use the Alpha Shape program for boundary recognition, and then perform grid division and grid mapping. Finally, the large deformation response of the soil layer with a fault dislocation amount of more than 1.5 m can be simulated, successfully breaking through the problems of grid distortion and calculation non-convergence caused by large deformation, as Figure 4Bas shown

[0076] The cross - fault site large - deformation response simulation method based on the particle finite - element method in this embodiment is as follows:

[0077] Modify the Python script parameters. Set the length of the two - dimensional model to 240 m, the width to 50 m, and apply a gravity load of 19.5 KN / m in the vertically downward direction 2 Conduct in - situ stress balance calculation. Apply fault dislocation displacement in the static analysis step. Apply a horizontal displacement of 0.5 m to the side boundary of the fault active disk, and simultaneously apply a horizontal displacement of 0.5 m and a vertical displacement of 0.5 m to the bottom. The initial step size is taken as 0.05. The material properties are shown in Table 1. The mesh division uses 6 - node modified quadratic plane strain triangular elements (CPE6M). Set *Restart, write, frequency = 5 in the ABAQUS input file, and submit the ABAQUS solver to perform the first small - deformation analysis calculation. The results are as Figure 7A as shown

[0078] Use the Python script program in Appendix 1 to access the ODB document generated by the ABAQUS solver, export the deformed node coordinates and write them into a file. Use the bisection method for iterative optimization, and use the Alpha Shape algorithm in Appendix 2 to calculate the boundary recognition results corresponding to α = 1, 2, 1.5, 1.75, 1.25, 1.125, 1.375 in sequence. Take α = 1.375 as the boundary recognition result when the vertical dislocation is 50 cm (as Figure 5 shown), and automatically output the boundary point coordinate list for secondary modeling in ABAQUS.

[0079] Based on the boundary recognition result, re - conduct mesh division in the new calculation domain through the built - in mesh generation method in ABAQUS. When performing mesh division, use global seeding, the global seeding size is 4, the mesh type is 6 - node modified quadratic plane strain triangular elements, and the mesh re - division results are as Figure 6 shown. Apply the fault dislocation displacement again. Apply a horizontal displacement of 0.5 m to the side boundary of the fault active disk, and simultaneously apply a horizontal displacement of 0.5 m and a vertical displacement of 0.5 m to the bottom. The initial step size is taken as 0.05.

[0080] Check that the maximum number of increment steps in the previous analysis calculation is 5. Specify the variable mapping command by modifying the input file *map solution, step = 1, inc = 5, and specify the new job and the old job in the ABAQUS command box to perform the second analysis calculation. The results are as Figure 10A as shown

[0081] Again, use the Python script program in Appendix 1 to access the ODB document generated by the ABAQUS solver, export the deformed node coordinates and write them to a file. Use the bisection method for iterative optimization, and use the Alpha Shape algorithm in Appendix 2 to calculate the boundary recognition results corresponding to α = 2, 3, 2.5, 2.75, 2.25, 2.125, and 2.375 in sequence. Take α = 2.125 as the boundary recognition result when the vertical dislocation is 100 cm (as Figure 8 shown), and automatically output the boundary point coordinate list for secondary modeling in ABAQUS.

[0082] Based on the boundary recognition results, re - mesh in the new calculation domain using the built - in mesh generation method in ABAQUS. When meshing, use global seeding, with a global seeding size of 3, and the mesh type is 6 - node modified quadratic plane strain triangular elements (as Figure 9 shown). Apply the fault dislocation displacement again. Apply a horizontal displacement of 0.5 m to the side boundary of the fault - active disk, and simultaneously apply a horizontal displacement of 0.5 m and a vertical displacement of 0.5 m to the bottom, with an initial step size of 0.05.

[0083] Check that the maximum number of increment steps in the previous analysis calculation is 23. Complete the variable mapping process by modifying the input file *map solution, step = 1, inc = 23. Specify the new job and the old job in the ABAQUS command box for the third analysis calculation, and the results are as Figure 11 shown. The mesh diagram after the last small - deformation calculation is as Figure 12 shown, and no mesh distortion problems occur.

[0084] Python script in Appendix 1 for exporting the node coordinates after model deformation:

[0085]

[0086]

[0087] Alpha Shape program in Appendix 2:

[0088]

[0089]

[0090]

[0091]

[0092] Map Solution script in Appendix 3:

[0093] *map solution, step = 1, inc = 10

[0094] **STEP: Step - 1

[0095] **

[0096] *Step, name = Step - 1, nlgeom = YES

[0097] *Static

[0098] 1., 1., 1e - 05, 1.

[0099] **

[0100] **BOUNDARY CONDITIONS

[0101] **

[0102] **Name: BC - 2 Type: Displacement / Rotation

[0103] *Boundary

[0104] _PickedSet131, 1, 1, 0.5

[0105] **Name: BC - 4 Type: Displacement / Rotation

[0106] *Boundary

[0107] _PickedSet133, 1, 1, 0.5

[0108] _PickedSet133, 2, 2, - 0.5

[0109] **

[0110] **OUTPUT REQUESTS

[0111] **

[0112] *Restart, write, frequency = 1

[0113] **

[0114] **FIELD OUTPUT: F - Output - 1

[0115] **

[0116] *Output, field

[0117] *Node Output

[0118] U,

[0119] *Element Output, directions=YES

[0120] PEEQ, S

[0121] **

[0122] **HISTORY OUTPUT: H-Output-1

[0123] **

[0124] *Output, history, variable=PRESELECT

[0125] *End Step

[0126] Although embodiments of the present invention have been shown and described, it will be understood by those of ordinary skill in the art that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the invention, and the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A simulation method for large deformation response of cross-fault sites based on the particle finite element method, characterized in that: It includes the following steps: S1. Use Python to establish a two-dimensional cross-fault site model, submit it to the ABAQUS solver, simulate the fault dislocation process by applying fault dislocation displacement on the two-dimensional model, and perform the first small deformation analysis and calculation; S2. Establish a mesh distortion warning mechanism, access the ODB file generated by the ABAQUS solver, and export the deformed node coordinates; S3. Use the Alpha Shape algorithm to identify the boundary of the deformed model, and output the boundary node coordinates to a list file; S4. Based on the obtained boundary node coordinate list file, create a deformed two-dimensional model; S5. Based on the boundary recognition result, re-mesh the new calculation domain; S6. Extract the mesh information, map the calculation information from the old mesh to the new mesh through variable mapping technology, continue to apply fault dislocation displacement on the deformed two-dimensional model, and perform the small deformation analysis and calculation again; S7. Repeat steps S2 to S6 until the preset fault dislocation amount is reached.

2. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 1, characterized in that: In step S1, by applying fault dislocation displacement on the two-dimensional model to simulate the fault dislocation process and submitting it to the ABAQUS solver for the first small deformation analysis and calculation, it includes: Before implementing the fault dislocation analysis, first complete the in-situ stress balance calculation: apply a gravity load in the vertically downward direction, create a static analysis step for stress balance calculation, and define the initial stress field that conforms to the in-situ stress distribution of the site; after completing the in-situ stress balance calculation, create a new static analysis step to apply fault dislocation displacement, with the initial step size taking 0.1% of the model height, using a certain amplitude of equivalent plastic strain penetration as the standard for soil rupture, dynamically adjusting the displacement increment step size according to the plastic zone expansion rate, and using 6-node modified quadratic plane strain triangular elements for the mesh element type. Complete the settings in the ABAQUS input file and submit it to the ABAQUS solver for the first small deformation analysis and calculation; Among them, the method of applying fault dislocation displacement is: apply horizontal constraints and no vertical constraints on the boundary of the passive side of the fault, constrain the horizontal and vertical degrees of freedom at the bottom, apply horizontal displacement and no vertical constraints on the boundary of the active side of the fault, and apply horizontal and vertical displacements simultaneously at the bottom. Simulate different bedrock fracture dips by controlling the ratio of the vertical displacement amount to the horizontal displacement amount.

3. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 1, characterized in that: In step S2, establishing a mesh distortion warning mechanism includes: adopting a dual-trigger mechanism for mesh distortion judgment. When the set increment step threshold is calculated or it is judged that the mesh is about to be distorted, use a script program written in the Python language to access the ODB document generated by the ABAQUS solver.

4. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 3, characterized in that: In step S2, access the ODB file generated by the ABAQUS solver and export the deformed nodal coordinates, including: S21, open the ODB file; S22, obtain the instances in the model; S23, obtain the analysis step and the last frame; S24, extract the original nodal coordinates and deformation data; S25, close the ODB file; S26, calculate the deformed nodal coordinates, where the deformed nodal coordinates = the original nodal coordinates + the deformation data; S27, close the ODB file and write the deformed nodal coordinates to a list file.

5. The large deformation response simulation method of cross-fault site based on the particle finite element method according to claim 4, characterized in that: In step S3, use the Alpha Shape algorithm to perform boundary recognition based on the exported deformed nodal coordinates, visualize the boundary points and output the boundary point coordinate list, including: S31, read the deformed nodal coordinate file exported in step S2 and change the format to the form of ID, X, Y; S32, identify the boundary and output the boundary points, including reading the point set, selecting the α value, extracting the boundary points, visualizing the results, plotting the original points and plotting the boundary; S33, output the boundary point coordinate list file and configure the file path.

6. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 5, characterized in that: In step S3, the process of selecting the α value is as follows: First, preliminarily determine the reasonable range of the α value, and then use the bisection method to iterate and optimize within this range. Calculate the boundary recognition results corresponding to the intermediate values in turn. Through the comparison error analysis with the boundary before deformation, gradually narrow the range of the α value until the accuracy requirement is met, and control the boundary recognition error within 5%.

7. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 1, characterized in that: In step S4, based on the obtained deformed boundary point coordinate list file, use the Python language to write a script program to create the deformed two-dimensional model, including: S41, read the nodal coordinate file and change the format to the form of X, Y; S42, create a sketch, draw points and connect lines in the sketch to create a new two-dimensional model.

8. The method for simulating large deformation response of cross-fault site based on particle finite element method according to claim 7, characterized in that: In step S5, based on the boundary recognition results, re-perform mesh generation in the new computational domain, including: Based on the boundary recognition results, re-perform mesh generation through the built-in mesh generation method in ABAQUS. Use global seeding when re-performing mesh generation, re-set the global seeding size, and the element type is a 6-node modified quadratic plane strain triangular element.

9. The method for simulating large deformation response of cross-fault site based on particle finite element method according to claim 8, characterized in that: In step S6, extract the mesh information, map the calculation information from the old mesh to the new mesh through the variable mapping technology, and continue to apply the fault dislocation displacement on the deformed two-dimensional model and perform small deformation analysis and calculation again, including: view the maximum number of increment steps of the previous analysis and calculation, call Map Solution, and modify the input file to specify the variable mapping command. Specify the new job and the old job again in the ABAQUS command box to perform small deformation analysis and calculation again.

10. The cross-fault site large deformation response simulation method based on the particle finite element method according to claim 1, characterized in that: The two-dimensional model includes the overlying soil layer thickness, fracture position, bedrock fracture dip angle, and fault dislocation amount. And in order to reduce the influence of the boundary on the deformation of the soil body center, the width of the two-dimensional model is set to 3-4 times the height.

Images