Dynamic relaxation numerical calculation method and simulation system

Abstract

The application discloses a dynamic relaxation numerical calculation method and a simulation system, wherein the method adopts a maximum iteration number and an iteration accuracy allowed value to control loading grading, forms a plurality of sub-loading steps, and adopts full-quantity displacement in the sub-loading steps to calculate strain and stress integration in a dynamic relaxation solving process of the sub-loading steps; a specific scheme for controlling loading grading is as follows: when ending iteration based on iteration accuracy: updating the starting position information based on the target displacement data; updating the starting stress state based on the target stress state; accumulating the total loading step based on the loading step increment; adjusting the loading step increment based on a preset adaptive rule; judging whether the calculation is completed based on the updated total loading step; when ending iteration based on the iteration number, adjusting the loading step increment based on the preset adaptive rule. The application can overcome the defects of iteration convergence difficulty of a conventional implicit static force solution and a false stress path of a conventional dynamic relaxation method, and is better applicable to a nonlinear complex rock-soil constitutive model.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geotechnical engineering, and in particular to a dynamic relaxation numerical calculation method and simulation system. BACKGROUND

[0002] Dynamic relaxation method is a numerical calculation method, which converts static problem into dynamic problem, and in the process of calculation, the iteration progress of displacement field is controlled by dynamic equation. If the displacement field does not satisfy the mechanical equilibrium condition, the unbalanced force existing in the system will force the displacement field to progress towards the position of equilibrium until convergence. In contrast, the implicit static method converts the static problem into the problem of solving nonlinear equations by forming the local stiffness matrix at the element level according to the material stiffness matrix of each element, and then assembling the global stiffness matrix. For linear problems or weakly nonlinear problems, the calculation efficiency of implicit static solution is usually higher, so it is often used in general finite element software (such as Abaqus and ANSYS). However, for strong nonlinear problems of geotechnical numerical analysis, due to the dramatic change of stiffness matrix, it is difficult to converge by using implicit static solution. In contrast, the explicit dynamic solution strategy of dynamic relaxation method adopts the strategy of "small step fast walk", which has great improvement in convergence in nonlinear problems. In addition, unlike the implicit static solution which requires the stiffness matrix of the constitutive model, dynamic relaxation method can also be used for material models that do not have stiffness matrix. For example, one branch of the constitutive model of geotechnological materials is the sub-plasticity model, which cannot give the incremental pseudo-linear form, and there is no stiffness matrix. If it is forced to be used for implicit static calculation, it will make the iteration convergence of static solution extremely difficult.

[0003] The main disadvantage of the conventional dynamic relaxation method is the problem of false stress path interference. SUMMARY

[0004] The purpose of the present application is to solve the defect of false stress path interference in the prior art.

[0005] In order to achieve the above-mentioned purpose of the application, the present application provides a dynamic relaxation numerical calculation method and simulation system;

[0006] A dynamic relaxation numerical calculation method is used for simulating and analyzing a target geotechnological constitutive model based on a target external force, comprising the following steps:

[0007] Obtain the total loading step and the loading step increment corresponding to the current;

[0008] Obtain the starting stress state and starting position information corresponding to the total loading step;

[0009] Obtain the virtual mass information corresponding to the total loading step;

[0010] Based on the total loading step and the loading step increment, corresponding action external force data is calculated and obtained;

[0011] Enter the sub-loading step, and based on the initial stress state, the virtual mass information and the action external force data, full-quantity displacement is used for dynamic relaxation iteration solution until a preset iteration end condition is reached to obtain corresponding target displacement data and target stress state; wherein, the iteration end condition is that the iteration accuracy reaches a preset iteration accuracy allowed value, or the iteration number reaches a preset maximum iteration number;

[0012] Based on the triggered iteration end condition, loading grading is performed, specifically;

[0013] When the iteration is ended based on the iteration accuracy:

[0014] The initial position information is updated based on the target displacement data;

[0015] The initial stress state is also updated based on the target stress state;

[0016] The total loading step is also accumulated based on the loading step increment;

[0017] The loading step increment is also adjusted based on a preset adaptive rule;

[0018] Whether the calculation is completed is also judged based on the updated total loading step;

[0019] When the iteration is ended based on the iteration number, the loading step increment is adjusted based on a preset adaptive rule.

[0020] As an implementable manner, the iteration process performed in the sub-loading step includes the following steps:

[0021] The historical unbalanced force, historical node speed and historical node displacement corresponding to each node in the last iteration step are obtained;

[0022] Based on the virtual mass information, the historical unbalanced force, historical node speed and historical node displacement, the node speed and node displacement corresponding to each node in the current iteration step are calculated;

[0023] For the node with displacement constraint boundary condition, the component of its node speed and node displacement in the constrained direction is set to 0;

[0024] Based on the node displacement, the element strain increment corresponding to each element is calculated;

[0025] Based on the element strain increment, element stress integration is performed with the initial value of the initial stress state to obtain the corresponding intermediate stress state;

[0026] Based on the intermediate stress state and the applied external force data, unbalanced force calculations are performed to obtain the unbalanced forces corresponding to each node in the current iteration step.

[0027] Determine whether the preset iteration end condition has been met. When it is determined that the iteration has ended, the node displacement is used as the target displacement data and the obtained intermediate stress state is used as the target stress state. When it is determined that the iteration calculation should continue, the unbalanced force is used as the historical unbalanced force, the node velocity is used as the historical node velocity, and the node displacement is used as the historical node displacement for use in the next iteration.

[0028] As one possible implementation method:

[0029] GPU acceleration is used for the calculation and updating of each unit variable and each node variable.

[0030] As one possible implementation method:

[0031] Stress state includes stress and state variables;

[0032] The formula for calculating the stress in the intermediate stress state is:

[0033] σ (k),(t) =σ (k),(T)) +△σ(σ (k),(T) ,s (k),(T) ,△ε (k),(t) )

[0034] in:

[0035] k represents the k-th tetrahedral element;

[0036] t represents the number of iterations corresponding to the current iteration step;

[0037] T represents the total loading step;

[0038] σ (k),(t) This represents the stress corresponding to element k in the current iteration step;

[0039] σ (k),(T) This represents the stress of element k corresponding to the total loading step T, that is, the stress corresponding to element k in the initial stress state;

[0040] s (k),(T) This represents the state variable of unit k corresponding to the total loading step T;

[0041] △ε (k),(t) This represents the element strain increment corresponding to element k in the current iteration step;

[0042] △σ(σ (k),(T) ,s(k),(T) (k),(t) represents the stress increment calculated by stress integration with initial state of σ (k),(T) and s (k),(T) , and strain path of △ε (k),(t) ;

[0043] The formula for calculating the state variable of the intermediate stress state is:

[0044] s (k),(t) = s (k),(T)) +△s(σ (k),(T) ,s (k),(T) ,△ε (k),(t) )

[0045] Wherein:

[0046] s (k),(t) represents the state variable corresponding to the element k in the current iteration step;

[0047] △s(σ (k),(T) ,s (k),(T) ,△ε (k),(t) ) represents the stress increment calculated by stress integration with initial state of σ (k),(T) and s (k),(T) , and strain path of △ε (k),(t) , and state increment △s obtained by stress integration calculation.

[0048] As an implementable manner, the formula for calculating the element strain increment is:

[0049] The formula for calculating the element strain increment is:

[0050]

[0051] Wherein:

[0052] △ε (k),(t) represents the element strain increment corresponding to the element k in the current iteration step;

[0053] x, y, z represent three dimensions;

[0054] The subscript i or j is used for the components of the second-order tensor △ε (k),(t) in the nine dimensions of xx, xy, xz, yx, yy, yz, zx, zy, zz;

[0055] (k, l) represents the surface opposite to the node l in the element k

[0056] ∑ l∈k represents the summation for each tetrahedral element k with respect to its corresponding nodes and corresponding surfaces; ​

[0057] V (k),(T) V represents the volume of the element k corresponding to the total loading step T;

[0058] u (l),(t) u represents the node displacement of the node l in the current iteration step;

[0059] n (k,l),(T) n represents the normal vector of the face (k, l) in the total loading step T;

[0060] S (k,l),(T) S represents the area of the face (k, l) in the total loading step T.

[0061] As an implementable manner, the specific steps of calculating the non-equilibrium force based on the intermediate stress state and the acting external force data are as follows:

[0062] The acting external force data includes the node acting external force of each node;

[0063] Based on the intermediate stress state, the force contribution applied to the node by the stress of each element in the current iteration step is calculated;

[0064] The force contributions are accumulated to obtain the node acting internal force of each node;

[0065] Based on the node acting external force and the node acting internal force, the non-equilibrium force corresponding to each node is calculated and obtained.

[0066] As an implementable manner:

[0067] The iteration precision is the ratio of the first non-equilibrium force maximum value and the second non-equilibrium force maximum value;

[0068] The first non-equilibrium force maximum value is the initial non-equilibrium force maximum value corresponding to the current sub-loading step;

[0069] The second non-equilibrium force maximum value is the maximum value of the non-equilibrium force of each node in the current iteration step.

[0070] As an implementable manner:

[0071] The acting external force data includes the node acting external force of each node;

[0072] The calculation formula of the node acting external force is:

[0073]

[0074] Wherein:

[0075] bs∈l represents that the face bs participating in accumulation is the face connected with the node l;

[0076] k∈l represents that the unit k participating in accumulation is the unit connected with the node l;

[0077] T represents the total loading step;

[0078] △T represents the loading step increment;

[0079] represents the node action external force corresponding to the node l in the total loading step T;

[0080] is the point load borne by the node l;

[0081] q (bs) is the surface load borne by the surface bs;

[0082] S (bs),(T) is the area of the surface bs;

[0083] ρ (k) is the density of the unit k, which is obtained based on the target rock-soil constitutive model;

[0084] g is the acceleration of gravity;

[0085] V (k),(T) is the volume of the unit k.

[0086] As an implementable manner, the adaptive rule for adjusting the loading step increment is specifically:

[0087] When the iteration is ended based on the iteration accuracy:

[0088] obtain the iteration number corresponding to the end of iteration;

[0089] generate an adjustment weight based on the iteration number and the maximum iteration number;

[0090] weight the loading step increment based on the adjustment weight to obtain a first candidate increment;

[0091] generate a second candidate increment based on the difference between the preset maximum total loading step and the updated total loading step;

[0092] take the smaller value between the first candidate increment and the second candidate increment as the new loading step increment;

[0093] When the iteration is ended based on the iteration number:

[0094] generate a third candidate increment based on the loading step increment, the third candidate increment being smaller than the loading step increment;

[0095] generate a fourth candidate increment based on the difference between the preset maximum total loading step and the total loading step;

[0096] taking the smaller value of the third candidate increment and the fourth candidate increment as a new loading step increment.

[0097] An emulation system for simulating and analyzing a target rock-soil constitutive model based on a target external force, comprising:

[0098] An acquisition module configured to acquire a total loading step and a loading step increment corresponding to a current target;

[0099] A calculation module configured to acquire starting stress state and starting position information corresponding to the total loading step, and to acquire virtual mass information corresponding to the total loading step, and to calculate corresponding acting external force data based on the total loading step and the loading step increment;

[0100] An iterative solution module configured to enter a sub-loading step, and to perform dynamic relaxation iterative solution based on the starting stress state, the virtual mass information and the acting external force data using full-quantity displacement until a preset iteration end condition is reached to obtain corresponding target displacement data and target stress state; wherein the iteration end condition is that an iteration accuracy reaches a preset iteration accuracy allowable value, or an iteration number reaches a preset maximum iteration number.

[0101] An updating module configured to perform loading grading based on the triggered iteration end condition, comprising a first updating unit and a second updating unit.

[0102] The first updating unit is configured to perform the following steps when iteration is ended based on iteration accuracy:

[0103] updating the starting position information based on the target displacement data;

[0104] updating the starting stress state based on the target stress state;

[0105] accumulating the total loading step based on the loading step increment;

[0106] adjusting the loading step increment based on a preset adaptive rule;

[0107] judging whether the calculation is completed based on the updated total loading step;

[0108] The second updating unit is configured to adjust the loading step increment based on a preset adaptive rule when iteration is ended based on iteration number.

[0109] Beneficial effects:

[0110] The application controls the loading stage by using the maximum iteration number and the iteration accuracy allowable value, forms several sub-loading steps, and uses the full displacement in the sub-loading step to calculate the strain and stress integral in the dynamic relaxation solving process of the sub-loading step, so as to overcome the defects of the false stress path of the conventional dynamic relaxation method and the convergence difficulty of the conventional implicit static force solving, and better apply to the nonlinear complex rock-soil constitutive model. BRIEF DESCRIPTION OF DRAWINGS

[0111] Figure 1 It is a work flow diagram of the dynamic relaxation numerical calculation method in embodiment 1.

[0112] Figure 2 It is a schematic diagram of the correspondence between the nodes and faces of the tetrahedron element.

[0113] Figure 3 It is a schematic diagram of the calculation results based on FLAC3D.

[0114] Figure 4 It is a schematic diagram of the calculation results based on the dynamic relaxation numerical calculation method in embodiment 1. DETAILED DESCRIPTION

[0115] The exemplary embodiments of the present disclosure are described below with reference to the accompanying drawings, including various details in order to facilitate understanding. They should be considered as merely exemplary. Therefore, those of ordinary skill in the art will recognize various changes and modifications of the embodiments described herein, without departing from the scope and spirit of the present disclosure. Also, in the following description, descriptions of well-known functions and structures are omitted in order to make the present disclosure clear and concise.

[0116] Embodiment 1, a dynamic relaxation numerical calculation method, for simulating and analyzing a target rock-soil constitutive model based on a target external force, referring to Figure 1 , comprising the following steps:

[0117] S100, data reading:

[0118] Obtain the total loading step T and the loading step increment △T corresponding to the current;

[0119] S200, data preparation:

[0120] S210, obtain the starting stress state and starting position information corresponding to the total loading step T;

[0121] The stress state includes stress and state variable, and the starting stress state includes the stress σ (k),(T) and state variable s (k),(T) of the total loading step T corresponding, wherein k represents the corresponding tetrahedron element, and T represents the corresponding total loading step.

[0122] The position information includes node space positions of the nodes, and the initial position information includes node space positions x of the nodes corresponding to the total loading step T (l),(T) Wherein, I represents the corresponding node;

[0123] Note:

[0124] The grid geometry information corresponding to the target rock-soil constitutive model is imported in advance, and the grid geometry information includes a plurality of tetrahedral elements and nodes corresponding to each tetrahedral element;

[0125] The initial stress state, i.e., the stress state of each node at T=0, is obtained based on the target rock-soil constitutive model at the beginning of dynamic relaxation numerical calculation; and the initial position information, i.e., the node space position of each node at T=0, is obtained based on the grid geometry information;

[0126] During the dynamic relaxation numerical calculation process, the initial stress state and the initial position information will be updated according to the iteration solution of the sub-loading step.

[0127] S220, calculating corresponding geometry information based on the total loading step T and the initial position information;

[0128] The geometry information includes:

[0129] The volume V of each tetrahedral element (k),(T) ;

[0130] The area S of each surface in the element (k,l),(T) , wherein (k, I) represents the surface opposite to the node I in the element k;

[0131] The normal vector n (k,l),(T) .

[0132] Taking a single tetrahedron as an example, the node and surface numbering conditions are shown in the accompanying Figure 2 .

[0133] S230, calculating the virtual mass corresponding to the current total loading step based on the geometry information;

[0134] Specifically:

[0135] S231, calculating the virtual mass contribution m (k,l),(T) of each element k to the connected node I based on the geometry information;

[0136] The calculation formula is:

[0137]

[0138] Wherein:

[0139] i is used to represent the component of the corresponding parameter in the x, y, z three dimensions, as in the formula, it represents the normal vector n (k ,l),(T) component in the x, y, z three dimensions;

[0140] a1 (k) It is the compressive modulus of unit k with side limits, which is obtained based on the target rock-soil constitutive model;

[0141] V (k),(T) It is the corresponding volume of unit k;

[0142] n (k,l),(T) It is the normal vector corresponding to surface (k, l);

[0143] S (k,l),(T) It is the area corresponding to surface (k, l).

[0144] S232, the mass contribution m (k,l),(T) is accumulated to obtain the virtual mass M (l),(T) of each node l.

[0145] The specific expression is:

[0146]

[0147] Wherein, k∈l indicates that the unit k participating in the accumulation calculation is the unit connected with the node l.

[0148] S240, based on the total loading step and the loading step increment, the corresponding action external force data is calculated and obtained;

[0149] The action external force data includes the node action external force corresponding to each node;

[0150] The node action external force is used to indicate the point load borne by the corresponding node and the force of the surface load and the body weight gravity load dispersed to the node, that is, the target external force includes the gravity load, the surface load and the point load applied externally.

[0151] The application adopts the maximum iteration number and the iteration accuracy allowed value to control the loading classification, and the loading classification

[0152] The calculation formula of the node action external force is:

[0153]

[0154] Wherein:

[0155] bs∈l indicates that the surface bs participating in the accumulation is connected with the node l;

[0156] k∈l indicates that the unit k participating in the accumulation is the unit connected with the node l;

[0157] T represents a total loading step;

[0158] ΔT represents a loading step increment;

[0159] represents a node action external force corresponding to node I in the total loading step T;

[0160] is a point load borne by node I;

[0161] q (bs) is a surface load borne by surface bs;

[0162] S (bs),(T) is an area of surface bs;

[0163] ρ (k) is a density of unit k, obtained based on a target rock-soil constitutive model;

[0164] g is a gravitational acceleration;

[0165] V (k),(T) is a volume of unit k.

[0166] S300, enter a sub-loading step, and perform dynamic relaxation iteration based on the initial stress state and the action external force data using total displacement until a preset iteration end condition is reached to obtain corresponding target displacement data and a target stress state;

[0167] The iteration end condition is that iteration accuracy reaches a preset iteration accuracy allowable value, or iteration times reaches a preset maximum iteration times.

[0168] The iteration process performed in the sub-loading step includes the following steps:

[0169] S310, obtain virtual mass information corresponding to the total loading step;

[0170] That is, extract the virtual mass M (l),(T) .

[0171] S320, obtain historical unbalanced forces, historical node velocities and historical node displacements corresponding to each node in the last iteration step;

[0172] When initial iteration calculation is performed, that is, the current iteration step t = 1, the initial unbalanced force corresponding to the moment t = 0 is calculated as the historical unbalanced force, and the node velocity and the node displacement of each node are initialized to 0, that is, each historical node velocity and each historical node displacement is initialized to 0.

[0173] S330, based on the virtual mass information, the historical unbalanced force, historical node speed and historical node displacement, calculating the node speed and node displacement corresponding to each node in the current iteration step;

[0174] S331, the calculation formula of the node speed corresponding to each node in the current iteration step is:

[0175]

[0176] Wherein:

[0177] v (l),(t) represents the node speed corresponding to node l in the current iteration step;

[0178] v (l),(t-1) represents the node speed corresponding to node l in the previous iteration step or the initial node speed, that is, the historical node speed;

[0179] represents the unbalanced force corresponding to node l in the previous iteration step or the initial unbalanced force, that is, the historical unbalanced force;

[0180] M (l),(T) represents the virtual mass corresponding to node l in the total loading step T;

[0181] c d represents the damping factor, and the damping factor c d is a constant between 0.3 and 0.8 in the embodiment.

[0182] S332, the calculation formula of the node displacement corresponding to each node in the current iteration step is:

[0183] u (l),(t) = u (l),(t-1) +v (l),(t)

[0184] Wherein:

[0185] u (l),(t) represents the node displacement corresponding to node l in the current iteration step;

[0186] u (l),(t-1) represents the node displacement corresponding to node l in the previous iteration step or the initial node displacement, that is, the historical node displacement.

[0187] S340, for the node with displacement constraint boundary condition, the component of the node speed and the node displacement in the constrained direction is set to 0;

[0188] S350, based on the node displacement, calculating the element strain increment corresponding to each element;

[0189] The formula for calculating the unit strain increment is:

[0190]

[0191] Wherein:

[0192] △ε (k),(t) represents the unit strain increment corresponding to the unit k in the current iteration step;

[0193] x, y, z represent three dimensions;

[0194] The subscript i or j is used for the second-order tensor △ε (k),(t) The components in the xx, xy, xz, yx, yy, yz, zx, zy, zz nine dimensions;

[0195] (k, l) represents the surface opposite to node l in unit k

[0196] ∑ l∈k represents that for each tetrahedral unit k, the corresponding nodes and the corresponding surfaces are accumulated;

[0197] V (k),(T) represents the volume of unit k corresponding to the total loading step T;

[0198] u (l),(t) represents the node displacement corresponding to node l in the current iteration step;

[0199] n (k,l),(T) represents the normal vector of face (k, l) in the total loading step T;

[0200] S (k,l),(T) represents the area of face (k, l) in the total loading step T.

[0201] S360, based on the unit strain increment, the unit stress integration is carried out with the initial value of the initial stress state, and the corresponding intermediate stress state is obtained;

[0202] S361, the formula for calculating the stress of the intermediate stress state is:

[0203] σ (k),(t) = σ (k),(T)) + △σ(σ (k),(T) , s (k),(T) , △ε (k),(t) )

[0204] Wherein:

[0205] k represents the kth tetrahedral unit;

[0206] t represents the iteration number corresponding to the current iteration step;

[0207] T represents the total loading step corresponding to;

[0208] σ (k),(t) represents the stress corresponding to element k in the current iteration step;

[0209] σ (k),(T) represents the stress of element k corresponding to the total loading step T, i.e. the stress corresponding to element k in the initial stress state;

[0210] s (k),(T) represents the state variable of element k corresponding to the total loading step T;

[0211] △ε (k),(t) represents the element strain increment corresponding to element k in the current iteration step;

[0212] △σ(σ (k),(T) ,s (k),(T) ,△ε (k),(t) ) represents the stress increment calculated by stress integration with σ (k),(T) and s (k),(T) as the initial state, and the strain path △ε (k),(t) ;

[0213] S362, the formula for calculating the state variable of the intermediate stress state is:

[0214] s (k),(t) = s (k),(T)) +△s(σ (k),(T) ,s (k),(T) ,△ε (k),(t) )

[0215] Wherein:

[0216] s (k),(t) represents the state variable corresponding to element k in the current iteration step;

[0217] △s(σ (k),(T) ,s (k),(T) ,△ε (k),(t) ) represents the state increment calculated by stress integration with σ (k),(T) and s (k),(T) as the initial state, and the strain path △ε (k),(t) ;

[0218] The present embodiment adopts full displacement for dynamic relaxation solution, i.e. strain stress integration is performed based on the initial stress state in each iteration of the sub-loading step. Since the starting point of each integration is fixed as the end point of the previous sub-loading step, the influence of the false stress path (usually uncontrollable cyclic loading) caused by iteration of the dynamic relaxation method on the soil stress and state variable can be effectively avoided.

[0219] S370, performing non-equilibrium force calculation based on the intermediate stress state and the external force data to obtain non-equilibrium forces corresponding to each node in the current iteration step;

[0220] The external force data includes the node external force of each node corresponding to the total loading step T calculated in advance

[0221] The specific steps are:

[0222] S371, calculating the force contribution of each node in the current iteration step based on the intermediate stress state

[0223] The calculation formula is:

[0224]

[0225] Wherein:

[0226] σ (k),(t) Indicates the stress corresponding to the unit k in the current iteration step;

[0227] n (k,l),(T) Indicates the normal vector of surface (k, l) in the total loading step T;

[0228] S (k,l),(T) Indicates the area of surface (k, l) in the total loading step T.

[0229] Note:

[0230] When calculating the initial non-equilibrium force, that is, when calculating the non-equilibrium force based on t = 0, σ (k),(t) Take σ (k),(T) .

[0231] S372, accumulating the force contribution to obtain the node internal force of each node

[0232] The calculation formula is:

[0233]

[0234] S373, calculating the non-equilibrium force corresponding to each node based on the node external force and the node internal force

[0235] The calculation formula is:

[0236]

[0237] S380, judging whether the preset iteration end condition is reached;

[0238] S381, The iteration accuracy reaches the preset allowable value for iteration accuracy;

[0239] In this embodiment, the iteration accuracy is the maximum value of the first unbalanced force. The second unbalanced force maximum value ratio when When the value is less than the allowable iteration precision TOL, the iteration is considered complete and the iteration solution is considered successful.

[0240] The first unbalanced force maximum value This is the maximum value among the initial unbalanced forces corresponding to each node in the current sub-loading step. That is, based on t=0, following the above calculation steps for unbalanced forces, the maximum value is selected from the calculated initial unbalanced forces to obtain the first unbalanced force maximum value.

[0241] The second maximum value of the unbalanced force is the maximum value of the unbalanced force at each node in the current iteration step.

[0242] S382, The number of iterations has reached the preset maximum number of iterations.

[0243] When the number of iterations t reaches the preset maximum value MAXITER, the iteration is considered complete and the iteration solution fails.

[0244] S383. Based on the judgment result, complete the sub-loading step or continue iterating;

[0245] When the determination iteration ends, the node displacement is taken as the target displacement data, and the obtained intermediate stress state is taken as the target stress state.

[0246] When it is determined to continue iterative calculation, the unbalanced force is taken as the historical unbalanced force, the node velocity is taken as the historical node velocity, and the node displacement is taken as the historical node displacement for use in the next iteration.

[0247] While adjusting the algorithm architecture and eliminating spurious stress paths, this embodiment still uses the "small steps, fast pace" strategy of the dynamic relaxation method in the core iterative process of its sub-loading steps, thus retaining the advantages of the conventional dynamic relaxation method for solving strongly nonlinear problems.

[0248] S400, Loading is graded based on the triggered iteration termination condition;

[0249] Specifically;

[0250] S410. When the iteration ends based on the iteration precision:

[0251] S411. Update the starting position information based on the target displacement data;

[0252] The node displacement u obtained by iterative calculation in the embodiment (l),(t) The displacement field u corresponding to the current total loading step T (l),(T) is accumulated to obtain the displacement field u corresponding to the next total loading step (l),(T+△T) , and the node displacement u (l),(t) The spatial position x of each node in the starting position data of the current total loading step T (l),(T) is updated to obtain the spatial position x of each node corresponding to the next loading step T+△T (l),(T+△T) , and the specific formula is:

[0253] u (l),(T+△T) =u (l),(T) +u (l),(t) ;

[0254] x (l),(T+△T) =x (l),(T) +u (l),(t) .

[0255] Note: After the iteration is completed, when the data is updated based on the data obtained by iterative calculation, the iteration step t is the iteration step at which the iteration is completed, that is, the last iteration step.

[0256] S412, updating the starting stress state based on the target stress state;

[0257] σ (k) , (T+△T) =σ (k),(t) ;

[0258] s (k),(T+△T) =s (k),(t) .

[0259] S413, accumulating the total loading step based on the loading step increment;

[0260] That is, let T+△T be the next total loading step.

[0261] S414, adjusting the loading step increment based on a preset adaptive rule;

[0262] The specific steps for adjusting the loading step increment when the iteration is completed based on the iteration accuracy are as follows:

[0263] Obtain the iteration number t corresponding to the end of iteration;

[0264] Generate an adjustment weight based on the iteration number t and the maximum iteration number MAXITER

[0265] Weighted calculation of the loading step increment △T based on the adjustment weight obtains a first candidate increment

[0266] based on the preset maximum total loading step T max and the difference between the updated total loading step and the preset maximum total loading step to generate a second candidate increment (T max -T), where T max is 1 in this embodiment.

[0267] the smaller one of the first candidate increment and the second candidate increment (T max -T) is taken as the new loading step increment.

[0268] S415, determining whether the calculation is completed based on the updated total loading step.

[0269] That is, when T≥T max , it is determined that the dynamic relaxation numerical calculation is completed, and the calculation is ended and the corresponding result data is output.

[0270] In this embodiment, the initial loading step increment and the preset maximum total loading step T max have the same value, that is, when the initial entering of the sub-loading step for iterative solving is successful, the calculation is directly completed, otherwise, the end determination is based on the updated total loading step after the successful solving.

[0271] S420, when the iteration is ended based on the iteration number, the loading step increment is adjusted based on a preset adaptive rule.

[0272] The specific steps of adjusting the loading step increment when the iteration is ended based on the iteration accuracy are as follows:

[0273] a third candidate increment is generated based on the loading step increment, and the third candidate increment is smaller than the loading step increment, and in this embodiment, the third candidate increment is ΔT / 2;

[0274] a fourth candidate increment (T max -T) is generated based on the preset maximum total loading step and the difference between the total loading step and the preset maximum total loading step;

[0275] the smaller one of the third candidate increment ΔT / 2 and the fourth candidate increment (T max -T) is taken as the new loading step increment.

[0276] In this embodiment, the hierarchical loading control is performed by the iteration accuracy allowable value and the maximum iteration number, that is, when the iteration is ended based on the iteration accuracy allowable value, it is determined that the solving is successful, at this time, the data is updated based on the data obtained by the iteration calculation, and the total loading step is updated based on the corresponding loading step increment; when the iteration is ended based on the maximum iteration number, it is determined that the solving fails, at this time, only the loading step increment is adjusted, and the sub-loading step is entered based on the total loading step, the starting position information and the starting stress state for iterative solving.

[0277] The main disadvantage of the conventional dynamic relaxation method is the existence of false stress path interference. Since the iterative solution is carried out in a dynamic manner, each calculation unit will experience multiple uncontrollable loading and unloading processes in the adjustment of the displacement field and the stress update calculation. The above loading and unloading processes are actually caused by iteration, rather than the real stress path of the material. The above false stress path has relatively small influence on the rock-soil constitutive model under the elastic-plastic framework, and the iteration process mainly occurs in the elastic region within the plastic yield surface. However, for complex rock-soil constitutive models that accumulate strain under cyclic loading, such as sub-plasticity models, cyclic loading boundary surface models and other strong stress history related models, the deformation of the soil body will be overestimated.

[0278] In the present application, this defect will be solved by splitting the loading step and using full displacement for strain stress calculation in the stress calculation part of the iteration cycle of the above sub-loading step, avoiding real-time update of stress and strain and thus introducing false stress path.

[0279] For the problem of large amount of calculation of the dynamic relaxation method, the present application uses the method of GPU high-performance numerical calculation to solve it. In the present embodiment, the calculation and update of each unit variable and each node variable are accelerated by GPU.

[0280] In the present embodiment, the node variables and the unit variables in each step are stored in the global variables of the GPU, so that in the core steps of data calculation and data update, each thread executes the same calculation instruction and operates on the continuous memory area for concurrent calculation.

[0281] The data calculation and data update, for example, include the following calculations:

[0282] 1) Calculation of geometric information in step S220, specifically, calculation of volume V (k),(T) , area S (k,l),(T) , and normal vector n (k ; ,l),(T)

[0283] 2) Calculation of node virtual mass contribution m (k,l),(T) in step S231;

[0284] 3) Calculation of node force contribution , including calculation in each iteration step and initial value calculation at t=0 moment;

[0285] 4) Calculation of non-equilibrium force of each node , including calculation in each iteration step and initial value calculation at t=0 moment;

[0286] 5) The nodal velocity v at time t=0 (l),(0) and nodal displacement u (l),(0) Initialize to 0;

[0287] 6) In step S330, the node velocity v (l),(t) and nodal displacement u (l),(t) Iterative calculation;

[0288] 7) The nodal velocity v of the displacement-constrained nodes in step S340 (l),(t) and nodal displacement u (l),(t) Zeroing out;

[0289] 8) In step S350, the element strain increment Δε (k),(t) Calculation;

[0290] 9) In step S360, the stress σ (k),(t) and state variable s (k),(t) Calculation;

[0291] 10) In step S410, the displacement field u (l),(T) Node spatial location x (l),(T) Stress σ (k),(T) State variable s (k),(T) Update.

[0292] In this embodiment, the atomic operations of the GPU are used in the process of accumulating and finding the maximum value to achieve concurrent operations of each thread on the same video memory address;

[0293] Calculating the maximum value of a sum includes, for example:

[0294] 1) In step S240, the surface load q is added using atomic addition. (bs) S (bs),(T) and gravity load ρ (k) gV (k),(T) Discretely accumulate the external forces at the connected nodes

[0295] 2) In step S232, the node virtual mass contribution m is calculated using atomic addition. (k,l),(T) Accumulated virtual mass M of computing nodes (l),(T) ;

[0296] 3) Use atomic addition operations to contribute nodal forces Accumulate the internal forces acting on the nodes This includes the calculation of nodal forces in each iteration step and the calculation of initial values ​​at time t=0;

[0297] 4) Using the atomic maximum value operation, find the maximum unbalanced force from the unbalanced forces at each node. or

[0298] To address the issue of the large computational load associated with the "small steps, fast pace" dynamic relaxation method, this embodiment utilizes the GPU approach for the core computations of the dynamic relaxation iteration process, taking advantage of its explicit iterative process and the characteristics of the GPU instruction set. This significantly improves the computational efficiency of the method.

[0299] The following is a detailed description of a dynamic relaxation numerical calculation method in this embodiment through a specific case.

[0300] In this case, the main variables that need to be allocated in GPU memory fall into two categories: ① node variables and ② cell variables. It should be noted that arrays in the GPU must be allocated as one-dimensional arrays, located in global memory, and arranged in column-major order. For example, let the total number of nodes be N. nd For node position pos[N nd *3] array, the addressing modes for accessing the x, y, and z coordinates of the l-th node are pos[l], pos[l+N], and pos[l+N], respectively. nd ] and pos[l+2*N nd This ensures that adjacent threads can strictly access contiguous address spaces during computation, thus improving access bandwidth through merged access.

[0301] Consider a total number of nodes of N nd The correspondence between the node variables that need to be applied for and the symbols in Example 1 is shown in the table below.

[0302] Table 1

[0303]

[0304]

[0305] The vast majority of the variables are vectors in three-dimensional space, therefore the array length is N. nd *3. Pay special attention to the array fixFlag[N] nd It contains the displacement boundary conditions used in the calculation, and is a boundary condition of length N. nd An integer array. The value range of this variable is 0 to 7, and its meaning is as follows:

[0306] 0: Indicates no constraint;

[0307] 1: Constrained only in the x-direction;

[0308] 2: Constrained only in the y-direction;

[0309] 3: Constrained only in the z-direction;

[0310] 4: constrained in both y and z directions;

[0311] 5: constrained in both z and x directions;

[0312] 6: constrained in both x and y directions;

[0313] 7: constrained in all x, y, z directions.

[0314] In the case where the displacement boundary condition is clear, the node velocity v (l),(t) and the node displacement u (l),(t) (i.e., the displacement within the sub-loading step) are set to 0 in the component of the constrained direction in step S340.

[0315] Consider that the total number of tetrahedral elements is N el , the required element variables to be applied correspond to the symbol correspondence in Embodiment 1 as shown in the following table.

[0316] Table 2

[0317]

[0318] Firstly, the connectivity relationship between nodes and elements needs to be considered, which is defined in the element layer in the program, i.e., the number of 4 nodes of each element is defined and stored in lkPts[N el *4];

[0319] The stress variable is a second-order symmetric tensor in three-dimensional space, which has 9 components of xx, xy, xz, yx, yy, yz, zx, zy, and zz; considering the symmetry, xy and yx, xz and zx, and yz and zy are equal to each other, so only 6 floating points are needed for the stress of each element (elSig and calElSig), i.e., stored according to xx, yy, zz, xy, xz, and yz. This storage form is commonly referred to as the Viogt form of the tensor.

[0320] The state variable and the material parameter are related to the selected constitutive model.

[0321] N sta in statev and calstatev is the number of internal variables defined in the model, and N prop in matProps is the number of material parameters defined in the model. For the geometry of the element surface and the external load, it is stored in 4 arrays according to face 1, 2, 3, and 4. For example, assuming that the area, normal vector, and surface load of the 2nd face of the kth element need to be accessed, the array elements to be accessed are:

[0322] Area: elSurf2[k]

[0323] Normal vector: (elSurf2[k+N el ], elSurf2[k+2*N el ], elSurf2[k+3*N el ])

[0324] Surface load: (qsSurf2[k], qsSurf2[k+N el ], qsSurf2[k+2*N el ]).

[0325] Most of the variables and constants used to control the calculation process in the algorithm are stored in the CPU, and the main loop steps and sub-loop steps are executed in the CPU.

[0326] The variables and constants used to control the calculation process are shown in the following table:

[0327] Table 3

[0328]

[0329] The core calculation process of the algorithm needs to be deployed in the GPU for execution.

[0330] This case gives a recommended deployment method, and the corresponding relationship between the functions of each kernel function and the steps and formulas in the application is as follows:

[0331] Kernel function gpu_calc_geom: This function is completely corresponding to the function of step S220, and the total number of threads is equal to the total number of tetrahedrons N el . Each thread obtains the node number of the tetrahedron it processes from the array lkPts in parallel, and obtains the spatial coordinates of the four nodes from the array loc according to the node number, calculates and updates the unit volume elVol array and the area and normal vector elSurf1-elSurf4 array.

[0332] Kernel function gpu_init_massfext: This function is used to initialize the dampMass array of virtual mass before executing the accumulation process of step S230, and also used to initialize the extForce array before executing the accumulation process of step S240, and the total number of threads is equal to the total number of nodes N nd . Each thread tid assigns dampMass[tid] to 0 in parallel. Initialize the node external force as the external force directly borne by the node, extForce[tid] = fPt[tid], extForce[tid+N nd ] = fPt[tid+N nd], extForce[tid + 2 * N nd ] = fPt[tid + 2 * N nd ].

[0333] Kernel function gpu_calc_massfext: This function is used to accumulate the dampMass array of virtual mass in step S230, and to accumulate the extForce array of external force in step S240. The total number of threads is equal to the total number of tetrahedrons N el . Each thread corresponds to each tetrahedron element, and calculates the external force contribution of the element to the node or calculates the mass contribution m (k,l),(T) . Then, the corresponding 4 nodes are read according to the lkPts array, and are accumulated to the node external force or virtual mass corresponding to the 4 nodes. The accumulation process must perform atomic operation atomic_add on the corresponding positions in the arrays extForce and dampMass to avoid writing of different threads to the same address.

[0334] Kernel function gpu_init_funbal: This function is used to clear the unBal array of unbalanced force before calculating the node action internal force (i.e., in the initial stage and steps S371 and S372 in each iteration step). The total number of threads is equal to the total number of nodes N nd .

[0335] Kernel function gpu_calc_fint: This function is used to accumulate the calculated element stress contribution and node action internal force (i.e., in the initial stage and steps S371 and S372 in each iteration step). The total number of threads is equal to the total number of tetrahedrons N el . Each thread corresponds to each tetrahedron element, and obtains the stress state in the current iteration step from calElsig, and calculates the element stress contribution . Then, the accumulation to the 4 nodes of the tetrahedron must perform atomic operation atomic_add on the corresponding positions in the array unBal to avoid writing of different threads to the same address. After the execution of this function, the node action internal force has been accumulated in unBal.

[0336] Kernel function gpu_calc_funbal: This function is used to perform node unbalanced force calculation (in the initial stage and corresponding step S373 in each iteration step). The total number of threads is equal to the total number of nodes N ndEach thread corresponds to each node, and adds the value in the extForce array to the corresponding value in the unBal array. After each thread finishes calculating unBal, it further calculates the absolute value of unBal, |unBal|, and compares it with the global maxUnBal. The atomic_max operation is used to write the larger value to the maxUnBal position.

[0337] Kernel function gpu_init_veldisp: This function is used to clear the velocity vel array and the displacement calDisp array in the sub-loading step before entering the sub-loading step, and the total number of threads is equal to the total number of nodes N nd .

[0338] Kernel function gpu_calc_veldisp: This function corresponds to the function of step S330 completely, and the total number of threads is equal to the total number of nodes N nd . Each thread corresponds to each node, and the velocity and displacement in the sub-loading step are calculated, and the vel array and the calDisp array are calculated.

[0339] Kernel function gpu_cal_fixvel: This function corresponds to the function of step S340 completely, and the total number of threads is equal to the total number of nodes N nd . Each thread corresponds to each node, reads the boundary condition constraint situation in the fixFlag array, and operates the vel array and the calDisp array according to different constraint situations.

[0340] Kernel function gpu_calc_calsig: This function corresponds to the functions of steps S350 and S360 completely, and the total number of threads is equal to the total number of tetrahedrons N el . Each thread corresponds to each tetrahedral unit, first obtains the 4 connected nodes from lkPts, reads the sub-cycle step displacement of the 4 nodes from calDisp, calculates the unit strain increment. Read the stress corresponding to the initial stress state and the state variable from elSig and statev, read the constitutive model parameters of the material from matProps, based on the calculated strain increment and the specific constitutive model, perform stress integration calculation, and write to calElSig and calstatev.

[0341] Kernel function gpu_update_veldisp: This function corresponds to the displacement field and node position update in step S411, and the total number of threads is equal to the total number of nodes N nd . Each thread corresponds to each node, reads the value in calDisp, and adds it to the disp array and the loc array.

[0342] Kernel function gpu_update_sig: This function corresponds to the stress and state variable update in step S412, and the total number of threads is equal to the total number of tetrahedrons N el Each thread corresponds to each tetrahedral unit, reads the values in calElSig and calStatev, and assigns them to the elSig and statev arrays.

[0343] The algorithm executes the control flow on the CPU and calls the kernel function on the GPU to perform the steps with large computational loads. In order to more clearly show the calling logic relationship of the functions, the core calculation loop in the dynamic relaxation numerical calculation method corresponding to embodiment 1 can be summarized as the following execution pseudo code. The following functions with the gpu prefix mean that the kernel function on the GPU is called by the CPU, and other control flows are completed on the CPU. The functions executed on the GPU are all floating point operations according to the formulas in the present application, and the specific code is not repeated.

[0344] The control code is as follows:

[0345]

[0346]

[0347] The present application is applied to the calculation of solving a certain slope gravity field, the material adopts the Drunker-Prager model, the slope has 10703 nodes and 9120 units. The calculation uses i7-7700, AMD R7 240 and Nvidia RTX 1050, and the calculation is compared with the commercial software FLAC3D to check the accuracy of the calculation and the acceleration effect of various computing devices. Figure 3 The calculation results of FLAC3D are, Figure 4 The results obtained based on the present application are visualized by using Paraview. The horizontal displacement at the slope angle is 1.19mm (FLAC3D) and 1.16mm (the present application), respectively, and the calculation results are basically consistent.

[0348] The calculation time is shown in the following table;

[0349] Table 4

[0350] Software type Commercial software The invention The invention The invention Computing platform CPU CPU GPU GPU Model i7-7700 i7-7700 AMD R7 240 Nvidia RTX 1050 Number of threads 4 4 384 640 Frequency 3.6 GHz 3.6 GHz 780 MHz 1354 MHz Computing time 9.8 min 2.25 min 0.67 min 0.69 min

[0351] Among them, even if all on the CPU, the calculation time of FLAC 3D (9.8min) is much larger than the self-developed program (2.25min). This is because FLAC 3D more considers the robustness of the program, and also needs to consider the needs of visual rendering, and the calculation efficiency is reduced.

[0352] In addition, the scheme proposed in the present application is significantly different from the CPU on the GPU platform. Even on a relatively low-end R240, more than 3 times acceleration effect (0.67 min) is obtained, and on an Nvidia RTX 1050, 3.26 times acceleration effect (0.69 min) is obtained.

[0353] The above results prove the effectiveness of the GPU-accelerated dynamic relaxation algorithm.

[0354] The above is only a preferred embodiment of the present application, and does not limit the present application in any form. Any simple modification or equivalent change made on the basis of the technical essence of the present application to the above embodiment falls within the protection scope of the present application.

[0355] Embodiment 2, a simulation system for simulating and analyzing a target rock-soil constitutive model based on a target external force, comprising:

[0356] The acquisition module is configured to acquire a total loading step corresponding to the current and a loading step increment;

[0357] The calculation module is configured to acquire starting stress state and starting position information corresponding to the total loading step, and to calculate corresponding action external force data based on the total loading step and the loading step increment;

[0358] The iterative solution module is configured to enter a sub-loading step, to perform dynamic relaxation iterative solution based on the starting stress state and the action external force data using full displacement, until a preset iteration end condition is reached, to obtain corresponding target displacement data and target stress state; wherein the iteration end condition is that the iteration accuracy reaches a preset iteration accuracy allowable value, or the iteration number reaches a preset maximum iteration number;

[0359] The update module is configured to perform loading grading based on the triggered iteration end condition, and comprises a first update unit and a second update unit;

[0360] The first update unit is configured to perform the following steps when the iteration is ended based on the iteration accuracy:

[0361] updating the starting position information based on the target displacement data;

[0362] updating the starting stress state based on the target stress state;

[0363] accumulating the total loading step based on the loading step increment;

[0364] adjusting the loading step increment based on a preset adaptive rule;

[0365] Further, it is judged whether the calculation is completed based on the updated total loading step;

[0366] The second updating unit is configured to adjust the loading step increment based on a preset adaptive rule when the iteration is ended based on the iteration number.

[0367] This embodiment is a device embodiment corresponding to the embodiment 1, and the details can be referred to the embodiment 1.

[0368] Those skilled in the art will understand that the embodiments of the present application can be provided as a method, device, or computer program product. Therefore, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present application can take the form of a computer program product implemented 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.

[0369] The present application is described with reference to flowcharts and / or block diagrams of the method, terminal device (system), and computer program product according to the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams, and the combination of flows and / or blocks in the flowcharts and / or block diagrams can be implemented by computer program instructions. These computer program instructions can be provided to a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing terminal device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing terminal device produce a device that implements the functions specified in the flowcharts and / or block diagrams. Figure 1 one or more flows and / or blocks Figure 1 an apparatus that performs the functions specified in one or more blocks.

[0370] These computer program instructions can also be stored in a computer-readable memory that can direct the computer or other programmable data processing terminal device to work in a specific manner, so that the instructions stored in the computer-readable memory produce a manufactured product including instruction apparatus, which implements the functions specified in the flowcharts and / or block diagrams. Figure 1 one or more flows and / or blocks Figure 1 an apparatus that performs the functions specified in one or more blocks.

[0371] These computer program instructions can also be loaded into a computer or other programmable data processing terminal device, so that a series of operation steps are performed on the computer or other programmable terminal device to produce a computer-implemented process, so that the instructions executed on the computer or other programmable terminal device provide a process for implementing the functions specified in the flowcharts and / or block diagrams. Figure 1 one or more flows and / or blocks Figure 1 an apparatus that performs the functions specified in one or more blocks.

Claims

1. A dynamic relaxation numerical calculation method for simulation analysis of a target geotechnical constitutive model based on a target external force, characterized in that, The method comprises the following steps: obtaining a total loading step corresponding to the current and a loading step increment; obtaining a starting stress state and starting position information corresponding to the total loading step; obtaining virtual mass information corresponding to the total loading step; based on the total loading step and the loading step increment, calculating corresponding external force data; entering a sub-loading step, based on the starting stress state, the virtual mass information and the external force data, using full-quantity displacement to perform dynamic relaxation iteration until a preset iteration end condition is reached, to obtain corresponding target displacement data and a target stress state; wherein the iteration end condition is that the iteration accuracy reaches a preset iteration accuracy allowable value, or the iteration number reaches a preset maximum iteration number; based on the triggered iteration end condition, performing loading grading, specifically: when the iteration ends based on the iteration accuracy: updating the starting position information based on the target displacement data; updating the starting stress state based on the target stress state; accumulating the total loading step based on the loading step increment; adjusting the loading step increment based on a preset adaptive rule; judging whether the calculation is completed based on the updated total loading step; when the iteration ends based on the iteration number, adjusting the loading step increment based on a preset adaptive rule; wherein the iteration process in the sub-loading step comprises the following steps: obtaining historical unbalanced forces, historical node velocities and historical node displacements corresponding to each node in the last iteration step; based on the virtual mass information, the historical unbalanced forces, the historical node velocities and the historical node displacements, calculating node velocities and node displacements corresponding to each node in the current iteration step; for nodes with displacement constraint boundary conditions, setting the components of the node velocities and the node displacements in the constrained directions to 0; based on the node displacements, calculating element strain increments corresponding to each element; based on the element strain increments, performing element stress integration with the starting stress state as the initial value to obtain a corresponding intermediate stress state; based on the intermediate stress state and the external force data, performing unbalanced force calculation to obtain unbalanced forces corresponding to each node in the current iteration step; judging whether the preset iteration end condition is reached, when determining that the iteration is ended, taking the node displacements as the target displacement data and taking the obtained intermediate stress state as the target stress state, when determining that the iteration is continued, taking the unbalanced forces as the historical unbalanced forces, taking the node velocities as the historical node velocities and taking the node displacements as the historical node displacements for use in the next iteration.

2. The dynamic relaxation numerical calculation method according to claim 1, wherein: GPU acceleration is used for calculation and updating of each element variable and each node variable.

3. The dynamic relaxation numerical calculation method according to claim 1, wherein: the stress state comprises stress and state variables; a formula for calculating the stress of the intermediate stress state is: σ (k),(t) =σ (k),(T)) +△σ(σ (k),(T) , s (k),(T) ,△ε (k),(t) ） ; wherein: k represents the 1 k st tetrahedral unit; t represents the iteration number corresponding to the current iteration step; T represents the total loading step to which it corresponds; Δε (k),(t) represents the cell in the current iteration step k corresponding stress; Δσ (k),(T) denotes the total loading step T the stress corresponding to the cell k the stress corresponding to the cell k in the initial stress state; s (k),(T) represents the total load step T the corresponding cell k state variable; Δε (k),(t) represents the cell in the current, iteration step k corresponding cell strain increment; Δσ (k),(T) , s (k),(T) a formula for calculating the state variable of the intermediate stress state is: (k),(t) ） representing the initial state Δs (k),(T) and s (k),(T) stress integral with initial state wherein: (k),(t) stress increment Δs obtained by stress integral calculation Δε s (k),(t) = s (k),(T)) + Δs (k),(T) , s (k),(T) Δε (k),(t) ）； ​ s (k),(t) represents the state variable corresponding to the cell in the current iteration step k in the current iteration step ​ (k),(T) , s (k),(T) ​ (k),(t) ） represents the initial state ​ (k),(T) and s (k),(T) the stress integral for the initial state ​ (k),(t) the state increment △s obtained by the stress integral calculation.

4. The dynamic relaxation numerical calculation method according to claim 3, characterized in that, The calculation formula of the unit strain increment is: The calculation formula of the unit strain increment is: ； Wherein: △ε (k),(t) represents the current iteration step unit k the corresponding unit strain increment; x, y, z represents three dimensions; subscript i or j for second order tensors △ε (k),(t) components in the xx, xy, xz, yx, yy, yz, zx, zy, zz nine dimensions; (k,l) representing unit k node l opposite surface represents the summation over each tetrahedron element k for its corresponding node and corresponding surface. V (k),(T) represents the total loading step T the corresponding unit k volume; u (l),(t) represents the node in the current iteration step l corresponding node displacement; n (k,l),(T) representing the total load step T normal vector of the median plane (k,l) of the median plane; S (k,l),(T) Indicates total loading steps T middle noodle (k,l) The area.

5. The dynamic relaxation numerical calculation method according to claim 1, characterized in that, The specific steps of non-equilibrium force calculation based on the intermediate stress state and the acting external force data are: The acting external force data includes the node acting external force of each node. Based on the intermediate stress state, the force contribution applied to the node by the stress of each unit in the current iteration step is calculated. The force contributions are accumulated to obtain the node acting internal force of each node. Based on the node acting external force and the node acting internal force, the non-equilibrium force corresponding to each node is calculated and obtained.

6. The dynamic relaxation numerical calculation method according to claim 5, wherein: The iteration accuracy is the ratio of the first non-equilibrium force maximum value to the second non-equilibrium force maximum value; The first non-equilibrium force maximum value is the initial non-equilibrium force maximum value corresponding to the current sub-loading step; The second non-equilibrium force maximum value is the maximum value of the non-equilibrium force of each node in the current iteration step.

7. The dynamic relaxation numerical calculation method according to any one of claims 1-6, wherein: The acting external force data includes the node acting external force of each node. The calculation formula of the node acting external force is: ； Wherein: bs∈l faces participating in the accumulation bs faces connected to the node l faces connected to the node k∈l a unit that indicates to participate in accumulation k a unit that is connected to a node l a unit that is connected to a node T represents the total loading step; △T represents a loading step increment; represents the total loading step T In the middle, the node l corresponding to the node role external force; node l point load q (bs) face bs face load S (bs),(T) for the area of bs the area; ρ (k) of the unit k density, obtained based on a target rock-soil constitutive model; g g is the gravitational acceleration; V (k),(T) is the volume of the cell k of the cell.

8. The dynamic relaxation numerical calculation method according to any one of claims 1 to 6, characterized in that, The adaptive rule for adjusting the loading step increment is specifically: When the iteration is ended based on the iteration accuracy: The iteration number corresponding to the end of iteration is obtained; An adjustment weight is generated based on the iteration number and the maximum iteration number; The loading step increment is weighted calculated based on the adjustment weight to obtain a first candidate increment; A second candidate increment is generated based on the difference between the preset maximum total loading step and the updated total loading step; The smaller value between the first candidate increment and the second candidate increment is taken as the new loading step increment; When the iteration is ended based on the iteration number: A third candidate increment is generated based on the loading step increment, and the third candidate increment is smaller than the loading step increment; A fourth candidate increment is generated based on the difference between the preset maximum total loading step and the total loading step; The smaller value between the third candidate increment and the fourth candidate increment is taken as the new loading step increment.

9. A simulation system for simulating and analyzing a target soil-rock constitutive model based on target external forces, characterized in that, It includes: An acquisition module is configured to acquire the total loading step and the loading step increment corresponding to the current iteration step; A calculation module is configured to acquire the starting stress state and the starting position information corresponding to the total loading step, and is further configured to acquire the virtual mass information corresponding to the total loading step, and is further configured to calculate the corresponding acting external force data based on the total loading step and the loading step increment; An iteration solving module is configured to enter a sub-loading step, perform dynamic relaxation iteration solving based on the starting stress state and the acting external force data using full-quantity displacement until a preset iteration end condition is reached, and obtain corresponding target displacement data and target stress state; wherein, the iteration end condition is that the iteration accuracy reaches a preset iteration accuracy allowed value, or the iteration number reaches a preset maximum iteration number; An update module is configured to perform loading grading based on the triggered iteration end condition, and includes a first update unit and a second update unit; The first update unit is configured to perform the following steps when the iteration is ended based on the iteration accuracy: updating the initial position information based on the target displacement data; updating the initial stress state based on the target stress state; accumulating the total loading step based on the loading step increment; adjusting the loading step increment based on a preset adaptive rule; judging whether the calculation is completed based on the updated total loading step; adjusting the loading step increment based on a preset adaptive rule when the iteration is ended based on the iteration number; wherein the iteration process in the sub-loading step comprises the following steps: obtaining historical unbalanced forces, historical node velocities and historical node displacements corresponding to each node in the last iteration step; calculating node velocities and node displacements corresponding to each node in the current iteration step based on the virtual mass information, the historical unbalanced forces, the historical node velocities and the historical node displacements; setting the components of the node velocities and the node displacements in the constrained direction to 0 for the nodes with displacement constraint boundary conditions; calculating element strain increments corresponding to each element based on the node displacements; performing element stress integration with the initial value of the initial stress state based on the element strain increments to obtain a corresponding intermediate stress state; calculating unbalanced forces corresponding to each node in the current iteration step based on the intermediate stress state and the external force data; judging whether a preset iteration end condition is reached, when determining that the iteration is ended, taking the node displacements as target displacement data and taking the obtained intermediate stress state as a target stress state, when determining that the iteration is continued, taking the unbalanced forces as historical unbalanced forces, taking the node velocities as historical node velocities and taking the node displacements as historical node displacements for the next iteration.

Images