Soil limit analysis adaptive mesh updating method, system, device and medium

Abstract

The application discloses a soil ultimate analysis adaptive grid updating method, system, equipment and medium. The method comprises the following steps: solving the kinematic upper limit and the static lower limit in the same grid; pairing the upper limit plastic dissipation and the lower limit stress to the internal work of the same upper limit plastic strain rate, taking the non-negative difference as the local original-dual power difference; obtaining the direction inconsistency index by unifying the plastic flow directions on both sides, and synthesizing the finite element unit comprehensive priority index to determine the adjacent graph edge cost; searching for a continuous potential instability channel connecting the source end and the target end boundary section, constructing anisotropic grid metrics according to the tangent and normal directions, and redrawing; projecting the obtained field to solve again only as the initial value; comparing the coincidence degree, the through relationship and the branch number of the adjacent two rounds of continuous potential instability channels in the fixed background space, and outputting when the channel index and the upper and lower limit difference and the target limit quantity change meet the preset stop condition, otherwise returning to pairing, so that the encryption area presents a boundary through form.

Description

Technical Field

[0001] This invention relates to the field of numerical calculation technology in geotechnical engineering, specifically to an adaptive grid update method, system, equipment, and medium for geotechnical limit analysis. Background Technology

[0002] In engineering scenarios such as slope stability evaluation, tunnel surrounding rock bearing capacity verification, determination of ultimate support pressure in underground caverns, and calculation of foundation ultimate bearing capacity, it is necessary to determine target limit quantities such as strength reduction coefficients, gravity load multipliers, or ultimate support pressures, given geometric contours, material zoning, and boundary conditions. Currently, the mainstream approach is based on the kinematic upper bound theorem and the static lower bound theorem of plasticity, transforming the solution of the target limit quantities into a convex optimization problem. The upper and lower bound solutions of the target limit quantities are obtained on a discrete finite element mesh, and the difference between the two is used to measure the accuracy level of the current discretization.

[0003] Within the aforementioned framework, the commonly used adaptive approach in this field is as follows: First, an initial finite element mesh is generated through preprocessing. Then, the upper bound kinematic model or lower bound static model is solved individually on the mesh to obtain the corresponding limit quantities, velocity fields, or stress fields. Subsequently, the error contribution of each element is calculated and ranked according to indices such as plastic dissipation, equivalent shear strain rate, stress gradient, yield proximity, or the change in limit quantities between two adjacent cycles. A certain proportion of the top-ranked elements are then selected for local subdivision. Another approach is to search for paths with higher local safety indices in the element adjacency relationships to identify possible failure surface directions. Yet another approach is to construct an anisotropy metric tensor based on the gradient, Hessian matrix, or principal direction of a certain field variable, and then perform directional subdivision accordingly. The typical iterative flow of the above approaches is as follows: initial subdivision - solution - calculation of element error indices - sorting and element selection - local subdivision - interpolation of field variables from the old mesh to the new mesh as initial values ​​for hot start - re-solution. This process is repeated until the global relative boundary difference or the change in the target limit quantity between two adjacent cycles is less than a given threshold, at which point the iteration terminates and the final mesh and limit quantity results are output.

[0004] The above approaches address the issues of local activity display, path connection, and directional subdivision, but still have some limitations under complex stratigraphic conditions.

[0005] First, local high values ​​in a single field can have overlapping sources: element distortion, integral caliber deviation, abrupt changes in material parameters, and residuals from unconverged solvers can all generate numerically local high values. Selecting elements based on these sources will direct limited degrees of freedom to regions not dominated by destructive mechanisms. Second, the refined regions formed by independent element-by-element selection are distributed in a discrete, patchy pattern, not necessarily connected to predetermined engineering boundaries, and tend to retain short, closed high-value loops within the computational domain, thus lacking engineering significance in terms of boundary continuity. Third, isotropic subdivision generates nodes without corresponding accuracy gains in the shear band extension direction, rapidly exhausting the element budget. Fourth, historical information is tied to element numbers; after mesh splitting, merging, or global reconstruction, the same physical location cannot be compared. Fifth, when numerical values ​​are used as the sole termination criterion, situations may arise where the target limit change is already small, while potentially unstable regions are still migrating, fracturing, or changing branches; in such cases, the iteration is judged to have converged.

[0006] The direct consequence of the above limitations is that the output limit does not match the damage mode presented, and the determination of support design and reinforcement scope loses a reliable basis; in order to avoid risks, it is often necessary to encrypt the entire system, which leads to a significant increase in the computational scale and time consumption. Summary of the Invention

[0007] (a) Technical problems to be solved

[0008] To address the shortcomings of existing technologies, this invention provides an adaptive mesh update method for geotechnical limit analysis. It synthesizes a comprehensive priority index for finite element units using the local original-dual power difference and direction inconsistency index to determine the edge cost of the element adjacency graph. This method searches for continuous potential instability channels connecting the source and target boundary segments, constructs anisotropic mesh metrics according to their tangential and normal directions, and re-divides the mesh. The field obtained by projecting the field on the mesh before re-divides is then used as initial values ​​and re-solved. In a fixed background space, the overlap, connectivity, and number of branches of two adjacent continuous potential instability channels are compared. Output is given when the channel index, upper and lower limit differences, and target limit changes all meet preset stopping conditions; otherwise, pairing is returned, resulting in a boundary-connected form in the encrypted region. This method solves the technical problems in the prior art.

[0009] (II) Technical Solution

[0010] To achieve the above objectives, the present invention provides the following technical solution:

[0011] An adaptive mesh update method for geotechnical limit analysis solves the upper kinematic limit and lower static limit on the same mesh separately. This includes pairing the internal work done by the upper limit plastic dissipation and the lower limit stress on the same upper limit plastic strain rate element by element, and taking the non-negative difference as the local original-dual power difference.

[0012] By grouping the upper limit plastic flow direction and the predicted plastic flow direction based on the lower limit static stress at the yield surface into the same representation space, a direction inconsistency index is obtained.

[0013] The finite element synthesis priority index is synthesized from the local original-dual power difference and direction inconsistency index, and the edge cost of the element adjacency graph is determined.

[0014] Search for continuous potential instability channels connecting the preset source and target boundary segments in the cell adjacency graph, and construct anisotropic mesh metrics and divide them according to the tangent and normal directions of the continuous potential instability channels.

[0015] The field obtained by projecting the field on the mesh before re-dividing is used as the initial value and then solved again on the mesh after re-dividing.

[0016] In a fixed background space that does not change with re-division, compare the overlap, connection relationship and number of branches of two adjacent consecutive potential unstable channels; when the overlap reaches the overlap threshold, the connection relationship and the number of branches meet the preset stability conditions, and the difference between the upper and lower limits and the change of the target limit quantity reach the corresponding threshold, output the result; otherwise, return to pairing.

[0017] Furthermore, the upper limit plastic dissipation includes the intra-element integral of the dissipation support function with respect to the same upper limit plastic strain rate, and the discontinuous velocity dissipation assigned to the common edge or common surface of the element.

[0018] Internal work includes the intra-unit double-point integral of the lower limit stress and the same upper limit plastic strain rate, and the work done by the lower limit interface force on the velocity discontinuity.

[0019] The non-negative difference is the larger of the difference between the upper limit of plastic dissipation and internal energy after adopting the same allocation rule and unified caliber, and zero.

[0020] The pairing process divides the contribution of the common edge or common face equally between two adjacent units, or according to the interface length, area, and integral weight conservation, and maintains the same allocation rule in all rounds.

[0021] An adaptive grid update system for geotechnical limit analysis includes:

[0022] The upper limit kinematics solution module and the lower limit statics solution module receive the same mesh and solve the upper limit kinematics and lower limit statics respectively on the same mesh, and output the upper limit plastic dissipation, upper limit plastic strain rate and lower limit stress to the dual mechanism analysis module;

[0023] The dual mechanism analysis module pairs the upper limit plastic dissipation and the lower limit stress with the internal work done by the upper limit plastic strain rate of the element for each element. The non-negative difference is taken as the local original-dual power difference. The upper limit plastic flow direction and the plastic flow direction predicted by the static lower limit determined by the lower limit stress at the yield surface are classified into the same representation space to obtain the direction inconsistency index. The local original-dual power difference and the direction inconsistency index are combined to form the finite element comprehensive priority index, and the finite element comprehensive priority index is output to the instability channel extraction module.

[0024] The instability channel extraction module determines the edge cost of the element adjacency graph using the finite element element comprehensive priority index, searches for continuous potential instability channels in the element adjacency graph that connect the preset source end and target end boundary segments, and outputs the continuous potential instability channels to the fixed background mapping module and the anisotropic mesh update module respectively.

[0025] The fixed background mapping module compares the overlap, connection relationship and number of branches of two adjacent rounds of continuous potential instability channels in a fixed background space that does not change with re-division, and outputs the comparison results to the constraint projection and stopping control module.

[0026] The anisotropic mesh update module constructs anisotropic mesh metrics based on the tangent and normal directions of continuous potential instability channels and performs re-meshing accordingly, then outputs the re-meshed mesh to the constraint projection and stop control module.

[0027] The constraint projection and stop control module takes the field obtained by projection from the field on the mesh before remapping as only the initial value, and schedules the kinematic upper limit solution module and the static lower limit solution module to resolve on the mesh after remapping; and triggers the result output module when the comparison result of the continuous potential instability channel meets the preset topological stability condition, the upper and lower limit difference between the kinematic upper limit and the static lower limit, and the change of the target limit quantity reaches the corresponding threshold; if any condition is not met, the result of this re-solution is sent to the dual mechanism analysis module to re-execute the pairing.

[0028] The results output module outputs the re-divided mesh, continuous potential instability channels, upper and lower limit differences, and stopping reasons.

[0029] An electronic device includes a memory and a processor. The memory stores a computer program, which, when executed by the processor, implements the adaptive mesh update method for geotechnical limit analysis.

[0030] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the aforementioned adaptive grid update method for geotechnical limit analysis.

[0031] (III) Beneficial Effects

[0032] This invention provides an adaptive mesh update method for geotechnical limit analysis, which has the following advantages:

[0033] By pairing the internal work done by the upper kinematic limit plastic dissipation and the lower static limit stress on the same upper kinematic limit plastic strain rate mechanism element by element, a local original-dual power difference with a unified physical power caliber is obtained. This distinguishes between spurious high values ​​caused by element distortion or solution residuals and the non-closure of the true mechanism. The directions of the upper kinematic limit plastic strain rate and the predicted plastic flow direction of the lower static limit are transformed into a unified plastic flow direction representation space. By comparing the directions along the unsigned axis, the inconsistency index at points where both sides have high values ​​but consistent directions is suppressed. The direction confidence coefficient further suppresses spurious directions.

[0034] The finite element element synthesis priority index simultaneously determines the element ordering and the edge cost of the element adjacency graph. Paths are searched between preset source and target boundary segments, transforming the densified region from discrete patches into continuous, potentially unstable channels with interconnected boundaries. Anisotropic mesh metrics are constructed using the tangential and normal directions of these continuous potentially unstable channels. The normal direction is reduced in size to distinguish cross-band abrupt changes, while the tangential direction retains a larger size, distributing computational degrees of freedom according to the mechanism direction. When direction confidence is insufficient, the mesh degenerates to isotropic meshes. On a fixed background partition where the geometric location and numbering do not change with re-subdivision, the overlap, geometric distance, connectivity, and number of major branches of adjacent finite-width unstable channel regions are compared, binding historical information to physical space rather than element numbers. Attached Figure Description

[0035] Figure 1 This is a schematic diagram of the overall computing architecture and closed-loop data flow of the present invention;

[0036] Figure 2 This is a schematic diagram illustrating the local original-dual power pairing and the unification of plastic flow direction in this invention;

[0037] Figure 3 This is a schematic diagram of the fixed background partitioning across the finite element mesh in this invention;

[0038] Figure 4 This is a schematic diagram of the unit adjacency graph and the extraction of continuous potential instability channels in this invention;

[0039] Figure 5 This is a schematic diagram of the finite width expansion of the continuous potential instability channel and the anisotropic mesh update of the present invention;

[0040] Figure 6 This is a schematic diagram of the velocity field constraint projection, stress field constraint projection, and combined stopping of the present invention;

[0041] Figure 7 These are schematic diagrams illustrating different geotechnical engineering calculation scenarios of the present invention, wherein, Figure 7a is a schematic diagram of the two-dimensional slope calculation domain. Figure 7 b is a schematic diagram of the calculation domain for a two-dimensional tunnel cross-section. Figure 7 c is a schematic diagram of the three-dimensional computational domain;

[0042] Figure 8 This is a schematic diagram of the virtual device, electronic device, and computer-readable storage medium of the present invention;

[0043] Figure 9 This is a schematic diagram of the abnormal state handling and loop rollback of the present invention;

[0044] Figure 10 This is a schematic diagram comparing the topological stability of adjacent cyclic continuous potential unstable channels according to the present invention;

[0045] Figure 11 This is a schematic diagram of the three-dimensional planar potential instability channel and anisotropic mesh of the present invention. Detailed Implementation

[0046] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0047] Please see Figures 1-11 This invention provides an adaptive mesh update method for geotechnical limit analysis.

[0048] The common practice in finite element limit analysis of geotechnical structures is to solve the upper kinematic limit model or lower static limit model separately on a pre-generated finite element mesh, and then select the elements to be refined based on plastic dissipation, shear strain rate, stress gradient, yield proximity, or changes in limit quantities between adjacent cycles. Another approach is to search for paths with higher local safety indices on the element adjacency relationships, or to construct an anisotropy metric tensor based on the gradient of a field variable, the Hessian matrix, or the principal direction. These approaches address the issues of local activity indication, path connectivity, and directional meshing, respectively. However, local high values ​​in a single field may originate from element distortion, integral bias, material mutations, or solution residuals, and the refined patches formed by element-by-element selection do not necessarily connect to the predetermined engineering boundaries. Even if the refinement location is determined using the difference between the upper kinematic limit and the lower static limit, if this difference is only used as a one-time element error index without further constraints on channel connectivity, mesh orientation, and stopping timing, situations may still occur where the limit quantity changes are already small, but continuously potentially unstable channels continue to migrate, break, or change branches.

[0049] The technical problem this implementation aims to solve is to ensure that each round of mesh update simultaneously responds to three types of information under limited computational degrees of freedom: where the current kinematic and static bearing mechanisms are not closed, whether this location has a substantial impact on the target limit, and whether the unclosed region can remain a continuous potential instability channel connecting the specified boundary after multiple rounds of repartitioning. Therefore, it is not simply a matter of sequentially piecing together the kinematic upper bound solution, the static lower bound solution, the graph path search, and the anisotropic repartitioning. Instead, a common data source, a cross-mesh spatial reference, and a feedback triggering relationship should be specified. This allows local indices such as the same local original-dual power difference, direction inconsistency index, and target limit sensitivity to jointly enter the finite element sorting, graph edge cost, and direction fitting. Furthermore, the same continuous potential instability channel should both generate anisotropic mesh metrics and participate in the iterative joint stopping judgment.

[0050] This implementation adopts the following technical approach: Under the same finite element mesh and consistent external work normalization caliber, the upper kinematic limit solution and the lower static limit solution are obtained respectively; the internal work done by the upper kinematic limit plastic dissipation and the lower static limit stress on the same upper kinematic limit plastic strain rate mechanism is paired element by element to obtain a non-negative local original-dual power difference; the direction of the upper kinematic limit plastic strain rate and the direction of the predicted plastic flow of the lower static limit are transformed into a unified plastic flow direction representation space; and the local original-dual power difference, direction inconsistency index, target limit quantity sensitivity, and fixed background partition persistence index are synthesized into a comprehensive priority index for finite element elements. Subsequently, based on the preset source-end boundary segment and target-end boundary segment, continuous potential instability channels are extracted from the cell adjacency graph. The discrete path is expanded into a finite-width instability channel region, and the tangential and normal directions are fitted to form an anisotropic mesh metric. After re-meshing, the velocity field and stress field are projected with feasibility constraints, but only the constraint projection results are used as the initial values ​​for the next round of solution. Finally, the kinematic upper limit and static lower limit difference, the change of the target limit quantity, the coincidence of the continuous potential instability channel region, the distance between the channel boundary or centerline, and the main branch relationship are checked simultaneously. The loop ends only after the conditions are continuously satisfied.

[0051] It should be noted that the local original-dual power difference, graph path based on physical field weights, fixed background partitioning, and anisotropic mesh metric can all be used as general numerical tools. This implementation does not take the isolated use of any one of these tools as the main technical focus. The difference lies in the fact that the same local original-dual power difference, along with the direction inconsistency index, target limit quantity sensitivity, and fixed background partitioning persistence index, collectively form the finite element element comprehensive priority index. The finite element unit comprehensive priority index The graph edge cost is directly entered; after the graph path undergoes validity checks and finite width expansion, the main metric direction is directly generated; the channel then returns to the same fixed background space to determine whether to stop. If any of the feedback relationships are deleted, such as refining only according to the local primitive-dual power difference, searching for the path only according to the kinematic upper limit plastic dissipation field, or stopping only according to the change of the target limit quantity, the closed-loop update mechanism based on the same data source described in the following embodiments cannot be formed.

[0052] Please refer to Figure 1 The electronic device 100 executing the above method receives the geometry, material, boundary, and target limit configurations submitted by the field operation terminal 300 through the communication interface 130. The processor 110 calls the virtual device 200 in the memory 120 via the communication bus 140. The virtual device 200 has an engineering data preprocessing module 210 with clear code boundaries, a kinematic upper limit solution module 220, a static lower limit solution module 230, a dual mechanism analysis module 240, a fixed background mapping module 250, an instability channel extraction module 260, an anisotropic mesh update module 270, a constraint projection and stopping control module 280, and a result output module 290. The kinematic upper bound solution module 220 and the static lower bound solution module 230 output the kinematic upper bound solution results and static lower bound solution results respectively on the same mesh version. The dual mechanism analysis module 240 completes the power pairing, direction unification, and target limit quantity sensitivity calculation based on the same kinematic upper bound plastic strain rate mechanism. The fixed background mapping module 250 saves the spatial history that does not change with the mesh number. The instability channel extraction module 260 converts the element index into a graph path from the source boundary segment to the target boundary segment. The anisotropic mesh update module 270 converts the path direction into a metric tensor and triggers re-meshing. The constraint projection and stopping control module 280 completes the field variable transfer, feasibility verification, and joint stopping judgment. The result output module 290 returns the final mesh, limit quantity, and process traceability record of each round to the field operation terminal 300. If the mesh version number, material version number, or cycle number received by any module is inconsistent, the result is not published to the subsequent level. Instead, the processor 110 restores the most recent complete checkpoint from the local memory 150, thereby preventing the mismatch and use of data from different cycles.

[0053] Figure 1In the attached figures, reference numerals 100, 110, 120, 130, 140, and 150 correspond to electronic devices, processors, memory, communication interfaces, communication buses, and local memory, respectively; reference numerals 200, 210, 220, 230, 240, 250, 260, 270, 280, and 290 correspond to virtual devices, engineering data preprocessing modules, kinematic upper limit solution modules, static lower limit solution modules, dual mechanism analysis modules, fixed background mapping modules, instability channel extraction modules, anisotropic mesh update modules, constraint projection and stopping control modules, and result output modules, respectively, and do not indicate additional components. Figure 1 The input arrows on the field operation terminal 300 indicate the data stream submitted for engineering configuration to the electronic device 100 via the communication interface 130. The hardware connections on the right indicate that the processor 110, memory 120, communication interface 130, and local memory 150 are connected via the communication bus 140. The full technical names mentioned above are used consistently throughout the text. The processor 110 executes the program instructions in the memory 120, and the field operation terminal 300 submits inputs via the communication interface 130 and receives the results returned by the result output module 290.

[0054] Figure 1 In the diagram, the forked arrows in the engineering data preprocessing module 210 indicate that the kinematic upper bound solution module 220 and the static lower bound solution module 230 receive the same version of the model, rather than having to work sequentially. The two branches converging into the dual mechanism analysis module 240 indicate that the dual mechanism analysis is only performed after the kinematic upper bound solution results and the static lower bound solution results have passed version and feasibility verification. The constraint projection and stopping control module 280 returns to the kinematic upper bound solution module 220 and the static lower bound solution module 230, indicating that if the joint stopping condition is not met, the solution is re-solved on the updated finite element mesh. The flow of the constraint projection and stopping control module 280 to the result output module 290 indicates that the results are published after the joint stopping condition is met or the resource stopping condition is reached.

[0055] To avoid ambiguity in parameter symbols and subscript / superscript references, the following conventions are used throughout the text: and For finite element indexes, For adaptive cycle numbering, For fixed background partition index; superscript Indicates the result of solving the upper limit of kinematics, superscript This represents the solution result for the lower limit of statics. This represents the common edge or common face interface of the unit; Represents the physical computational domain. Indicates the first Each finite element unit region; superscript Indicates supply to the first The initial values ​​of the constraint projections used in the re-solution are... Unless otherwise specified, ... The quantity corresponding to the subscript A finite element element, with The quantity corresponding to the superscript or subscript is the first... Wheel adaptive loop. Composite superscript Indicates the work done by the lower limit static stress on the upper limit kinematic plastic strain rate mechanism; superscript The colon ":" indicates the transpose of a matrix or vector; the colon between second-order tensors indicates a double dot product. (Slope height) Hausdorff distance index with cyclic subscript The two symbols are different; the former represents the length quantity, while the latter, after being normalized by the characteristic length, is a dimensionless topological distance index.

[0056] In this implementation, the model consistency check value is a consistency check result formed by encoding geometric data, material partitions and material parameters, boundary conditions, target limit quantities, unit systems, and version fields in a predetermined order. It is used to confirm whether data generated by different functional modules correspond to the same input model and does not participate in solving the mechanical equations of the finite element limit analysis of geotechnical structures. The task mutual exclusion state is a mutually exclusive execution state set for the same task identifier to prevent the same task from being repeatedly started or concurrently written. The checkpoint verification summary is used to verify whether the loop number, model consistency check value, mesh version, and submitted result status are consistent.

[0057] In this implementation, geometric tolerance, mapping tolerance, solution residual tolerance, path validity threshold, mesh quality threshold, and joint stopping threshold are all configured parameters for the computational task. These parameters are determined before the task begins by the solver's numerical accuracy, the fixed background partition resolution, and the required engineering accuracy. The same set of decision criteria is used within the same adaptive phase. When it is necessary to switch phase thresholds or weights, the switch is only performed after the current loop is completed and the checkpoint is saved, to avoid changes in the decision criteria within the same loop. None of the above thresholds require undisclosed fixed constants as a necessary condition.

[0058] Combination Figure 1 In a single adaptive cycle, the kinematic upper bound model and the static lower bound model use a discrete carrier with consistent geometric domains, material partitions, structural interfaces, and mesh versions. The two models can employ different interpolation spaces or unknown arrangements that meet their respective theoretical requirements, but a traceable spatial correspondence must be established for each element and its common edge or common surface. This constraint ensures that local power pairing occurs within the same physical region, preventing the misinterpretation of high values ​​on different partitions as dual differences after simple nearest-neighbor concatenation.

[0059] refer to Figure 2Strict kinematic upper bound and strict static lower bound refer to results that satisfy the premises of the kinematic upper bound theorem and the static lower bound theorem of limit analysis, respectively. For example, the material adopts the assumptions of rigid-plastic or ideal plasticity, a convex yield criterion, and a related flow rule adapted to the kinematic upper bound theorem. The kinematic upper bound discrete field satisfies kinematic permission and external work normalization, while the static lower bound discrete field satisfies equilibrium, stress boundary, and yield feasibility. When the material adopts non-related flow, a non-convex yield surface, or a constitutive relation approximated by smoothing, and cannot provide strict theoretical guarantees, the following process can still be used for finite element mesh updates, but the corresponding results are called kinematic upper bound type results and static lower bound type results.

[0060] Combination Figure 2 The local primitive-dual power difference is the nonnegated power difference between the plastic dissipation required by the kinematic upper limit velocity mechanism on a single element and its distribution interface, and the internal work done by the static lower limit stress state on the same kinematic upper limit plastic strain rate tensor and the same interface velocity discontinuity. It is not a typical stress peak, nor is it a misnomer of the static lower limit stress as static lower limit dissipation; its engineering boundary includes the element integral contribution, the interface contribution of the conserved distribution, and calibration terms set for a unified integral caliber. It should be verified by order of magnitude and conservation against the global boundary difference obtained from the kinematic upper limit limit and the static lower limit limit.

[0061] refer to Figure 2 The unified plastic flow direction representation space converts the plastic flow direction information obtained from solving the kinematic upper limit and the static lower limit predicted plastic flow direction information obtained from solving the static lower limit into a directly comparable representation space. The kinematic upper limit side is obtained by normalizing the kinematic upper limit plastic strain rate tensor; the static lower limit side first projects the static lower limit stress state to the nearest yield surface point, and then obtains the static lower limit predicted plastic flow direction from the yield function gradient or the applicable plastic potential function gradient. The two directions are compared according to the unsigned axial relationship. When the directions are opposite but represent the same shear zone axis, the axial sign itself does not increase the directional difference, and the directional inconsistency index is suppressed. The direction confidence coefficient is used to suppress unstable directional information at positions where plastic activity is close to 0, the static lower limit stress is far from the yield surface, or the principal eigenvalues ​​are not sufficiently separated. The value range is 0 to 1, where 0 indicates that the directional information at that point is unreliable, and 1 indicates that the directional information at that point has the highest confidence.

[0062] Combination Figure 2The target limit quantity sensitivity is the dimensionless response of the target ultimate load multiplier, safety factor, strength reduction factor, ultimate support pressure, or other specified limit quantities to the local strength parameters or yield constraint scale of the element. This response is preferably obtained from the original variables, dual variables, and partial derivatives of the Lagrangian function at the optimal solution; if the solver cannot provide stable multipliers, it can be obtained from two-sided perturbations of local parameters or adjoint calculations. Sensitivity is only used to allocate mesh resources and construct path costs, and does not change the material parameters input by the user.

[0063] Figure 2 The plastic strain rate tensor shown refers to the plastic strain rate tensor output by the upper limit kinematics solution module 220, rather than the cumulative plastic strain; the upper limit kinematics plastic dissipation contribution, the lower limit static stress state, the local original-dual power difference, the direction inconsistency index, and the target limit sensitivity are all data or calculation results and do not constitute new hardware modules. Figure 2 In the process, the plastic strain rate tensor is directly input into the dual mechanism analysis module 240 to perform power pairing and plastic flow direction normalization based on the same upper limit plastic strain rate mechanism; the lower limit static stress state is directly input into the dual mechanism analysis module 240 to perform internal work on the same upper limit plastic strain rate tensor and interface velocity discontinuities; to avoid line congestion, Figure 2 Not all connections between the data items were drawn. The direction confidence coefficient is... Figure 2 The independent box represents the direction inconsistency index, which is constrained by the direction confidence coefficient within the dual mechanism analysis module 240.

[0064] refer to Figure 3 A fixed background partition is an auxiliary pixel, voxel, triangle, tetrahedron, tetrapod leaf region, or octagonal leaf region that covers the physical computation domain before the adaptive loop begins and whose geometric position and number remain unchanged in subsequent re-meshment. It does not participate in the equation assembly of the kinematic upper bound model or the static lower bound model; it only receives the area or volume overlap ratio of the current candidate unstable region, saves the fixed background partition persistence index updated exponentially, and back-maps this index to the updated finite element mesh. Because historical information is bound to the physical space rather than the finite element element number, the same position can still be compared after elements are split, merged, or reconnected.

[0065] Figure 3 The fixed background partition and its persistence index are shown. These correspond to fixed spatial partitions and the cross-cycle persistent states stored in those partitions, respectively; candidate unstable regions are mapped to fixed background partitions according to their area or volume coverage ratio for updating. Then, the finite element persistence index is formed by inverse mapping the overlap weights of the updated finite element mesh and the fixed background partition. The three thick solid lines in the figure represent the spatial extent of the same candidate region in different discrete carriers, and do not represent material interfaces or the final actual slip surface. The polygonal mesh on the right only summarizes the adjacency changes after re-subdivision; the actual finite element elements still use the triangular, quadrilateral, or hybrid meshes determined in step S101 and its alternative implementations.

[0066] Combination Figure 4 A continuous potential instability channel is a path or set of paths in the element adjacency graph that leads from a virtual source node to a virtual target node, checked for boundary connectivity, minimum length, prohibited boundary crossings, and the proportion of finite element elements with high overall priority. This term is not equivalent to the final true slip surface, but rather represents a continuous candidate region where the current kinematic upper limit discretization and static lower limit discretization still exhibit inconsistencies in their contribution to failure and direction of plastic flow, and which influence the target limiting quantity. The path passes through... Figure 5 Only after the centerline is smoothed and the finite half-width is expanded to form a finite-width unstable channel region can it be used as the spatial input for actual mesh remapping.

[0067] Figure 4 The small hollow circles and thin solid lines in the diagram represent element graph nodes and candidate adjacent graph edges, respectively. The large circles on the left and right represent virtual source nodes and virtual target nodes that do not correspond to actual finite element elements. Only representative virtual connections are shown in the diagram; in reality, all qualified candidate finite element elements for the source and target boundary segments are connected. A cross indicates that the corresponding candidate graph edge has been excluded from the path search set. The thick solid arrows indicate the traversal order of the main channel from the source to the target. Dashed lines only highlight non-overlapping sections of branch channels relative to the main channel; shared sections are still connected by traversable graph edges. The text consistently uses "target end" and "branch channel" as technical terms.

[0068] refer to Figure 5 Anisotropic mesh metric is a symmetric positive definite tensor composed of eigenvectors and eigenvalues ​​representing the tangential, normal, and corresponding target dimensions of a continuous potential instability channel. Two-dimensional continuous potential instability channels are described by one tangential and one normal; three-dimensional planar potential instability channels are described by two in-plane tangentials and one normal. This metric stipulates that the mesh generator uses smaller dimensions in the normal direction to resolve cross-band abrupt changes, and retains larger dimensions in the tangential direction to avoid invalid degrees of freedom. When there is insufficient confidence in the direction, excessive curvature, proximity to geometric sharp corners, or substandard element quality, the anisotropy ratio is suppressed until it degenerates into isotropy.

[0069] Figure 5From left to right, the diagram shows the discrete path, smooth centerline, finite-width unstable channel region, its surrounding transition layer, and the updated finite element mesh. Tangential and normal arrows indicate the unsigned principal axes of the anisotropic mesh metric; reference numeral 270 corresponds to the anisotropic mesh update module 270. The dashed lines between the anisotropic mesh update module 270 and the centerline, and between the updated finite element mesh and the received channel information and the updated finite element mesh that has passed the mesh quality check, respectively, do not indicate a re-search for the path. The text consistently uses "finite-width unstable channel region" and "transition layer" as their corresponding technical names.

[0070] Combination Figure 6 Constraint projection is an initial value transfer operation between the old and new meshes, rather than replacing the updated finite element mesh solution with the interpolated result. Velocity field constraint projection minimizes the difference between the old and new velocities while satisfying fixed, symmetric, and defined velocity boundaries; stress field constraint projection minimizes the difference between the old and new stresses while satisfying discrete equilibrium, defined stress boundaries, and yield inequalities. After projection, the kinematic upper limit solution module 220 and the static lower limit solution module 230 must still resolve on the updated finite element mesh. Only results that again meet the corresponding feasibility requirements will proceed to the next mechanism judgment.

[0071] refer to Figure 6 and Figure 10 The topological stability of a continuous potential instability channel refers to the stability of adjacent adaptive cycle continuous potential instability channels in a common background space, simultaneously satisfying the following conditions: overlap of finite-width instability channel regions, geometric distance between channel boundaries or centerlines, continuity from the source to the target, number of major branches, and stability of the boundary segment to which the main endpoint belongs. A small change in the target limit alone is insufficient to determine the topological stability of a continuous potential instability channel; even with a high region overlap ratio, joint stopping is not triggered if the main continuous potential instability channel breaks or the number of major branches changes. Therefore, the joint stopping condition simultaneously constrains both the numerical boundary difference and the morphology of the continuous potential instability channel.

[0072] Figure 6 In the process, the updated finite element mesh is subjected to velocity field constraint projection and stress field constraint projection to form two types of initial constraint projections. After merging, these projections are used for complete re-solution on the updated finite element mesh, topological stability judgment, and joint stopping judgment. Reference numeral 280 corresponds to constraint projection and stopping control module 280. The dashed line in the figure indicates that the joint stopping judgment is a control function of constraint projection and stopping control module 280. When all necessary criteria are continuously satisfied, the process proceeds to result output module 290; if any necessary criterion is not satisfied, the process proceeds to the next adaptive loop. The results of velocity field constraint projection and stress field constraint projection are always used only as initial values ​​for re-solution.

[0073] Combination Figure 1The internal calculation loop shown in this main embodiment uses the limit stability analysis of a two-dimensional slope as an example, but the following closed-loop update mechanism does not depend on specific slope height, slope angle, or material values. Electronic device 100 can be a standalone workstation, a computing server, or a computing cluster composed of multiple computing nodes; its processor 110 is connected to memory 120, local memory 150, and communication interface 130 via communication bus 140. Field operation terminal 300 submits geometric files, material partitions, weak interlayers, free surfaces, displacement constraints, load boundaries, target limit quantities, and calculation tolerances to engineering data preprocessing module 210. Engineering data preprocessing module 210 converts the input into a model object with a version number; kinematic upper bound solution module 220 and static lower bound solution module 230 each receive the same model object; dual mechanism analysis module 240 only accepts kinematic upper bound solution results and static lower bound solution results where the mesh version and material version are simultaneously consistent. The fixed background mapping module 250, the instability channel extraction module 260, the anisotropic mesh update module 270, and the constraint projection and stopping control module 280 sequentially return the fixed background mapping verification status, the continuous potential instability channel effective status, the updated finite element mesh quality status, and the joint stopping status, forming a feedback chain of calculation results - continuous potential instability channels - anisotropic mesh re-partitioning - re-solution. If any status fails, the result output module 290 only records the failure stage and the recovery point, and does not display the intermediate fields that failed the verification as the final continuous potential instability channels.

[0074] Figure 7 include Figure 7 a, Figure 7 b and Figure 7 c, where Figure 7 'a' represents the two-dimensional slope calculation domain. Figure 7 b represents the computational domain of the two-dimensional tunnel cross-section. Figure 7 c represents the three-dimensional computational domain. Figure 7 a and Figure 7 The black dots in b are only used as schematic anchor points for the source or target boundary section. The actual calculation uses the nodes or finite element set within the corresponding boundary section. Figure 7 In section c, the source and target boundary surface segments correspond to the endpoint constraints of the 3D scene. The regular mesh is only used to illustrate partitions and relative positions, and does not limit the finite element type or mesh density.

[0075] like Figure 7 As shown in Figure a, the exemplary slope model uses slope height... As a characteristic length, the bottom edge is set with a fixed constraint, and the side edge is set with a normal constraint that matches the limit analysis discretization used. The slope is a free surface. The material parameters are directly adopted from the user-input unit weight, cohesion, internal friction angle and necessary weak interlayer parameters. Values ​​not provided are not assumed in the program.

[0076] The main embodiment can employ a convex Mohr-Coulomb yield criterion and associated flow rule, using the strength reduction factor or gravity load multiplier as the target limit quantity; if other convex yield criteria are used, only the dissipative support function, yield surface projection, and corresponding constraints need to be replaced. The boundary outside the computational domain is kept at a distance sufficient to reduce the truncation effect from the slope; this distance is verified by the user based on boundary sensitivity and is not claimed to be a fixed multiple in this embodiment. The above settings provide a reproducible boundary and constitutive call relationship while avoiding the fabrication of specific engineering dimensions and measured material data. Among these, The slope height is the slope height of the two-dimensional slope model, and in this main embodiment, it is used as a characteristic length for spatial dimension normalization.

[0077] In step S101, as Figure 1 As shown, the engineering data preprocessing module 210 obtains the slope calculation domain, the location of the slope top and toe, the material interface, and the boundary conditions, with the slope height as the basis. Or computational domain reference length The spatial dimensions are unified, and an initial triangular mesh is generated. The initial mesh maintains geometric consistency for slope surfaces, slope crest angles, slope toe angles, weak interlayers, and material interfaces. The element size near the weak layer is no larger than the input thickness of that layer to prevent the first round of discretization from directly crossing the structure. The engineering data preprocessing module 210 simultaneously generates a geometric template for a fixed background partition, but this template does not participate in the finite element solution.

[0078] For example, the background partition feature size can be set to 1 to 2 times the expected minimum grid size, and this value is only used as the program initialization range. If the background partition cannot cover the entire computational domain or duplicate numbers appear, the processor 110 locks step S102 and requests the engineering data preprocessing module 210 to re-partition. This coverage check ensures that subsequent historical quantities have unique spatial locations. The reference length is predetermined based on the geometric scale of the computational domain and remains unchanged within the same computational task; for two-dimensional slope scenarios, slope height is preferred. In scenarios where the slope height is not clearly defined, the following approach is adopted. .

[0079] Acceptance of step S101 and Figure 1In an alternative implementation of the preprocessing relationship shown, the initial mesh can be a quadrilateral or a triangle-quadrilateral hybrid mesh, and the fixed background partition can be a regular pixel, a quadrilateral leaf area, or an independent triangular partition independent of the finite element mesh. When the background partition crosses the material interface, the engineering data preprocessing module 210 saves the coverage share for different materials separately, without directly averaging the persistence index of the fixed background partitions on both sides of the material interface; when the boundary projection error exceeds the geometric tolerance, the engineering data preprocessing module 210 retains the original geometric edge and reduces the size of the adjacent background cells. This processing prevents the erroneous merging of adjacent but physically different continuous potential instability channels due to an excessively coarse common background.

[0080] In step S102, combined Figure 1 The parallel solution relationship shown and Figure 2 The data chain on the left, the kinematic upper limit solution module 220, uses nodal velocities, element velocity parameters, or common edge velocities as unknowns on the initial mesh, applies essential velocity boundaries and external work normalization, and obtains the kinematic upper limit limit. Velocity field , No. Plastic strain rate of each element And element and interface dissipation. The statics lower limit solution module 230 uses element stress, interface force, or equilibrium variables as unknowns on the same mesh, applies discrete equilibrium, stress boundary, and yield feasibility conditions, and obtains the statics lower limit limit. Stress field Yield margin and optimization multipliers. Both modules can call interfaces for linear programming, nonlinear convex programming, second-order cone programming, or semidefinite programming, as long as the returned solutions satisfy the corresponding license conditions; there are no restrictions on specific commercial software or processor models.

[0081] After the feasibility check is passed for the solution results on both sides of step S102, the constraint projection and stop control module 280 first calculates the global relative boundary difference.

[0082]

[0083] Combination Figure 2 In the formula, For dimensionless global relative boundary difference, and These are the upper kinematic limit and lower static limit on the current grid, respectively. To prevent positive numbers with a denominator of zero, their values ​​should be higher than the zero-value noise level caused by numerical precision, and significantly lower than the typical absolute value of the target limit in normal calculations.

[0084] If the target limit is allowed to be negative, the denominator uses its absolute value or a preset dimensionless reference quantity. The constraint projection and stop control module 280 also reads the kinematically permissible residuals, static equilibrium residuals, boundary residuals, and maximum yield exceedance; if any residual exceeds the allowable value in the solver configuration, the dual mechanism analysis module 240 does not calculate subsequent indices, but instead requires the corresponding solver module to increase iterations, reduce the step size, or restart from the most recent feasible checkpoint. This avoids using local high values ​​where the solver has not yet converged as the basis for mesh refinement. and They have the same units of measurement to ensure global relative boundary differences. It is a dimensionless quantity.

[0085] In step S103, combined Figure 2 The relationship between the upper limit of kinematic plastic dissipation and the lower limit of static stress state information shown in the dual mechanism analysis module 240 is incorporated into the dual mechanism analysis module 240. Calculate the upper limit of kinematic plastic dissipation contribution for each unit:

[0086]

[0087] And the internal work contribution of the static lower limit stress to the same kinematic upper limit plastic strain rate tensor and the corresponding interface velocity discontinuity:

[0088]

[0089] Combination Figure 2 , For the first Each unit area For the dissipative support function corresponding to the convex yield criterion, For the discontinuous velocity dissipation of the common edge or common surface allocated to this cell, This refers to the work done by the interface force at the lower limit of statics on the same velocity discontinuously.

[0090] The contribution of the common interface can be allocated 1 / 2 to each of the two adjacent units, or it can be allocated according to the interface length, area, or integral weight conservation; once selected, the same rule is maintained throughout the entire cycle. This power pairing based on the same kinematic upper limit plastic strain rate mechanism places the kinematic upper limit solution and the static lower limit solution in the same physical power caliber, thus enabling the differentiation between regions where the plastic activity is large on the kinematic upper limit side but the static lower limit side can still meet the bearing conditions, and regions where there is indeed a local inconsistency between the original-dual power on the kinematic upper limit side and the static lower limit side. and All have the dimension of power; Indicates the first A common edge or common surface interface between a finite element and its adjacent elements.

[0091] refer to Figure 2 The dual mechanism analysis module 240 further calculates:

[0092]

[0093] Combination Figure 2 , The difference between the non-negative local primitive and dual power. This is a calibration correction term used to unify the caliber of the integral contribution and the interface allocation contribution within the cell under the selected common interface conservation allocation rule. Its value corresponds to the calibration correction amount allocated to the cell when the local original-dual power difference of each cell is calibrated proportionally according to the global difference. When the selected allocation rule has already made the two types of contribution calibers consistent and no calibration is required, this calibration correction term is 0. To normalize the local original-dual power difference index, To regularize positive numbers.

[0094] Duality analysis module 240 with by Comparison of converted global power differences; when the difference is within the integral and rounding tolerance, it can be compared according to the global difference. The calibration is performed proportionally. If the difference exceeds the limit, the normalization of external work, the order of unit integral, and the interface contribution are rechecked for duplication. The calibration only changes the allocation caliber and does not change the upper limit of kinematics or the lower limit of statics. This conservation check prevents the continued execution of loops where the normalized local original-dual power difference index and the global limit are contradictory. and Having the same power dimensions, Therefore, it is a dimensionless normalized index. The summation index in the normalized local primitive-dual power difference index formula traverses all finite element elements in the current cycle.

[0095] Combination Figure 2 In cases where the solver cannot output complete interface forces or replacement paths for dual variables, the dual mechanism analysis module 240 can normalize the upper kinematic limit plastic dissipative activity and the lower static limit yield proximity, respectively, and then calculate the difference between the two as a simplified local primitive-dual difference index. However, the processor 110 writes the simplified pairing state into the result metadata and does not interpret this index as a strict local primitive-dual power difference. If only a small amount of interface data is missing, the corresponding interface share is recovered first by conservation of adjacent elements; if the global power verification still fails after recovery, the anisotropic mesh update of the loop is locked. This provides a degradation method compatible with the old solver interface while maintaining the semantic boundaries between strict pairing results and replacement results.

[0096] In step S104, combined Figure 2As shown in the yield projection and flow direction link, the dual mechanism analysis module 240 normalizes the upper limit plastic strain rate to:

[0097]

[0098] Combination Figure 2 , Let be the dimensionless upper kinematic upper limit plastic flow direction tensor. For the Frobenius norm, To prevent division of positive numbers by zero by zero active units.

[0099] For the static lower limit stress state, the dual mechanism analysis module 240 sets the stress according to the preset stress distance. Projected to the nearest yield point and calculate

[0100]

[0101] refer to Figure 2 , To predict the direction of plastic flow using the dimensionless lower limit of statics. Input the yield function of the material. For regularization positive numbers. In the two-dimensional case, the dual mechanism analysis module 240 uses the tensor from the kinematic upper limit plastic flow direction. And the lower limit of statics for predicting the direction of plastic flow Candidate zone normals are extracted and orthogonal tangents are constructed; in three dimensions, tensor eigenvectors, maximum shear direction, or local principal component directions are used. The projection algorithm can be analytical projection, constraint minimum distance solution, or lookup interpolation on the yield surface parameters; its output is only used for direction comparison and does not replace the static lower limit feasible stress state. In the formula... Let be the independent variable of the stress tensor corresponding to the yield function. It has the same dimensions as the Frobenius norm of the yield function gradient, and thus only serves as a regularization function to prevent the denominator from being zero.

[0102] in, Indicates the first The nearest yield surface point of the static lower limit stress state of a finite element element under the adopted stress distance metric; It has the same dimensions as the plastic strain rate norm.

[0103] Combination Figure 2 The direction inconsistency index output shown can be adjusted according to the direction confidence coefficient. Construct, in which, To normalize the plastic activity, To normalize the yield proximity, The primary eigenvalue separation degree, and all three are limited to 0-1, therefore It is also limited to 0-1. The indicator for inconsistency in direction is...

[0104]

[0105] Combination Figure 2 , The value is not less than zero, and the sign of the axial direction is eliminated for the squared term; When approaching zero, even if the numerical eigenvectors are rotated, the inconsistency index of direction is suppressed. For example, when the plastic activity or yield proximity is below the robust quantile threshold of the current cycle, the corresponding normalized plastic activity or normalized yield proximity can be directly set to 0; this threshold is derived from the field distribution of the current cycle rather than a fixed material constant. When using the non-associated flow rule, the gradient of the yield function is replaced by the gradient of the plastic potential function, and the result is labeled as a mechanism for comparing kinematic upper bound results with static lower bound results. This confidence constraint prevents arbitrary characteristic directions in the stable region from being amplified by the mesh generator.

[0106] Normalized plasticity, normalized yield proximity, and principal eigenvalue separation are all mapped to 0 to 1 based on the original quantities corresponding to the current cycle according to the pre-fixed robust quantiles, maximum values, or equivalent normalization scales; the same normalization rule is maintained within the same computational task.

[0107] In step S105, the dual mechanism analysis module 240 selects the mechanism that can characterize the first... Local parameters of the unit strength or yield constraint scale Among them, the dimensionless Lagrangian function obtained by normalizing the objective limit quantity through the corresponding optimization problem is... Calculate the original sensitivity:

[0108]

[0109] Combination Figure 2 , The dimensionless absolute sensitivity before normalization. It can correspond to dimensionless transformations of cohesion, internal friction angle, tensile strength, or Hawke-Brown parameters, but the selection rules remain consistent within the same embodiment. The dual mechanism analysis module 240 normalizes it according to robust quantile values, maximum values, or sums. .

[0110] If the dual multiplier is unstable, then the dual multiplier is unstable. Apply exemplary small perturbations of equal positive and negative proportions and resolve separately, using two-sided differencing for verification; the perturbation amplitude must be greater than the fluctuation of the target limit caused by the solver tolerance while still remaining within the local linear range. To control computational load, differencing is only applied to... Higher-ranking candidate elements are processed, while the remaining finite element elements use the available analytical sensitivity; if analytical sensitivity cannot be stably obtained, it is set to 0. This step preserves connected elements whose local activity is not prominent but which control the continuity of potential instability channels and affect the target limit. For dimensional local strength parameters, they are first dimensionlessly proportional to a pre-determined input benchmark value within the same material partition that remains unchanged throughout all adaptive cycles, and then used as the basis for further processing. When the parameter itself is dimensionless, its input value is used directly. The sensitivity normalization rule remains consistent across the same computational task; when the normalization denominator is 0, the corresponding... Take 0.

[0111] In step S106, as Figure 3 As shown, the fixed background mapping module 250 calculates the first... Round of candidate instability regions and the first Background partitions Area or volume overlap ratio And recursively:

[0112]

[0113] Combination Figure 3 , For a persistent indicator of a fixed background partition, Take values ​​from 0 to 1. Historical weights ranging from 0 to 1. For example, during initial program debugging, the weights can be... This indicates that current coverage accounts for 30% and historical coverage for 70%; this value is only used to illustrate the recursive relationship, and the formal calculation is calibrated according to the degree of grid disturbance and the expected migration speed of the mechanism. First round of orders The fixed background mapping module 250 then presses the current unit and each... The overlap ratio, Weighted mapping to finite element element persistence index .

[0114] If the overlap weight and deviation exceed the mapping tolerance, the fixed background mapping module 250 re-executes the geometric intersection without using the nearest element number, thus preserving the history of the same physical region after complete re-partitioning. The fixed background partitioning persistence index is then inversely mapped to the finite element persistence index. When the current finite element element overlaps with each fixed background partition in area or volume ratio, the non-negative weight is used, and all weights corresponding to the finite element element are normalized to a weight sum of 1.

[0115] like Figure 3As shown, when an overly coarse background partition causes two spatially separated, potentially unstable channels to fall into the same partition, the fixed background mapping module 250 can enable quadtree or octree encryption for that local area, but the parent-child partition identifiers remain unchanged in subsequent loops; in memory-constrained replacement paths, a sparse background key-value structure can be used, saving only the partitions covered by the candidate region. When communication is interrupted, the fixed background mapping module 250... The model consistency check value, mesh version, and cycle number are written to local memory 150. After recovery, recursion continues only if the model consistency check value is consistent. This persistence and consistency check suppresses single-round noise and prevents breakpoint recovery from bringing the history of another model into the current calculation.

[0116] Combination Figure 3 No results yet in the first round At that time, the dual mechanism analysis module 240 is available. , and Temporary candidate scores are generated based on the normalized local primitive-dual power difference index. Inconsistent direction indicators and target limit sensitivity The high quantile intervals of the temporary candidate scores are selected, or candidate regions are selected based on the cumulative local primitive-dual power contribution; from the second round onwards, the mapping obtained in the previous round is used. participate Then from the current Update candidate coverage.

[0117] When using the quantile rule, at least a preset proportion of temporary candidate finite element elements with higher scores are retained in each round; when using the cumulative contribution rule, selection is made from high to low until... The preset ratio is achieved. Minimum retention rules can also be set near the source and target ends, but these only take effect within the configured boundary buffer. This clarifies the data sequence between the first and subsequent loops, preventing the fixed background partition persistence metric from being referenced in a loop before it has been generated.

[0118] In step S107, the dual mechanism analysis module 240 utilizes Figure 3 Inverse mapping and to , , and Perform quantile truncation and 0-1 normalization respectively, and calculate the finite element synthesis priority index.

[0119]

[0120] Combination Figure 3 , , , and All are non-negative and sum to 1, controlling the contributions of the local original-dual power difference, direction inconsistency index, target limit sensitivity, and fixed background partition persistence index, respectively. For example, the early stages of non-first rounds can use... First, locate the non-closed region using the positioning mechanism; when Once the threshold is reached, it can be switched to This directs more grid resources towards the target's critical area. The values ​​above are non-limiting initialization examples and do not represent the optimal weights measured in practice.

[0121] The duality analysis module 240 only performs duality analysis within a neighborhood that does not cross material interfaces or structural discontinuities. Smoothing is performed 1 to several times, preserving the local peak positions, thereby reducing noise without smoothing out the narrow band region that truly controls the formation of continuous potential unstable channels. , , and All are dimensionless weights, therefore the finite element element comprehensive priority index It is a dimensionless quantity.

[0122] like Figure 7 As shown in Figure a, in the two-dimensional slope of this main embodiment, the field operation terminal 300 configures the preset free surface section near the slope crest and the preset free surface section near the slope toe as the source boundary section and target boundary section of the continuous potential instability channel, respectively. If the input model contains a weak interlayer, the boundary area between the weak interlayer and the upper and lower free surfaces or the boundary of the computational domain can be changed to the endpoint boundary section. This configuration is part of the engineering boundary conditions, rather than being determined by the instability channel extraction module 260 after the fact based on the high value location. For symmetrical slopes or slopes with multiple exits, multiple combinations of source boundary section and target boundary section can be configured and searched separately, and finally filtered by a common path validity rule. Clearly defining the endpoint boundary section can exclude short closed high value loops within the computational domain, giving the continuous potential instability channel an engineering meaning of boundary continuity.

[0123] In step S108, as Figure 4 As shown, the instability channel extraction module 260 converts each finite element element into a graph node, and establishes shared edges, shared surfaces, or two elements that satisfy a preset proximity distance as graph edges. Each node carries... The system includes the tangential direction, direction confidence, target limit sensitivity, material type, and boundary type of candidate continuous potential instability channels. Each edge carries the center distance between two adjacent finite element elements, the direction of the center connection, the cross-boundary type, and the passability status. When encountering a complete lining, an insurmountable structural boundary, or a user-specified blocking interface, the instability channel extraction module 260 deletes the corresponding graph edge. When encountering a weak interlayer or a slip-allowed interface that is already clearly defined in the input, it only reduces the cost along the interface extension direction, without automatically creating a weak layer based on calculated high values. This structure binds discrete field values ​​to geometric connectivity conditions, preventing paths from crossing physically blocking entity boundaries.

[0124] Combination Figure 4 When the mesh has suspended nodes, mismatched interfaces, or mixed elements, adjacency relationships can be determined by common geometric measures rather than node numbers; if the common edge length or common surface area of ​​two finite element elements is lower than the geometric tolerance, a direct connection is not established. The instability channel extraction module 260 performs verification on the number of graph nodes, graph edges, and connected components, and returns isolated nodes to the engineering data preprocessing module 210 for mesh topology checking. In an alternative implementation, a graph can be built on a fixed background partition, and then continuous potential instability channels can be back-mapped to finite element elements; however, the graph node weights must still originate from the same local original-dual power difference, direction inconsistency index, target limit sensitivity, and fixed background partition persistence index data link. This alternative method is suitable for scenarios where finite element meshes undergo frequent global reconstruction.

[0125] In step S109, combined Figure 4 The diagram illustrates the connection relationships between virtual source nodes, virtual target nodes, and candidate elements. The instability channel extraction module 260 establishes virtual source nodes and virtual target nodes in the diagram, connecting the virtual source nodes to all candidate finite element elements in the source boundary segments, and connecting all candidate finite element elements in the target boundary segments to the virtual target nodes. Virtual connection edges use zero cost or a uniform non-negative small cost, ensuring the search process selects the starting and ending points with lower actual costs. Candidate finite element elements in the endpoint boundary segments must simultaneously be within the buffer distance of the configured segment and connected to effective elements in the computational domain; the buffer distance can be given according to the local initial mesh size. If the set of candidate finite element elements in the source boundary segments or the set of candidate finite element elements in the target boundary segments is empty, the instability channel extraction module 260 does not arbitrarily assign a comprehensive priority index to the finite element elements in the entire domain. Instead of changing the highest finite element to the endpoint, return the endpoint configuration error to prevent the obtained path from losing the predetermined engineering boundary constraints.

[0126] In step S110, for Figure 4 The adjacent units shown and The unstable channel extraction module 260 employs an implementable non-negative edge cost.

[0127]

[0128] Combination Figure 4 , For the cost of graph edges, The distance between the centers of two adjacent finite element elements. To prevent positive numbers with a denominator of zero, The directional confidence combination value is obtained by combining the directional confidence coefficients of two adjacent finite element elements according to a pre-fixed combination rule and is limited to 0 to 1. The direction of the line connecting the centers, The predicted instability tangent is formed by unifying the signs of the upper limit plastic flow direction and the lower limit predicted plastic flow direction for each of two adjacent finite element elements. Then, the signs of the two predicted instability tangents are unified, and an average direction is formed according to a pre-fixed weighting rule. This weighting rule remains unchanged throughout all adaptive loops of the same computational task. The direction penalty coefficient, For material or structural boundary penalties, The corresponding weights. and The larger the value, the lower the basic path cost; when the connection deviates from the high-confidence prediction direction, the sine square term increases; when the direction confidence is insufficient, the direction term weakens on its own; when crossing an interface that is not suitable for crossing, the multiplicative boundary term increases.

[0129] From the same Simultaneously controlling the candidate unit sorting and path cost can prevent graph search from deviating from the preceding dual mechanism data and establishing a separate set of evaluation criteria. , , and All are treated as dimensionless quantities. For dimensionless regularized positive numbers; and It is a non-negative dimensionless adjustment coefficient. The directional confidence coefficients of two adjacent finite element elements are obtained according to a pre-fixed combination rule and are limited to 0 to 1. The combination rule remains unchanged throughout all adaptive cycles of the same computational task. Within the same material region or at a permissible interface, a value of 0 is acceptable. Non-negative values ​​are used at material or structural boundaries requiring penalty. For boundaries explicitly prohibited from crossing, the corresponding graph edges are directly deleted without assigning a finite penalty. It has the dimension of length. It has a length-based cost dimension and is only used for comparing path costs within the same graph. and Using the same angular measurement caliber, the difference between the two is only used for the sinusoidal square direction penalty term, therefore this direction penalty term is a dimensionless quantity.

[0130] refer to Figure 4 ,all When non-negative, the unstable channel extraction module 260 can use Dijkstra's algorithm to obtain the main continuous potential unstable channel; when multiple branches are required, the previous method can be used. Perform deduplication on the shortest path, minimum cost flow, or paths with similar costs. A positive integer configured for the program does not imply that a continuous potential instability channel must have a fixed number of branches.

[0131] Another alternative implementation is to use graph cut or region growing with orientation constraints, but it still needs to meet three conditions: the source end must be connected to the target end, the region cost must be affected by the local primitive-dual power difference, the orientation inconsistency index, the target limit sensitivity, and the fixed background partition persistence index, and the output must be able to be converted into a finite-width unstable channel region. If a general directed weight is used and non-negativity cannot be guaranteed, a path algorithm adapted to the corresponding weight conditions can be used instead, and this scheme is not limited by the algorithm name.

[0132] like Figure 4 As shown, the instability channel extraction module 260 checks each candidate path: the path connects the specified source and target ends; the physical length is not less than the minimum value set according to the slope characteristic length; the length ratio of finite element elements with a high comprehensive priority index reaches a preset threshold; the path does not cross the prohibited boundary; and there are no sudden jumps in the corners of adjacent polyline segments that cannot be explained by local curvature. For example, the current round... The top 20% are considered as candidates for higher finite element unit (FEM) synthesis priority indicators, and their length in the path must account for no less than 60%. This value is only used to illustrate the judgment action. Locally closed loops, isolated short paths, and cross-boundary paths that fail to pass the test are eliminated. Multiple paths with similar costs and spatial separation are retained as potential instability channels with continuous branches, while highly overlapping paths are merged according to common segments. This path validity screening further restricts the mathematical shortest path to one that is engineering-accessible.

[0133] In step S111, as Figure 5 As shown, the unstable channel extraction module 260 connects the center of the unit according to the path node sequence to form a polyline, and then moves within the sliding window according to... Principal component fitting or spline fitting is performed using directional confidence weights to obtain a smooth centerline. The half-width of the finite-width instability channel region is determined by the local normal target size, the current mesh size, and the minimum number of refinement layers, covering at least one layer of candidate elements on each side of the centerline. Elements intersecting the centerline buffer region form the finite-width instability channel region, with 1-3 layers of transition elements added around it to gradually bring the size back to the stable region scale. The half-width represents the numerically resolved influence zone and is not interpreted as the actual shear band thickness in an ideal plasticity model without material internal scale. The conversion from path to finite area ensures that the subsequent mesh generator receives an operable region rather than a zero-width polyline.

[0134] Combination Figure 5 At path bifurcation points, the instability channel extraction module 260 fits the centerlines of each branch and uses a smooth transition in the common segment. When the angle or spacing between two branches is below the geometric tolerance, they are merged into a single finite-width instability channel region. When two branches remain spatially separated and both pass the endpoint check, they are retained separately. When a continuous potential instability channel is near a slope toe angle, material cusp, or load discontinuity, the window must not cross the discontinuity to avoid the curve smoothing pulling the centerline of the continuous potential instability channel out of the effective material domain. In an alternative implementation, the path can be projected onto a fixed background partition first, and then the centerline can be obtained using a connected region skeleton; the obtained centerline is still extended back into the finite element physical coordinates. This local constraint maintains the true geometric transformation while suppressing the pseudo-rotation caused by mesh jaggedness.

[0135] In step S112, combined Figure 5 As shown, for the centerline and local principal axis, the anisotropic mesh update module 270 performs a weighted covariance decomposition on the path point coordinates within the neighborhood of each centerline, and takes the maximum principal direction as the local tangential unit vector. Take its two-dimensional orthogonal direction as the normal unit vector. Adjacent tangential directions first undergo sign uniformity, making... Then smooth it within the range that does not cross material interfaces and geometric angles.

[0136] If the separation of the principal eigenvalues ​​is low or the local curvature of the centerline is too large, the anisotropic mesh update module 270 reduces the direction confidence and shortens the fitting window; if a stable direction still cannot be obtained, the location is marked as an isotropic transition region. This degradation rule prevents the crossover of slender cells caused by the reversal of positive and negative eigenvectors or the non-uniqueness of directions near bifurcation.

[0137] In step S113, as Figure 5 As shown, the anisotropic mesh update module 270 updates the mesh based on the finite-width unstable channel region. and Determine the normal target size

[0138]

[0139] Combination Figure 5 , Based on the size, and It is a non-negative adjustment coefficient. and As the upper and lower limits of the allowable size, the above calculations of the maximum and minimum values ​​of the inner and outer layers are used to limit the normal target size within a certain range. and between; and These are the target dimensions in the normal and tangential directions, respectively. The anisotropy ratio is not less than 1 and does not exceed the upper limit allowed by the mesh generator. When the direction confidence is low, the curvature of continuous potential instability channels is large, the mesh is near a bend, or the current mesh quality is poor, the anisotropy mesh update module 270 will adjust the anisotropy ratio. Compression towards 1; an increase in the anisotropy ratio is only permitted when the curvature of the continuous potential instability channel is small and the direction is stable. This size allocation increases resolution along the normal to the continuous potential instability channel, while reducing resolution at nodes where there is no corresponding error gain along the tangential to the continuous potential instability channel. , , , and All have the dimension of length; , and All of these are dimensionless parameters.

[0140] refer to Figure 5 A two-dimensional metric tensor is written as

[0141]

[0142] Combination Figure 5 , It is a symmetric positive definite metric tensor. The orthogonal basis is formed by tangential and normal directions. The anisotropic mesh update module 270 updates the mesh outside the finite-width instability channel region based on a lower finite element synthesis priority index. and lower target limit sensitivity Maintain the original size or restrict coarsening, and apply adjacent size gradient constraints to the entire metric field.

[0143] If the eigenvalue is not positive, the adjacent size ratio exceeds the limit, or the direction change exceeds the mesh generator's capacity, the anisotropic mesh update module 270 first truncates the eigenvalue and reduces its value. Instead of transmitting the disease attitude quantity to the repartitioning procedure, the transition layer range is expanded. This pre-screening reduces the likelihood of repartitioning producing inverted or minimal corner units. The eigenvalues ​​have the dimension of the reciprocal of the square of the length, and are used to convert the side length of physical space into the side length of metric space; It is the orthogonal directional basis matrix.

[0144] In step S114, combined Figure 5 The updated finite element mesh shown can be generated by the anisotropic mesh update module 270, which can utilize anisotropic Delaunay meshing, leading edge advancement, constrained quadtree, directional bisection, or background metric-based re-meshmentation. During generation, the positions of material boundaries, weak interlayers, free surfaces, load boundaries, and structural interfaces are maintained. New nodes on curved boundaries are projected back to their original geometry; interfaces where topological changes are not permitted are either kept in node matching or discretized using predetermined interfaces. The anisotropic mesh update module 270 then checks the minimum angle, aspect ratio, Jacobian determinant, metric space side length, and boundary deviation. If the check fails, the local anisotropy ratio is sequentially reduced, a transition layer is added, and the process is retried; if it still fails, the updated finite element mesh is revoked, the previous valid finite element mesh is retained, and a failure status is returned to the constraint projection and stop control module 280. This closed loop does not use the generated file as a substitute for mesh availability judgment.

[0145] Combination Figure 5 In the local re-meshing replacement method, only finite element elements in the finite-width unstable channel region and transition layer are deleted, while the far-field nodes remain unchanged, thereby reducing the projection range. In the global re-meshing method, the mesh is regenerated using the entire metric field, thereby improving the overall element quality. Both methods must preserve the overlap mapping between the old finite element mesh elements, the updated finite element mesh elements, and the fixed background partition. If computational resources reach the preset limit, the anisotropic mesh update module 270 first checks the finite element synthesis priority index. Lower and target limit sensitivity Lower-range coarsening is performed to free up the unit budget, followed by normal refinement of the finite-width unstable channel region; if the upper limit is still exceeded, refinement is stopped and the reason for resource termination is recorded. This budget constraint aligns refinement priority with the target limit sensitivity and the mechanism of continuous potential unstable channels.

[0146] In step S115, as Figure 6 As shown, the constraint projection and stop control module 280 uses least squares with boundary constraints or Projection yields updated initial values ​​of finite element mesh velocity.

[0147]

[0148] Its boundary constraints are Combining Figure 6 , For the physical computation domain, The definition is in the physical computation domain. The second norm on, For the candidate velocity field on the updated finite element mesh, For the first The velocity field obtained by solving the upper limit of wheel kinematics. This represents the sampling mapping of the velocity field from the old mesh to the updated finite element mesh. For discrete boundary constraint operators or matrices corresponding to fixed, symmetric, and specified velocity boundaries. This is the specified value vector corresponding to the boundary constraint. Only indicates the first The initial values ​​for the constraint projection are used to solve the upper limit of the wheel's kinematics. When using velocity discontinuous elements, the constraint projection and stop control module 280 also maps the old interface discontinuous conservation to the neighborhood of the newly generated finite-width unstable channel region, or recalculates the velocities from the adjacent finite element elements after projection. If the boundary constraints conflict with the sampled values, the boundary constraints take precedence; if the projection problem is unsolvable, a warm start is discarded and a cold start is performed using the zero field or analytically permissible field that satisfies the boundary conditions. This process utilizes historical direction information without compromising kinematic permissibility. , and Use the same velocity dimension.

[0149] In step S116, combined Figure 6 The stress field constraint projection branch shown, constraint projection and stop control module 280 solve the stress constraint projection.

[0150]

[0151] Simultaneously apply discrete static equilibrium, specified stress boundaries, continuity of interface normal forces, and Combining Figure 6 , For the candidate stress field on the updated finite element mesh, For the first The stress field obtained by solving the lower limit of wheel statics. This is a sampling mapping from the old stress field to the updated finite element mesh. For the first Initial values ​​of constraint projections for solving the lower limit of wheel statics. This represents the material yield function value. When the directly sampled value is outside the yield surface, the constraint projection and stop control module 280 pulls it back into the yield feasible region in the constraint problem; when the constraint projection is infeasible, the yield or equilibrium conditions are not relaxed, but the static lower bound model is cold-started, and if it still fails, it is rolled back and re-partitioned. This priority ensures that the initial value transfer used to accelerate the solution belongs to the static lower bound theoretical conditions. , and The same stress axis is used.

[0152] In step S117, as Figure 6 As shown, the processor 110 reassembles the kinematic upper bound solution module 220 and the static lower bound solution module 230 on the updated finite element mesh to solve the optimization problem, respectively using... and A warm start allows for the complete acquisition of new [products / services]. , The velocity field, stress field, and dual variables are considered. Only when the new feasible residual passes the same threshold as in step S102 is the result registered as a valid input for the mechanism analysis in the next cycle. The continuous potential instability channels in the current cycle compared in step S118 are the continuous potential instability channels that have been formed and saved in steps S108 to S112 of this cycle and used to drive the current anisotropic mesh update. Figure 6 The connection from completely resolving the updated finite element mesh to determining topological stability serves as a general representation of feasibility verification, state registration, and topological comparison. Initial constraint projection values ​​and new solution results are saved separately, and the result output module 290 is not allowed to mistakenly label initial value fields as final fields. Therefore, after each mesh change, the upper kinematic limit and lower static limit properties are re-established based on the complete constraints on the updated finite element mesh.

[0153] In step S118, when the current cycle is the first round or the previous round's consecutive potential instability channels have been marked as historical references and prohibited from participating in comparisons, only the current round's consecutive potential instability channels are saved as subsequent comparison benchmarks, and topological stability comparisons between adjacent rounds are not performed; starting from the cycle where there are valid consecutive potential instability channels in the previous round that can participate in comparisons, combined with Figure 6 Topological stability determination process and Figure 10 The common background comparison relationship shown is used by the fixed background mapping module 250 to map the first background. Wheel and the first The finite-width unstable channel regions of the wheel are mapped to a common background space, and the following calculations are performed:

[0154]

[0155] Combination Figure 10 , For the intersection and union ratio of regions, For the first The size of the unstable channel region of the wheel with finite width is taken as area or number of background partitions in two dimensions, and volume or number of voxels in three dimensions. The fixed background mapping module 250 also converts the centerlines or boundaries of the two wheels into a common coordinate point set, and calculates... or Normalized Hausdorff distance The comparison process involves checking whether the source and target ends remain connected, whether the number of main branches is consistent, whether the main endpoints belong to the same boundary segment, and whether the average direction change exceeds the limit. Spatial comparison using a common background eliminates the interference of different grid numbers and cell shapes on the determination of topological stability. Indicates the first Limited width instability channel area for wheels; via slope height Or computational domain reference length After normalization, This is a dimensionless geometric distance index. When any cycle to be compared does not form a finite-width unstable channel region that passes the validity check, the region crossover ratio is not calculated, and the cycle is determined not to meet the topological stability condition for a continuous potential unstable channel.

[0156] Figure 10 The solid line outline and solid line center line in the text represent the first... The finite-width instability channel region of the wheel and its centerline, the dashed outline and the dashed centerline represent the first... The diagram shows the finite-width instability channel region of one wheel and its centerline. The diagonal line indicates the intersection of the finite-width instability channel regions of two wheels. The bidirectional distance arrows are only used to indicate the geometric separation between the centerlines or boundaries of two wheels; the actual Hausdorff distance is not shown. Calculated based on the centerline point set or boundary point set in a common background space. The topological stability of a continuous potentially unstable channel is determined by the following factors: region overlap, geometric distance, connection between the source and target ends, number of major branches, and the boundary segment to which the main endpoint belongs. If any one of these necessary conditions is not met, the channel is deemed topologically unstable. Figure 10 The unstable channel (finite width) of the current cycle and the unstable channel (finite width) of the previous cycle marked in the text correspond to the finite width unstable channel regions of the current and previous cycles, respectively. The finite width unstable channel region is used as the technical name in the main text.

[0157] Combination Figure 10In one embodiment, the field operation terminal 300 is equipped with a display interface. After the constraint projection and stop control module 280 completes the topology stability judgment for the current round, the result output module 290 sends the comparison data to the field operation terminal 300. The display interface overlays the current and previous rounds of continuous potential instability channels with a finite width of half-width along the normal direction, the intersection area of ​​the two, the bidirectional geometric distance between the centerlines of the two rounds of continuous potential instability channels, and the current values ​​of overlap, continuity, and number of branches. When the overlap reaches the overlap threshold, the continuity and number of branches meet the preset stability conditions, and the difference between the upper and lower limits and the change of the target limit quantity reach the corresponding thresholds, the display interface displays the stop reason text of the output result; when any condition is not met, the display interface displays the name of the unmet item and its value, and displays the return pair status indicator, allowing the user to review the specific reason why the current round failed to stop. The above display does not change the joint stop judgment itself, but only visualizes the topology comparison result and the joint stop basis.

[0158] In step S119, if there is no valid target limit from the previous round in the first round, only the target limit of the current round is saved as the subsequent comparison benchmark, and the rate of change of adjacent cycles is not calculated; starting from the cycle where there is a valid solution result from the previous round, if... Figure 6 As shown, the constraint projection and stop control module 280 calculates the rate of change of the target limit quantity in adjacent cycles:

[0159]

[0160] Combination Figure 6 , In one embodiment, the first is fixedly selected. The upper limit of kinematics, the lower limit of statics, and the arithmetic mean of the upper limit of kinematics and the lower limit of statics of a wheel, or a specified safety factor. Indicates the first The target limit quantity selected for wheels with the same diameter To regularize positive numbers. When Not greater than the threshold difference Not greater than the target limit change threshold Not less than the overlap threshold Only when the distance threshold is not greater than the distance threshold, and the connection between the source and target ends, the number of main branches, and the boundary segment to which the main endpoint belongs are all stable, is this round counted as a single satisfaction.

[0161] In some embodiments, the average direction change not exceeding a direction threshold may also be used as an additional criterion. A non-limiting configuration for illustrating the judgment logic may be... , , , The formal threshold is determined by the solver tolerance, background resolution and engineering accuracy, and must not be set below the achievable numerical noise level. and Having the same dimensions, therefore It is a dimensionless rate of change; , and normalized All of these are dimensionless criteria.

[0162] like Figure 6 As shown, if the target limit has stabilized and If the branching relationship is unstable, the constraint projection and stop control module 280 continues to cycle and prioritizes maintaining the weight of the inconsistent direction index and the weight of the fixed background partition persistence index; if the continuous potentially unstable channel is stable but Still large, continuing in the normalized local primitive-dual power difference index Higher regional refinement; if all mechanism indicators are stable but the unit budget has reached its limit, then stop with resource conditions and explicitly record the unmet criteria, rather than pseudo-joint convergence. In an alternative implementation, the average symmetric surface distance can be used instead of the Hausdorff distance, or background connected component labels can be used instead of explicit branch counts, but the three types of information—region overlap, geometric distance, and connectivity—must not be completely omitted.

[0163] In step S120, Figure 1 The result output module 290 outputs the final finite element mesh to the field operation terminal 300. , Specify target limit quantity The parameters include: velocity field, stress field, centerline of continuous potential instability channels, finite-width instability channel region, anisotropic mesh metric field and stopping reasons, and the mesh version, model consistency check value, solution tolerance, and feasibility residual for each round. , , , Image search results , The branch status is written to local memory 150.

[0164] The output file includes a cycle number and a checkpoint verification summary. If the field terminal 300 fails to read the file, it can be retransmitted without resolving. This process traceability chain allows users to verify which local primitive-dual power difference, continuous potential instability channels, and stopping states triggered a particular grid update.

[0165] like Figure 7As shown in b, the first parallel embodiment is used for the ultimate rock limit analysis of a two-dimensional tunnel cross-section. The engineering data preprocessing module 210 reads the tunnel outline, surrounding rock zones, ground stress or load boundaries, and support boundaries. Based on the failure mode to be analyzed, it sets the boundary sections near the arch crown or abutment as source-end boundary sections and the sections near the sidewalls or invert as target-end boundary sections. When comparing left and right failures, two sets of endpoint boundary sections are established separately for the left and right sides. The kinematic upper limit solution module 220 and the static lower limit solution module 230 still operate on the same cross-sectional grid. Figure 2 The calculation methods for the local original-dual power difference, the static lower limit prediction of plastic flow direction, and the target limit sensitivity remain unchanged. Therefore, this embodiment alters the engineering geometry and the endpoint constraints of the continuous potential instability channel, instead of setting a separate refinement basis for the single high stress value at the tunnel shoulder.

[0166] Combination Figure 7 b and Figure 4 The preferred embodiment of the tunnel is the front The shortest path is used to extract potentially unstable channels in the surrounding rock that may separate to the left and right. The unstable channel extraction module 260 first verifies whether each path connects the specified arch and sidewall or invert arch sections, and then compares the common segments and spatial spacing of the paths. Paths with similar costs and spatial separation are retained as two spatially separated branches of potentially unstable channels. When the paths highly overlap, their common center lines are merged. The fixed background mapping module 250 saves the persistence index of the fixed background partition corresponding to each branch with regular pixels or independent triangular backgrounds, and simultaneously saves the continuous potential unstable channel attribution label for each background partition to prevent left and right branches from being mistakenly counted as the same branch in the overlapping projection. If the symmetric model only calculates half of the image, the fixed background mapping module 250 restores the correspondence according to the symmetric boundary rules when comparing the number of branches, rather than judging the full-image topology change based on the number of half-image paths.

[0167] Figure 7 The dashed lines in a only represent paths or centerlines. The anisotropic mesh update still uses the finite-width unstable channel region formed in step S111 as input. Figure 7 If no path or centerline is drawn in b, whether it is continuous depends on... Figure 4 The adjacency graph and validity check shown are the criteria. Figure 7 b. Linings or supports not shown are given by the engineering input and are discretized for inspection according to prohibited crossing edges or given interfaces.

[0168] refer to Figure 7 b and Figure 5 At the arch shoulder bends, lining ends, and near material interfaces, the anisotropic mesh update module 270 reduces the curvature or geometric angle of the mesh based on the larger curvature or geometric angle. If necessary, this local area is designated as an isotropic transition zone; larger target dimensions are retained tangentially in continuous potential instability channel sections at the gentle arch and sidewalls. Anisotropic mesh updates can replace only the finite-width instability channel region and its transition layer to preserve the far-field mesh; alternatively, the entire cross-section can be re-divided according to the metric field. Regardless of the path taken, when the lining or complete support boundary is configured as a no-crossing zone, the instability channel extraction module 260 deletes the graph edges crossing that boundary, and the anisotropic mesh update module 270 also maintains the interface geometry. The direct effect is to allow the mesh orientation to respond to continuous potential instability channels in the surrounding rock that have passed the validity check, rather than arbitrarily stretching along the tunnel curve boundary.

[0169] Combination Figure 7 b、 Figure 6 and Figure 10 The tunnel implementation still performs velocity field constraint projection and stress field constraint projection separately, and completely resolves on the updated finite element mesh. Joint stopping judgment, except... , , and In addition, the number of left and right branches, and the arch and sidewall sections connected to each branch are compared. If the change in ultimate support pressure is less than the threshold but a branch with a continuous potential instability channel migrates from the sidewall to the invert arch, the constraint projection and stop control module 280 does not stop; if a branch with a continuous potential instability channel temporarily disappears in one round, the persistent history of the fixed background partition corresponding to that branch is retained but the current coverage is reduced, and the number of stable branches is updated only after the continuous cycle no longer occurs. This method avoids frequent switching of branch states caused by accidental graph search results.

[0170] like Figure 7 As shown in c, the second parallel embodiment is used for three-dimensional tunnel faces, underground chambers, or three-dimensional slopes. The initial mesh can be tetrahedral, hexahedral, or a hybrid mesh, and the fixed background uses regular voxels, tetrahedral background regions, or octagonal leaf regions. The kinematic upper limit solution module 220 and the static lower limit solution module 230 output three-dimensional velocity, plastic strain rate, and stress tensor. The dual mechanism analysis module 240 still uses the same kinematic upper limit plastic strain rate mechanism for power pairing. The fixed background mapping module 250 updates the fixed background partition persistence index according to the volume overlap ratio. The source-end boundary surface segment and the target-end boundary surface segment are configured by the actual engineering boundary surface, such as between the near-field region of the tunnel face and the overlying free surface or a specified far-field segment. The algorithm does not subsequently change the endpoints based on the obtained high values.

[0171] Combination Figure 7 c and Figure 4The 3D element adjacency graph uses finite element elements as nodes and common surfaces as edges. Potential 3D planar instability channels can be obtained through multiple representative paths and minimum cost flows, or by using graph cuts constrained by the source-target endpoint to form connected candidate instability regions, and then reconstructing the 3D planar potential instability channels from multiple path point sets. When using graph cuts, the region term is still determined by... The boundary terms are still constrained by direction confidence and prohibited crossing information. When using multi-path reconstruction, the unstable channel extraction module 260 removes non-connected, excessively short, and boundary-crossing paths before fitting the common surface of the three-dimensional planar potential unstable channels. Thus, the three-dimensional implementation changes the overall route of topology extraction, while the local original-dual power difference, direction inconsistency index, target limit quantity sensitivity, and fixed background partition persistence index remain common inputs.

[0172] refer to Figure 7 c. Figure 5 and Figure 11 The anisotropic mesh update module 270 performs weighted covariance decomposition on the local point cloud of the three-dimensional planar potentially unstable channel, obtaining two in-plane tangents. , and 1 normal direction ,structure

[0173]

[0174] Combination Figure 11 , , and These are the target dimensions in the tangential direction of two planes. This is the normal dimension. When a linear potentially unstable channel cannot stably fit the surface, one direction is set as the path tangent, and the other two directions are determined by the local plastic activity gradient or isotropic dimensions. If the three principal eigenvalues ​​are not sufficiently separated, the anisotropic mesh update module 270 reduces the neighborhood of the voxel to isotropic to avoid the arrangement of three-dimensional slender tetrahedrons around an uncertain axis. , and All have the dimension of length. It is an orthogonal directional basis matrix consisting of two in-plane tangent unit vectors and one normal unit vector.

[0175] Figure 11 In the diagram, the source-end boundary surface segment, the target-end boundary surface segment, and the black dot represent the source-end constraint, target-end constraint, and their schematic anchor points in the 3D scene, respectively, without limiting the endpoints to a single node. The in-plane tangent, in-plane tangent, and normal arrows are the principal axes of anisotropic mesh measurement; the two tangential target dimensions correspond to... and Normal target size corresponds The three together form The updated finite element mesh in the figure represents a schematic diagram of an oriented volumetric mesh in the finite-thickness neighborhood of a three-dimensional planar potential instability channel, rather than an independent surface detached from the surrounding volumetric mesh.

[0176] Combination Figure 7 c. Figure 6 , Figure 10 and Figure 11 In the 3D implementation, the mesh can be globally re-partitioned to maintain volume mesh quality, or the finite-thickness regions on both sides of the 3D planar potential instability channel can be locally re-partitioned; both re-partitioning methods preserve voxel mapping and perform 3D equilibrium and yield constraint projection. When performing topological stability assessment of continuous potential instability channels, Calculated by volume The calculation is performed based on the point set of the three-dimensional planar potential instability channel, while simultaneously comparing whether the three-dimensional planar potential instability channel connects the specified boundary surface segment, the number of major connected pieces, and the change in surface normal. Therefore, even if the three-dimensional planar potential instability channel changes from a single surface to two spatially separated surfaces, and the change in the limiting quantity is very small, the constraint projection and stop control module 280 continues to loop. This embodiment illustrates that this scheme is not limited to the two-dimensional shortest path, but covers the three-dimensional connected surface and its oriented volume mesh formed by the aforementioned local original-dual power difference, direction inconsistency index, target limiting quantity sensitivity, and fixed background partition persistence index.

[0177] like Figure 9 As shown, the processor 110 treats an adaptive loop as a computational transaction to be submitted, and sequentially performs input and version checks, kinematic and static feasibility checks, continuous potential instability channel validity checks, updated finite element mesh quality checks, velocity field constraint projection and stress field constraint projection checks, and communication and version checks.

[0178] If any check fails, the following steps are executed: resubmit input, cold start, restricted isotropic refinement, reduction of local anisotropy ratio, undoing the updated finite element mesh, or restoring the checkpoint, and then entering a rollback loop. Undoing unsubmitted updated finite element meshes, initial values ​​for constraint projections, temporary states of fixed background partitions, and result records, and restoring the most recent complete checkpoint. The current loop is only committed after all checks pass; if the joint stopping condition is not met, the loop proceeds to the next loop with this committed state. Figure 9 This is a summary diagram of abnormal state control, without changing the order of S101 to S120.

[0179] Combination Figure 9 Input check failure branch and Figure 1 and Figure 7a) When the geometric boundary is not closed, material regions overlap, material parameters are missing, unit systems are inconsistent, or the target limit quantity is undefined, the engineering data preprocessing module 210 sets the model state to invalid input, does not call the kinematic upper limit solution module 220 and the static lower limit solution module 230, and does not use preset constants to replace missing parameters. After the field operation terminal 300 resubmits the corrected configuration, the engineering data preprocessing module 210 recalculates the model consistency verification value; only after the model consistency verification value, version information, and geometric verification all pass will the task mutual exclusion state be released. This blocking prevents incomplete engineering input from being packaged into a numerical instability mechanism in the future.

[0180] Combination Figure 9 Feasibility check failed branches and Figure 1 and Figure 6 When step S102 or S117 exceeds the solution time limit, the iteration stops but the feasibility residual exceeds the limit. If the numerical tolerance is exceeded, or the yield limit cannot be eliminated, the constraint projection and stop control module 280 marks the current solution as invalid and prohibits the fixed background mapping module 250 from updating the fixed background partition persistence index. The processor 110 first cancels the warm start and resolves using the cold start initial values ​​that satisfy the boundary conditions; if no feasible solution is still obtained, it restores the mesh, background state, and results of the previous complete cycle, outputting the reason for the solution's non-convergence without outputting new kinematic upper and static lower bound conclusions. After restoration, if the user chooses to continue, the numerical tolerance or mesh quality configuration can be adjusted, but it is not allowed to force the solution to pass by abandoning the kinematic or static permission conditions.

[0181] Combination Figure 9 Channel check failed branches and Figure 4 When step S110 fails to obtain a valid through path from the source end to the target end, the unstable channel extraction module 260 first checks the candidate finite element set of the source end boundary section, the candidate finite element set of the target end boundary section, and the prohibited crossing boundary, and then expands the endpoint buffer and adds the front within a preset range. Candidate paths or switching graph cuts, minimum cost flow, and region growing with directional constraints are considered. If the alternative algorithm still fails the path validity check, no validly connected potential unstable channels are declared in this round; priority indices can be integrated in the finite element element. Higher-order regions undergo restricted isotropic refinement and are re-searched in the next round. Alternatively, potential instability channels from the previous round can be retained but marked as historical references and prohibited from participating in joint stop judgment. This downgrade maintains computational progress while prioritizing the finite element synthesis index. A clear distinction is made between higher candidate regions and effectively connected continuous potential instability channels.

[0182] Combination Figure 9Quality inspection failure branches and Figure 5 When the mesh generated in step S114 has inverted elements, boundary offsets, aspect ratio exceeding limits, or unqualified metric space quality, the anisotropic mesh update module 270 first reduces the local [quality]. Then, a transition layer is added and the mesh is re-divided; if consecutive retries still fail, it degenerates into a locally isotropic size. The degenerate mesh must still pass the same quality threshold; otherwise, the update is canceled and the original mesh is restored. The anisotropic mesh update module 270 writes the geometric location, metric eigenvalue, and number of retries for each failure to the checkpoint, ensuring that failed combinations are not reused in the next round. This mechanism prioritizes mesh quality constraints over pursuing a larger anisotropy ratio.

[0183] Combination Figure 9 Projection check failed branches and Figure 6 When the velocity field constraint projection in step S115 conflicts with the essential boundary, or the stress field constraint projection in step S116 cannot simultaneously satisfy equilibrium, stress boundary, and yield constraints, the constraint projection and stop control module 280 discards the direct interpolation value and cold-starts the corresponding solver; if cold-starting is still not feasible, the updated finite element mesh is canceled. Projection failure will not modify the persistence history of the fixed background partition, because the persistence index of the fixed background partition is only updated after the complete kinematic upper bound solution and static lower bound solution on the updated finite element mesh pass. When the projected field fails the check, the previously stored potential unstable channels in the fixed background space remain unchanged and are not replaced by channels that failed the check in this round. This transactional sequence avoids inconsistencies where the mesh has been updated but the mechanism history still comes from old feasible solutions.

[0184] Combination Figure 9 Communication check failed branches and Figure 1 When communication between the field operation terminal 300 and the electronic device 100 is interrupted, the output receipt times out, the same task is repeatedly started, or the checkpoint verification fails, the processor 110 refuses to write repeatedly using the task identifier, loop number, and task mutual exclusion status; the currently executing module writes the current complete loop to the local memory 150 and then stops publishing externally. When communication is restored, the field operation terminal 300 first compares the model consistency verification value with the most recently completed loop number, and then continues from that loop; if the verification is inconsistent, a new task is created without splicing the old history. Thus, data transmission anomalies only affect task scheduling and do not change the computational state that has already passed feasibility and topology verification.

[0185] Figure 9When communication checks fail, verifiable complete results are first saved as unsubmitted; transmission resumes when communication is restored and the loop number, model consistency check value, and checkpoint check summary are consistent. The most recent checkpoint is only restored when the identifier is inconsistent, the data is incomplete, or the checkpoint fails; when resubmitting input causes the model consistency check value to change, a new task is created, and the old history is not appended.

[0186] like Figure 8 As shown, the virtual device 200 can be implemented as a function library in the same process, several classes with clear interfaces, multiple local processes, or services connected through a task queue. The engineering data preprocessing module 210, corresponding to step S101, receives engineering input and outputs an initial model with a version; the kinematic upper bound solution module 220 and the static lower bound solution module 230, corresponding to steps S102 and S117, receive the same finite element mesh and output the kinematic upper bound solution results that satisfy the kinematic permission conditions and the static lower bound solution results that satisfy the static feasibility conditions, respectively; the dual mechanism analysis module 240, corresponding to steps S103 to S105 and S107, receives the field variables obtained from the kinematic upper bound solution, the field variables obtained from the static lower bound solution, and the optimization multipliers, and outputs... , , and If any input version is inconsistent, the module terminates the call chain through an exception return value and does not produce a default result.

[0187] Combination Figure 8 The fixed background mapping module 250, corresponding to steps S106 and S118, receives candidate unstable regions, old finite element mesh geometry, updated finite element mesh geometry, and fixed background partition status. It then performs overlap mapping, fixed background partition persistence index recursion, and common background space topology comparison, and... , and The data is sent to the dual mechanism analysis module 240 or the constraint projection and stopping control module 280. The instability channel extraction module 260, corresponding to steps S108 to S111, receives the finite element synthesis priority index. The anisotropic mesh update module 270, corresponding to steps S112 to S114, receives the centerline of the continuous potential unstable channel, the finite width unstable channel region, the mesh target size limit, and the original engineering geometry. It then fits the local orientation, constructs metrics, and outputs the updated finite element mesh that passes quality checks or a clearly failed state. If the unstable channel extraction module 260 has no valid path, it does not send an empty continuous potential unstable channel to the anisotropic mesh update module 270 as a success result. The module also performs mapping, searching, verification, and outputs the centerline and finite width unstable channel region.

[0188] refer to Figure 8 The constraint projection and stop control module 280, corresponding to steps S115 to S119, receives the old finite element mesh, the updated finite element mesh, the historical velocity field and historical stress field, the feasibility residual, and the topological stability index of the continuous potential instability channel, and performs velocity field constraint projection, stress field constraint projection, re-solution scheduling, and joint stop judgment, respectively. The result output module 290, corresponding to step S120, receives the status confirmed by stop or abnormal termination, and generates the final result package and the round-by-round process traceability record. The above module names only indicate functional boundaries and do not require them to physically occupy independent processors. Two or more functional modules can be merged into the same program component, and one functional module can also be executed in parallel by multiple computing processes, as long as the input, output, version verification, and abnormal status maintain the above relationships.

[0189] like Figure 8 As shown, the electronic device 100 includes a processor 110, a memory 120, a communication interface 130, a communication bus 140, and a local memory 150. The processor 110 can be one or more general-purpose processing units, graphics processing units, or computational units suitable for numerical optimization; brand and model are not limited. The memory 120 stores the running mesh, matrix, field variables, graph structure, and program instructions. The local memory 150 stores project inputs, checkpoints, and result files. The communication interface 130 connects to the field operation terminal 300 or other computing nodes. The communication bus 140 transmits instructions, data addresses, and status between various hardware components. After reading instructions from memory 120, processor 110 sequentially calls the engineering data preprocessing module 210, kinematic upper limit solution module 220, static lower limit solution module 230, dual mechanism analysis module 240, fixed background mapping module 250, unstable channel extraction module 260, anisotropic mesh update module 270, constraint projection and stop control module 280, and result output module 290 in virtual device 200 according to the abnormal return value, and submits, rolls back, or stops a loop.

[0190] Figure 8In the accompanying drawings, reference numerals 100, 130, 140, 150, 300, and 400 correspond to electronic devices, communication interfaces, communication buses, local memory, field operation terminals, and computer-readable storage media, respectively; the full technical names mentioned above are used consistently in the text. The virtual device 200 is located within memory 120, indicating that program instructions and runtime data can be saved and loaded, without limiting each functional module to independent hardware or residing in the same process. Lines without arrows in the figure indicate direct or via communication bus 140 communication coupling; the bidirectional arrows between the field operation terminal 300 and communication interface 130 indicate configuration submission and result return; the arrows between communication bus 140 and computer-readable storage media 400 indicate accessibility, with read-only media only performing read operations. The module calling order remains as per steps S101 to S120.

[0191] Combination Figure 8 In the distributed implementation, the kinematic upper bound solution module 220 and the static lower bound solution module 230 can be executed by different computational processes. The processor 110 or the scheduling process only triggers the dual mechanism analysis module 240 after the model consistency verification values, mesh versions, and tolerance configurations returned by both are consistent. When communication fails, each process writes the complete results to the local memory 150 and returns to an uncommitted state. This parallel configuration only changes the task execution location and does not change the pairing on the same finite element mesh or the requirement for resolving the updated finite element mesh. In the single-machine implementation, all modules share the memory 120 and are called sequentially, which eliminates network transmission but still retains the version field and checkpoints.

[0192] like Figure 8 As shown, the computer-readable storage medium 400 can be a non-volatile solid-state memory, flash memory, read-only memory, erasable and rewritable memory, magnetic storage medium, or other tangible medium capable of storing program instructions. It can also be used in conjunction with volatile memory serving as operating space. When the program stored in the medium 400 is executed by the processor 110, the electronic device 100 completes steps S101 to S120, and retains the most recent complete checkpoint in case of program errors. The medium embodiment protects the instructions and data structures that can execute the above methods, and does not treat the pure propagation signal itself as a tangible storage medium.

[0193] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. An adaptive grid update method for limit analysis of soil and rock, wherein the upper kinematic limit and the lower static limit are solved separately on the same grid, characterized in that, include: For each element, the upper limit plastic dissipation and the lower limit stress are paired with the internal work done by the lower limit stress on the same upper limit plastic strain rate, and the non-negative difference is taken as the local original-dual power difference; By grouping the upper limit plastic flow direction and the predicted plastic flow direction based on the lower limit static stress at the yield surface into the same representation space, a direction inconsistency index is obtained. The finite element synthesis priority index is synthesized from the local original-dual power difference and direction inconsistency index, and the edge cost of the element adjacency graph is determined. Search for continuous potential instability channels connecting the preset source and target boundary segments in the cell adjacency graph, and construct anisotropic mesh metrics and divide them according to the tangent and normal directions of the continuous potential instability channels. The field obtained by projecting the field on the mesh before re-dividing is used as the initial value and then solved again on the mesh after re-dividing. In a fixed background space that does not change with re-division, compare the overlap, connection relationship and number of branches of two adjacent consecutive potential unstable channels; when the overlap reaches the overlap threshold, the connection relationship and the number of branches meet the preset stability conditions, and the difference between the upper and lower limits and the change of the target limit quantity reach the corresponding threshold, output the result; otherwise, return to pairing.

2. The adaptive mesh update method for geotechnical limit analysis according to claim 1, characterized in that: The upper limit plastic dissipation includes the intra-element integral of the dissipation support function with respect to the same upper limit plastic strain rate, and the discontinuous velocity dissipation assigned to the common edge or common surface of the element. Internal work includes the intra-unit double-point integral of the lower limit stress and the same upper limit plastic strain rate, and the work done by the lower limit interface force on the velocity discontinuity. The non-negative difference is the larger of the difference between the upper limit of plastic dissipation and internal energy after adopting the same allocation rule and unified caliber, and zero. The pairing process divides the contribution of the common edge or common face equally between two adjacent units, or according to the interface length, area, and integral weight conservation, and maintains the same allocation rule in all rounds.

3. The adaptive grid update method for geotechnical limit analysis according to claim 2, characterized in that: Before pairing, read the kinematic allowable residual, static equilibrium residual and maximum yield exceedance. If any one of them exceeds its allowable value, pairing is not performed and the solution is restarted from the most recent checkpoint. If none of them exceed their allowable values, compare the sum of the local original-dual power differences of each element with the global power difference converted from the upper and lower limit differences. When the difference between the two is within the integral and rounding tolerance, the local original-dual power difference is calibrated proportionally to the global power difference, and the step of classifying the upper limit plastic flow direction and the static lower limit predicted plastic flow direction into the same representation space continues; when it exceeds the tolerance, it is checked whether the external work normalization and the contribution of the common edge or common surface are repeatedly counted, and no repartitioning is performed.

4. The adaptive mesh update method for geotechnical limit analysis according to claim 2, characterized in that: The upper limit plastic flow direction is obtained by dividing the same upper limit plastic strain rate by the sum of its Frobenius norm and regularized positive number; the lower limit predicted plastic flow direction is obtained by projecting the lower limit stress to the nearest yield surface point and normalizing the gradient of the yield function at that point in the same way; the two directions are compared in the same representation space according to the unsigned axial relationship.

5. The adaptive mesh update method for geotechnical limit analysis according to claim 1, characterized in that: The fixed background space consists of fixed background partitions whose numbers remain unchanged during the repartitioning. The persistence index of each fixed background partition is updated recursively according to the overlap ratio between the candidate unstable region in this round and the fixed background partition, based on historical weights. The persistence index of the fixed background partition is mapped inversely to the persistence index of the finite element element by normalizing the weights based on the overlap ratio between the current finite element element and each fixed background partition. The finite element element comprehensive priority index is a normalized weighted sum of the local original-dual power difference, direction inconsistency index, sensitivity of the target limit quantity to the finite element element strength, and finite element element sustainability index; each weight is non-negative and the sum is 1, the weight of the finite element element sustainability index in the first round is 0, and the weight of the index in subsequent rounds is non-negative.

6. The adaptive mesh update method for geotechnical limit analysis according to claim 5, characterized in that: The edge cost of the unit adjacency graph is the product of the basic term, the direction penalty term, and the boundary penalty term; the basic term has the center distance between two adjacent finite element units as the numerator and the sum of the regularization positive number and the comprehensive priority index of the two finite element units as the denominator. The directional penalty term is 1 plus a directional penalty coefficient, a directional confidence combination value obtained from the directional confidence coefficients of two adjacent finite element elements according to a pre-fixed combination rule and limited to 0 to 1, the direction of the line connecting the centers of the two finite element elements, and the predicted instability tangent formed by unifying the signs of the upper limit plastic flow direction and the lower limit predicted plastic flow direction of the two adjacent finite element elements, and the product of the square of the sine of the difference between the two predicted instability tangents and the average direction formed by unifying the signs of the two predicted instability tangents according to a pre-fixed weighting rule. The boundary penalty is set to 1 for allowed passage interfaces and greater than 1 for required penalty boundaries; for prohibited passage interfaces, the corresponding graph edges are deleted without assigning a finite cost.

7. The adaptive mesh update method for geotechnical limit analysis according to claim 1, characterized in that: The search employs at least one of the following: shortest path, top K shortest paths, minimum cost flow, graph cut, and region growing. The results connect the boundary segments of the source and target ends, the cost is affected by the comprehensive priority index of the finite element elements, and can be converted into a finite-width unstable channel region that expands along the normal direction to a preset half-width. The candidate path is checked to ensure that its length is not less than the preset minimum value, the proportion of finite element elements with high priority in the current round is not less than the preset proportion, and it does not cross the prohibited crossing boundary. Those that pass the inspection are retained as continuous potential instability channels, and those that fail are removed; multiple channels with similar costs and spatial separation are retained as branch continuous potential instability channels, and those with spatial overlap are merged according to common segments.

8. The adaptive mesh update method for geotechnical limit analysis according to claim 7, characterized in that: If no through path is found connecting the source and target boundary segments, check the candidate cell set and prohibited boundary crossings of the boundary segment in sequence, expand the boundary buffer distance of the candidate cell set, increase the number of candidate paths and switch the search mode; if a through path is found, use it to determine the continuous potential instability channels and continue to perform the construction of anisotropic mesh measurement. If the result is still not obtained, restricted isotropic refinement is performed in the region with a high priority index of finite element synthesis, and the search is re-executed in the next round. The previous round of continuous potential instability channels are marked as historical references and prohibited from participating in the comparison.

9. The adaptive mesh update method for geotechnical limit analysis according to claim 8, characterized in that: When there is no complete lower bound interface force and dual variable in the pairing, the activity of the upper bound plastic dissipation and the proximity of the lower bound stress to the yield surface are normalized respectively, and the difference between the two is calculated as a simplified local original-dual difference index, and written into the simplified pairing state in the metadata. When common edge or common surface data is missing, the share is first restored by adjacent units according to the interface length or area; after restoration, if the sum of the local original-dual power difference passes the global power verification, repartitioning is performed; if it fails, the current round of repartitioning is locked. In the case of unrelated flow, the gradient of the plastic potential function at the yield surface is used to determine the lower limit of statics to predict the direction of plastic flow, and the results are labeled as upper limit kinematic results and lower limit static results.

10. The adaptive mesh update method for geotechnical limit analysis according to claim 5 or 7, characterized in that: The anisotropy ratio is the ratio of the target size along the tangential to the normal direction of the anisotropic mesh metric. The repartition is divided into two types: local repartition, which only deletes the finite-width unstable channel region and transition layer unit that extends along the normal direction of the continuous potential unstable channel, and global repartition, which is regenerated according to the anisotropic mesh metric. Both types preserve the mapping between the mesh before and after repartition and the fixed background space. After re-meshing, the mesh is checked for minimum angle and aspect ratio. If it passes, the solution is re-solved. If it fails, the local anisotropy ratio is reduced and a transition layer is added before regeneration. If it still fails, the re-meshing is canceled and the previous mesh is retained.

11. The adaptive mesh update method for geotechnical limit analysis according to claim 10, characterized in that: The output results include the re-meshed mesh, the limiting quantities of the upper kinematic limit and the lower static limit, the difference between the upper and lower limits, the centerline of the continuous potential instability channel and its finite-width instability channel region that extends along the normal direction to a preset half-width, the field distribution of the anisotropic mesh metric and the reason for stopping; the mesh version number, model consistency check value, local original-dual power difference, directional inconsistency index, overlap, connection relationship and number of branches of each round are written into the process traceability record; Both the output results and the process traceability records contain a cycle number and a checkpoint verification summary; if the receiving end fails to read, it will retransmit according to the cycle number, and if the reading is successful, the solution will not be executed repeatedly.

12. The adaptive mesh update method for geotechnical limit analysis according to claim 11, characterized in that: The same display interface overlays the current and previous rounds of continuous potential instability channels, which extend along the normal direction to a finite width instability channel region with a preset half width, as well as the intersection region of the two, the bidirectional geometric distance markings between the center lines of the two rounds of continuous potential instability channels, and the current values ​​of overlap, connection relationship and number of branches. When the overlap reaches the overlap threshold, the connection relationship and the number of branches meet the preset stability conditions, and the difference between the upper and lower limits and the change of the target limit quantity reach the corresponding threshold, the stop reason text of the output result is displayed on the interface. If any condition is not met, the interface will display the name of the unmet item and its value, as well as a status indicator indicating that a pairing has been returned.

13. The adaptive grid update method for geotechnical limit analysis according to claim 12, characterized in that: The finite-width instability channel region covers at least one layer of units on each side of the centerline, and one to three layers of transition units are set around it. The target size of the anisotropic mesh metric along the normal is truncated by preset upper and lower limits; the anisotropy ratio is not less than 1, and is compressed to 1 when the curvature of the centerline of a continuous potential instability channel is higher than the upper limit of curvature, or when there are adjacent geometric sharp corners or unqualified element quality; when none of the above conditions are met and the curvature of the centerline is small and the direction is stable, the anisotropy ratio is allowed to increase, and the anisotropy ratio does not exceed the upper limit allowed by the mesh generator; The output result is executed when the overlap reaches the overlap threshold, the connection relationship and the number of branches meet the preset stability conditions, the difference between the upper and lower limits and the change of the target limit quantity reach the corresponding threshold, and the above conditions are met for at least 2 consecutive rounds; otherwise, the pairing is returned.

14. The adaptive mesh update method for geotechnical limit analysis according to claim 10, characterized in that: Treat an update as a pending transaction and perform input and version checks, permission and feasibility checks, validity checks of continuous potential instability channels, mesh quality checks after repartitioning, field checks of projected results, and communication checks in sequence. If all comparisons pass, submit the current round and execute the output based on the comparison results or proceed to the next round; If the input and version check fails, resubmit the input; if the permission and feasibility check fails, perform a cold start or restore the most recent checkpoint; if the validity check of a continuous potential instability channel fails, perform restricted isotropic refinement; if the mesh quality check after remapping fails, reduce the anisotropy ratio or cancel the remapping; if the field check of the projected result fails, perform a cold start or cancel the remapping; if the communication check fails, restore the most recent checkpoint or continue transmitting the saved complete result. If the projected field fails to pass, the previous round of potential unstable channels stored in the fixed background space will not be updated.

15. An adaptive grid update system for geotechnical limit analysis, characterized in that, include: The upper limit kinematics solution module and the lower limit statics solution module receive the same mesh and solve the upper limit kinematics and lower limit statics respectively on the same mesh, and output the upper limit plastic dissipation, upper limit plastic strain rate and lower limit stress to the dual mechanism analysis module; The dual mechanism analysis module pairs the upper limit plastic dissipation and the lower limit stress with the internal work done by the upper limit plastic strain rate of the element for each element. The non-negative difference is taken as the local original-dual power difference. The upper limit plastic flow direction and the plastic flow direction predicted by the static lower limit determined by the lower limit stress at the yield surface are classified into the same representation space to obtain the direction inconsistency index. The local original-dual power difference and the direction inconsistency index are combined to form the finite element comprehensive priority index, and the finite element comprehensive priority index is output to the instability channel extraction module. The instability channel extraction module determines the edge cost of the element adjacency graph using the finite element element comprehensive priority index, searches for continuous potential instability channels in the element adjacency graph that connect the preset source end and target end boundary segments, and outputs the continuous potential instability channels to the fixed background mapping module and the anisotropic mesh update module respectively. The fixed background mapping module compares the overlap, connection relationship and number of branches of two adjacent rounds of continuous potential instability channels in a fixed background space that does not change with re-division, and outputs the comparison results to the constraint projection and stopping control module. The anisotropic mesh update module constructs anisotropic mesh metrics based on the tangent and normal directions of continuous potential instability channels and performs re-meshing accordingly, then outputs the re-meshed mesh to the constraint projection and stop control module. The constraint projection and stop control module takes the field obtained by projection from the field on the mesh before remapping as only the initial value, and schedules the kinematic upper limit solution module and the static lower limit solution module to resolve on the mesh after remapping; and triggers the result output module when the comparison result of the continuous potential instability channel meets the preset topological stability condition, the upper and lower limit difference between the kinematic upper limit and the static lower limit, and the change of the target limit quantity reaches the corresponding threshold; if any condition is not met, the result of this re-solution is sent to the dual mechanism analysis module to re-execute the pairing. The results output module outputs the re-divided mesh, continuous potential instability channels, upper and lower limit differences, and stopping reasons.

16. An electronic device, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, it implements the adaptive grid update method for geotechnical limit analysis as described in any one of claims 1 to 14.

17. A computer-readable storage medium, characterized in that, It stores a computer program, which, when executed by a processor, implements the adaptive grid update method for geotechnical limit analysis as described in any one of claims 1 to 14.

Images