Method for rapidly and stably solving buckling path after multiple critical points based on generalized displacement control method

Abstract

The invention discloses a method for rapidly and stably solving a buckling path after multiple critical points based on a generalized displacement control method. The method comprises the steps that S1, a state variable and an anchoring type secant predictor are initialized; s2, predicting a reference displacement, a calibration displacement and a load increment factor based on the state variable and the anchored secant predictor; s3, calculating the displacement increment of the current iteration step based on the reference displacement, the calibration displacement and the load increment factor, and updating the total displacement and the total load of the structure; s4, updating the geometric state of the structure and the internal force of the unit based on the displacement increment of the current iteration step; s5, calculating the difference value between the total internal force and the external load of the structure based on the updated geometric state of the structure and the unit internal force to obtain the unbalanced force of the structure, and performing error checking based on the unbalanced force of the structure; s6, executing the steps S1 to S5 circularly until the unbalanced force of the structure meets the balance condition, and starting the next increment step; and incremental analysis is executed until user requirements are met, and a buckling path after multiple critical points is output.

Description

Technical Field

[0001] The present invention belongs to the technical field of nonlinear finite element and structural nonlinear analysis, and in particular relates to a method for quickly and stably solving post-buckling paths of multiple critical points based on a generalized displacement control method. Background Art

[0002] Over the past few decades, with advancements in material preparation technology, engineers have placed their hopes on integrating high-performance materials into structures, achieving lightweight structures while simultaneously improving their load-bearing capacity and durability. However, when such thin-walled components are subjected to abnormal or extreme loads, they often experience considerable nonlinear deformation due to structural instability. To meet the current demands for thin-walled structural design in various engineering projects, such as high-rise buildings, aerospace, mechanical manufacturing, and bridge structures, many design codes require stability analysis as part of structural design practice.

[0003] Traditionally, stability analysis has only been applied to some simple structures, and for complex structures, approximate techniques are mostly used. In order to make the designed structure reliable, the design specifications require a higher safety factor, which greatly increases the consumption of materials. Today, with the development of finite element technology, engineers are able to understand the mechanical response of the structure more realistically. As my country's infrastructure construction is shifting from the incremental construction stage to the high-quality development stage, the requirements for structural safety are constantly increasing, and structural analysis is becoming more refined. Therefore, it is urgent to study more robust, accurate, and efficient stability analysis methods. Summary of the Invention

[0004] In order to solve the above technical problems, the present invention proposes a method for quickly and stably solving the post-buckling path of multiple critical points based on the generalized displacement control method to solve the problems existing in the above-mentioned prior art.

[0005] To achieve the above objectives, the present invention provides a method for quickly and stably solving post-buckling paths at multiple critical points based on a generalized displacement control method, comprising:

[0006] S1, initialize state variables and anchored secant predictor;

[0007] S2. Predicting a reference displacement, a calibration displacement, and a load increment factor based on the state variable and the anchored secant predictor;

[0008] S3. Calculating the displacement increment of the current iteration step based on the reference displacement, calibration displacement and load increment factor, and updating the total displacement and total load of the structure;

[0009] S4. updating the structural geometric state and unit internal force based on the displacement increment of the current iteration step;

[0010] S5. Calculate the difference between the total internal force of the structure and the external load based on the updated structural geometric state and the unit internal force to obtain a structural unbalanced force, and perform an error check based on the structural unbalanced force;

[0011] S6. Execute S1-S5 cyclically until the structural unbalanced force satisfies the equilibrium condition and the next incremental step begins. When the incremental analysis is performed until the structural load or displacement meets the user's requirements, the multiple critical point post-buckling path is output.

[0012] Optionally, it is characterized in that the state variables include: node coordinates, unit coordinate basis, node displacement, unit internal force, external load, load increment factor and unit stiffness matrix.

[0013] Optionally, the process of predicting reference displacements, calibration displacements, and load increment factors includes:

[0014] Predicting a displacement increment of a current incremental step based on the state variable and the anchored secant predictor, wherein the initial displacement increment includes: a reference displacement and a calibration displacement;

[0015] The load increment factor of the current iteration step is calculated using a generalized displacement control method based on the reference displacement and the calibration displacement.

[0016] Optionally, the process of predicting the displacement increment for the current increment step includes:

[0017] If the current iteration step is 1, the reference displacement and calibration displacement are calculated based on the anchored secant predictor; wherein the expression for calculating the reference displacement and calibration displacement based on the anchored secant predictor is:

[0018]

[0019] If the current iteration step is greater than 1, the reference displacement and calibration displacement are calculated based on the definition; wherein the expression for calculating the reference displacement and calibration displacement based on the definition is:

[0020]

[0021] Where, Indicates the reference displacement when the current iteration step is 1, Indicates the calibration displacement when the current iteration step is 1, {ΔU i-1} represents the displacement increment of each iteration of the i-1th incremental step The cumulative value of λ i-1 Indicates the incremental load factor for each iteration of the i-1th incremental step The cumulative value of Indicates the reference displacement when the current iteration step is greater than 1, Indicates the calibration displacement when the current iteration step is greater than 1, represents the structural stiffness matrix at the beginning of the jth iteration step of the i-th increment, and They represent the user-preset reference load and the unbalanced force calculated in the previous iteration step respectively.

[0022] Alternatively, the generalized displacement control method is used to calculate the load increment factor for the current iteration step:

[0023] If the current iteration step is 1:

[0024]

[0025] Where, It is the initial load increment factor preset by the user, and GSP is the generalized stiffness parameter;

[0026] If the current iteration step is greater than 1:

[0027]

[0028] Where, represents the load increment factor, Represents the reference displacement of the 1st iteration of the initial increment (i-1).

[0029] Optionally, the process of calculating the displacement increment of the current iteration step and updating the total displacement and total load of the structure includes:

[0030] Based on the superposition principle, the displacement increment caused by the external load is superimposed on the calibration displacement to obtain the displacement increment of the current iteration step; wherein the displacement increment caused by the external load is the product of the load increment factor and the reference displacement;

[0031] The external load increment is obtained based on the product of the load increment factor and the reference load;

[0032] The total displacement and total load of the structure are updated based on the external load increment and the displacement increment of the current iteration step.

[0033] Optionally, the expressions for updating the total displacement and total load of the structure based on the external load increment and the displacement increment of the current iteration step are:

[0034]

[0035] Where, and are the total loads applied to the structure at the j-1th iteration step and the jth iteration step of the i-th increment, respectively. and are the total displacement responses of the structure from the initial state to the i-th incremental step, the j-1th iteration step, and the j-th iteration step, respectively. represents the external load increment, Indicates the displacement increment of the current iteration step.

[0036] Optionally, the process of updating the structural geometric state and the unit internal force based on the displacement increment of the current iteration step includes:

[0037] updating the structural geometric state based on the updated total displacement and total load of the structure and the state variables to obtain updated node coordinates and an updated unit coordinate basis;

[0038] Based on the updated node coordinates and the updated unit coordinate basis, the elastic stiffness matrix and the geometric stiffness matrix of each unit are calculated using the state variables to obtain the updated unit elastic stiffness matrix and the geometric stiffness matrix;

[0039] The updated unit internal force is obtained based on the updated node coordinates, the updated unit coordinate basis, the updated unit elastic stiffness matrix and the updated geometric stiffness matrix.

[0040] Optionally, the equilibrium condition is:

[0041]

[0042] Where, represents the structural unbalanced force, ε r represents the relative error limit, represents the total load applied to the structure at the jth iteration step.

[0043] Compared with the prior art, the present invention has the following advantages and technical effects:

[0044] The method of the present invention is based on nonlinear finite element theory and establishes an efficient, stable and accurate structural nonlinear and post-buckling analysis method. It can effectively identify the post-buckling path and multiple critical points of the structure, which is of great significance for predicting the bearing capacity or ultimate strength of the structure. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] The accompanying drawings, which constitute part of this application, are intended to provide a further understanding of this application. The exemplary embodiments and descriptions of this application are intended to explain this application and do not constitute an improper limitation on this application. In the accompanying drawings:

[0046] Figure 1 This is a flow chart of a method for quickly and stably solving post-buckling paths with multiple critical points based on a generalized displacement control method according to an embodiment of the present invention;

[0047] Figure 2 The embodiment of the present invention is based on an incremental-iterative mechanism for updating the Lagrangian equation;

[0048] Figure 3The prediction model proposed in the embodiment of the present invention: anchored secant predictor;

[0049] Figure 4 Schematic diagram of a two-dimensional hinged deep arch structure in an embodiment of the present invention;

[0050] Figure 5 The post-buckling response of multiple critical points under symmetrical load conditions calculated based on the method of the present invention in an embodiment of the present invention;

[0051] Figure 6 is the post-buckling response of multiple critical points under eccentric load conditions calculated based on the method of the present invention in an embodiment of the present invention;

[0052] Figure 7 Schematic diagram of a three-dimensional open spherical shell structure in an embodiment of the present invention;

[0053] Figure 8 This is the post-buckling response of multiple critical points calculated based on the method of the present invention in an embodiment of the present invention. DETAILED DESCRIPTION

[0054] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments in this application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0055] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0056] Example 1

[0057] like Figure 1 As shown, this embodiment provides a method for quickly and stably solving the post-buckling path of multiple critical points based on the generalized displacement control method. It is a method for solving structural nonlinearity and post-buckling problems based on the nonlinear finite element method, which is particularly beneficial for and not limited to stability analysis using only elastic stiffness.

[0058] The present invention novelly proposes an "anchored secant predictor" to predict the first iteration step of each incremental step, that is, the secant stiffness of the previous incremental step is used to approximate the equivalent tangent stiffness of the current first iteration step. This can complete the prediction stage of the first iteration step without the need for a stiffness matrix, greatly improving the computational efficiency.

[0059] The present invention combines the "anchored secant predictor" with the "generalized displacement control method" to achieve a fast, stable and accurate solution to the post-buckling path of multiple critical points.

[0060] The present invention allows users to perform structural post-buckling analysis using only elastic stiffness in moderate nonlinear problems, while overcoming the problem of reduced efficiency in the prior art when only elastic stiffness matrix is ​​used for buckling analysis.

[0061] When solving highly nonlinear problems, the present invention strengthens the orthogonality condition of the generalized displacement control method by adding a rigid body qualified geometric stiffness matrix in subsequent iterations to ensure that the iteration can converge stably near the post-buckling path with large curvature.

[0062] In addition to the necessary initialization definitions (including the initial setting of state variables and the initial condition definition of the anchored secant predictor), the solution process can be divided into three stages: prediction stage, calibration stage, and error checking stage. To maintain the generality of the method, the calculation process of the jth iteration step of the i-th increment step is described here:

[0063] Prediction stage:

[0064] Step 1: Inherit the state variables from the previous iteration, including node coordinates, element coordinate basis, node displacements, element internal forces, external loads, initial load increment factors, and element stiffness matrices. If the current iteration j is greater than 1, the coordinate transformation matrix [T] is further calculated based on the state variables, and the element stiffness matrices are assembled into the structural stiffness matrix [K]. The structural stiffness matrix is ​​also modified based on the boundary conditions.

[0065] Step 2: Predict structural reference displacements and calibration displacement If the current iteration step j = 1, the reference displacement With calibration displacement Will be calculated by the proposed "anchored secant predictor":

[0066]

[0067] If the current iteration step j>1, the reference displacement With calibration displacement Then it is calculated by the definition:

[0068]

[0069] In the above formula {ΔU i-1} is the displacement increment of each iteration of the i-1th incremental step The cumulative value of λ i-1 The incremental load factor for each iteration of the i-1th incremental step The cumulative value of . is the structural stiffness matrix at the beginning of the jth iteration step in the i-th incremental step, which is assembled from the element stiffness matrices calculated in the j-1th iteration step. and Respectively represent the user preset reference load and the unbalanced force calculated in the previous iteration step. (If the current iteration step j = 1, the structure is in equilibrium, and there is an unbalanced force Therefore, the calibration displacement calculated from the unbalanced force is ).

[0070] Step 3: Calculate the load increment factor according to the generalized displacement control method using the calculated reference displacement and calibration displacement For the first iteration the load increment factor is:

[0071]

[0072] In the above formula The initial load increment factor is preset by the user, which determines the length of the incremental step. GSP is the generalized stiffness parameter, which is defined as the ratio of the initial displacement increment of the structure to the current displacement increment:

[0073]

[0074] in and are the reference displacements of the first iteration of the initial increment, the first iteration of the i-1th increment, and the first iteration of the i-th increment. For the jth iteration (j>1), the load increment factor is:

[0075]

[0076] and are the reference displacement and calibration displacement of the jth iteration step of the i-th incremental step, respectively.

[0077] Step 4: Using the superposition principle, the displacement increment caused by the external load (expressed as the load increment factor With reference displacement The product of the unbalanced force and the calibration displacement caused by Superposition, to obtain the displacement increment of the jth iteration step of the i-th increment step The corresponding external load increment Expressed as load increment factor With reference load The product of The total displacement and total load of the structure in the current iteration step are updated as follows:

[0078]

[0079] In the above formula and are the total loads applied to the structure at the j-1th iteration step and the jth iteration step of the i-th increment, respectively; and are the total displacements of the structure from the initial state to the i-th incremental step, the j-1th iteration step, and the j-th iteration step, respectively.

[0080] Calibration phase:

[0081] Step 5: Update the structural geometry (including node coordinates and unit coordinate basis) based on the predicted structural displacement increments in step 4. Based on the updated geometric state, calculate the coordinate transformation matrix [T] of each unit in turn, and use the displacement increments obtained in step 4 Obtain the displacement increment of each unit through the coordinate transformation matrix [T]

[0082] Step 6: Traverse each unit and calculate the elastic stiffness matrix [k e ] and the geometric stiffness matrix [k g ]. (For moderate nonlinear problems, the geometric stiffness matrix can be omitted, but when solving highly nonlinear problems with large curvature at the rebound point in the post-buckling path, the geometric stiffness matrix of the rigid body must be included. Here, for the versatility of the program, the geometric stiffness matrix [k g ] steps are listed)

[0083] Step 7: Based on the rigid body criterion and the theory of small deformation increments, traverse all elements and update the internal forces of the elements in the deformed configuration:

[0084]

[0085] In the above formula is the element internal force inherited from step 1 at the j-1th iteration without updating the element geometry, and is referenced to the element coordinate basis of the unupdated geometry; and The internal force of the unit after the structural geometry is updated in step 5, with the unit coordinate basis after the structural geometry is updated as the reference. (Here the symbol The upper left script "1" indicates that the variable belongs to the h-1th iteration step, while the lower left script "1" indicates that the variable is based on the coordinate basis of the j-1th iteration step; and the symbol Similarly, it represents the jth iteration step)

[0086] Error checking stage:

[0087] Step 8: Traverse all units and update the internal forces of each unit Stored so that it can be called when the element internal force is updated in the next iteration step, that is, as the element internal force in the next iteration step Afterwards, the internal force of each unit The total internal force of the structure is obtained by converting the coordinate transformation matrix [T] to the structural coordinate basis and superimposing all unit internal forces.

[0088] Step 9: Calculate the structural unbalanced force based on the difference between the total external load and the total internal force And determine whether the structure has reached equilibrium conditions (Usually the relative error limit ε r It can be set to 0.0001) If the equilibrium condition is not met, return to step 1 and start the j+1th iteration step; if the equilibrium condition is met, the i-th incremental step is considered completed, and return to step 1 and start the 1st iteration step of the i+1th incremental step.

[0089] Furthermore, the theoretical basis of the technical solution includes:

[0090] In continuum mechanics, the motion and deformation of a structure can be described by three configurations: the initial configuration C 0 , the previous step knows the configuration C i-1 and the currently unknown configuration C i According to the updated Lagrangian formula, select the known configuration C in the previous step i-1 As a reference configuration to describe the structure from configuration C i-1 In place shape C i Incremental behavior between . Figure 3 The incremental-iterative mechanism based on the Lagrangian series is demonstrated, such as Figure 2 As shown, it serves as the theoretical basis of the present invention (the upper left subscript of each variable in the figure and the subsequent theory indicates the incremental step when it appears, and the lower right subscript indicates the iterative step.):

[0091] (1) The structure is in a known configuration C i-1 Time (at Figure 3 Point a), under load {P i-1} is in equilibrium, and its corresponding displacement response is {U i-1 Apply a small load increment {ΔP i} Make the load from {P i-1} changes to {P i}.

[0092] (2) Prediction stage: According to the stiffness equation The load increment {ΔP i} displacement increment {ΔU i}(At this time, along the trajectory ab reaches Figure 3 b point), and at the same time, the displacement increment {Δu i}.

[0093] (3) Calibration stage: According to the required displacement increment {ΔU i}Update structure geometry (from Figure 3 From point a to point d), and calculate the total internal force of the unit (by Figure 3 The calculation of the total internal force of the unit includes two parts: the change of the internal force of the unit caused by the rotation of the rigid body (represented by the line segment ce). Figure 3 Line segment oa to line segment cd) and the unit force increment caused by natural deformation (by Figure 3 (represented by line segment de).

[0094] (4) Error checking stage: Calculate and superimpose the internal forces of each unit according to the coordinate transformation to obtain the total internal force of the structure {F i}, and then the structural unbalanced force is obtained by calculating the total external load and the total unit internal force, that is, {R i}={P i}-{F i}(Depend on Figure 3 (represented by line segment eb)

[0095] (5) Next iteration: The unbalanced force is applied to the structure as an external load, and a new iteration step is started to calculate the structural displacement and unit force; if the unbalanced force meets the preset allowable conditions, the structure is considered to be in equilibrium, the calculation of the i-th incremental step is terminated, and the next i+1-th incremental step is started.

[0096] According to the prediction stage in the above incremental-iterative mechanism, the incremental stiffness equation at the i-th incremental step and the j-th iteration step can be expressed in the general form:

[0097]

[0098] in is the structural stiffness matrix at the beginning of the jth iteration step of the i-th increment, is the displacement increment to be determined, is the reference load set, is the load increment factor to be determined, is the unbalanced force of the structure in the j-1th iteration step. Usually, the displacement increment to be calculated can be divided into and Two parts. It is called the reference displacement and is determined by the reference load caused by; It is called the calibration displacement, which is caused by the unbalanced force Caused by:

[0099]

[0100] Equation (9) contains a total of N equations, but contains N+1 unknowns (including the displacement increment of N degrees of freedom With constant ). Therefore, in order to solve equation (9), a constraint equation needs to be added, and its generalized form is:

[0101]

[0102] In the above formula, k, {C} and These are all constraint parameters that need to be assumed in advance. Different constraint parameters can distinguish different path tracking methods. and It is subsequently determined by the initial value conditions. Combining equations (9) and (13), the load increment factor can be obtained as:

[0103]

[0104] At the beginning of the first iteration of each incremental step, the structure is in equilibrium and the displacement is calibrated. Formula (14) is transformed into:

[0105]

[0106] In the generalized displacement control method, it is assumed that In each incremental step, it is kept as a constant value, that is, in the first iteration step, j = 1, In subsequent iterations j≥2, So for the second iteration:

[0107]

[0108] Substitute initial conditions From this solution we can get Substituting it back into formula (15) we can obtain:

[0109]

[0110] GSP is the generalized stiffness parameter, which is defined as the ratio of the initial displacement increment to the current displacement increment of the structure. It is used in the generalized displacement control method to reflect the change in structural stiffness:

[0111]

[0112] The generalized stiffness parameter (GSP) is a control parameter of the generalized displacement control method. Its characteristic is that it is negative only in the first increment after passing the limit point. At this time, the direction of the load increment should be opposite to the previous increment. In simple terms, when GSP>0 (the reference displacement of the i-1th increment is negative), and the reference displacement of step i acute angle), the load increment factor for the i-th increment step The sign of the load increment in step i-1 should be the same as that of the load increment in step i-1. and the reference displacement of step i obtuse angle), the direction of the load increment factor should be opposite to the sign of the load increment factor in step i-1.

[0113] Traditionally, for the prediction phase of each iteration, the reference displacement and calibration displacement They can be obtained from the definitions (10) and (11) respectively. Substituting them into the equations (17) and (16) respectively, we can obtain the load increment factors (Iteration 1) and (j>1 for the jth iteration). However, for the first iteration, if the reference displacement According to the proposed anchored secant predictor calculation, the calculation time can be greatly saved while ensuring the prediction accuracy.

[0114] Furthermore, the anchored secant prediction principle includes:

[0115] The concept of secant stiffness is extended from the previous increment to predict the first iteration of the i-th increment, as Figure 3 As shown. For loads from {P i-1} to {P i}, the line segments ac and ab are the first iteration steps calculated using the secant stiffness and tangent stiffness respectively. Their trajectories will pass through the subsequent iteration points m and f, and will eventually converge to the vicinity of point s. When the incremental step size is relatively small, that is, point n is very close to point a. At this time, the secant stiffness anchored at the i-1th incremental step (the slope of the line segment na or ac) will also be very close to the tangent stiffness (the slope of the line segment ab) of the first iteration step of the current i-th incremental step, that is, It is worth noting that Figure 5 Points c and b in the figure are very close, but for the sake of convenience, they are exaggerated here. Figure 5 The secant stiffness equation for the i-1th incremental step can be written as:

[0116]

[0117] In the above formula is the secant stiffness of the i-1th incremental step; {ΔU i-1} and λ i-1 are the cumulative displacement increment and cumulative load increment factor of the i-1th incremental step respectively. According to the definition, the cumulative load {ΔPi-1} can be expressed as a reference load and cumulative load increment factor λ i-1 The product of .

[0118] According to the reference displacement of the first iteration step The definition of calculation (10) uses the property that the secant stiffness and the tangent stiffness are approximately equal in the sense of incremental step. The calculation of the reference displacement of the first iteration step can be written as:

[0119]

[0120] Combining equations (19) and (20), we can simplify to:

[0121]

[0122] Equation (21) is called the anchored secant predictor. It is worth noting that in order to start the anchored secant predictor, the definition equation (10) should be used to calculate the first iteration of the first increment. The characteristic of the anchored secant predictor is that once the i-1th increment is completed, the concept of secant stiffness can be directly used to predict the first iteration of the i-th increment. Note that Equation (21) does not require any information related to the stiffness equation, and only the data of the previous increment needs to be used for assignment operations. Therefore, the steps of assembling the structural stiffness matrix and solving the linear equations in the first iteration can be omitted, which can greatly save computing time.

[0123] Furthermore, the orthogonal properties of the generalized displacement control method include:

[0124] The purpose of subsequent iterations is to eliminate the error caused by the approximation of the incremental method in the first iteration. Based on the generalized displacement control method, the constraint equation for the j-th iteration j>1 is:

[0125]

[0126] In the above formula is the displacement increment to be determined at the current j-th iteration step, is the reference displacement increment of the first iteration step of the i-1th increment, which can be regarded as the guiding vector of the iterative process. The constraint equation (22) can be regarded as the displacement to be positioned. With guidance vector The orthogonality condition between them, that is, the subsequent iteration steps and etc. always with the guide displacement direction is perpendicular to .

[0127] Furthermore, the unit internal force calibration based on the rigid body criterion includes:

[0128] The calibration of unit internal forces is performed at the unit level and is divided into two aspects: (1) recovering the unit force increment {Δf} from the displacement increment; (2) updating the total unit internal force acting on the unit. After the internal force calibration phase is completed, the total internal force of the element is used to calculate the structural unbalanced force {R}. It is conceivable that if the total internal force of the element is If the update is inaccurate, the calculated unbalanced forces will naturally be ill-conditioned, and subsequent iterations will naturally fail to converge to the correct position. Therefore, the update of the total internal forces of the element during the element internal force calibration phase determines the accuracy of the obtained solution.

[0129] According to physical objectivity, any finite element displacement {Δu} can be divided into two parts: rigid body displacement {Δu} r and natural deformation {Δu} n When the structure is unstable, the rigid body rotation of the unit {Δu} r It accounts for the vast majority of the unit displacement increment {Δu}; in comparison, the natural deformation of the unit is much smaller. This objective phenomenon requires that the change in the unit internal force caused by the rigid body displacement must be accurately considered in the incremental-iterative analysis, because it represents the main part of the unit displacement increment; on the contrary, the unit force increment caused by the natural deformation can be regarded as the unit force increment within the small deformation range and can be calculated only using the elastic stiffness:

[0130] {Δf}=[k e ]{Δu} n (twenty three)

[0131] Since the elastic stiffness matrix needs to meet the requirements of the piecewise test, no unit force increment will be generated under the rigid body displacement. Therefore, the above formula (23) can also be rewritten as:

[0132] {Δf}=[k e ]{Δu} n =[k e ]({Δu}-{Δu} r )=[k e ]{Δu} (24)

[0133] On the other hand, the changes in the internal forces of the unit caused by the displacement of the rigid body can be accurately handled according to the concept of the rigid body criterion. The rigid body criterion is a simple and intuitive physical criterion and is a universally valid principle for testing the quality of finite elements subjected to initial forces. For a finite element subjected to a set of balanced forces before the rigid body rotates, the internal forces acting on the unit during the rigid body rotation will rotate along with the rigid body rotation of the finite element. After the rigid body rotation, the unit is still in equilibrium, and the forces acting on the unit have only changed in direction compared to before the rotation, while the magnitude remains unchanged. The rigid body criterion can be expressed mathematically as:

[0134]

[0135] The change of unit internal force caused by the rigid body rotation can be considered completely accurately through equation (25). The program is also very simple. No additional coordinate changes are required. It is only necessary to change the unit internal force of the previous iteration step. Numerically, it is directly regarded as the internal force of the unit after the rigid body rotates By superimposing the changes in unit internal force caused by rigid body displacement and the increments in unit internal force caused by natural deformation, the total internal force of the unit on the updated geometry can be obtained:

[0136]

[0137] The characteristic of the above formula is that only the elastic stiffness matrix is ​​needed to complete the update of the unit internal force, without being affected by the approximation of the geometric stiffness derivation.

[0138] The method of the present invention is based on nonlinear finite element theory and establishes an efficient, stable and accurate structural nonlinear and post-buckling analysis method. It can effectively identify the post-buckling path and multiple critical points of the structure, which is of great significance for predicting the bearing capacity or ultimate strength of the structure.

[0139] Example 2

[0140] As a specific implementation of this embodiment, further numerical verification of the technical solution of the present invention is performed.

[0141] Example 1 parameters:

[0142] The geometry of the two-dimensional articulated deep arch structure is as follows Figure 4 As shown, under the action of vertical load P, the arch span L = 254 cm, the interface moment of inertia Iz = 41.62 cm 4 , cross-sectional area A=64.52cm 2 , elastic modulus E = 1378kPa. In the finite element analysis, the deep circular arch is first divided into 25 equal two-dimensional beam elements; then, the dome element is divided into two elements to provide a central vertex. Therefore, a total of 26 two-dimensional beam elements are used to simulate this problem. In this embodiment, two loading conditions will be considered: one is symmetrical loading, in which the vertical load P will be applied to the center node of the dome; and the other is eccentric loading, in which the vertical load P will be applied to the adjacent nodes that deviate from the center node. Set the reference load in the incremental analysis Initial load increment factor It should be added that the numerical simulations of all the embodiments in the present invention were performed on the 13th Gen Core TMi5-13600KF, equipped with a 3.50GHz CPU and 32.0GB memory, running under Windows 11, 64-bit operating system.

[0143] To verify the improvement in computational efficiency of the anchored secant predictor proposed in this invention, we conducted numerical simulations and compared the computational efficiency with that of the traditional predictor, which is defined as Equation (10). Figure 5 (a) is the load-deflection curve of two predicted sub-deep circular arches under symmetrical load conditions. Figure 5 (b) is the GSP-deflection curve of the two predictors under symmetrical load conditions. Figure 6 (a) is the load-deflection curve of two predicted sub-deep circular arches under eccentric load conditions. Figure 6 (b) is the GSP-deflection curve of the two predictors under symmetrical load conditions.

[0144] from Figure 5 (a) and Figure 6 (a) It can be seen that both predictors can track the post-buckling path of the structure and identify multiple critical points well; Figure 5 (b) and Figure 6 (b) also shows that the two predictors are basically the same in predicting the generalized stiffness parameter GSP, and both can identify the zero points on the GSP curve (which correspond to the extreme points in the post-buckling path). The computational efficiency comparison results of the anchored secant predictor and the traditional predictor under the two loading conditions are shown in Tables 1 and 2:

[0145] Table 1

[0146]

[0147]

[0148] Table 2

[0149]

[0150] The calculation results show that under both central and eccentric load conditions, the incremental steps and average iterations used by the two predictors are comparable. However, it is clear that the anchored secant predictor takes less time to compute than the traditional predictor, demonstrating the proposed anchored secant predictor's superiority in improving computational efficiency.

[0151] Example 2 parameters:

[0152] The geometric shape of the three-dimensional open spherical shell structure is as follows Figure 7As shown, the shell edges are hinged and the center is subject to a vertical load P. The shell side length 2a = 1569.8mm, thickness t = 99.45mm, radius of curvature R = 2540mm, Poisson's ratio v = 0.3, and elastic modulus E = 6895kPa. Due to the axisymmetric nature of the shell, only one-quarter of the shell is modeled as an 8×8 grid for analysis during finite element modeling. Setting the reference load in incremental analysis Initial load increment factor

[0153] To verify the improvement of the computational efficiency of the method of the present invention, we conducted numerical simulations and compared the computational efficiency with two traditional predictors. Both of these traditional predictors are completed by the definition formula (10), where the stiffness matrix of the traditional predictor 1 is It also contains the elastic stiffness matrix [K e ] and the geometric stiffness matrix [K g ]; while the traditional predictor 2 only contains the elastic stiffness matrix [K e ]. Figure 8 (a) is the load-deflection curve of three predicted sub-spherical shells under central load conditions. Figure 8 (b) is the GSP-deflection curve of three predicted sub-spherical shells under central load conditions.

[0154] from Figure 8 (a) It can be seen that the three predictors can track the post-buckling path of the structure and identify multiple critical points well; however, it should be noted that Figure 8 The difference between the GSP curve of the traditional predictor 2 in (b) and the other two cannot be ignored. This is because a large error occurs when the predictor using only elastic stiffness predicts the displacement of the first iteration step, resulting in an error in the GSP calculation. This also explains why the problem of reduced computational efficiency occurs when only elastic stiffness is used for prediction, because the error in the first iteration step requires more subsequent iteration steps to eliminate, resulting in a waste of computation. However, it can be seen that the anchored secant predictor proposed in the present invention can cleverly solve the problem of reduced efficiency when only elastic stiffness is used. Further, the computational efficiency results of the three predictors are shown in Table 3:

[0155] Table 3

[0156]

[0157] The efficiency improvement of the anchored secant predictor in Example 2 is even more significant, with the computational efficiency nearly double that of the traditional predictor. Furthermore, the anchored secant predictor can overcome the low computational efficiency problem reported in the literature when using only elastic stiffness.

[0158] The above are merely preferred embodiments of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in this application should be included in the scope of protection of the present application. Therefore, the scope of protection of the present application should be based on the scope of protection of the claims.

Claims

1. A method for quickly and stably solving post-buckling paths with multiple critical points based on the generalized displacement control method, characterized in that: The following steps are involved: S1, initialize state variables and anchored secant predictor; S2. Predicting a reference displacement, a calibration displacement, and a load increment factor based on the state variable and the anchored secant predictor; S3. Calculating the displacement increment of the current iteration step based on the reference displacement, calibration displacement and load increment factor, and updating the total displacement and total load of the structure; S4. updating the structural geometric state and unit internal force based on the displacement increment of the current iteration step; S5. Calculate the difference between the total internal force of the structure and the external load based on the updated structural geometric state and the unit internal force to obtain a structural unbalanced force, and perform an error check based on the structural unbalanced force; S6, looping through S1-S5 until the structural unbalanced force satisfies the equilibrium condition, and starting the next incremental step; Incremental analysis is performed until the structural load or displacement meets the user's requirements, and multiple critical point post-buckling paths are output.

2. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 1, characterized in that: The state variables include: node coordinates, unit coordinate basis, node displacement, unit internal force, external load, initial load increment factor and unit stiffness matrix.

3. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 1, characterized in that: The process of predicting the reference displacement, calibration displacement, and load increment factor involves: Predicting a displacement increment of a current incremental step based on the state variable and the anchored secant predictor, wherein the initial displacement increment includes: a reference displacement and a calibration displacement; The load increment factor of the current iteration step is calculated using a generalized displacement control method based on the reference displacement and the calibration displacement.

4. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 3, characterized in that: The process of predicting the displacement increment for the current increment step includes: If the current iteration step is 1, the reference displacement and calibration displacement are calculated based on the anchored secant predictor; wherein the expression for calculating the reference displacement and calibration displacement based on the anchored secant predictor is: If the current iteration step is greater than 1, the reference displacement and calibration displacement are calculated based on the definition; wherein the expression for calculating the reference displacement and calibration displacement based on the definition is: Where, Indicates the reference displacement when the current iteration step is 1, Indicates the calibration displacement when the current iteration step is 1, {ΔU i-1 } represents the displacement increment of each iteration of the i-1th incremental step The cumulative value of λ i-1 Indicates the incremental load factor for each iteration of the i-1th incremental step The cumulative value of Indicates the reference displacement when the current iteration step is greater than 1, Indicates the calibration displacement when the current iteration step is greater than 1, represents the structural stiffness matrix at the beginning of the jth iteration step of the i-th increment, and They represent the user-preset reference load and the unbalanced force calculated in the previous iteration step respectively.

5. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 4, characterized in that: The expression for calculating the load increment factor of the current iteration step using the generalized displacement control method is: If the current iteration step is 1: Where, It is the initial load increment factor preset by the user, and GSP is the generalized stiffness parameter; If the current iteration step is greater than 1: Where, represents the load increment factor, Represents the reference displacement of the 1st iteration of the initial increment (i-1).

6. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 5, characterized in that: The process of calculating the displacement increment of the current iteration step and updating the total displacement and total load of the structure includes: Based on the superposition principle, the displacement increment caused by the external load is superimposed on the calibration displacement to obtain the displacement increment of the current iteration step; wherein the displacement increment caused by the external load is the product of the load increment factor and the reference displacement; The external load increment is obtained based on the product of the load increment factor and the reference load; The total displacement and total load of the structure are updated based on the external load increment and the displacement increment of the current iteration step.

7. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 6, characterized in that: The expressions for updating the total displacement and total load of the structure based on the external load increment and the displacement increment of the current iteration step are: Where, and are the total loads applied to the structure at the j-1th iteration step and the jth iteration step of the i-th increment, respectively. and are the total displacement responses of the structure from the initial state to the i-th incremental step, the j-1th iteration step, and the j-th iteration step, respectively. represents the external load increment, Indicates the displacement increment of the current iteration step.

8. The method for rapidly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 1, characterized in that: The process of updating the structural geometric state and element internal forces based on the displacement increment of the current iteration step includes: updating the structural geometric state based on the updated total displacement and total load of the structure and the state variables to obtain updated node coordinates and an updated unit coordinate basis; Based on the updated node coordinates and the updated unit coordinate basis, the elastic stiffness matrix and the geometric stiffness matrix of each unit are calculated using the state variables to obtain the updated unit elastic stiffness matrix and the geometric stiffness matrix; The updated unit internal force is obtained based on the updated node coordinates, the updated unit coordinate basis, the updated unit elastic stiffness matrix and the updated geometric stiffness matrix.

9. The method for quickly and stably solving post-buckling paths at multiple critical points based on the generalized displacement control method according to claim 1, characterized in that: The equilibrium condition is: Where, represents the structural unbalanced force, ε r represents the relative error limit, represents the total load applied to the structure at the jth iteration step.

Images