An elasto-plastic numerical simulation method

Abstract

This invention discloses an elastoplastic numerical simulation method, comprising: acquiring input data of a numerical model of the geological body to be simulated; constructing elastoplastic incremental control equations based on the input data, and transforming the elastoplastic incremental control equations into a convex conical complementary problem model; solving the convex conical complementary problem model using a model basis matrix splitting iterative algorithm, which uses optimal iteration parameters for iterative calculation to obtain the displacement field, stress field, and plastic state field of the geological body to be simulated; and outputting the displacement field, stress field, and plastic state field as the elastoplastic numerical simulation results. By uniformly transforming the elastoplastic problem into a convex conical complementary model and theoretically deriving the optimal iteration parameters, the convergence efficiency and computational scale of the numerical solution are improved, solving the technical bottlenecks of poor convergence, high computational complexity, and difficulty in parameter tuning in traditional methods for ultra-large-scale elastoplastic simulations.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology, and in particular to a method for elastoplastic numerical simulation. Background Technology

[0002] In cutting-edge research on crustal dynamics, analyzing the stress evolution of complex fault systems is crucial for revealing the mechanisms of earthquake gestation and occurrence. This core problem ultimately boils down to solving a large-scale, highly nonlinear elastoplastic boundary value problem. Theoretically, this problem can be formulated as a quadratic programming problem with inequality constraints, but its theoretical analysis and numerical implementation are extremely challenging and have long remained unsolved in the field of computational geodynamics.

[0003] Currently, widely used traditional numerical methods are mainly based on the finite element framework, including the Goodman element method with penalty functions, Newton's method, and non-smooth Lagrange multiplier methods. These methods have clear mechanical concepts, relatively simple procedures, and are easy to implement. However, they generally suffer from inherent defects in convergence, especially for highly nonlinear problems, where convergence is often difficult to guarantee; their convergence accuracy also fails to meet the modern requirements for high-precision simulation of stress evolution in complex fault systems; furthermore, these methods are difficult to implement for ultra-large-scale computations. For example, the non-smooth Lagrange multiplier method involves multiplier calculations in a globally full matrix, resulting in a geometric increase in memory requirements, and convergence becomes extremely difficult as the problem size increases. Therefore, the computational geodynamics and computational mathematics communities have been exploring more efficient numerical algorithms.

[0004] From a computational mathematics perspective, the complementarity method offers a novel modeling approach with a rigorous mathematical foundation for such problems. This method expresses conditions such as contact and plastic flow as complementary relationships. Its algorithm boasts provable convergence and controllable accuracy, offering hope for overcoming the shortcomings of traditional methods. However, when applying this model to practical engineering problems, obtaining the system matrix and vectors required for the standard complementary form necessitates multiple inversions of the overall stiffness matrix, resulting in a computationally high complexity. O ( n ³). This enormous computational cost has long severely hampered the widespread application of complementary models in large-scale geophysical problems. Summary of the Invention

[0005] The purpose of this invention is to provide a numerical simulation method for elastoplasticity. By uniformly transforming the elastoplastic problem into a complementary convex cone model and theoretically deriving the optimal iterative parameters, the convergence efficiency and computational scale of numerical solutions are improved. This solves the technical bottlenecks of traditional methods in ultra-large-scale elastoplastic simulations, such as poor convergence, high computational complexity, and difficulty in parameter tuning.

[0006] To address the aforementioned technical problems, a first aspect of this invention provides a method for numerical simulation of elasticity and plasticity, comprising the following steps: Obtain the input data for the numerical model of the geological body to be simulated; Based on the input data, an elastoplastic incremental control equation is constructed, and the elastoplastic incremental control equation is transformed into a convex cone complementary problem model. The convex cone complementary problem model is a mathematical model that expresses the elastoplastic constitutive relation, yield criterion and equilibrium equation as complementary conditions under convex cone constraints. The model of the complementary convex cone problem is solved by the model basis matrix splitting iterative algorithm. The model basis matrix splitting iterative algorithm uses the optimal iteration parameters to perform iterative calculations and obtain the displacement field, stress field and plastic state field of the geological body to be simulated. The optimal iteration parameters are determined by making the implicit iteration equation of the model basis matrix splitting iterative algorithm explicit and minimizing the norm of the resulting explicit iteration matrix. The displacement field, stress field, and plastic state field are output as the results of the elastoplastic numerical simulation.

[0007] Furthermore, the input data includes discretized mesh data, mechanical boundary condition data, material elastoplastic constitutive parameters, and initial stress field and displacement field data; The step of constructing an elastoplastic incremental control equation based on the input data, and transforming the elastoplastic incremental control equation into a complementary convex cone problem model, includes: Based on the discretized mesh data and the material's elastoplastic constitutive parameters, an elastoplastic constitutive complementary system that couples the yield function, the plastic multiplier increment, and their complementary conditions is constructed within the incremental step. According to the yield criterion, the constraints related to the yield function in the elastoplastic constitutive complementary system are equivalently mapped to second-order cone constraints, which are convex cone sets defined on the vector space composed of stress tensor invariants. Based on the discretized grid data, the mechanical boundary condition data, and the initial stress field and displacement field data, the overall equilibrium equation is assembled, and the overall equilibrium equation, the second-order cone constraint, and the elastoplastic constitutive complementary system are combined to construct the convex cone complementary problem model.

[0008] Further, the step of equivalently mapping the constraints related to the yield function in the elastoplastic constitutive complementary system to second-order cone constraints according to the yield criterion includes: Based on the yield criterion, determine the coefficients of the linear term related to hydrostatic pressure and the norm term related to deviatoric stress in the yield function; Based on the linear term coefficients and the norm term coefficients, a linear inequality is constructed consisting of the hydrostatic pressure and the Euclidean norm of the deviatoric stress tensor. The linear inequality is identified as a second-order cone constraint with the hydrostatic pressure and deviatoric stress tensor components as variables.

[0009] Furthermore, the assembly of the overall equilibrium equations, and the simultaneous establishment of the overall equilibrium equations, the second-order cone constraint, and the elastoplastic constitutive complementary system, to construct the convex cone complementarity problem model, includes: Substituting the constitutive relations in the elastic-plastic constitutive complementary system into the global equilibrium equations, we obtain a set of governing equations with nodal displacement increments and plastic multiplier increments as the basic unknowns. By combining the governing equations, the second-order cone constraint, and the complementary conditions in the elastoplastic constitutive complementary system, a hybrid complementary equation set with the nodal displacement increment and the plastic multiplier increment as unknown variables is constructed as the model for the convex cone complementary problem.

[0010] Furthermore, the step of solving the complementary convex cone problem model using the modular basis matrix splitting iterative algorithm includes: Initialize the iteration parameters, initial iteration vector, and convergence tolerance of the modular basis matrix splitting iteration algorithm, and set its coefficient matrix and right-hand vector based on the convex cone complementarity problem model, and set the initial iteration vector as the current iteration vector; In each iteration, the absolute value vector of the current iteration vector is calculated. The absolute value vector is a vector formed by calculating the magnitude of each sub-vector of the current iteration vector corresponding to the second-order cone constraint. Based on the iteration parameters and the absolute value vector, construct intermediate auxiliary variable pairs that satisfy the convex cone constraint and the complementary condition; Determine whether the norm of the residual vector formed by the intermediate auxiliary variable pair is less than or equal to the convergence tolerance, wherein the residual vector is used to measure the degree to which the intermediate auxiliary variable pair deviates from the linear relationship defined by the convex cone complementarity problem model; If the norm of the residual vector is greater than the convergence tolerance, then based on the splitting form of the coefficient matrix, the iteration parameters, the current iteration vector, the absolute value vector, and the right-hand vector, a system of linear equations is constructed and solved. The solution vector obtained is used as the new current iteration vector, and the step of calculating the absolute value vector of the current iteration vector is returned. Otherwise, the intermediate auxiliary variable pair is output as the solution to the convex cone complementarity problem model.

[0011] Furthermore, the construction and solution of the linear equation system includes: Based on the inherent properties of the coefficient matrix of the complementary convex cone problem model, a matrix splitting method is selected so that the system matrix of the constructed linear equation system is a diagonal matrix or a triangular matrix. The inherent properties of the coefficient matrix are determined by the input data and the mathematical structure of the complementary convex cone problem model. Taking advantage of the structural characteristic that the system matrix is ​​a diagonal or triangular matrix, the linear equation system is solved using the forward substitution method or the backward substitution method.

[0012] Furthermore, the selection of matrix splitting methods based on the inherent properties of the coefficient matrix of the complementary convex cone problem model includes: Identify the sparse structure and non-zero element distribution pattern of the coefficient matrix; Based on the sparse structure and distribution pattern, the coefficient matrix is ​​split into a diagonal part and a non-diagonal part, or into a lower triangular part, a diagonal part and an upper triangular part, so that the system matrix is ​​a diagonal matrix or a triangular matrix.

[0013] Furthermore, the construction and solution of the linear equation system includes: The coefficient matrix of the model for the complementary convex cone problem is identified as a symmetric positive definite matrix. Based on the symmetric positive definiteness, the coefficient matrix is ​​split into a matrix splitting method of symmetric positive definite part and remainder part to construct the linear equation system; Solve the system of linear equations to obtain the solution vector required for updating the iteration vector.

[0014] Furthermore, before performing iterative calculations using the optimal iteration parameters, the modular basis matrix splitting iterative algorithm further includes: For the current iteration vector, a mapping matrix is ​​determined according to a predefined rule based on the geometric region in which it is located within the second-order cone constraint. The predefined rule is as follows: when the current iteration vector is located inside or on the boundary of the second-order cone constraint, the mapping matrix is ​​an identity matrix or a variant thereof; when the current iteration vector is located outside the second-order cone constraint, the mapping matrix is ​​a negative identity matrix or a variant thereof. Using the mapping matrix, the absolute value vector in the modular basis matrix splitting iteration algorithm is represented as the product of the mapping matrix and the current iteration vector, thus transforming the implicit iteration equation into an explicit iteration form; Based on the explicit iteration form, an explicit iteration matrix is ​​constructed with the iteration parameters as variables; The optimal iteration parameters are obtained by minimizing the matrix norm of the explicit iteration matrix.

[0015] Further, the step of obtaining the optimal iteration parameters by minimizing the matrix norm of the explicit iteration matrix includes: The square of the matrix norm of the explicit iteration matrix is ​​expressed as a quadratic function of the iteration parameters; By solving the condition for the quadratic function to reach its minimum value, a closed-form expression for the coefficient matrix and the trace of the mapping matrix of the optimal iterative parameters with respect to the complementary convex cone problem model is obtained. The optimal iteration parameters are calculated based on the closed-form expression.

[0016] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described elastoplastic numerical simulation method.

[0017] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described elastoplastic numerical simulation method.

[0018] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By systematically transforming the elastoplastic incremental governing equations into a convex cone complementary problem model, a unified mathematical expression for multiphysics constraints is achieved. The traditionally separate elastoplastic constitutive relations, yield criteria, and global equilibrium equations are organically integrated into a rigorous mathematical framework based on convex cones and complementary conditions. This modeling approach not only fundamentally avoids the convergence difficulties and accuracy losses commonly encountered in traditional finite element methods when dealing with inequality constraints, but also provides a high-fidelity numerical description of the nonlinear mechanical behavior of complex geological bodies, significantly improving the physical reliability and numerical stability of the simulation results. 2. By theoretically deriving and explicitly determining the optimal iterative parameters of the modular basis matrix splitting iterative algorithm, the huge computational burden of parameter tuning in traditional trial-and-error methods is completely eliminated. The implicit iterative equation is made explicit by introducing a mapping function, and the closed-form expression of the optimal parameters is directly derived by minimizing the norm of the explicit iterative matrix. This enables the algorithm to efficiently and automatically obtain high-performance parameters before each solution, avoiding the problem of slow convergence or even failure due to improper parameters. The parameter optimization cost, which may have originally been dominant, is reduced to almost zero, which greatly improves the efficiency and reliability of the overall solution process. 3. By designing targeted solution strategies for ultra-large-scale sparse systems, a fundamental reduction in computational complexity has been achieved, from cubic to quadratic and even linear levels. Fully utilizing the sparsity and structural characteristics of the system matrix after discretization of the elasto-plastic convex conical complementary model, and combining efficient numerical techniques such as matrix splitting, back-splitting, and Krylov subspace iteration, the solution process can fully leverage the advantages of parallel computing and significantly reduce memory consumption. This makes it possible to perform full-size, high-resolution elasto-plastic simulations of geological bodies containing massive numbers of elements and nodes, effectively solving the long-standing bottleneck of traditional methods being unable to be applied to practical large-scale engineering and scientific problems due to computational resource limitations. Attached Figure Description

[0019] Figure 1 This is a flowchart of the elastoplastic numerical simulation method provided in the embodiments of the present invention; Figure 2 This is a schematic diagram of the elastoplastic numerical simulation logic provided in an embodiment of the present invention; Figure 3a This is a schematic diagram illustrating the time cost of using the cone complementarity algorithm provided in this embodiment of the invention to solve the frictional contact problem; Figure 3b This is a schematic diagram illustrating the time cost of using the cone complementarity algorithm provided in this embodiment of the invention to solve elastoplastic problems; Figure 4a This is a schematic diagram illustrating the number of iteration steps when using the cone complementarity algorithm provided in this embodiment of the invention to solve the frictional contact problem; Figure 4b This is a schematic diagram showing the number of iteration steps when using the cone complementarity algorithm provided in this embodiment of the invention to solve elastoplastic problems. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0021] Please refer to Figure 1 The first aspect of this invention provides a method for numerical simulation of elastic-plastic properties, comprising the following steps: Step 100: Obtain the input data for the numerical model of the geological body to be simulated.

[0022] For regional stress evolution simulations involving complex fault systems, the input data mainly includes: discretized grid data describing the geometry and internal structure of geological bodies, usually obtained through 3D geological modeling and mesh generation; mechanical boundary condition data reflecting tectonic loads and constraints, such as displacement constraints or force loads applied to the model boundaries; material elastoplastic constitutive parameters characterizing the mechanical response of soil and rock media, such as Young's model system, Poisson's ratio, and cohesion and internal friction angle in the Drucker-Prager or Mohr-Coulomb yield criteria determined based on laboratory tests; and initial stress field and displacement field data simulating the initial stress state, which are often provided through geostress inversion or previous simulation results.

[0023] Step 200: Based on the input data, construct the elastoplastic incremental control equation and transform it into a convex cone complementary problem model. The convex cone complementary problem model is a mathematical model that unifies the elastoplastic constitutive relation, yield criterion and equilibrium equation into complementary conditions under convex cone constraints.

[0024] Based on the above input data, the governing equations of the elastoplastic mechanics problem are established within the time increment step and transformed into a convex conical complementary form. First, by combining the discretized mesh and constitutive parameters, an elastoplastic constitutive complementary system is constructed, coupling the yield function, the plasticity multiplier increment, and the complementarity conditions. Then, according to the selected yield criterion, the yield-related constraints in the system are equivalently transformed into second-order conical constraints—constraints that define a convex cone in the vector space spanned by the hydrostatic pressure and the deviatoric stress tensor norm. Finally, equilibrium equations are introduced, and the overall equilibrium relationship is obtained through finite element integration. This relationship is then combined with the constitutive complementary system and the second-order conical constraints to form a unified convex conical complementary problem model. This model transforms the elastoplastic boundary value problem with inequality constraints into a mathematical programming problem that can be rigorously described within the convex conical framework.

[0025] Step 300: The model of the complementary convex cone problem is solved by using the model basis matrix splitting iterative algorithm. The model basis matrix splitting iterative algorithm uses the optimal iteration parameters to perform iterative calculations and obtain the displacement field, stress field and plastic state field of the geological body to be simulated. The optimal iteration parameters are determined by making the implicit iteration equation of the model basis matrix splitting iterative algorithm explicit and minimizing the norm of the explicit iteration matrix obtained therefrom.

[0026] Before the iteration begins, based on the mathematical structure of the convex cone complementarity problem model, the implicit iterative equations of the algorithm are made explicit by introducing a mapping function. The norm of the iteration matrix is ​​minimized based on this explicit form, thus theoretically deriving the optimal iterative parameters and avoiding the computational cost of traditional trial-and-error parameter tuning. During the solution process, the algorithm initializes the iteration vector and convergence tolerance. In each iteration, the absolute value vector of the current vector is calculated, auxiliary variable pairs satisfying the cone constraint and complementarity conditions are constructed, and their linear relationship as defined in the complementarity problem is checked. If not, a system of linear equations is constructed and solved based on the splitting form of the coefficient matrix, the optimal parameters, and other iterative information. The iteration vector is updated until the residuals satisfy the convergence condition, and finally, the displacement field, stress field, and plastic state field are output.

[0027] Step 400: Output the displacement field, stress field, and plastic state field as the results of the elastoplastic numerical simulation.

[0028] The system outputs key mechanical field variables obtained from numerical simulations, including displacement, stress, and plastic state fields that indicate plastic deformation. These results are output in structured data format and can be directly used for subsequent analyses. For example, visualization tools can be used to generate deformation contour maps, stress distribution maps, and plastic zone expansion maps of geological bodies, providing quantitative data for assessing fault activity, regional stability, or seismic hazard. The output data can also serve as input for other coupled simulations (such as seepage-stress coupling), supporting multiphysics integrated analysis.

[0029] This invention establishes a convex conical complementary unified model for elastoplastic problems, incorporating complex mechanical constraints into a rigorous mathematical programming framework. This significantly improves the convergence and numerical stability of traditional methods when dealing with nonlinear and inequality constraints. By theoretically deriving the optimal iterative parameters, the invention avoids the blindness and computational burden of algorithm parameter tuning, thereby improving solution efficiency and reliability. Furthermore, by employing matrix splitting and sparse solution techniques, the invention drastically reduces computational and storage complexity, making large-scale and even ultra-large-scale elastoplastic dynamic simulation of geological bodies feasible. This provides an effective tool for refined numerical analysis in the fields of geophysics and engineering geology.

[0030] Specifically, in one embodiment of the present invention, the input data includes discretized mesh data, mechanical boundary condition data, material elastoplastic constitutive parameters, and initial stress field and displacement field data.

[0031] The aforementioned discretized mesh data describes the geometry and spatial discrete structure of the geological body to be simulated, typically generated based on 3D geological modeling and exploration data (such as seismic profiles and borehole data). Taking the simulation of a regional stress field containing faults as an example, the mesh needs to accurately characterize the geometry, attitude, and contact relationship with the surrounding rock of the fault plane, and the model domain is partitioned using hexahedral or tetrahedral elements. This data provides the computational basis for finite element analysis, determining the element node coordinates, element topology, and material grouping information, and is the foundation for subsequent structural equilibrium equations and numerical integration. The density and quality of the mesh directly affect the computational accuracy and efficiency; local refinement is required in key areas (such as near faults). Mechanical boundary condition data defines the kinematic or dynamic constraints of the model boundaries, reflecting the effect of the regional tectonic environment on the simulated area. For example, when studying the stress accumulation process of a strike-slip fault system, boundary conditions can be based on the relative plate motion velocity or GPS observation data, applying uniform displacement constraints to the lateral boundaries of the model, applying fixed normal constraints to the bottom boundaries, and treating the surface as a free surface. These data transform macroscopic structural drivers into computable boundary value problems, serving as external excitation sources controlling the mechanical response of the simulated system. Their proper setting is crucial to the realism of the simulation results. Material elastoplastic constitutive parameters quantitatively describe the mechanical behavior of geological media (such as rocks and soils), typically obtained through indoor triaxial tests, true triaxial tests, or field testing. Taking the simulation of brittle-plastic rock behavior as an example, the parameters include Young's model system and Poisson's ratio in the elastic stage, and cohesion, internal friction angle, and dilatation angle in the plastic stage determined according to the Drucker-Prager or Mohr-Coulomb yield criteria. These parameters are assigned to the corresponding material groups in the mesh to calculate stress-strain relationships, determine the yield state, and update plastic internal variables within incremental steps, representing the intrinsic physical properties controlling the nonlinear mechanical response of the medium. Initial stress and displacement field data define the mechanical state of the geological body at the start of the simulation, serving as the historical starting point for incremental loading or evolution analysis. For regional-scale simulations, the initial stress field can be derived from interpolation of geostress measurement data, theoretical calculations based on self-weight and tectonic background, or the output of previous tectonic evolution simulations; the initial displacement field is often set to zero or inherited from previous simulations. This data provides the initial conditions of the mechanical system, ensuring that the simulation evolves from a physically reasonable equilibrium or quasi-equilibrium state, and has a significant impact on the accurate simulation of stress redistribution and plastic development processes.

[0032] Accordingly, step S200 involves constructing an elastoplastic incremental control equation based on the input data, and transforming the elastoplastic incremental control equation into a convex cone complementary problem model, including: Step S210: Based on the discretized mesh data and the material's elastoplastic constitutive parameters, construct an elastoplastic constitutive complementary system that couples the yield function, the plastic multiplier increment, and the complementary conditions of the two within the incremental step.

[0033] Within each time increment step, based on the elements and integration points defined by the discretized mesh data, and combined with the material's elastoplastic constitutive parameters (such as the internal friction angle and cohesion in the Drucker-Prager criterion), a constitutive complementary system describing the elastoplastic behavior of material points is constructed. This system couples the yield function, the plasticity multiplier increment, and the complementary conditions between them (i.e., the yield function is less than or equal to zero, the plasticity multiplier increment is greater than or equal to zero, and their product is zero). For each material point, based on the incremental stress-strain relationship, the stress update is expressed as the sum of the elastic prediction and the plasticity correction, where the direction of the plasticity correction is determined by the gradient of the yield function. Through element shape functions and numerical integration, the local complementary systems of all material points are integrated into a global elastoplastic constitutive complementary system, which fully describes the nonlinear relationship between stress, strain, and plasticity intrinsic variables in the discrete model and implicitly satisfies the plasticity irreversibility condition.

[0034] Step S220: According to the yield criterion, the constraints related to the yield function in the elastoplastic constitutive complementary system are equivalently mapped to second-order cone constraints. Second-order cone constraints are a set of convex cones defined on the vector space composed of stress tensor invariants.

[0035] Based on the selected yield criterion (e.g., the Drucker-Prager criterion), the constraints related to the yield function within the elastoplastic constitutive complementary system in step S210 are transformed mathematically to an equivalent mathematical form. The yield function is typically expressed as a combination of the first invariant (hydrostatic pressure) and the second invariant (deviatoric stress tensor norm) of the stress tensor. Through algebraic operations, the yield inequality is reformulated as a linear inequality between hydrostatic pressure and the deviatoric stress norm. This inequality defines a second-order cone (also known as a Lorentz cone or ice cream cone) in the vector space composed of the components of the hydrostatic pressure and deviatoric stress tensors, in the form of a standard second-order cone constraint: the first component of a point within the cone is not less than the Euclidean norm formed by the subsequent components. This mapping transforms the complex yield surface constraints into standard cone constraints in convex optimization, providing a unified mathematical form for applying convex cone complementarity theory.

[0036] Step S230: Based on the discretized grid data, mechanical boundary condition data, and initial stress field and displacement field data, assemble the overall equilibrium equations, and combine the overall equilibrium equations, second-order cone constraints, and elastoplastic constitutive complementary system to construct a convex cone complementary problem model.

[0037] Based on discretized mesh data, mechanical boundary condition data, and initial stress and displacement field data, the global equilibrium equations are assembled using the finite element method. The specific process includes: calculating the element stiffness matrix and equivalent nodal force vectors according to the element type and shape function; and integrating the global stiffness matrix and global load vector according to the node number to form a system of linear equations with nodal displacement increments as the basic unknowns. Subsequently, the elastoplastic constitutive complementary system established in step S210 (where stress is expressed as a function of strain through constitutive relations and thus related to displacement increments) is substituted into the equilibrium equations, while the second-order cone constraint defined in step S220 is introduced. The global equilibrium equations, the second-order cone constraint, and the complementary conditions in the constitutive complementary system are combined to form a hybrid complementary problem with nodal displacement increments and plastic multiplier increments as unknowns. This problem has a complementary structure under convex cone constraints, which is the constructed convex cone complementary problem model, transforming the original elastoplastic boundary value problem into a standard form that can be efficiently solved within a mathematical programming framework.

[0038] Further, in step S220, according to the yield criterion, the constraints related to the yield function in the elastoplastic constitutive complementary system are equivalently mapped to second-order cone constraints, including: Step S221: Determine the coefficients of the linear term related to hydrostatic pressure and the norm term related to deviatoric stress in the yield function according to the yield criterion.

[0039] For a selected yield criterion (e.g., the Drucker-Prager criterion applicable to soil and rock materials), the linear term coefficients related to the material's internal friction properties and the norm term coefficients related to the material's shear strength are determined based on its mathematical expression. Taking the Drucker-Prager criterion as an example, its yield function can be expressed as a linear combination of hydrostatic pressure and the deviatoric stress norm. The linear term coefficients (reflecting the weight of hydrostatic pressure on yield) and the norm term coefficients (usually 1) can be calculated from the internal friction angle and cohesion parameters obtained through material testing. When simulating the plastic behavior of rocks in fault zones, these coefficients reflect the strengthening effect of confining pressure on yield strength and the critical condition for shear failure, serving as a crucial bridge connecting the material's physical properties and mathematical constraints.

[0040] Step S222: Based on the coefficients of the linear term and the coefficients of the norm term, construct a linear inequality consisting of the Euclidean norm of the hydrostatic pressure and the deviatoric stress tensor.

[0041] Using the linear term coefficients and norm term coefficients determined in step S221, a linear inequality is constructed with hydrostatic pressure and the Euclidean norm of the deviatoric stress tensor as variables. The hydrostatic pressure is calculated from the first invariant of the stress tensor, characterizing the spherical stress state; the deviatoric stress tensor norm is calculated from the second invariant of the stress deviator, characterizing the shear stress intensity. The constructed linear inequality is in the form that the linear term coefficient multiplied by the hydrostatic pressure plus the deviatoric stress norm is less than or equal to a constant related to the material strength. This inequality defines a half-space in the stress invariant space, the boundary of which is the yield surface. When the stress point is located within this half-space, the material is in an elastic or rigid state; when it is located at the boundary, plastic flow begins.

[0042] Step S223 identifies the linear inequality as a second-order cone constraint with hydrostatic pressure and deviatoric stress tensor components as variables.

[0043] The linear inequality constructed in step S222 is reformulated as a standard second-order cone constraint. Specifically, the inequality is rewritten as the deviatoric stress tensor norm being less than or equal to a linear expression concerning hydrostatic pressure. By introducing an auxiliary variable, this inequality can be equivalently expressed as: the Euclidean norm of a vector composed of the components of the deviatoric stress tensor does not exceed another scalar variable, and this scalar variable is linearly related to the hydrostatic pressure. This precisely satisfies the definition of a second-order cone—the first component of a point is not less than the magnitude of the vector composed of its subsequent components. Therefore, the original yield inequality is mapped to a second-order cone constraint in an extended vector space composed of hydrostatic pressure and deviatoric stress components, thus transforming the physical yield condition into a standard geometric constraint in convex optimization.

[0044] Further, in step S230, the overall equilibrium equations are assembled, and the overall equilibrium equations, the second-order cone constraint, and the elastoplastic constitutive complementary system are combined to construct a convex cone complementarity problem model, including: Step S231: Substitute the constitutive relations in the elastoplastic constitutive complementary system into the global equilibrium equations to obtain a set of governing equations with nodal displacement increments and plastic multiplier increments as the basic unknowns.

[0045] Substituting the incremental constitutive relations of the elastoplastic constitutive complementary system into the assembled global equilibrium equations yields a set of coupled governing equations with the total nodal displacement increment and the plastic multiplier increment at all integration points as the basic unknowns. Specifically, the global equilibrium equations are essentially a discretized virtual work principle, expressing that the work done by internal and external forces at any virtual displacement is equal. Internal forces are obtained by integrating the element stress field, and through the constitutive complementary system, the stress increment at each integration point can be expressed as a function of the strain increment (obtained from the nodal displacement increment through the strain-displacement matrix) and the plastic multiplier increment. Substituting this relationship into the internal force calculation formula, the original equilibrium equation with displacement as the single unknown is extended into a set of algebraic equations simultaneously containing the nodal displacement increment and the plastic multiplier increment. This set of governing equations comprehensively reflects the force balance relationship of the discrete system and the elastoplastic constitutive response at all material points.

[0046] Step S232: Combine the governing equations, the second-order cone constraint, and the complementary conditions in the elastoplastic constitutive complementary system to construct a hybrid complementary equation set with nodal displacement increments and plastic multiplier increments as unknown variables, which serves as the model for the convex cone complementary problem.

[0047] By combining the coupled control equations obtained in step S231, the second-order cone constraint defined in step S220, and the inherent complementarity conditions in the elastoplastic constitutive complementary system, the final convex cone complementarity problem model is constructed. The control equations provide the linear equality relationship that must be satisfied between the system variables (displacement increment and plastic multiplier increment). The second-order cone constraint, in the form of inequalities, restricts the stress state (which, in turn, relates displacement and plastic multiplier through constitutive relations) to be within the allowable domain defined by the yield criterion. The constitutive complementarity condition requires that, at each integration point, the plastic multiplier increment and the yield function value satisfy the complementary relationship of "non-negativity, complementarity, and orthogonality". Combining these three types of mathematical objects—equality, cone inequality, and complementarity condition—as a whole forms a standard hybrid complementarity problem. The unknowns of this problem are the global nodal displacement increment vector and the plastic multiplier increment vector. The constraints consist of linear equality, second-order cone constraint, and point-to-point complementarity condition, thus completing the construction of a convex cone complementarity problem model that can be directly solved within the framework of mathematical programming, transforming a physical mechanics problem into one that can be directly solved within the framework of mathematical programming.

[0048] Specifically, step S300 involves solving the complementary convex cone problem model using the modular basis matrix splitting iterative algorithm, including: Step S310: Initialize the iteration parameters, initial iteration vector, and convergence tolerance of the modular basis matrix splitting iteration algorithm, and set its coefficient matrix and right-hand vector based on the convex cone complementary problem model, and set the initial iteration vector as the current iteration vector.

[0049] First, the key operating parameters of the algorithm are set: iteration parameters (usually including relaxation and scaling factors), initial iteration vector (often set as the zero vector or an initial solution based on physical predictions), and convergence tolerance (a preset, minimal positive number used to determine computational accuracy). Simultaneously, the coefficient matrix and right-hand side vector corresponding to the standard form of the constructed convex conical complementarity problem model (i.e., finding vector pairs that satisfy conical constraints, linear relationships, and orthogonality) are extracted. Finally, the initial iteration vector is designated as the current iteration vector, preparing for the iteration loop. This initialization process sets a definite starting point and termination criterion for the algorithm and establishes a direct connection with the mathematical model to be solved.

[0050] Step S320: In each iteration, calculate the absolute value vector of the current iteration vector. The absolute value vector is the vector formed by calculating the magnitude of each sub-vector in the current iteration vector corresponding to the second-order cone constraint.

[0051] In each iteration, a crucial preprocessing step is performed on the current iteration vector: calculating its absolute value vector. This operation is not simply taking the absolute value of each component of the vector, but rather performing a specialized calculation on the block sub-vectors in the current iteration vector corresponding to each second-order cone constraint. For each sub-vector, based on its spectral decomposition in the corresponding second-order cone, the absolute value of its spectral value is calculated and reassembled according to the original spectral decomposition structure. The resulting vector is called the modulus representation of that sub-vector. Concatenating the modulus representations of all sub-vectors in sequence yields the complete absolute value vector. This calculation utilizes Jordan's algebra on the second-order cone and is one of the core operations of the modular basis algorithm that implicitly encodes cone constraints and complementarity conditions into the iterative format.

[0052] Step S330: Based on the iteration parameters and the absolute value vector, construct intermediate auxiliary variable pairs that satisfy the convex cone constraint and the complementary condition.

[0053] Using the current iteration parameters (denoted as γ and ω) and the absolute value vector calculated in step S320, a pair of intermediate auxiliary variables (denoted as z and w) are constructed. The construction formulas are: z = (absolute value vector + current iteration vector) / γ, w = ω * (absolute value vector - current iteration vector) / γ. Based on the mathematical properties of second-order cones and complementarity theory, the variable pair (z, w) constructed in this way has the following characteristics: z and w are both automatically located within their respective second-order cones and satisfy the positive complementary condition (i.e., the dot product of z and w is zero). Therefore, this variable pair has satisfied the two core requirements of the convex cone complementarity problem model regarding "within the cone" and "mutually orthogonal complementarity". The difference between it and the final solution lies only in whether the linear relationship w = A defined in the model is satisfied. z + b.

[0054] Step S340: Determine whether the norm of the residual vector formed by the intermediate auxiliary variable pairs is less than or equal to the convergence tolerance. The residual vector is used to measure the degree to which the intermediate auxiliary variable pairs deviate from the linear relationship defined by the complementary convex cone problem model.

[0055] Calculate the residual vector formed by the intermediate auxiliary variable pair (z, w) constructed in step S330: r = w - (A Let (z, w) = (z + b), where A and b are the coefficient matrix and right-hand side vector of the convex cone complementarity problem model, respectively. The residual vector r quantitatively characterizes the degree to which the current auxiliary variable deviates from the linear relationship required by the model. A norm (such as the 2-norm) of the residual vector r is calculated and compared with a pre-defined convergence tolerance. If the norm is less than or equal to the convergence tolerance, the currently constructed (z, w) is considered to sufficiently satisfy all the conditions of the convex cone complementarity problem, and the iteration converges; otherwise, the iteration needs to continue.

[0056] Step S350: If the norm of the residual vector is greater than the convergence tolerance, then construct and solve a system of linear equations based on the split form of the coefficient matrix, the iteration parameters, the current iteration vector, the absolute value vector, and the right-hand side vector. Use the solution vector obtained as the new current iteration vector and return to the step of calculating the absolute value vector of the current iteration vector. Otherwise, output the intermediate auxiliary variable pair as the solution to the convex cone complementarity problem model.

[0057] If the convergence condition is not met, iterative updates are performed. Based on a pre-defined split of the coefficient matrix A (e.g., A = M - N, where M is an invertible matrix), the iteration parameter ω, the current iteration vector x, the absolute value vector, and the right-hand vector b, a system of linear equations is constructed according to the fixed format of the modular basis algorithm: (M + ωI). x new = N x + (ωI - A) * (absolute value vector) - γ b, where I represents the identity matrix of the same order as the coefficient matrix A. The system matrix (M + ωI) of this system of equations usually possesses favorable properties (such as symmetric positive definiteness, diagonal or triangular dominance) and can be solved using efficient linear solvers (such as the forward-backward substitution method and the conjugate gradient method). The obtained solution vector x new Use this as the current iteration vector for the next iteration and return to step S320 to start a new round of iteration. If the convergence condition has been met, directly output the current intermediate auxiliary variable pair (z, w) as the final solution of the convex cone complementary problem model, from which the required displacement field, stress field and other information can be extracted.

[0058] In one embodiment of the present invention, the construction and solution of the linear equation system in step S350 includes: Step S3511: Based on the inherent properties of the coefficient matrix of the complementary convex conic problem model, select a matrix splitting method so that the system matrix of the constructed linear equation system is a diagonal matrix or a triangular matrix. The inherent properties of the coefficient matrix are determined by the input data and the mathematical structure of the complementary convex conic problem model.

[0059] Based on the inherent mathematical properties of the coefficient matrix of the convex conic complementarity problem model, a matrix splitting method is selected to give the system matrix of the final constructed iterative linear equation system a special structure of diagonal or triangular matrices. The inherent properties of the coefficient matrix, such as its symmetry, positive definiteness, sparsity pattern (e.g., banded or blocky structure), and diagonal dominance, are jointly determined by the original input data (e.g., mesh connectivity, material parameter distribution) and the mathematical construction of the convex conic complementarity problem model itself (e.g., the coupling form of equilibrium equations and conic constraints). When simulating large-scale geological bodies, model discretization often produces large sparse matrices. Based on the analysis of the matrix's non-zero element distribution and eigenvalue distribution, classic splitting strategies such as Jacobi splitting (making the split matrix M a diagonal matrix) or Gauss-Seidel splitting (making M a lower triangular matrix) can be selected. This selection process decomposes the original coefficient matrix A into two parts, M and N (A = M - N), with the goal of ensuring that the system matrix of the iterative equation (M + ωI)x new =..., the system matrix (M + ωI) also presents a diagonal or triangular form due to the structure of M, thus laying the foundation for the subsequent highly efficient direct solution.

[0060] Step S3512: Utilize the structural characteristic that the system matrix is ​​a diagonal or triangular matrix, and solve the linear equation system using the forward substitution method or the backward substitution method.

[0061] For the system matrix with a diagonal or triangular structure obtained in step S3511, a highly optimized forward substitution method or backward substitution method is used for solution. If the system matrix is ​​a diagonal matrix, its solution degenerates into a simple operation of dividing each component sequentially by the diagonal element. If it is a (lower) triangular matrix, the forward substitution method (for lower triangular matrices) or backward substitution method (for upper triangular matrices) is used, and the solution can be completed through a series of sequential scalar operations involving only a single unknown, without the need for complex matrix decomposition or iteration. These methods fully utilize the special characteristics of the matrix structure, reducing the computational complexity of the solution from that of a general linear solver. O ( n ³) Significantly reduced to O ( n )or O ( n ²). In implementation, this process is highly vectorizable and has limited but very regular parallelism, making it particularly suitable for modern processor architectures. It can complete the core solution tasks in the iterative steps with extremely low computational overhead, and is a key link in achieving efficient computation for dealing with ultra-large-scale problems.

[0062] Furthermore, the inherent properties of the coefficient matrix in step S3511 based on the complementary convex cone problem model, and the selection of the matrix splitting method, include: Step S35111: Identify the sparse structure and non-zero element distribution pattern of the coefficient matrix.

[0063] This paper identifies and analyzes the sparse structure and non-zero element distribution patterns of the coefficient matrix in a model of complementary convex cones. The sparse structure of the coefficient matrix reflects the connectivity between nodes in the discretized mesh. For example, in finite element discretization, non-zero elements typically appear at positions representing coupling between adjacent nodes, forming specific sparse patterns (such as banded, blocky, or more complex irregular structures). The distribution patterns of non-zero elements further describe the aggregation characteristics of these non-zero elements in the matrix, such as whether they exhibit a significant diagonal dominance, or whether they have a blocky diagonal or blocky triangular structure. This identification process can be completed by analyzing the graphical representation of the matrix (such as an adjacency graph) or by directly scanning the non-zero element indices, aiming to provide an objective basis for subsequently selecting the most suitable splitting method. When simulating large-scale geological bodies, due to the large mesh size and potential local refinement, the coefficient matrix usually has the characteristics of high sparseness but uneven distribution of non-zero elements. Accurately identifying these characteristics is a prerequisite for achieving efficient splitting.

[0064] Step S35112: Based on the sparse structure and distribution pattern, the coefficient matrix is ​​split into a diagonal part and a non-diagonal part, or into a matrix splitting method of a lower triangular part, a diagonal part and an upper triangular part, so that the system matrix is ​​a diagonal matrix or a triangular matrix.

[0065] Based on the sparse structure and distribution pattern identified in step S35111, a specific matrix splitting operation is performed. If the non-zero elements of the coefficient matrix are mainly concentrated near the diagonal and the off-diagonal elements are relatively dispersed, the matrix is ​​split into a diagonal part and an off-diagonal part, i.e., Jacobi splitting. In this method, all diagonal elements of the original matrix are extracted to form a diagonal matrix in the diagonal part, while the off-diagonal part contains all off-diagonal elements and is negative. If the coefficient matrix has a significant lower triangular (or upper triangular) clustering feature of non-zero elements, it is split into a lower triangular part, a diagonal part, and an upper triangular part, i.e., Gaussian-Seidel splitting. Specifically, the original matrix is ​​decomposed into a combination of a lower triangular part (containing elements at and below the diagonal), a diagonal part (containing only diagonal elements), and an upper triangular part (containing elements above the diagonal). Depending on the needs of the iterative format, the combination of the lower triangular part (or upper triangular part) and the diagonal part is usually used as the main matrix after splitting. Both of the above splitting methods can ensure that the system matrix (i.e., the sum of the splitting matrix and the parameter matrix) in the subsequently constructed iterative linear equation system presents a simple structure of a diagonal or triangular matrix, thus creating conditions for efficient solution using the forward substitution method or the backward substitution method.

[0066] In another embodiment of the present invention, the construction and solution of the linear equation system in step S350 includes: Step S3521: Identify that the coefficient matrix of the complementary convex cone problem model is a symmetric positive definite matrix.

[0067] Identifying the mathematical properties of the coefficient matrix in a convex conical complementary problem model hinges on determining whether it is a symmetric positive definite matrix. Symmetric positive definiteness means the matrix is ​​equal to its transpose, and all eigenvalues ​​are positive. After a elastoplastic problem is discretized using the finite element method and transformed into a convex conical complementary form, the symmetry of the coefficient matrix often stems from the system's self-adjoint properties (such as the symmetry of the elastic stiffness matrix), while positive definiteness is typically guaranteed by the material's stable constitutive response (such as the positive definiteness of the elastic tensor) and appropriate boundary conditions. Identifying this property can be done by examining the matrix's symmetric structure and using numerical methods (such as calculating the minimum eigenvalues ​​or attempting Cholesky decomposition). Confirming the symmetric positive definiteness of the coefficient matrix provides a crucial basis for subsequently selecting a specific matrix splitting strategy with theoretical convergence guarantees.

[0068] Step S3522, splitting method, to construct a system of linear equations.

[0069] Based on the property that the coefficient matrix is ​​identified as a symmetric positive definite matrix, this step selects a matrix splitting method that can utilize and preserve the advantage of this property to construct the linear equation system required for iteration. Typical splitting strategies include decomposing the original symmetric positive definite matrix A into the difference of two matrices A = M - N, where the splitting matrix M is deliberately constructed to also be symmetric positive definite and is generally easier to invert or solve than A. For example, M can be chosen as the diagonal part of A (when A is diagonally dominant), a block diagonal approximation of A, or the product of a sparse approximation lower triangular matrix obtained through incomplete Cholesky decomposition and its transpose. This splitting method aims to ensure that the subsequently constructed iterative linear equation system (of the form M * x) new The system matrix M in ( = N * x + ...) is not only symmetric and positive definite, but may also have more favorable sparsity, condition number or structure, thus supporting an efficient and stable solution process.

[0070] Step S3523: Solve the system of linear equations to obtain the solution vector required for updating the iteration vector.

[0071] The linear equations constructed in step S3522 are numerically solved to obtain the solution vector needed to update the iteration vector. Given that the system matrix M is designed to be symmetric positive definite, efficient solvers optimized for such matrices can be invoked. Common methods include direct methods such as Cholesky decomposition (especially when the sparse structure of M allows for stable decomposition) or iterative methods such as the preprocessed conjugate gradient method (when M is extremely large or very sparse). The solution process fully utilizes the symmetric positive definite structure of the matrix to obtain the solution vector with relatively low computational and storage requirements. This solution vector is then used to update the current iteration vector in the modular basis iterative algorithm, driving the entire iterative process toward a solution to the convex cone complementarity problem.

[0072] Furthermore, before the modular basis matrix splitting iterative algorithm in step S300 performs iterative calculations using the optimal iterative parameters, it also includes: Step S301: For the current iteration vector, based on its geometric region within the second-order cone constraint, determine a mapping matrix according to predefined rules. The predefined rules are: when the current iteration vector is located inside or on the boundary of the second-order cone constraint, the mapping matrix is ​​the identity matrix or a variant thereof. When the current iteration vector is located outside the second-order cone constraint, the mapping matrix is ​​the negative identity matrix or a variant thereof.

[0073] For the current iteration vector, a corresponding mapping matrix is ​​determined based on its geometric region (i.e., inside, boundary, or outside) within the second-order cone constraint, according to a predefined mathematical rule. Specifically, if the current iteration vector is inside or on the boundary of the second-order cone constraint, the mapping matrix is ​​the identity matrix or a variant thereof (e.g., a constant multiple of the identity matrix); if the current iteration vector is outside the second-order cone constraint, the mapping matrix is ​​the negative identity matrix or a variant thereof. This rule is based on the spectral decomposition property of absolute value operations on a second-order cone: for a vector inside or on the cone surface, both spectral values ​​in its spectral decomposition are non-negative, and its modulus vector is equal to the vector itself or a specific linear transformation; while for a vector outside the cone, one spectral value is positive and the other negative, and its modulus vector is equal to the inverse of the vector or a corresponding linear transformation. By examining the relationship between the Euclidean norms of the first component and its subsequent components, the region to which the vector belongs can be automatically determined, and a mapping matrix can be generated according to the rule.

[0074] Step S302: Using the mapping matrix, the absolute value vector in the modular basis matrix splitting iteration algorithm is expressed as the product of the mapping matrix and the current iteration vector, thus transforming the implicit iteration equation into an explicit iteration form.

[0075] Using the mapping matrix determined in step S301, the key nonlinear term (absolute value vector) in the modular basis matrix splitting iterative algorithm is expressed as the product of the mapping matrix and the current iteration vector. Specifically, for the current iteration vector x, its absolute value vector abs(x) can be written precisely or approximately as the product of the mapping matrix φ(x) and x, i.e., abs(x) = φ(x) * x. Substituting this relationship into the original implicit iterative equation of the modular basis algorithm (usually in the form of...) In this process, the nonlinear term abs(x) can be eliminated, yielding a result that depends only on the current iteration vector x and the mapping matrix. The explicit iterative equation. This transformation converts the original iterative scheme, which inherently contains nonlinearity due to absolute value operations, into a formally linear equation, so that the update of each iteration can be viewed as the current iteration vector passing through a... The linear transformation that determines the outcome.

[0076] Step S303: Based on the explicit iteration form, construct an explicit iteration matrix with iteration parameters as variables.

[0077] Based on the explicit iterative form obtained in step S302, an explicit iterative matrix is ​​constructed with the iterative parameter (usually the relaxation factor ω) as the variable. Specifically, the explicit iterative equation is rewritten as x new = G(ω) * x + c, where G(ω) is the explicit iteration matrix, which consists of the splitting matrices M and N of the coefficient matrix A, the parameter ω, and the mapping matrix. Together they constitute, and the specific expression is as follows: This matrix explicitly depends on the iteration parameter ω, while other quantities are known given the current iteration vector and model coefficients. By constructing an explicit iteration matrix, the convergence speed and stability of the iterative algorithm are determined by the spectral radius (or norm) of this matrix, thus transforming the parameter selection problem into an optimization problem of the spectral properties or norm of the matrix G(ω).

[0078] Step S304: The optimal iteration parameters are obtained by minimizing the matrix norm of the explicit iteration matrix.

[0079] The optimal iterative parameter ω is obtained by minimizing a certain matrix norm (such as the Frobenius norm or 2-norm) of the explicit iterative matrix G(ω) constructed in step S303. opt Minimizing the norm of the iteration matrix aims to maximize the contractility of the iteration operator, thus theoretically guaranteeing the fastest convergence speed. Specifically, the matrix norm ||G(ω)|| is expressed as a function of the parameter ω, and matrix calculus or numerical optimization techniques are used to find the minimum value of this function. In the specific settings of this application, due to the mapping matrix... With its special structure (taking only variants of the identity matrix or negative identity matrix), this optimization problem can often be further simplified, and even a closed-form analytical expression for the optimal parameters can be derived, such as ω. opt The optimal parameters can be obtained by calculating the trace of the coefficient matrix A and the spectral radius of the splitting matrix M. This eliminates the need for time-consuming trial and error or numerical search, allowing the algorithm to be directly used to initialize the model basis iterative algorithm, ensuring near-optimal efficiency.

[0080] Furthermore, in step S304, the optimal iteration parameters are obtained by minimizing the matrix norm of the explicit iteration matrix, including: Step S3041: Express the square of the matrix norm of the explicit iteration matrix as a quadratic function of the iteration parameters.

[0081] The square of a chosen matrix norm (usually the Frobenius norm) of the explicit iteration matrix G(ω) is systematically expressed as a quadratic function of the iteration parameter ω. Specifically, the structure of matrix G(ω) can be represented as follows: ,in Let G(ω) be the mapping matrix. The square of its Frobenius norm is ∥G(ω)∥ F 2 It equals the sum of squares of all elements of the matrix, which can be further expanded into the operation of the trace of the matrix: ∥G(ω)∥ F 2 =trace(G(ω) T By substituting into the expression for G(ω) and utilizing the linear and cyclic properties of the matrix trace, the expression can be rearranged into a quadratic form with respect to the parameter ω: ∥G(ω)∥ F 2 =aω 2 +bω+c. Where the coefficients a, b, c are given by the matrices M, N, A, The trace and its combination operations determine the value. This transformation reduces the infinite-dimensional optimization problem of minimizing the matrix norm to the problem of minimizing a single-variable quadratic function.

[0082] Step S3042: By solving the condition for the quadratic function to obtain the minimum value, a closed-form expression for the trace of the coefficient matrix and the mapping matrix of the optimal iterative parameters with respect to the complementary convex cone problem model is obtained.

[0083] By solving the quadratic function f(ω) = aω 2 The condition for +bω+c to reach a minimum value yields the optimal iterative parameter ω. opt The closed-form expression is derived from the coefficient matrix A of the complementary convex cone problem model, its splitting matrices M and NN, and the mapping matrix. The trace operation is determined. Specifically, the optimal parameters are given by the following formula: ; in, Represents the trace of a matrix. Represents the mapping matrix The square of, express The transpose of the molecule. When the condition is met, the above formula represents the optimal parameter; otherwise, the optimal iterative parameter should be a positive value that is as small as possible.

[0084] Step S3043: Calculate the optimal iteration parameters based on the closed-form expression.

[0085] The optimal iterative parameter ω is derived based on step S3042. opt The closed-form expression is used to perform specific numerical calculations to obtain its value. The calculation process involves processing known matrices A, M, N and the mapping matrix corresponding to the current iteration. Perform the specified trace operation. The trace can be efficiently calculated by summing the diagonal elements of the matrix. Since the expression is analytic, the calculation only requires basic matrix operations (such as extracting diagonal elements and matrix multiplication and summation), without any iterative optimization or trial-and-error process. The resulting ω opt It can be directly used as the parameter input for the modular basis matrix splitting iterative algorithm. During the algorithm iteration process, if the mapping matrix... As the iteration vector changes, steps S3041 to S3043 can be re-executed periodically or at each step to update the parameters, thereby achieving dynamic optimal parameter tuning.

[0086] like Figure 2 As shown, Figure 2 This paper presents the overall approach to transforming the three-dimensional elastoplastic initial-boundary value problem into a standard conical complementary problem within a time-incremental framework. First, the elastoplastic governing equations are implicitly discretized in time, transforming the continuous evolution problem into an incremental problem from t_n to t_{n+1}. Then, the yield condition, flow rule, and consistency condition are uniformly expressed as complementary constraints. Based on the mathematical structure of different yield criteria, they are represented as second-order conical (SOC) or rotated second-order conical (RSOC) constraints, thus achieving a unified description of yield models such as Mohr–Coulomb, Drucker–Prager, and Von Mises / J2 within the same convex optimization framework. Building upon this, extended variables are introduced, and a global algebraic system containing equilibrium equations, constitutive relations, and plastic evolution equations is constructed. Finally, the incremental elastoplastic problem is standardized into an optimization problem combining linear mapping and conical complementary constraints, which can be directly solved using an efficient primordial-dual conical algorithm. Specifically, Figure 2For a detailed explanation of the symbols used, please refer to Table 1.

[0087] Table 1

[0088] Based on the above, this invention focuses on the three-dimensional elasto-plastic initial-boundary value problem. Addressing the challenges of strong material nonlinearity, diverse yield criteria, and insufficient convergence and uniformity of traditional return-mapping algorithms in complex fault systems, it proposes a unified numerical modeling and solution framework based on conical complementarity theory. Through time-increment discretization, yield-flow condition complementarity, and the convex conical representation of the yield criterion, the three-dimensional viscoelastic-plastic problem is systematically transformed into a standard second-order conical complementarity problem, thereby achieving efficient and stable solutions for multiple yield models within the same mathematical framework. The process of this invention is further summarized below: (1) Elastic-plastic initial-boundary value problem and incremental discretization Starting with the elastoplastic initial-boundary value problem of a continuous medium, the governing equations include equilibrium equations, strain-displacement relations, constitutive relations, and corresponding boundary and initial conditions. To characterize the path-dependent behavior of the material, an implicit time integral scheme is used to discretize the governing equations, transforming the problem into an incremental elastoplastic problem from t to tn+1. Within this incremental framework, unknowns such as displacement increment, strain increment, and plastic strain increment are introduced, allowing stress updates, constitutive relations, and plastic evolution to be expressed in incremental form. This transforms the evolution problem into a series of nonlinear but structurally clear static subproblems.

[0089] (2) Formalization of the complementary system of yield condition and flow law Within each time step, plastic behavior is described by the yield function and the flow rule. The yield condition, the nonnegativity of the plasticity multiplier, and the consistency condition can be uniformly written as the Kuhn–Tucker complementarity relation, that is, the yield function and the plasticity multiplier satisfy nonnegativity and positive interactive complementarity constraints. Through this treatment, the "yield-loading-unloading" discrimination process in traditional plasticity theory is systematically transformed into a nonlinear complementary problem, laying the foundation for the subsequent introduction of convex optimization tools.

[0090] (3) Correspondence between yield criterion and cone type Depending on the specific yield criterion used, the complementary constraints are further convex conicalized: for the Von Mises (J2) yield criterion, its yield condition can be naturally expressed as a standard second-order conical constraint; for the Drucker–Prager yield criterion, by introducing the coupling relationship between the pressure invariant and the deviatoric stress norm, the yield condition is transformed into a rotating second-order conical constraint; for the Mohr–Coulomb yield criterion, by stress space transformation or convex approximation, its non-smooth yield surface is embedded in the second-order conical constraint system.

[0091] Therefore, different rock mechanics yielding models can be described within a unified second-order cone or rotating second-order cone framework.

[0092] (4) Introduction of extended variables and construction of the global algebraic system To achieve a unified expression of the yield condition, constitutive relation, and equilibrium equations, extended variables such as plasticity multipliers and stress invariants are introduced, transforming the nonlinear relationship into a combination of linear mapping and cone constraints. After spatial discretization, the incremental problem can be represented as a coupled algebraic system consisting of three parts: discrete equilibrium equations, stress-strain constitutive relations, and complementary conditions for plastic evolution and yield. All unknowns are uniformly assembled into an extended variable vector, thus forming a well-structured overall algebraic system.

[0093] (5) Form and numerical solution of the standard conical complement problem Through the above steps, the original three-dimensional elastoplastic problem is ultimately transformed into a standard form of cone complementarity problem: under the action of a linear operator, the unknown variables simultaneously satisfy the second-order cone (or rotated second-order cone) constraint and the complementary orthogonality condition. This standard form can be directly solved using mature cone optimization and primal-dual algorithms, avoiding the complex yield determination and local nonlinear iteration processes in traditional return mapping methods, significantly improving the robustness and parallel computing potential of the algorithm.

[0094] like Figure 3a and Figure 3b As shown, the present invention (ours) can significantly reduce computation time, reducing computation time by about 90% compared to traditional methods, effectively alleviating the bottleneck of system matrix construction and solution; especially when the scale increases, the cost of solving the system matrix will become increasingly large, becoming a bottleneck restricting the development of this model.

[0095] like Figure 4a and Figure 4b As shown, the optimal parameters proposed in this invention significantly reduce the number of iterations when solving this type of problem, greatly improving the automation and robustness of the algorithm, and further illustrating the importance of studying optimal parameter selection strategies.

[0096] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor and a memory connected to the at least one processor. The memory stores instructions executable by the at least one processor, which, when executed by the at least one processor, cause the at least one processor to perform the aforementioned elastoplastic numerical simulation method.

[0097] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described elastoplastic numerical simulation method.

[0098] The embodiments of the present invention aim to protect an elastic-plastic numerical simulation method, which has the following effects: 1. By systematically transforming the elastoplastic incremental governing equations into a convex cone complementary problem model, a unified mathematical expression for multiphysics constraints is achieved. The traditionally separate elastoplastic constitutive relations, yield criteria, and global equilibrium equations are organically integrated into a rigorous mathematical framework based on convex cones and complementary conditions. This modeling approach not only fundamentally avoids the convergence difficulties and accuracy losses commonly encountered in traditional finite element methods when dealing with inequality constraints, but also provides a high-fidelity numerical description of the nonlinear mechanical behavior of complex geological bodies, significantly improving the physical reliability and numerical stability of the simulation results. 2. By theoretically deriving and explicitly determining the optimal iterative parameters of the modular basis matrix splitting iterative algorithm, the huge computational burden of parameter tuning in traditional trial-and-error methods is completely eliminated. The implicit iterative equation is made explicit by introducing a mapping function, and the closed-form expression of the optimal parameters is directly derived by minimizing the norm of the explicit iterative matrix. This enables the algorithm to efficiently and automatically obtain high-performance parameters before each solution, avoiding the problem of slow convergence or even failure due to improper parameters. The parameter optimization cost, which may have originally been dominant, is reduced to almost zero, which greatly improves the efficiency and reliability of the overall solution process. 3. By designing targeted solution strategies for ultra-large-scale sparse systems, a fundamental reduction in computational complexity has been achieved, from cubic to quadratic and even linear levels. Fully utilizing the sparsity and structural characteristics of the system matrix after discretization of the elasto-plastic convex conical complementary model, and combining efficient numerical techniques such as matrix splitting, back-splitting, and Krylov subspace iteration, the solution process can fully leverage the advantages of parallel computing and significantly reduce memory consumption. This makes it possible to perform full-size, high-resolution elasto-plastic simulations of geological bodies containing massive numbers of elements and nodes, effectively solving the long-standing bottleneck of traditional methods being unable to be applied to practical large-scale engineering and scientific problems due to computational resource limitations.

[0099] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0100] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0101] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0102] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0103] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A numerical simulation method for elastoplasticity, characterized in that, Includes the following steps: Obtain the input data for the numerical model of the geological body to be simulated; Based on the input data, an elastoplastic incremental control equation is constructed, and the elastoplastic incremental control equation is transformed into a convex cone complementary problem model. The convex cone complementary problem model is a mathematical model that expresses the elastoplastic constitutive relation, yield criterion and equilibrium equation as complementary conditions under convex cone constraints. The model of the complementary convex cone problem is solved by the model basis matrix splitting iterative algorithm. The model basis matrix splitting iterative algorithm uses the optimal iteration parameters to perform iterative calculations and obtain the displacement field, stress field and plastic state field of the geological body to be simulated. The optimal iteration parameters are determined by making the implicit iteration equation of the model basis matrix splitting iterative algorithm explicit and minimizing the norm of the resulting explicit iteration matrix. The displacement field, stress field, and plastic state field are output as the results of the elastoplastic numerical simulation.

2. The elastoplastic numerical simulation method according to claim 1, characterized in that, The input data includes discretized mesh data, mechanical boundary condition data, material elastoplastic constitutive parameters, and initial stress field and displacement field data. The step of constructing an elastoplastic incremental control equation based on the input data, and transforming the elastoplastic incremental control equation into a complementary convex cone problem model, includes: Based on the discretized mesh data and the material's elastoplastic constitutive parameters, an elastoplastic constitutive complementary system that couples the yield function, the plastic multiplier increment, and their complementary conditions is constructed within the incremental step. According to the yield criterion, the constraints related to the yield function in the elastoplastic constitutive complementary system are equivalently mapped to second-order cone constraints, which are convex cone sets defined on the vector space composed of stress tensor invariants. Based on the discretized grid data, the mechanical boundary condition data, and the initial stress field and displacement field data, the overall equilibrium equation is assembled, and the overall equilibrium equation, the second-order cone constraint, and the elastoplastic constitutive complementary system are combined to construct the convex cone complementary problem model.

3. The elastoplastic numerical simulation method according to claim 2, characterized in that, The step of mapping the constraints related to the yield function in the elastoplastic constitutive complementary system to second-order cone constraints according to the yield criterion includes: Based on the yield criterion, determine the coefficients of the linear term related to hydrostatic pressure and the norm term related to deviatoric stress in the yield function; Based on the linear term coefficients and the norm term coefficients, a linear inequality is constructed consisting of the Euclidean norm of the hydrostatic pressure and the deviatoric stress tensor. The linear inequality is identified as a second-order cone constraint with the hydrostatic pressure and deviatoric stress tensor components as variables.

4. The elastoplastic numerical simulation method according to claim 2, characterized in that, The assembly of the overall equilibrium equations, combined with the second-order cone constraint and the elastoplastic constitutive complementary system, constructs the model of the convex cone complementarity problem, including: Substituting the constitutive relations in the elastoplastic constitutive complementary system into the global equilibrium equations, we obtain a set of governing equations with nodal displacement increments and plastic multiplier increments as the basic unknowns. By combining the governing equations, the second-order cone constraint, and the complementary conditions in the elastoplastic constitutive complementary system, a hybrid complementary equation set with the nodal displacement increment and the plastic multiplier increment as unknown variables is constructed as the model for the convex cone complementary problem.

5. The elastoplastic numerical simulation method according to claim 2, characterized in that, The method of solving the complementary convex cone problem model using the modular basis matrix splitting iterative algorithm includes: Initialize the iteration parameters, initial iteration vector, and convergence tolerance of the modular basis matrix splitting iteration algorithm, and set its coefficient matrix and right-hand vector based on the convex cone complementarity problem model, and set the initial iteration vector as the current iteration vector; In each iteration, the absolute value vector of the current iteration vector is calculated. The absolute value vector is a vector formed by calculating the magnitude of each sub-vector of the current iteration vector corresponding to the second-order cone constraint. Based on the iteration parameters and the absolute value vector, construct intermediate auxiliary variable pairs that satisfy the convex cone constraint and the complementary condition; Determine whether the norm of the residual vector formed by the intermediate auxiliary variable pair is less than or equal to the convergence tolerance, wherein the residual vector is used to measure the degree to which the intermediate auxiliary variable pair deviates from the linear relationship defined by the convex cone complementarity problem model; If the norm of the residual vector is greater than the convergence tolerance, then based on the splitting form of the coefficient matrix, the iteration parameters, the current iteration vector, the absolute value vector, and the right-hand vector, a system of linear equations is constructed and solved. The solution vector obtained is used as the new current iteration vector, and the step of calculating the absolute value vector of the current iteration vector is returned. Otherwise, the intermediate auxiliary variable pair is output as the solution to the convex cone complementarity problem model.

6. The elastoplastic numerical simulation method according to claim 5, characterized in that, The construction and solution of the linear equation system includes: Based on the inherent properties of the coefficient matrix of the complementary convex cone problem model, a matrix splitting method is selected so that the system matrix of the constructed linear equation system is a diagonal matrix or a triangular matrix. The inherent properties of the coefficient matrix are determined by the input data and the mathematical structure of the complementary convex cone problem model. Taking advantage of the structural characteristic that the system matrix is ​​a diagonal or triangular matrix, the linear equation system is solved using the forward substitution method or the backward substitution method.

7. The elastoplastic numerical simulation method according to claim 6, characterized in that, The selection of matrix splitting methods based on the inherent properties of the coefficient matrix of the complementary convex cone problem model includes: Identify the sparse structure and non-zero element distribution pattern of the coefficient matrix; Based on the sparse structure and distribution pattern, the coefficient matrix is ​​split into a diagonal part and a non-diagonal part, or into a lower triangular part, a diagonal part and an upper triangular part, so that the system matrix is ​​a diagonal matrix or a triangular matrix.

8. The elastoplastic numerical simulation method according to claim 5, characterized in that, The construction and solution of the linear equation system includes: The coefficient matrix of the model for the complementary convex cone problem is identified as a symmetric positive definite matrix. Based on the symmetric positive definiteness, the coefficient matrix is ​​split into a matrix splitting method of symmetric positive definite part and remainder part to construct the linear equation system; Solve the system of linear equations to obtain the solution vector required for updating the iteration vector.

9. The elastoplastic numerical simulation method according to any one of claims 5-8, characterized in that, Before performing iterative calculations using the optimal iteration parameters, the modular basis matrix splitting iterative algorithm also includes: For the current iteration vector, a mapping matrix is ​​determined according to a predefined rule based on the geometric region in which it is located within the second-order cone constraint. The predefined rule is as follows: when the current iteration vector is located inside or on the boundary of the second-order cone constraint, the mapping matrix is ​​an identity matrix or a variant thereof; when the current iteration vector is located outside the second-order cone constraint, the mapping matrix is ​​a negative identity matrix or a variant thereof. Using the mapping matrix, the absolute value vector in the modular basis matrix splitting iteration algorithm is represented as the product of the mapping matrix and the current iteration vector, thus transforming the implicit iteration equation into an explicit iteration form; Based on the explicit iteration form, an explicit iteration matrix is ​​constructed with the iteration parameters as variables; The optimal iteration parameters are obtained by minimizing the matrix norm of the explicit iteration matrix.

10. The elastoplastic numerical simulation method according to claim 9, characterized in that, The process of minimizing the matrix norm of the explicit iteration matrix to obtain the optimal iteration parameters includes: The square of the matrix norm of the explicit iteration matrix is ​​expressed as a quadratic function of the iteration parameters; By solving the condition for the quadratic function to reach its minimum value, a closed-form expression for the coefficient matrix and the trace of the mapping matrix of the optimal iterative parameters with respect to the complementary convex cone problem model is obtained. The optimal iteration parameters are calculated based on the closed-form expression.

Images