Ground subsidence simulation risk method and system based on finite elements

Abstract

The application is specifically a ground subsidence simulation risk method and system based on finite elements, and relates to the technical field of geological disaster simulation and numerical calculation, comprising: applying control equations, material parameters and boundary conditions on the current grid, calling a finite element solver for numerical solution, and obtaining displacement field values and stress field values of each grid element. In the application, Kelly error estimators are used to calculate local error indicator values of each grid element, and based on this, the grid elements are automatically subjected to refinement, coarsening or maintaining operations, thereby solving the technical problems that the traditional uniform grid method is difficult to balance between precision and calculation efficiency and is low in grid efficiency depending on manual experience adjustment.

Description

Technical Field

[0001] This invention relates to the field of geological disaster simulation and numerical calculation technology, and in particular to a method and system for simulating ground subsidence risk based on the finite element method. Background Technology

[0002] Ground subsidence is a common geological hazard that has a significant impact on urban infrastructure and building safety. The finite element method (FEM) is the primary numerical method for ground subsidence simulation and analysis, and one of its core steps is meshing the geological model.

[0003] In existing technologies, finite element ground settlement simulation typically uses a uniform grid to discretize the computational domain and relies on the human experience of engineers to adjust the grid density in local areas.

[0004] Therefore, the following defects exist that need to be improved: When using a coarser uniform grid, the solution accuracy in stress concentration areas or areas of abrupt geological changes is insufficient, causing the simulation results to deviate from the actual deformation characteristics. Using a finer, uniform mesh increases the overall computational load, resulting in a waste of computational resources. On the other hand, relying on manual experience to adjust the local mesh is time-consuming and it is difficult to strike a balance between solution accuracy and computational efficiency, which limits both the reliability and efficiency of simulation analysis. Summary of the Invention

[0005] The purpose of this invention is to propose a finite element method and system for simulating the risks of ground subsidence in order to solve the above-mentioned problems.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: Finite element method-based risk simulation methods for ground settlement include: The computational domain is discretized based on the geometric information of the geological model to generate an initial mesh composed of multiple mesh elements. Governing equations, material parameters, and boundary conditions are applied to the current mesh, and a finite element solver is invoked to perform numerical solutions, obtaining the displacement and stress field values ​​for each mesh element. Based on the solution results for each mesh element, an error estimator is used to calculate the local error index value corresponding to each mesh element. The local error index value of each mesh element is compared with a refinement threshold and a coarsening threshold, respectively. Mesh elements with local error index values ​​greater than the refinement threshold are marked as refined, those with local error index values ​​less than the coarsening threshold are marked as coarsened, and those with local error index values ​​between the coarsening and refinement thresholds are marked as... To maintain the state, an adaptive adjustment operation is performed on the current mesh based on the marking results. This involves refining the mesh elements in the refined state and coarsening the mesh elements in the coarsened state, generating an adaptive mesh. The solution results from the old mesh are mapped to the adaptive mesh, and the finite element method is re-solved based on the adaptive mesh to obtain updated solution results. A global error index is calculated based on the updated solution results, and the global error index is compared with a preset global error tolerance. If the global error index is greater than the global error tolerance, the adaptive mesh is used as the current mesh, and the process returns to the step of calculating the local error index value to continue iterating. If the global error index is not greater than the global error tolerance, the current solution result is output as the final solution result.

[0007] Further, the step of using an error estimator to calculate the local error index value corresponding to each grid cell includes: using a Kelly error estimator, taking the displacement field values ​​of each grid cell as input, obtaining the feature size of the grid cell and the set of boundary surfaces of the grid cell for each grid cell; calculating the jump amount of the normal gradient of the displacement field numerical solution between adjacent cells on both sides of the boundary surface for each boundary surface of the grid cell, wherein the jump amount is the difference of the normal gradient values ​​of adjacent cells on both sides of the boundary surface; performing area integration on the boundary surface on the square of the jump amount, summing the integration results of all boundary surfaces, multiplying by the ratio of the feature size of the grid cell to a preset constant, and taking the square root of the resulting product to obtain the local error index value corresponding to the grid cell.

[0008] Furthermore, the refinement threshold and the coarsening threshold are determined as follows: the local error index values ​​of all grid cells are sorted in descending order of value; the local error index value corresponding to the first proportional parameter position after sorting is used as the refinement threshold; and the local error index value corresponding to the second proportional parameter position after sorting is used as the coarsening threshold. The sum of the first proportional parameter and the second proportional parameter is less than one, and the coarsening threshold is less than the refinement threshold.

[0009] Furthermore, before performing adaptive adjustment operations on the current mesh based on the marking results, the method further includes: performing smooth constraint processing on the marking results, traversing all mesh cells marked as refined, checking the marking status of their adjacent mesh cells, and when a mesh cell marked as refined is directly adjacent to a mesh cell marked as coarse, modifying the marking of the mesh cell marked as coarse to a maintain state, so that the size transition between adjacent mesh cells meets the preset grade constraint conditions.

[0010] Furthermore, the refinement and coarsening operations include: in the two-dimensional case, the refinement operation is to split a triangular mesh cell into four sub-triangular cells along the lines connecting the midpoints of each side; in the three-dimensional case, the refinement operation is to split a tetrahedral mesh cell into eight sub-tetrahedral cells; the coarsening operation is the reverse process of the refinement operation, which re-merges multiple sub-cells generated by refining the same parent cell back into the same parent cell.

[0011] Further, mapping the solution results on the old grid to the adaptive grid includes: performing the mapping using a shape function-based interpolation method; for the target node in the adaptive grid, locating the grid cell where the target node is located in the old grid, and obtaining the local coordinates of the target node in the located grid cell; calculating the shape function values ​​corresponding to each node in the located grid cell based on the local coordinates; multiplying the known field quantity values ​​at each node by the corresponding shape function values ​​and summing them to obtain the mapped value at the target node; wherein the shape function satisfies the interpolation condition of taking a value of one at the corresponding node and a value of zero at the other nodes.

[0012] Furthermore, locating the grid cell where the target node is located in the old grid includes: traversing all grid cells in the old grid, converting the global coordinates of the target node into the local coordinates of each grid cell, determining whether the local coordinates satisfy the reference domain range constraint of the cell type, and if so, confirming that the target node is located in the currently traversed grid cell, and using the local coordinates as the input for subsequent shape function interpolation.

[0013] Furthermore, before resolving the finite element problem after mapping the solution results from the old mesh to the adaptive mesh, the method further includes: performing an equilibrium check on the mapped initial field values; calculating the residual force vector of the mapped stress field on the adaptive mesh; when the norm of the residual force vector is greater than a preset equilibrium tolerance, performing an equilibrium iterative correction on the mapped initial field values ​​until the norm of the residual force vector is not greater than the equilibrium tolerance; obtaining the equilibrium-corrected initial field values; and using the equilibrium-corrected initial field values ​​as the initial conditions for resolving the problem.

[0014] Further, the global error index is calculated as follows: the local error index values ​​of all grid cells in the adaptive mesh are squared and summed, and the square root of the sum is taken. The convergence condition also includes an iteration limit. When the number of iterations reaches the preset maximum number of iterations, the iteration loop is terminated and the current solution result is output as the final solution result. After the global error index meets the global error tolerance, the calculation further includes: extracting the local error index values ​​of each grid cell in the preset region of interest, squaring the local error index values ​​of all grid cells in the region of interest, summing them, and taking the square root of the sum to obtain the local error index; comparing the local error index with the preset local error tolerance, and when the local error index is greater than the local error tolerance, re-performing the local error index value calculation, adaptive labeling, adaptive adjustment, and mapping solution operations for the grid cells in the region of interest until the local error index is not greater than the local error tolerance; wherein, the local error tolerance is not greater than the global error tolerance.

[0015] This invention discloses a multi-scale adaptive mesh generation system, comprising: an initial mesh generation module, used to discretize the computational domain based on the geometric information of a geological model to generate an initial mesh composed of multiple mesh elements; a finite element solution module, used to apply governing equations, material parameters, and boundary conditions to the current mesh, call a finite element solver to perform numerical solutions, and obtain the displacement field and stress field values ​​of each mesh element; an error estimation module, used to calculate the local error index value corresponding to each mesh element based on the solution results of each mesh element using an error estimator; an adaptive marking module, used to compare the local error index value of each mesh element with a refinement threshold and a coarsening threshold respectively, and mark each mesh element as a refinement state, a coarsening state, or a maintainable state according to the comparison results; and a mesh adjustment module, used to adjust the mesh according to the marked values. The algorithm performs an adaptive adjustment operation on the current mesh, refines the mesh elements in the refined state, and coarsens the mesh elements in the coarsened state to generate an adaptive mesh. A solution result mapping module maps the solution results from the old mesh to the adaptive mesh and calls the finite element solution module to re-solve the finite element problem based on the adaptive mesh to obtain updated solution results. A convergence judgment module calculates a global error index based on the updated solution results, compares the global error index with a preset global error tolerance, and if the global error index is greater than the global error tolerance, returns the adaptive mesh as the current mesh to the error estimation module for further iteration; if the global error index is not greater than the global error tolerance, outputs the current solution result as the final solution result.

[0016] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are: 1. This invention solves the technical problems of traditional uniform grid methods, such as the difficulty in balancing accuracy and computational efficiency, and the low efficiency of relying on manual experience to adjust the grid, by using the Kelly error estimator to calculate the local error index value of each grid cell and automatically performing refinement, coarsening or maintenance operations on the grid cells accordingly.

[0017] 2. This invention reduces the number of redundant grid cells in low-error regions by allocating computing resources to regions with larger errors, thereby reducing the overall computational scale. The grid adjustment process is automatically executed based on error estimation results, without the need for human experience, thus improving the automation level of grid adjustment. By mapping the solution results of the old grid to the adaptive grid as the initial condition, the computational overhead of resolving after each grid adjustment is reduced. Attached Figure Description

[0018] Further details, features, and advantages of this application are disclosed in the following description of exemplary embodiments in conjunction with the accompanying drawings, in which: Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0019] Several embodiments of this application will now be described in more detail with reference to the accompanying drawings to enable those skilled in the art to implement this application. This application may be embodied in many different forms and for various purposes and should not be limited to the embodiments set forth herein. These embodiments are provided to make this application thorough and complete, and to fully convey the scope of this application to those skilled in the art. The embodiments described do not limit this application.

[0020] Unless otherwise defined, all terms used herein (including technical and scientific terms) shall have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains. It will be further understood that terms such as those defined in commonly used dictionaries shall be interpreted as having a meaning consistent with their meaning in the relevant field and / or the context of this specification, and shall not be interpreted in an idealized or overly formal sense unless expressly defined herein.

[0021] Example 1 Its specific implementation method is combined with the appendix Figure 1 Please provide a detailed explanation.

[0022] In this embodiment, it includes: Step 1: Generate the initial mesh based on the geometric information of the geological model. Geometric information of the geological model involved in the ground subsidence simulation is obtained, including the boundary contour of the computational domain, the location of the stratigraphic interface, and fault distribution data; based on the geometric information, the computational domain is discretized to generate an initial grid covering the entire computational domain, the initial grid consisting of multiple grid cells.

[0023] The initial mesh is a relatively coarse, uniform mesh, meaning that the size of each mesh cell is basically the same. The cell size is determined based on the overall range of the computational domain and the preset number of initial cells. The type of mesh cell can be triangular or tetrahedral, or quadrilateral or hexahedral, depending on the dimensions and geometric complexity of the geological model.

[0024] Step 2: Perform finite element analysis on the current mesh and obtain the solution results for each mesh element. The control equations, material parameters, and boundary conditions for the ground settlement simulation are applied to the current mesh. The finite element solver is then called to perform numerical solutions on the current mesh, and the solution results for each mesh element are obtained. The solution results include the displacement field values ​​and stress field values ​​for each mesh element.

[0025] Step 3: Based on the solution results of each grid cell, calculate the local error index value of each grid cell. The Kelly error estimator is used, taking the displacement and stress field values ​​of each grid cell as input, to calculate the local error index value corresponding to each grid cell. The local error index value is used to quantitatively characterize the solution accuracy of the corresponding mesh element. For mesh elements... Its local error index value The calculation formula is: ; in, For grid cells Feature size, For grid cells The set of boundary surfaces, For the boundary surface, For boundary surface The outer normal vector, For the numerical solution obtained on the current grid, For numerical solutions gradient, The normal gradient of the numerical solution at the interface of adjacent elements The jump amount is specifically defined as the difference in normal gradient values ​​between adjacent cells on both sides of the interface, i.e. ,in and These are the shared interfaces. Two adjacent grid cells, For the interface The area of ​​the infinitesimal element is defined as follows: The jump amount reflects the degree of discontinuity in the solution results between adjacent elements; the larger the jump amount, the greater the discretization error in that region.

[0026] Furthermore, in ground settlement simulation, the numerical solution... The displacement field value is taken as the input quantity for error estimation, that is, the normal gradient jump of the displacement field at the interface between adjacent elements is used as the input quantity for error estimation, thereby making the local error index value... It can reflect the level of discrete error of the displacement field at that grid cell.

[0027] Step 4: Based on the local error index value of each grid cell, adaptively label each grid cell and generate the labeling results. The local error index value of each grid cell is compared with the refinement threshold. and coarsening threshold A comparison is performed, and based on the comparison results, a corresponding operation tag is assigned to each grid cell, generating a tagging result. In the tagging result, each grid cell is marked as one of the following three states: when the local error index value... Greater than the refinement threshold When the local error index value is [value missing], mark the grid cell as refined; Less than the coarsening threshold When the local error index value is... Between the coarsening threshold With refined threshold When in between, mark the grid cell as in a hold state.

[0028] The refinement threshold and coarsening threshold The determination method is as follows: First, sort the local error index values ​​of all grid cells in descending order of value, and then select the cells that are ranked first... The local error index value corresponding to the proportional position is used as the refinement threshold. After sorting, the position will be at the end. The local error index value corresponding to the proportional position is used as the coarsening threshold. ,in and The preset scaling parameters, and That is, the one with the largest local error index value. The scaled-down mesh cells will be marked as refined, with the smallest local error index value being the last one. The proportional mesh cells will be marked as coarsened.

[0029] Furthermore, the refined threshold With coarsening threshold Between The constraint relationships are used to ensure the local error index value of any grid cell. It can be uniquely classified into one of the three states: refinement, preservation, or coarsening, thus avoiding ambiguity in the marked state.

[0030] In this embodiment of the application, in order to avoid numerical instability during the solution process due to excessive size differences between adjacent mesh cells, after adaptively marking each mesh cell, the following steps are also included: smoothing constraint processing of the marking results, that is, traversing all mesh cells marked as refined, checking the marking status of their adjacent mesh cells, and when a mesh cell marked as refined is directly adjacent to a mesh cell marked as coarse, modifying the marking of the mesh cell marked as coarse to the maintain state, thereby ensuring that the size transition between adjacent mesh cells meets the preset grade constraint conditions.

[0031] Step 5: Based on the marking results, perform an adaptive adjustment operation on the current mesh to generate an adaptive mesh. Based on the marking results obtained in step 4, corresponding adjustment operations are performed on each mesh cell in the current mesh to generate an adaptive mesh. Specifically, a refinement operation is performed on mesh cells marked as refined, splitting the mesh cell into multiple smaller sub-cells; a coarsening operation is performed on mesh cells marked as coarsened, merging multiple adjacent small-sized mesh cells into a larger-sized mesh cell; and no adjustment operation is performed on mesh cells marked as maintained.

[0032] One specific form of the refinement operation is as follows: in two dimensions, a triangular mesh cell is split into four sub-triangular cells along the lines connecting the midpoints of its sides; in three dimensions, a tetrahedral mesh cell is split into eight sub-tetrahedral cells. The coarsening operation is the reverse process of the refinement operation, that is, multiple sub-cells generated by refining from the same parent cell are merged back into the same parent cell.

[0033] Step 6: Map the solution results from the old mesh to the adaptive mesh, and re-solve the finite element method based on the adaptive mesh to obtain the updated solution results. Using the displacement and stress field values ​​from the old mesh as input, a shape function-based interpolation method is employed to map the solution results from the nodes and elements of the old mesh to the corresponding nodes and elements of the adaptive mesh, thereby obtaining the initial field values ​​on the adaptive mesh. The specific process of the shape function-based interpolation method is as follows: For the target node in the adaptive mesh, the mesh element containing the target node is first located in the old mesh. Then, the shape function of this mesh element is used to perform weighted interpolation on the node values ​​in the old mesh to obtain the mapped values ​​at the target node. The specific method of the shape function-weighted interpolation is as follows: Let the local coordinates of the target node in the mesh element of the old mesh be... The grid cell has a total of The node, the The shape function value corresponding to each node is , No. The known field quantities at each node are: Then the mapping value at the target node Calculated by the following formula: ; in, For standard finite element shape functions, satisfying the condition in the first... The value at each node is [value]. The value at the remaining nodes is... Interpolation conditions, The total number of nodes in the located grid cell. The node's sequence number. Position the first grid cell in the old grid The known field quantities at each node are summed, with the summation covering all nodes of that grid cell. from Received .

[0034] Furthermore, the process of locating the target node within the old grid cell is as follows: traversing all grid cells in the old grid, and converting the global coordinates of the target node into the local coordinates of each grid cell. It determines whether the local coordinates satisfy the reference domain range constraint of the cell type. If they do, it confirms that the target node is located within the currently traversed mesh cell and uses the local coordinates. As input for subsequent shape function interpolation.

[0035] The initial field values ​​are used as initial conditions on the adaptive mesh. The finite element solver is then called to perform a new numerical solution based on the adaptive mesh to obtain the updated solution results.

[0036] In this embodiment of the application, in order to reduce the impact of mapping errors introduced by mesh adjustment on subsequent solutions, after mapping the solution results on the old mesh to the adaptive mesh and before resolving the finite element problem, the following steps are also included: performing balance verification on the mapped initial field values, that is, calculating the residual force vector of the mapped stress field on the adaptive mesh. When the norm of the residual force vector is greater than the preset balance tolerance, the mapped initial field values ​​are iteratively corrected for balance until the norm of the residual force vector is less than or equal to the balance tolerance, thereby obtaining the initial field values ​​after balance correction, and using the initial field values ​​after balance correction as the initial conditions for resolving the problem.

[0037] Step 7: Determine if the convergence condition is met. If not, return to Step 3 to continue iterating. If the condition is met, output the final solution. Based on the updated solution results, calculate the global error index for the current iteration. The global error index The calculation formula is: ; in, For each grid cell, calculate and obtain the local error index value in the current iteration round. For each grid cell, the summation range is all grid cells in the adaptive grid. The global error index is used. Compared with the preset global error tolerance Comparison: When When the convergence condition is met, the solution result on the current adaptive mesh is output as the final solution result; when... If the convergence condition is not met, the current adaptive mesh is used as the new current mesh, and the process returns to step 3. Based on the updated solution results, the local error index values ​​of each mesh element are recalculated, and subsequent marking, adjustment, mapping, and solving operations are performed until the convergence condition is met.

[0038] Furthermore, the global error tolerance Dimensions and local error index values The dimensions are consistent with those of the numerical solution. The displacement fields they represent have the same dimensions, thus making the global error index... With global error tolerance The comparison between them has a clear physical meaning, that is, when When this occurs, it indicates that the overall discrete error level of all grid cells on the adaptive grid has been reduced to an acceptable range.

[0039] The convergence condition may also include an iteration count limit, i.e., when the number of iterations reaches a preset maximum number of iterations. Even if the global error index has not yet met the global error tolerance requirement, the iteration loop is terminated and the current solution result is output as the final solution result.

[0040] In this embodiment of the application, based on the convergence judgment, in order to obtain higher solution accuracy in specific high-risk areas, the following steps are also included: [The text abruptly ends here, likely due to an incomplete sentence or missing information.] Global error tolerance has been met. Next, the local error index values ​​of each grid cell within the preset area of ​​interest are extracted, and the local error index of the area of ​​interest is calculated. : ; in, For the preset area of ​​interest, To focus on the grid cells within the region, To determine the local error index values ​​for each grid cell within the region of interest, the summation range covers all grid cells within that region. The local error index... Compared with the preset local error tolerance When comparing, At this point, only the mesh cells within the region of interest are re-executed through steps 3 to 6, resulting in targeted mesh refinement and re-solution of the region of interest, until... The area of ​​interest can be pre-defined according to the analysis requirements of ground settlement simulation, such as the stratum directly below the building foundation or the area with drastic changes in groundwater level.

[0041] Furthermore, the local error tolerance With global error tolerance Between The constraints are used to ensure that the solution accuracy within the region of interest is not lower than the global accuracy requirement, thereby enabling targeted and encrypted solutions for high-risk regions.

[0042] The multi-scale adaptive mesh generation method provided in this embodiment obtains local error index values ​​by calculating the error of the solution results of each mesh element during the finite element solution process. Based on these local error index values, the mesh elements are adaptively marked and adjusted. Regions with larger errors are automatically refined into denser meshes, while regions with smaller errors are automatically coarsened into sparser meshes. Therefore, computational resources are allocated to regions with insufficient solution accuracy, rather than being uniformly distributed across the entire computational domain. This reduces the number of mesh elements in unnecessary regions and lowers the overall computational scale while meeting preset accuracy requirements. Furthermore, since the mesh adjustment process is automatically executed based on the error estimation results and iteratively optimizes the mesh distribution until convergence conditions are met, no human experience is required to intervene in the mesh density adjustment decision process. Therefore, this overcomes the inefficiency and difficulty in achieving a stable balance between accuracy and efficiency inherent in traditional methods that rely on human experience to adjust the mesh. Furthermore, by mapping the solution results on the old mesh to the adaptive mesh after mesh adjustment, the solution information from the previous round is used to provide initial conditions for the new mesh, avoiding the need to start from scratch after each mesh adjustment, which helps to reduce the computational cost of finite element solution in each iteration round.

[0043] Application scenarios During the construction of an urban rail transit project, finite element simulation analysis of ground settlement in the surrounding strata of the tunnel is required. The computational domain covers a depth of 60 meters below the surface and extends horizontally for 120 meters, involving three geological layers: an upper fill layer (0 to 8 meters), a middle silty clay layer (8 to 35 meters), and a lower gravel layer (35 to 60 meters). A near-vertical fault exists at a depth of 22 meters. The tunnel has an outer diameter of 6.2 meters, and its centerline is located 18 meters below the surface, belonging to a typical area of ​​abrupt geological changes. The initial conditions are: the groundwater level is 5 meters below the surface, the release of ground pressure caused by tunnel excavation is applied as the load, and the global error tolerance is... Meters, Area of ​​Focus The geological zone is defined as a stratum 20 meters wide and 10 to 26 meters deep above the tunnel centerline, with a local error tolerance. rice.

[0044] Application Examples Step 1: Generate the initial mesh based on the geometric information of the geological model. The geometric information of the tunnel engineering geological model was obtained, including the boundary outline of a 120m × 60m rectangular computational domain, the locations of the three strata interfaces (horizontal interfaces at depths of 8m and 35m), and the fault line at a depth of 22m. Based on the above geometric information, the computational domain was discretized to generate an initial uniform triangular mesh. The total area of ​​the computational domain was 7200 square meters, with a preset initial number of 480 elements, and each element had the following characteristic dimensions. It is approximately 5.5 meters.

[0045] Table 1 Basic parameters of the initial mesh Step 2: Perform finite element analysis on the current mesh and obtain the solution results for each mesh element. The tunnel excavation ground pressure release load, the parameters of three geological layers (elastic modulus of 12 MPa for the fill layer, 28 MPa for the silty clay layer, and 65 MPa for the gravel layer), and the bottom fixed boundary condition were applied to the initial mesh. The finite element solver was then used to complete the first round of calculations. The results showed that the maximum settlement occurred at the surface node directly above the tunnel, and stress concentration areas were concentrated near the fault and in the tunnel arch region. The results of the first round of calculations for six representative mesh elements are listed below.

[0046] Table 2. Solution results for typical grid elements in the first round. Step 3: Based on the solution results of each grid cell, calculate the local error index value of each grid cell. The Kelly error estimator is used, with the displacement field values ​​of each grid cell as the numerical solution. The input is used to calculate the local error index value of each unit. For cell 047, its feature size is... The normal gradient jump between adjacent elements in the vault region is relatively large, and the cumulative value of the jump integral along the four boundary surfaces is... m² / m, substituting into the formula: ; For element 412 at the edge of the computational domain, the displacement change between adjacent elements is gradual, and the cumulative value of the jump integral is only... m² / m, substituting into the formula: ; The local error index values ​​for each typical unit are summarized below.

[0047] Table 3. Local error index values ​​of typical units in the first round. The global error index for round 1 is: ; This value is much greater than The meter needs further iteration.

[0048] Step 4: Based on the local error index value of each grid cell, adaptively label each grid cell and generate the labeling results. Sort the local error index values ​​of the 480 units from largest to smallest, and set the scaling parameter. , After sorting, the index value corresponding to the top 25% (120th position) is the refinement threshold. The value of the index corresponding to the meter, located in the bottom 30% (position 337), is the coarsening threshold. Rice, to satisfy .

[0049] Judgment of the marking of each typical unit: Local error index value of unit 051 Meters greater than the refinement threshold Marked as refined; local error index value of unit 412 Meters less than the coarsening threshold Marked as coarsened; Local error index value of unit 089 Meters between the two thresholds are marked as held.

[0050] After smoothing constraint processing, it was found that some elements marked as coarsened were directly adjacent to refined elements (such as the transition zone on the side of the tunnel). Their markings were changed to be preserved. The final marking results are as follows.

[0051] Table 4. Statistics of Adaptive Tagging Results in Round 1 Step 5: Based on the marking results, perform an adaptive adjustment operation on the current mesh to generate an adaptive mesh. Perform a quadripartite splitting operation along the midpoints of each side on 120 refined elements, splitting each element into 4 sub-triangles. The number of elements was reduced by approximately 95 (126 coarsened elements were merged, recombining four sub-triangles generated from the same parent element into a single parent element, thus reducing the number of elements by about 95; some coarsened elements were in the initial layer and had no parent elements to merge with); 234 elements were left unchanged. After adaptive mesh generation, the total number of elements increased from 480 to approximately 745. The mesh feature size in the fault region and tunnel arch region was refined from 5.5 meters to approximately 1.4 meters, and the size of the computational domain edge region was coarsened from 5.5 meters to approximately 9.8 meters.

[0052] Step 6: Map the solution results from the old mesh to the adaptive mesh, and then perform the finite element solution again based on the adaptive mesh. Using the displacement and stress field values ​​on the old mesh as input, shape function-based interpolation mapping is performed on each target node in the adaptive mesh (745 elements).

[0053] Taking the newly added target node NODE-T23 in the adaptive mesh as an example, this node is located in the refined region of the tunnel arch, and its global coordinates are transformed to the local coordinates of the old mesh cell 047. The old mesh element 047 is a triangular element with 3 nodes. The vertical displacement values ​​at the three nodes are as follows: mm mm mm, and the corresponding shape function values ​​are respectively , , The mapped value is: ; After completing the mapping of all nodes, the norm of the residual force vector is calculated for the mapped stress field. kN, less than the preset balance tolerance kN, without the need for balancing iterative correction, directly uses the mapping result as the initial condition, and calls the finite element solver to complete the second round of solution based on the adaptive mesh.

[0054] Step 7: Determine if the convergence condition is met, and iterate until the final solution is output. After multiple iterations, the changes in the global error index and the local error index of the area of ​​interest in each iteration are as follows.

[0055] Table 5 Error Indicators and Convergence Judgments for Each Iteration Round Round 3 Global Error Indicators Meters less than or equal to the global error tolerance The value of the sample is 1 meter, satisfying the global convergence condition. Further extraction of the region of interest is needed. The local error index of each unit within the stratum (20 meters wide and 10 to 26 meters deep above the tunnel) was calculated. Meters greater than the local error tolerance The meter did not meet the local accuracy requirements. Steps 3 to 6 were performed individually on the units within the area of ​​interest. After another round of targeted refinement, the local error index was reduced to Meters, less than or equal to the local error tolerance The results were obtained in meters, meeting the local error tolerance requirements. The final output was the ground settlement solution on the 4th round adaptive grid (1247 elements). The maximum ground settlement above the tunnel arch was 23.47 mm, and the maximum principal stress at the fault was 476.3 kPa.

[0056] The foregoing has only described certain exemplary embodiments of the present invention by way of illustration. Undoubtedly, those skilled in the art can modify the described embodiments in various ways without departing from the spirit and scope of the present invention. Therefore, the foregoing drawings and descriptions are illustrative in nature and should not be construed as limiting the scope of protection of the claims of the present invention.

[0057] It should be noted that, in this document, the use of relational terms such as "first" and "second" is merely for distinguishing one entity or operation from another, and does not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes the element.

[0058] It should be understood that in the various embodiments of this application, the sequence number of each process does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0059] In addition, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

Claims

1. A finite element method for simulating the risk of ground settlement, characterized in that, include: The computational domain is discretized based on the geometric information of the geological model to generate an initial grid composed of multiple grid cells; Apply governing equations, material parameters, and boundary conditions to the current mesh, call the finite element solver to perform numerical solutions, and obtain the displacement field and stress field values ​​of each mesh element. Based on the solution results of each grid cell, an error estimator is used to calculate the local error index value corresponding to each grid cell; The local error index value of each grid cell is compared with the refinement threshold and the coarsening threshold respectively. Grid cells with local error index values ​​greater than the refinement threshold are marked as refined, grid cells with local error index values ​​less than the coarsening threshold are marked as coarsening, and grid cells with local error index values ​​between the coarsening threshold and the refinement threshold are marked as maintained. Based on the marking results, perform an adaptive adjustment operation on the current mesh, refine the mesh cells in the refined state, and coarse the mesh cells in the coarse state to generate an adaptive mesh; The solution results on the old mesh are mapped to the adaptive mesh, and the finite element solution is re-performed based on the adaptive mesh to obtain the updated solution results. The global error index is calculated based on the updated solution results. The global error index is compared with the preset global error tolerance. If the global error index is greater than the global error tolerance, the adaptive mesh is used as the current mesh and the step of calculating the local error index value is returned to continue the iteration. If the global error index is not greater than the global error tolerance, the current solution result is output as the final solution result. The local error index value corresponding to each grid cell is calculated using an error estimator, including: The Kelly error estimator is used, with the displacement field values ​​of each grid cell as input. For each grid cell, the characteristic size and the set of boundary surfaces of the grid cell are obtained. For each boundary surface of the grid cell, the jump amount of the normal gradient of the displacement field numerical solution between the adjacent cells on both sides of the boundary surface is calculated. The jump amount is the difference of the normal gradient values ​​of the adjacent cells on both sides of the boundary surface. The square of the jump amount is integrated on the boundary surface, and the integration results of all boundary surfaces are accumulated and multiplied by the ratio of the characteristic size of the grid cell to a preset constant. The square root of the product is taken to obtain the local error index value corresponding to the grid cell. After mapping the solution results from the old mesh to the adaptive mesh and before resolving the finite element method, the process also includes: The mapped initial field values ​​are balanced and verified. The residual force vector of the mapped stress field on the adaptive grid is calculated. When the norm of the residual force vector is greater than the preset balance tolerance, the mapped initial field values ​​are balanced and iteratively corrected until the norm of the residual force vector is not greater than the balance tolerance. The balanced and corrected initial field values ​​are obtained and used as the initial conditions for resolving the problem.

2. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, The refinement threshold and coarsening threshold are determined as follows: The local error index values ​​of all grid cells are sorted in descending order of value. The local error index value corresponding to the first proportional parameter position after sorting is used as the refinement threshold, and the local error index value corresponding to the second proportional parameter position after sorting is used as the coarsening threshold. The sum of the first proportional parameter and the second proportional parameter is less than one, and the coarsening threshold is less than the refinement threshold.

3. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, Before performing adaptive adjustments to the current mesh based on the labeling results, the following steps are also included: The marking results are smoothed by iterating through all the mesh cells marked as refined and checking the marking status of their adjacent mesh cells. When a mesh cell marked as refined is directly adjacent to a mesh cell marked as coarse, the marking of the mesh cell marked as coarse is changed to the maintain state, so that the size transition between adjacent mesh cells meets the preset grade constraint conditions.

4. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, The refinement operation and the coarsening operation include: In the two-dimensional case, the refinement operation splits a triangular mesh cell into four sub-triangular cells along the lines connecting the midpoints of each side; in the three-dimensional case, the refinement operation splits a tetrahedral mesh cell into eight sub-tetrahedral cells; the coarsening operation is the reverse process of the refinement operation, which merges multiple sub-cells generated by refining the same parent cell back into the same parent cell.

5. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, Mapping the solution results from the old mesh to the adaptive mesh includes: Mapping is performed using a shape function-based interpolation method. For a target node in the adaptive mesh, the mesh cell in which the target node is located is located in the old mesh, and the local coordinates of the target node in the located mesh cell are obtained. Based on the local coordinates, the shape function values ​​corresponding to each node in the located mesh cell are calculated. The known field values ​​at each node are multiplied by the corresponding shape function values ​​and then summed to obtain the mapped value at the target node. The shape function satisfies the interpolation condition that it takes a value of one at the corresponding node and a value of zero at the other nodes.

6. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, Locating the target node's location within the existing mesh includes: Traverse all grid cells in the old grid, convert the global coordinates of the target node into the local coordinates of each grid cell, determine whether the local coordinates satisfy the reference domain range constraint of the cell type, if they do, confirm that the target node is located in the currently traversed grid cell, and use the local coordinates as the input for subsequent shape function interpolation.

7. The finite element method for simulating ground settlement risk according to claim 1, characterized in that, The global error index is calculated as follows: the local error index values ​​of all grid cells in the adaptive grid are squared and summed, and then the square root of the sum is taken. The convergence condition also includes an iteration count limit. When the number of iterations reaches the preset maximum number of iterations, the iteration loop is terminated and the current solution result is output as the final solution result. After the global error index meets the global error tolerance, the method further includes: extracting the local error index values ​​of each grid cell within a preset area of ​​interest; squaring the local error index values ​​of all grid cells within the area of ​​interest and summing them; then taking the square root of the summation result to obtain the local error index; comparing the local error index with the preset local error tolerance; if the local error index is greater than the local error tolerance, re-performing the local error index value calculation, adaptive labeling, adaptive adjustment, and mapping solution operations on the grid cells within the area of ​​interest until the local error index is not greater than the local error tolerance; wherein the local error tolerance is not greater than the global error tolerance.

8. A finite element method-based ground settlement simulation risk system, as described in any one of claims 1-7, characterized in that, include: The initial mesh generation module is used to discretize the computational domain based on the geometric information of the geological model and generate an initial mesh composed of multiple mesh cells. The finite element solution module is used to apply governing equations, material parameters and boundary conditions to the current mesh, call the finite element solver to perform numerical solutions, and obtain the displacement field and stress field values ​​of each mesh element. The error estimation module is used to calculate the local error index value corresponding to each grid cell based on the solution results of each grid cell using an error estimator. The adaptive labeling module is used to compare the local error index value of each grid cell with the refinement threshold and the coarsening threshold respectively, and to label each grid cell as a refinement state, a coarsening state or a hold-at-all state based on the comparison result. The mesh adjustment module is used to perform adaptive adjustment operations on the current mesh based on the marking results, perform refinement operations on mesh cells in the refinement state, and perform coarsening operations on mesh cells in the coarsening state to generate an adaptive mesh. The solution result mapping module is used to map the solution results on the old mesh to the adaptive mesh, and call the finite element solution module to re-solve the finite element problem based on the adaptive mesh to obtain the updated solution results; The convergence judgment module is used to calculate the global error index based on the updated solution result, compare the global error index with the preset global error tolerance, and if the global error index is greater than the global error tolerance, return the adaptive mesh as the current mesh to the error estimation module to continue iterating. If the global error index is not greater than the global error tolerance, output the current solution result as the final solution result.

Images