Finite element approximate calculation method for thin-wall arch structure and lightweight system

Abstract

The invention relates to the field of dynamic jump behavior simulation analysis of a thin-wall arch structure, in particular to a thin-wall arch structure finite element approximate calculation method and a lightweight system. The invention particularly provides a thin-wall arch structure finite element approximate calculation method and a lightweight system, and the method comprises the steps: carrying out one-time quasi-static equilibrium path analysis, obtaining a characteristic value closest to zero of each increment step, recognizing characteristic values before and after sign change, carrying out two-time perturbation finite element analysis, and carrying out two-time perturbation finite element analysis; the method comprises the following steps of: obtaining an approximate transient dynamic jump parameter and a feature vector, obtaining an approximate jump parameter and a feature vector at a saddle point based on a feature value interpolation closest to zero, finally obtaining an approximate transient jump response through a Gaussian hyper-geometric function, and proposing a Python / C-based hybrid programming method to realize the lightweight of a finite element analysis system. And the calculation cost and the analysis difficulty of the flexible shallow arch saddle point transient jump analysis can be effectively reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of microelectromechanical intelligent applications of flexible structures and finite element simulations, and particularly relates to a finite element approximate calculation method and a lightweight system for a thin-walled arched structure. Background Art

[0002] A thin-walled arched structure is an arched structure with an arch height much smaller than its length. Typical geometric shapes include circular arches, parabolic arches, etc. A flexible thin-walled arched structure refers to a thin-walled arched structure with a small material stiffness and is prone to elastic deformation. Such structures are prone to deformation, and this deformation does not represent structural failure but is what is required artificially and can be used to design various microelectromechanical control components. There is a saddle-point instability phenomenon in thin-walled arches. Figure 1 An example of the saddle-point instability phenomenon is given, that is, the load-deformation equilibrium path diagram of the thin-walled arched structure under an external load force and the corresponding configuration diagram of the thin-walled arch. In this example, the original thin-walled arched configuration protrudes upward. At the initial stage of loading, as the load slowly increases, a downward displacement is generated at the center of the thin-walled arch, and the system is still in a balanced state. The deformation amount gradually increases. When the load reaches the saddle point (i.e., the highest point of the load), the system is in a critical equilibrium state (corresponding to the saddle-point configuration). When the load further increases, the thin-walled arch cannot maintain the equilibrium state and will undergo a dynamic shape jump. The configuration after the jump is quite different from the saddle-point configuration. The jump of the flexible thin-walled arch makes it easy to fabricate the thin-walled arch into a switch-type component, which is widely used in the microelectromechanical field. For example, for a thin-walled arched structure made of a metal sheet, after passing a direct current, under the action of a uniform magnetic field, the Ampere force can be used as the driving load, and the jump of the thin-walled arch can be controlled by changing the current intensity or the magnetic field intensity. If the thin-walled arch is locally charged, the thin-walled arch can be placed between capacitors, and the uniform electric field generated by the capacitors drives the thin-walled arch to jump under the action of the Coulomb force. The transient jump dynamic characteristics (such as the time required for the jump and the parameters related to the relationship between the jump deformation and time) are important design parameters. In the analysis of the transient jump process, if conventional general-purpose finite element software such as Abaqus and Ansys is used, to obtain the response of the dynamic jump, explicit or implicit dynamic finite element analysis needs to be carried out, and calculations need to be carried out for multiple time increment steps in sequence. The operation is complex, the calculation cost is high, the analysis efficiency is low, and it poses challenges to the parameter optimization of related structures.

[0003] Aiming at the deficiencies of the above technologies, the present invention proposes a perturbation finite element approximate calculation method for the transient response of the saddle-point type jump of a flexible thin-walled arch and a related lightweight finite element system development method. Without conventional dynamic analysis, an approximate transient dynamic jump response can be obtained only by quasi-static equilibrium path analysis and two perturbation finite element analyses, shortening the analysis time and accelerating the design optimization in related fields. Summary of the Invention

[0004] The object of the present invention is to overcome the deficiencies of the above-mentioned prior art, and to provide a finite element approximate calculation method and a lightweight system for a thin-walled arched structure, aiming to solve the technical problem of too long analysis time in the prior art.

[0005] To achieve the above object, the present invention proposes a finite element approximate calculation method for a thin-walled arched structure, which at least includes the following steps: S1. Generation of a geometric model and finite element mesh of the thin-walled arch; S2. Tracking of the equilibrium path of the thin-walled arch and extraction of eigenvalues based on the Riks method; S3. Finite element calculation of saddle point perturbation of the thin-walled arch and parameter interpolation based on two eigenvalues; S4. Approximate calculation of the transient response based on the Gaussian hypergeometric function.

[0006] Preferably, in step S1, according to the set parameters of the thin-walled arch, a three-dimensional geometric model of the thin-walled arch is established, and a solid finite element mesh of the eight-node hexahedron type is divided, and the node and element connection information is output.

[0007] Preferably, S2 includes the following steps: a. Define the interpolation function of the global field variable, and the specific method is as follows: (1) Use eight-node hexahedron elements and adopt the following Lagrange interpolation method: , , , , , where is the Lagrangian coordinate of the standard element, is the coordinate at the field variable, is the numerical value of the field variable at the eight nodes, is the interpolation function of each node inside the element; (2) Let the three-dimensional coordinate vector before deformation be , and define of the interpolation form as ; From the one-to-one correspondence relationship between and defined by Equation [2], regard as a function of ; ; (3) Uniformly number all the nodes , where represents the number of nodes of the finite element model of the thin-walled arch structure. The global field variables in the final thin-walled arch structure are determined by the following interpolation formula: , where represents the value of the global field variable at , is the global field interpolation function, is the value of the field variable at the -th node; b. Calculate the equilibrium path under the external load. The specific calculation method is as follows: Adopt the Green-Lagrange strain , where is the deformation gradient tensor, represents the transpose of the tensor, represents the unit tensor; Based on the principle of virtual work, establish the variational equation: , where represents the degrees of freedom in the X, Y, and Z directions, represents the node number, is the region occupied by the thin-walled arch structure before deformation, represents the variational operation, is the fourth-order tensor of isotropic linear elastic stiffness. ":" represents the contraction operation of tensors, represents the Lagrange multiplier that constrains the -th degree of freedom of the -th node, represents the -th degree of freedom of the -th node, is the value of the degrees of freedom defined by the boundary conditions, represents the load force value applied in the direction of the -th degree of freedom of the -th node, is the load factor that characterizes the change in the load force value, represents the set of node numbers and degree-of-freedom pairs that are constrained in the boundary conditions, represents the set of node numbers and degree-of-freedom pairs where the load force is applied; Based on the variational principle, establish the discrete equilibrium equation: , where represents the tensor product operation. When , When , When , , take ; The constraint conditions are defined; the load factor and the displacement degree of freedom are used as unknowns at the same time, and the arc-length method is used to solve the discrete equilibrium equations to obtain the equilibrium path under the external load ; c. Extract the eigenvalue closest to zero for each arc-length method increment step. The specific method is as follows: Based on and linearization, a tensor equation in the following form is obtained (for all repeated indices, the Einstein summation convention is used, is the Kronecker symbol ( ), represents the differential): When , , ; When , , After flattening and reshaping the tensor, it is transformed into matrix form: where represents the 1D column vector obtained by flattening the constrained , with the dimension set to , represents the 1D column vector obtained by flattening the remaining , with the dimension set to , represents the 1D column vector obtained by flattening , with the dimension equal to , is the identity matrix, is a square matrix of order , is a matrix of order , is a matrix of order , is a matrix of order , is a zero matrix block, is a 1D column vector related to the loading force; Set the large gain coefficient , for each increment step, solve the eigenvalues of the modified matrix , eliminate the interference of the minimum non-zero eigenvalue caused by the Lagrange multiplier term, and extract the eigenvalue closest to zero at this increment step.

[0008] Preferably, S3 includes the following steps: a. Calculate the approximate transient response of the saddle point jump of the thin-walled arched structure, and the specific method is as follows: Based on D'Alembert's principle, establish the variational equation: , where is the acceleration vector field, is the equilibrium load of the saddle point, is the additional load coefficient, represents the additional load vector, that is, the additional load is and product, is the material density, represents the determinant of the deformation gradient tensor; Assume that the displacement vector field satisfies the perturbation expansion condition: , Assume that the Lagrange multiplier satisfies the perturbation expansion: , where is the flattened vector of all , and denote , is the flattened vector of the Lagrange multiplier at the saddle point, is the coefficient vector of each term of the perturbation expansion, is the physical real time and the perturbation parameter scaled time parameter obtained by the product, is the transient displacement vector field, is the displacement vector field at the saddle point, is the coefficient of the perturbation expansion; Based on order expansion variational equation, obtain the corresponding equations of each order: (1) order equation: and ; (2) order equation: , where , is the variable component, representing the unit vectors in the X, Y, and Z directions; Fs is the deformation gradient tensor at the saddle point; Convert the equation of order to the eigenvalue form, and for all repeated indices, sum according to the Einstein summation convention: , when , ; Based on the equation of order, after node numbering, tensor reshaping, and flattening, it is transformed into an eigenvalue problem of a matrix to obtain the eigenvector field corresponding to the zero eigenvalue , that is , where, is the eigenvector of the eigenvalue problem, is the eigenvector field; after normalizing the eigenvector , denote the normalized eigenvector field as ; (3) The equation of order is: , , where, is the determinant of the deformation gradient tensor at the saddle point; Decompose into ; Based on the solvability condition of the equation of order , obtain the ordinary differential equation about : , where, , ; b. Based on the interpolation of the eigenvalue closest to zero approximate calculation Approximately obtain at the saddle point based on the interpolation of the positive and negative eigenvalues closest to zero , and the specific implementation method is as follows: (1) After calculating the quasi-static equilibrium path, extract the eigenvalue closest to zero for each increment step, arrange them in sequence, and denote the eigenvalues before and after the sign change in the sequence as , ; (2) Respectively take and corresponding increment steps. Assume that the results of these two increment steps are both in the saddle point state, and carry out two perturbation finite element calculations to respectively obtain the corresponding , where and are obtained by using the same normalization method; (3) Use the following linear interpolation formula: , , ; Thus, obtain the value of the saddle point approximation.

[0009] As an optimization, the specific calculation method in S4 is as follows: Assume the initial condition , establish the relationship between and : , where the function is the Gaussian hypergeometric function; Take the limit , obtain the approximate jump time of the perturbation finite element for the thin-walled arch structure; Further, from the approximate expression , establish the relationship between the transient displacement field and time , and finally establish it from the relationship between the transient displacement field and the real physical time , to realize the calculation of the approximate transient response of the saddle point jump of the thin-walled arch structure.

[0010] As an optimization, obtain the eigenvector field corresponding to the zero eigenvalue. In the calculation, take the eigenvector corresponding to the eigenvalue closest to zero as the approximation of the eigenvector field .

[0011] As an optimization, for the normalization of the eigenvector , use the maximum value normalization, that is, take the maximum value of each component , and divide the eigenvector by this maximum value to obtain the normalized eigenvector.

[0012] To achieve the above objectives, the present invention proposes a lightweight system, which uses the above-mentioned finite element approximate calculation method, is based on Python / C hybrid programming, and includes the following modules: input file parsing module, arc length method analysis module, saddle point perturbation finite element analysis module, implicit dynamic analysis module and result visualization module.

[0013] Preferably, the input file parsing module is responsible for inputting the nodes and degrees of freedom directions of load loading, the nodes and degrees of freedom directions of geometric constraints, and the finite element node coordinate table and the node list contained in the unit after the input file is imported into the abaqus platform for modeling and meshing, and is implemented in Python language; the arc length method analysis module is responsible for the Riks analysis step proposed in the present invention, and through the C / Python API method, a C language module with high execution efficiency is compiled and encapsulated into a Python module for Python calling.

[0014] As a preferred method, the saddle point perturbation finite element analysis module first uses C language to calculate the calculation part of each unit involved in the equation, uses Python's own eigsh function to obtain the eigenvector V, and finally uses Python's own library function to calculate the Gaussian hypergeometric function; the implicit dynamic analysis module directly calculates the transient dynamic process of the thin-walled arch structure jump, which is mainly used to verify the correctness of the results of the saddle point finite element module; the result visualization module uses the mature Abaqus python development platform to realize the visualization of the results of each analysis module.

[0015] Compared with the prior art, the finite element approximate calculation method for a thin-walled arch structure and a lightweight system provided by the present invention have the following beneficial effects:

[0016] First, a one-time quasi-static equilibrium path analysis is carried out, and a gain coefficient is introduced to eliminate the interference of extremely small non-zero eigenvalues ​​caused by Lagrange multiplier terms, to obtain the eigenvalue closest to zero for each incremental step, to identify the eigenvalues ​​before and after the change of positive and negative signs, and to carry out two perturbation finite element analyses to obtain approximate transient dynamic jump parameters and eigenvectors. Based on the interpolation of the eigenvalues ​​closest to zero, the approximate jump parameters and eigenvectors at the saddle point are obtained. Finally, the approximate transient jump response is obtained by the Gaussian hypergeometric function, which can effectively reduce the computational cost and difficulty of transient jump analysis of flexible thin-walled arch saddle points. A hybrid programming method based on Python / C is proposed, which can effectively improve the computational efficiency of computationally intensive parts such as the unit stiffness matrix. Matrix calculations are realized through the Python scientific computing library, and the system results are visualized using the Abaqus platform, which provides a certain reference value for the lightweight of related finite element analysis systems.

[0017] When the load is greater than the saddle-point equilibrium load, the flexible shallow arch undergoes dynamic jumps. Direct dynamic jump analysis requires setting multiple time increment steps and performing iterative calculations for each time increment step, resulting in a large computational workload. In contrast, the calculation of the equilibrium path of the shallow arch only requires quasi-static analysis, with lower computational costs and difficulty.

[0018] The features and advantages of the present invention will be described in detail through embodiments in conjunction with the accompanying drawings. Description of the Drawings

[0019] Figure 1 . Schematic diagrams of the equilibrium path, saddle point, and dynamic jump of the shallow arch.

[0020] Figure 2 . Overall step framework diagram of the present invention.

[0021] Figure 3 . Framework diagram of the finite element analysis system for saddle-point perturbation of the shallow arch proposed by the present invention.

[0022] Figure 4 . Schematic diagram of an eight-node hexahedron standard element.

[0023] Figure 5 . Schematic diagram and finite element mesh diagram of a circular shallow arch with hinged ends under a central concentrated force.

[0024] Figure 6 . Geometric parameter diagram of the shallow arch.

[0025] Figure 7 . Node diagram of constraints in symmetric boundary conditions.

[0026] Figure 8 . Node diagram of constraints in hinged boundary conditions.

[0027] Figure 9 . Loading nodes of the concentrated force (node numbers 805, 604, 403, 202, 1) and the central node number 2413 of the central end face.

[0028] Figure 10 . Load factor in the equilibrium path Variation relationship with the total arc length.

[0029] Figure 11 . Eigenvalue distribution when K = 1 in the case of no gain.

[0030] Figure 12 . Gain coefficient K = 10 4 Eigenvalue distribution at this time.

[0031] Figure 13 . Comparison of calculation results between perturbation finite element and implicit finite element.

[0032] Figure 14 . The deformation process of the shallow arch obtained by implicit dynamics calculation.

[0033] Figure 15 . Comparison of the displacement-time curves between the perturbation finite element method of the present invention and the direct dynamics implicit finite element calculation.

[0034] Figure 16 . c0 values for different increment steps.

[0035] Figure 17 . c1 values corresponding to different increment steps.

[0036] Figure 18 . Relationship diagram of the load factor and the total arc length after increasing the arc length increment of the increment step (the arc length increment for each step is 10).

[0037] Figure 19 . Eigenvalues and parameter c0 at increment steps 7, 8, 9, and 10.

[0038] Figure 20 . Eigenvalues and parameter c1 at increment steps 7, 8, 9, and 10.

[0039] Figure 21 . Equilibrium path curve of the shallow arch with non-uniform wall thickness (relationship between the load factor and the total arc length).

[0040] Figure 22 . Corresponding relationship between the increment step, the eigenvalue, and parameter c0.

[0041] Figure 23 . Corresponding relationship between the increment step, the eigenvalue, and parameter c1.

[0042] Figure 24 . Initial shape of the shallow arch, and the shapes of the shallow arch at increment step 8 and increment step 9 after superimposing the scaled eigenvectors. Detailed implementation manners

[0043] To make the objectives, technical solutions, and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. However, it should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the scope of the present invention. In addition, in the following description, the descriptions of well-known structures and technologies are omitted to avoid unnecessarily confusing the concepts of the present invention.

[0044] An embodiment of the present invention provides a lightweight system, a shallow arch saddle point perturbation finite element analysis system based on Python / C hybrid programming. Among them, a finite element approximate calculation method for a flexible thin-walled arched structure is adopted, and the thin-walled arched structure can also be called a shallow arch structure.

[0045] Refer toFigure 2 , the finite element approximation calculation method for the flexible thin-walled arch structure at least includes the following steps: 【1】Generation of the geometric model and finite element mesh of the flexible shallow arch; 【2】Tracking of the equilibrium path of the flexible shallow arch and extraction of the eigenvalue closest to zero based on the Riks method; 【3】Finite element calculation of the saddle point perturbation of the flexible shallow arch and parameter interpolation based on two eigenvalues; 【4】Approximate calculation of the transient response based on the Gaussian hypergeometric function.

[0046] Step [1]

[0047] According to the geometric parameters of the flexible shallow arch (such as wall thickness, width, shape, etc.), establish a geometric model of the shallow arch in the Abaqus platform, and divide the finite element mesh of the eight-node hexahedron type (such as C3D8, C3D8R, C3D8I), output the input file of the Abaqus platform containing node coordinates and element connection information. Finally, import the element node coordinates, the node list contained in the element, the loading force information, and the constraint node information into the perturbation finite element analysis system of the present invention.

[0048] Step [2] specifically includes step 201, step 202, and step 203.

[0049] Step 201. Interpolation of field variables based on eight-node hexahedron elements

[0050] Adopt eight-node hexahedron elements, such as Figure 4 The schematic diagram of the standard element of a regular hexahedron with a side length of 2 is given, where g, h, and r respectively represent the coordinate values in three orthogonal directions with the center as the origin. The number of element integration points is 8, and standard Gaussian numerical integration is used.

[0051] Adopt the following Lagrangian interpolation formula 【1】 , , , ,

[0052] Among them, is the Lagrangian coordinate of the standard element, is the coordinate at the field variable, is the numerical value of the field variable at 8 nodes, is the interpolation function of each node inside the element.

[0053] For non-standard cells, let the three-dimensional coordinate vector before deformation be , in the same interpolation form as the field variable:

[0054] Within each cell, Equation 【2】defines a and one-to-one correspondence. Thus, the field variable and interpolation function within the cell can both be regarded as functions, that is: .

[0055] By uniformly numbering all the nodes , where represents the number of nodes in the model, the global field variable within the shallow arch can be obtained from the following interpolation formula , 【3】 For each cell, each node has internal numbers 1, 2, 3, 4, 5, 6, 7, 8 within this cell.

[0056] However, the finite element model consists of multiple cells, and all the nodes have global numbers .

[0057] For each node , find all the cells that contain this node (denote the cell numbers as ).

[0058] Within each cell internally, the internal number corresponding to the node in the cell is a number among 1, 2, 3, 4, 5, 6, 7, 8. Therefore, this node's internal number in cell is a definite function of and , which may be denoted as . Obviously, .

[0059] Then within cell internally, take (this function is determined by the aforementioned internal interpolation function of the cell). Note that the domain of this function is within cell . It can extend the domain to the entire finite element model by zero-value extension, that is, when is within cell , define , and when is not within cell internally, ; In this way is a function defined on the entire finite element model; Define the node global interpolation function , that is, traverse all the elements that contain the node (denoted by the number ) and accumulate the corresponding global function . Obviously is a function defined on the entire finite element model and related to the global node number . For example, assume the global node , and the numbers of the elements that contain this node are (that is, element 2 and element 3 include this node), then define ;

[0060] Among them, represents the value of the global field variable at , is the global field interpolation function (when restricted to each element, it is the interpolation function), is the value of the field variable at the th node.

[0061] Step 202. Calculation of the equilibrium path under external load

[0062] Each node (number) has three translational degrees of freedom , respectively representing the degrees of freedom in the X, Y, and Z directions.

[0063] Adopt Green-Lagrange strain , among which, is the deformation gradient tensor, represents the transpose of the tensor, represents the unit tensor. Obtained from the principle of virtual work: , 【4】

[0064] Among them, is the area occupied by the shallow arch before deformation, represents the variational operation, is the isotropic linear elastic stiffness tensor (uniquely determined by the elastic modulus and Poisson's ratio), ":" represents the contraction operation of the tensors, represents the Lagrange multiplier that restricts the th degree of freedom of the th node, represents the th degree of freedom of the Number of displacement degrees of freedom, The numerical value of the displacement degrees of freedom defined by the boundary conditions, Indicates the load applied to the th node in the direction of the degree of freedom (unchanged in the equilibrium path), Indicates the load factor of the concentrated force (changes in the equilibrium path, characterizing the change in the magnitude of the concentrated force), Indicates the set of node numbers and degree-of-freedom pairs corresponding to the constraints in the boundary conditions, Indicates the set of node numbers and degree-of-freedom pairs where the load is applied. The calculation of the equilibrium path requires determining the deformations of the shallow arch under static equilibrium for different .

[0065] From Equation [4] and the variational principle, we obtain: ,

[0066] where, is an arbitrary three-dimensional vector field determined by finite element discretization, is at the node in the direction of the degree of freedom.

[0067] From the global interpolation function, , where, represents the unit vectors in the X, Y, and Z directions.

[0068] From the variational principle, we obtain the equation: , [5]

[0069] where, represents the tensor product operation. When , ; when , . Similarly, when , , when, take ;

[0070] In addition, due to the arbitrariness of , we obtain the equation corresponding to the boundary conditions, that is, when , . [6]

[0071] Taking the load factor and the displacement field as unknowns, using the arc-length method (Riks) to obtain the solutions of Equation [5] and Equation [6], and obtaining the external load The equilibrium path below

[0072] Step 203. Extraction of the eigenvalue closest to zero

[0073] Based on the saddle point being a singular point, the overall matrix (i.e., the matrix composed of 9 block matrices in Equation

[10] ) is singular, that is, there is an eigenvalue of zero at the saddle point.

[0074] From Equations 【5】 and 【6】, based on and the linearization of obtains a tensor equation in the following form (for all repeated indices, the Einstein summation convention is adopted; is the Kronecker symbol ( ), that is, it takes 1 when and takes 0 when );

[0075] When , 【7】 ; 【8】

[0076] When : 【9】

[0077] Write the above equation in the form of the following matrix:

[10]

[0078] Where, represents the flattened 1D column vector of the constrained (assuming the dimension is ), represents the flattened 1D column vector of the remaining (assuming the dimension is , represents the flattened 1D column vector of (the dimension is equal to ), is the identity matrix, is a square matrix ( order), is a order matrix, is a order matrix, is a order matrix, is a zero matrix block, is the flattened one-dimensional column vector related to the loading force.

[0079] The finite element model consists of There are nodes, and each node has three degrees of freedom. For the sake of convenience of expression, therefore, the degrees of freedom involved in the above equations , that is, there are degrees of freedom, which are unknowns; similarly is also an unknown, and the number is determined by the number of constraints in the constraint set B; all the equations involved are non-homogeneous linear equations about (theoretically, non-homogeneous linear equations can all be written in matrix form ); For the need of solving, the equations involved need to be transformed into matrix form, that is, in the form of equations (where is a matrix, is a column vector, is a column vector), so that the standard matrix eigenvalue algorithm or matrix equation solving algorithm (such as Gaussian elimination method) can be used to solve; therefore, all , need to be written in the form of a column vector (flattened); At the same time, the coefficients related to , in the equation can be written in the form of a matrix ; the non-homogeneous terms in the equation that are not related to , can also be flattened into a column vector form. The following explains the relevant steps.

[0080] The specific steps are as follows: Number the elements of set B as follows: where #B represents the number of elements in the set; (2) Define the set of all degrees of freedom as , take , and “-” is the set difference operation; Number the elements in as follows, , where #D represents the number of elements in set D.

[0081] (3) Define , represents transpose, is a column vector of dimension #D; Similarly, , is a column vector of dimension #B; similarly, write the Lagrange multipliers as a column vector , so, is a column vector of #B dimensions; (4) That is, The column vector formed by arranging all elements in sequence, that is, .

[0082] (5) Take , , (including #B zero elements), , that is, is The column vector obtained by multiplying the formed column vector by .

[0083] (6) The elements of correspond to the equations [7] (see the specific implementation part) in sequence, The elements of correspond to the equations [9] in sequence, The elements of correspond to the equations [8] in sequence; After arranging the above equations [7], [8], [9] in the corresponding order, extract the corresponding coefficients of the elements about in each equation (in the order of elements), to form the matrix , Therefore, each equation corresponds to a row of a matrix , and each variable corresponds to column of.

[0084] (7) The obtained by the above method has the form of the matrix in equation

[10] , and equation

[10] is the so-called equation .

[0085] When directly calculating the eigenvalues of the matrix in equation

[10] , due to the existence of the Lagrange multiplier term, there may be multiple eigenvalues that are extremely small and close to zero. These extremely small and non-zero eigenvalues make it difficult to track the eigenvalue with the smallest absolute value (refer to Figure 11 ), Therefore, the present invention proposes a method for tracking the eigenvalue closest to zero, specifically as follows:

[0086] Set a large gain coefficient , and solve the eigenvalues of the following matrix for each increment step, which can effectively eliminate the interference of the extremely small non-zero eigenvalues caused by the Lagrange multiplier term (refer to Figure 12 ), and achieve the tracking of the eigenvalue closest to zero.

[0087]

[0088] Step [3] includes Step 301 and Step 302

[0089] Step 301. Calculation of the approximate transient response of the shallow arch saddle point jump

[0090] When the structural load exceeds the equilibrium load of the saddle point, the shallow arch undergoes dynamic jumping. Based on D'Alembert's principle, the following variational equation is established ,

[11]

[0091] where is the acceleration vector field, is the equilibrium load of the saddle point, is the additional load coefficient, represents the additional load vector, that is, the additional load is and the product of, is the material density, is the determinant of the deformation gradient tensor.

[0092] Assume that the displacement vector field satisfies the perturbation expansion condition:

[12]

[0093] Assume that the Lagrange multiplier satisfies the perturbation expansion:

[13]

[0094] where is the flattened vector of all and denote , is the flattened vector of the Lagrange multiplier at the saddle point, is the coefficient vector of the perturbation expansion, is the physical real time and the perturbation parameter the scaled time parameter obtained by the product of, is the transient displacement vector field, is the displacement vector field at the saddle point, is the coefficient of the perturbation expansion.

[0095] Substitute the perturbation expansions (Equation

[12] and Equation

[13] ) into the variational equation

[11] , based on order, to obtain:

[0096] (1) For order, 【14-1】

[0097] and ; 【14-2】

[0098] (2) For order, [15 - 1]

[0099] [15 - 2]

[0100] Among them, in equation [15 - 1] , is an arbitrary value, represents the deformation gradient tensor field at the saddle point.

[0101] Using arbitrariness, write equation [15 - 1] in the form of an eigenvalue problem (for all repeated indices, sum according to the Einstein summation convention): [15 - 3] , where, when , when at this time, .

[0102] Equation [15 - 3] and equation [15 - 2] define a zero eigenvalue problem. After node numbering and flattening, it is transformed into an eigenvalue problem of a matrix to obtain the eigenvector field corresponding to the zero eigenvalue , that is , where, is the eigenvector of the eigenvalue problem, is the eigenvector field. In actual numerical calculations, due to errors, the eigenvalue is not exactly zero, and it is replaced by the eigenvector corresponding to the eigenvalue with the smallest absolute value.

[0103] After normalizing the eigenvector (as a preference, use maximum normalization, that is, take the maximum value of each component , and divide the eigenvector by this maximum value to obtain the normalized eigenvector), denote the normalized eigenvector field as .

[0104] (3) For order, [16 - 1] , , [16 - 2]

[0105] Among them, is the determinant of the deformation gradient tensor at the saddle point, which characterizes the volume ratio of the saddle point state and the state before deformation.

[0106] Satisfying equation [15 - 3], it can be decomposed into ;

[0107] Based on the solvability conditions of equations [16 - 1] and [16 - 2], an ordinary differential equation about is obtained:

[0108] wherein, 【17 - 1】 【17 - 2】 【17 - 3】

[0109] Step 302. Based on the interpolation of the eigenvalue closest to zero, and are approximately calculated

[0110] As can be seen from the process of step 301, it is necessary to extract the parameters and at the saddle point. The approximate position of the saddle point is calculated from the equilibrium path of step [2]. However, in actual operation, there is no way to ensure that the equilibrium state of the incremental step calculation is exactly at the saddle point position, unless an extremely small incremental step size is used. However, reducing the step size will significantly increase the total number of incremental steps, resulting in a substantial increase in the calculation cost.

[0111] The present invention proposes a method for approximately obtaining at the saddle point based on the interpolation of positive and negative eigenvalues closest to zero, and the specific method is as follows:

[0112] (1)In the calculation of the equilibrium path Riks step, extract the eigenvalues closest to zero in each incremental step to form a sequence, and denote the eigenvalues before and after the sign change as , respectively;

[0113] (2)Respectively take the incremental steps corresponding to and . Approximately assume that the results of these two incremental steps are both in the saddle point state, and perform two perturbation finite element calculations through the method of step 301 to obtain the corresponding , wherein, and are obtained by using the same normalization method;

[0114] (3)Adopt the following interpolation formula: , , .

[0115] Thus, the value of the saddle point approximation is obtained. Value.

[0116] Step [4]

[0117] Obtain the approximate transient response of the shallow arch saddle point jump from this step [4].

[0118] The specific method is as follows: Assume the initial conditions , establish the relationship between and :

[0119] , [17 - 4] , where The function is the Gaussian hypergeometric function. Take the limit , obtain the theoretical shallow arch jump time of the perturbation finite element

[0120] . [17 - 5]

[0121] Thus, the parameters can be extracted from Equation [17 - 1], Equation [17 - 2] and Equation [17 - 3], and the relationship between and is established by the Gaussian hypergeometric function [17 - 4]. Furthermore, the relationship between the transient displacement field and time is established by the approximate expression . Finally, the relationship between the transient displacement field and the real physical time is established by to realize the calculation of the approximate transient response of the shallow arch saddle point jump.

[0122] The present invention proposes a lightweight system for shallow arch saddle point perturbation finite element analysis based on Python / C hybrid programming, which includes five modules: (1) Parameter input and input file parsing module based on the abaqus platform; (2) Riks step analysis module; (3) Saddle point perturbation finite element analysis module; (4) Implicit dynamic analysis module; (5) Result visualization module based on the abaqus platform.

[0123] Module (1) is responsible for inputting the nodes and degrees of freedom directions of load loading, the nodes and degrees of freedom directions of geometric constraints, the finite element node coordinate table after importing the input file into the abaqus platform for modeling and mesh generation, and the node list included in the elements. This module is fully implemented in the Python language;

[0124] Module (2) is responsible for the Riks analysis step proposed in the present invention. Through the C / Python API method, a C language module with high execution efficiency is compiled and encapsulated into a Python module for Python to call. The C language part is mainly responsible for calculating the Jacobian stiffness matrix and the internal force residual vector of each element (corresponding to Equation 【5】). Python combines the Jacobian matrix and the internal force residual vector of the element into a global Jacobian matrix and a global internal force residual vector, converts the matrix into a sparse matrix form, and uses the spsolve function of the scipy library function in Python to solve the matrix equation, and uses the eigsh function of the sparse module of the scipy library in Python to calculate the eigenvalues of the corresponding matrix (the matrix corresponding to Equation

[10] ); Module (3) is responsible for the saddle point perturbation finite element analysis step proposed in the present invention. First, the C language is used to calculate the calculation part of each element involved in Equation 【15-3】, and the eigsh function in Python is used to obtain the eigenvector V, and the C language is used to calculate the part of each element involved in Equations 【17-2】 and 【17-3】. Finally, the Gaussian hypergeometric function is calculated using the hyp2f1 library function of scipy.special in Python; Module (4) is a module proposed in the present invention for directly calculating the shallow arch jump transient dynamics process, mainly used to verify the correctness of the results of the saddle point finite element module proposed in the present invention, and is divided into an initial acceleration calculation step and an implicit dynamics analysis step based on the GN22 algorithm; Module (5) is used for visualizing the results of each analysis module. In order to avoid repeated development, the visualization of the results of the perturbation finite element analysis system involved in the present invention is realized by means of the mature Abaqus python development platform.

[0125] Taking the visualization of displacement or deformation in the calculation results as an example, the specific implementation method of Module (5) is described as follows:

[0126] From the aforementioned input file containing node and element information, add an empty static analysis load step (i.e., without adding any loads), and submit it to the abaqus background for calculation to generate the corresponding result odb file;

[0127] Create a new analysis step step (implemented by the odb.Step command of abaqus python) in the odb file, and create an increment step (implemented by the Frame command of abaqus python) under this analysis step. Create a field variable under this increment step (implemented by the FieldOutput function of the abaqus python command)

[0128] (3) Organize the displacement fields obtained by each analysis module into a node number table and a displacement vector table, and finally import them into the field variables by the addData function of abaqus's python;

[0129] (4) After saving the odb file, the visualization function of abaqus can be used to view information such as the deformation nephogram of the shallow arch.

[0130] Example 1

[0131] Take Figure 5 the shown circular shallow arch with hinged ends as an example. The width of the shallow arch is 10 mm, the thickness is 1 mm, the radius of the upper surface is 1 m, and the radius of the lower surface is 999 mm (as Figure 6 shown), the shallow arch bears a central concentrated force, the shape of the shallow arch is symmetric about the YZ plane, and the line connecting the end face and the center of the circle makes a 10-degree angle with the Y-axis. Therefore, the shallow arch includes an angular range of 20 degrees; assume that the shallow arch is a linearly elastic material with an elastic modulus of 2.2 GPa, a Poisson's ratio of 0.3, and a density of .

[0132] This example aims to realize the deformation calculation of the shallow arch during the gradual increase of the central concentrated force, the accurate positioning of the saddle point based on eigenvalues, and the calculation of the dynamic parameters near the saddle point by the perturbation finite element method, and give the comparison and verification between the shallow arch jump transient response given by the perturbation finite element method and the direct finite element dynamic analysis.

[0133] (1) Establishment of the finite element mesh model of the shallow arch

[0134] Step 101. Finite element mesh division, boundary conditions, and loading conditions

[0135] In the Abaqus platform, due to the structural symmetry, establish half of the finite element mesh model as shown in Figure 5 the figure. Use C3D8R elements to divide regular hexahedral meshes. Among them, 4 meshes are evenly divided in the width direction, 4 meshes are evenly divided in the thickness direction, and 200 meshes are evenly divided in the circumferential direction. Therefore, the total number of meshes is 3200. Export the node coordinates and the node list included in the element through the input file of Abaqus. An example of the node information and element connection relationship information included in the input file is as follows:

[0136] Node list, each line contains the node number and the three-dimensional coordinates of X, Y, and Z before deformation.

[0137] 1, 0.173474535, 0.983822942, 0.00999999978 2, 0.17261605, 0.98397547, 0.00999999978 The list of unit connection relationships includes the unit node numbers and the numbers of the 8 nodes it contains. For example: 1, 1006, 1007, 1208, 1207, 1, 2, 203, 202 2, 1007, 1008, 1209, 1208, 2, 3, 204, 203

[0138] Figure 7 The nodes constrained in the central end face symmetric boundary condition are given, that is, the displacement of the nodes in the X direction needs to be guaranteed to be zero; Figure 8 The positions of the nodes constrained under the hinged boundary condition are given, that is, 4 nodes corresponding to the inner radius of the end face are taken, and the displacements of each node in the X, Y, and Z directions are set to zero; A downward concentrated force (-Y) is applied to 5 nodes on the upper surface of the central end face, and the force value on each node is -0.001 N, as Figure 9 .

[0139] (2) Calculation of the equilibrium path during the increase of the concentrated force

[0140] Each node (number) has three translational degrees of freedom , which respectively represent the degrees of freedom in three directions.

[0141] The Green-Lagrange strain is adopted , where is the deformation gradient tensor, represents the transpose of the tensor, represents the unit tensor.

[0142] It is obtained from the principle of virtual work: ,

[0143] where is the region occupied by the shallow arch before deformation, that is, the volume region of the entire finite element model, represents the variational operation, is the isotropic linear elastic stiffness tensor (determined by the elastic modulus and Poisson's ratio ), and its component form is: ," :" represents the contraction operation of tensors, represents the Lagrange multiplier that constrains the th node's th degree of freedom, Indicates the th degree of freedom of the th node, which is the degree of freedom defined by the boundary conditions, Indicates the reference concentrated force (unchanged in the equilibrium path) applied in the direction of the th degree of freedom of the th node, that is ( , see Figure 9 ), Indicates the load factor of the concentrated force (changing in the equilibrium path, characterizing the change in the magnitude of the concentrated force), Indicates the set of node numbers and degree-of-freedom pairs corresponding to the constraints in the boundary conditions, determined by the Figure 7 and Figure 8 nodes and related degrees of freedom shown in Indicates the set of node numbers and degree-of-freedom pairs where the concentrated force is applied, that is .

[0144] Derived from the above equations: ,

[0145] where, is an arbitrary three-dimensional vector field allowed by the finite element discretization, is at the node in the direction of the th degree of freedom.

[0146] From the global interpolation function, , where, Indicates the unit vectors in the X, Y, Z directions, that is , , .

[0147] Based on the arbitrariness, a tensor equation is obtained:

[0148] , where, Indicates the tensor product operation. When , , when , , and similarly, when , , when, take .

[0149] In addition, due to the arbitrariness, the equation corresponding to the boundary conditions is obtained, that is, when When .

[0150] The arc-length method (Riks) is adopted to solve the above tensor equation to obtain the equilibrium path. Figure 10 The corresponding relationship between the load factor and the total arc length in the equilibrium path of this example is given. The total arc length is a quantity without a clear physical meaning obtained in the Riks method, which characterizes the "curve length" of the load factor-displacement path from the start of loading to the current state. The arc length increment at each step is taken as a fixed value of 6.

[0151] From Figure 10 it can be seen that the equilibrium path of the shallow arch includes a load rising section and a load falling section, and the saddle point of the shallow arch corresponds to the maximum load point of the equilibrium path.

[0152] (3) Precise saddle point localization method based on matrix eigenvalues

[0153] According to the method proposed in the present invention for tracking the eigenvalue closest to 0, the specific implementation is as follows:

[0154] Figure 11 gives Figure 10 the distribution of the eigenvalues of the global matrix in each increment step in Figure 10 (taking 50 load steps, corresponding to the data points in

[0155] From Figure 11 it can be seen that when no treatment is performed, that is , there are multiple non-zero but extremely small eigenvalues in each gain step. These extremely small non-zero eigenvalues are caused by the Lagrange multiplier term. After taking a larger gain coefficient , the eigenvalue diagrams of each increment step are as shown in Figure 12 .

[0156] From Figure 12 it can be seen that when , the eigenvalue closest to 0 is no longer disturbed and can be clearly identified and extracted. Moreover, the corresponding variation law is as follows: First, the smallest positive eigenvalue decreases and decreases to 3.28e-5 at the 14th increment step. Subsequently, at the 15th increment step, it significantly decreases to a negative value Figure 10 marks the smallest eigenvalues of some gain steps obtained therefrom. This embodiment shows that by setting a larger gain , the eigenvalue closest to zero can be effectively obtained to obtain the approximate position of the saddle point.

[0157] The numerical value of the load factor at the saddle point position of this embodiment is approximately (corresponding to Figure 10in the 14th increment step), that is, when a is applied in the Y-axis direction of the central end face, the saddle point position is reached, and the shallow arch is in a critical equilibrium state at this time. At this time, the three-dimensional displacement vector of the central node of the central end face is calculated as , that is Figure 9 the node No. 2413 in Figure 14 . The original undeformed configuration is as shown in Configuration 1 in Figure 14 , and the deformed configuration at the saddle point is as shown in Configuration 2 in

[0158] (4) Extraction of saddle point jump parameters based on perturbation finite element

[0159] When the load at the saddle point changes from to , where is the additional load magnitude coefficient relative to the saddle point, is a fixed quantity related to the direction of the additional load ( represents the vector of the additional load). In this example, , that is ( , see Figure 9 ), and ( ). Obviously, determines the magnitude of the additional load, and here is taken.

[0160] Assume that the vector field satisfies the perturbation expansion condition: ,

[0161] Lagrange multiplier perturbation expansion: ,

[0162] where is the flattened vector of all , and denote , is the flattened vector of the Lagrange multiplier at the saddle point, is the coefficient vector of each term of the perturbation expansion, is the physical real time and the perturbation parameter the scaled time parameter obtained by multiplying, is the displacement vector field, is the displacement vector field at the saddle point, is the coefficient related to the perturbation expansion and time .

[0163] In this embodiment, to illustrate the execution steps of the perturbed finite element method, the saddle point is approximately taken at the increment step 14. Then, in the above equations, is the displacement field at the end of the increment step 14, is the vector composed of all Lagrange multipliers at the increment step 14.

[0164] Further calculate the matrix eigenvalue problem determined by Equation [15-3] and Equation [15-2], and then obtain the eigenvector field , that is , where is the eigenvector of the matrix eigenvalue problem, is the eigenvector;

[0165] The magnitude of the eigenvector is arbitrary. Therefore, normalization is required. In this embodiment, the maximum value normalization is adopted, that is, take the maximum value of each component , and divide the eigenvector by this maximum value to obtain the normalized eigenvector, denoted as . Since, in this embodiment, among the eigenvectors, the absolute value of the displacement degree of freedom in the Y direction of the nodes at the central end face is the largest. Therefore, in the normalized eigenvector, the Y direction degree of freedom of the central end face is 1.

[0166] Finally, carry out the perturbed finite element analysis step, and obtain the parameters , . To clarify the practicality of the calculation, carry out a direct implicit dynamics finite element analysis, using the GN22 algorithm, that is, the relationships of displacement, velocity, and acceleration satisfy Equation [18-1] and Equation [18-2].

[0167] 【18-1】 【18-2】

[0168] Among them, take , represents time the displacement field, velocity field, and acceleration field at time represents the displacement field, velocity field, and acceleration field at time is the time step.

[0169] Thus, only need to calculate the unknown quantity , and then the displacement and velocity at time can be obtained.

[0170] Due to the additional load Resulting in an initial acceleration, to accurately calculate the dynamic response, the initial acceleration needs to be correctly calculated. Therefore, first calculate the initial acceleration after loading, and the specific method is as follows:

[0171] (1) As can be seen from Equation

[11] , at the moment when the load increases at the saddle point, the following variational equation is satisfied:

[19]

[0172] and the acceleration constraint condition When , where the acceleration constraint condition is caused by the displacement constraint condition ( ). By using the finite element interpolation formula for acceleration,

[0173] where is the th node's th acceleration component, is the Lagrange multiplier at the saddle point. Through the variational principle, a linear equation system about the unknowns and is obtained. Furthermore, after changing the variables' shape and flattening, by solving a matrix equation, the initial acceleration field can be obtained.

[0174] (2) The initial acceleration field obtained by the above method is . Set the zero velocity condition , (the displacement field is the same as that at the saddle point), then according to Equation [18-1] and Equation [18-2], combined with Equation

[11] , the acceleration field at each moment can be solved, and then the displacement field and velocity field at each moment can be obtained.

[0175] Figure 13 respectively gives the comparison diagram of the results of the perturbed finite element calculation and the direct implicit dynamics calculation of the present invention for the radial displacement at the center of the shallow arch ( Figure 9 the displacement increment in the -Y direction of the 2413 node relative to the saddle point state in ), where, , , is the load at the saddle point ( ), thus, the load during the jump increases by compared with the saddle point load.

[0176] Figure 14 gives the shape of the global deformation obtained by the direct implicit dynamics finite element calculation, corresponding to the configurations 2, 3, 4, 5, 6 marked in Figure 13 .

[0177] As can be seen from the implicit dynamics finite element calculation in Figure 13 : at time 0 s, it corresponds to the saddle point; as time increases, dynamic jumps occur and the displacement increases significantly. At , the displacement is the largest; subsequently, the displacement oscillates reciprocally. Since no damping is added to the implicit dynamics finite element model, it can be predicted that the structure will oscillate significantly infinitely. In an actual structure, due to the existence of damping, the structure usually gradually stabilizes after slight vibration after the displacement reaches the maximum. Therefore, the time interval from the start to the first time the displacement reaches the maximum is the most physically meaningful. During this time interval, the perturbation finite element method proposed by the present invention can accurately predict the displacement-time response. At the same time, the jump time predicted by the perturbation finite element and the time when the displacement first reaches the maximum are in good agreement.

[0178] The present invention is applicable to different , and there is no need to carry out dynamic finite element analysis for each case. Only perturbation finite element analysis needs to be carried out to obtain the corresponding results of different : As shown in Figure 15 , the displacement-time path at is given (i.e., the -Y direction displacement of node 2413 relative to the saddle point state). Thus, the perturbation finite element results and the direct dynamic finite element results are in good agreement within the time interval from the start to the first time the displacement reaches the maximum. At the same time, the jump time obtained by the perturbation finite element is close to the corresponding time when the displacement first reaches the maximum, indicating that the perturbation finite element method of the present invention can approximately obtain the dynamic response.

[0179] It should be noted that a series of results can be obtained by perturbation finite element calculation, and there is no need to carry out perturbation finite element analysis multiple times for . For direct dynamic finite element analysis, calculations need to be carried out for each . At the same time, in each calculation, iterative calculations for multiple time steps are required, which is time-consuming and cumbersome, posing an efficiency challenge to the optimal design of thin-walled arched structures.

[0180] This embodiment verifies the effectiveness and feasibility of the perturbation finite element method proposed by the present invention.

[0181] Example 2

[0182] In Embodiment 1, the fixed arc length increment for each step is taken as 6, and there is exactly a very small eigenvalue of 3.28e-5 in the equilibrium path calculation. Therefore, this point can be used as an approximate saddle point to carry out the perturbation finite element analysis calculation. However, in actual calculation, since it is not easy to exactly obtain a point with a very small eigenvalue in the equilibrium path calculation step, there is a difficulty in saddle point positioning, unless the arc length increment taken for each increment step is extremely small. At this time, a large number of increment steps and high computational costs are required.

[0183] In order to further improve the accuracy of the saddle point, the present invention proposes a method of linear interpolation of eigenvalues, which is as follows:

[0184] It is only necessary to perform two perturbation finite element calculations on the eigenvalues (denoted as and ) before and after the change of the eigenvalue sign, and respectively obtain the normalized (using the same maximum value normalization method as in Embodiment 1) and the corresponding parameters and . Assume that the at the saddle point is determined by the following formula: , , .

[0185] Using the same parameters and equilibrium path calculation as in Embodiment 1, for the sake of illustration, a perturbation finite element analysis is carried out for each increment step to obtain the parameter .

[0186] Figure 16 and Figure 17 illustrate the rationality of the interpolation method of the present invention:

[0187] It can be seen that when the eigenvalue is near 0, both approximately change linearly with the eigenvalue. Extract the results of 6 increment steps (increment steps 10, 11, 12, 13, 14, 15) in Figure 10 . It can be seen that the last positive eigenvalue is , and the first negative eigenvalue is (increment step 15), . Then, according to the interpolation formula, the

[0188]

[0189] Since the eigenvalue of increment step 14 relative to the eigenvalue of increment step 15 is very close to zero. Therefore, the numerical value of the true saddle point is approximated to that of increment step 14 and , which demonstrates the effectiveness of the interpolation method proposed by the present invention.

[0190] Embodiment 3

[0191] To further illustrate the effectiveness of the eigenvalue interpolation method based on the closest-to-zero eigenvalue proposed by the present invention, the same parameters and model as in Embodiment 1 are adopted, but the arc length of each equilibrium path calculation step is increased (the increment step size is larger), making it artificially difficult to accurately obtain the saddle point. It is hoped that through this embodiment, it can be shown that even with a larger increment step size, the interpolation proposed by the present invention can still accurately obtain the value at the saddle point.

[0192] Figure 18 shows the relationship diagram between the load factor and the total arc length when the fixed arc length increment per step is 10. Compared with Embodiment 1 (arc length increment is 6), the arc length increment is significantly increased, indicating that the step size of each increment step is significantly increased. It can be seen that after increasing the arc length increment, the number of gain steps reaching the highest point decreases, shortening the calculation time. The numerical value corresponding to each increment step in the figure represents the eigenvalue closest to zero: as the loading progresses, this eigenvalue continuously decreases, from 3.24e-4 to the negative number -3.89e-4; due to the larger arc length increment, compared with Embodiment 2, the eigenvalue is larger, and using any calculation point as an approximate saddle point will inevitably bring a larger error.

[0193] According to the interpolation method based on the eigenvalue closest to zero proposed by the present invention, only the information of two points with eigenvalues of 3.24e-4 and -3.89e-4 needs to be calculated to approximately obtain the accurate result. Similarly, for the convenience of explanation, the perturbation finite element analysis is still carried out for each increment step, and the corresponding parameters are extracted. The results are as shown in Figure 17 and Figure 18 .

[0194] From the eigenvalue closest to zero, it can be seen that the sign change occurs between increment step 8 and increment step 9: In increment step 8, the eigenvalue is , , ; In increment step 9, the eigenvalue is , , .

[0195] Thus, from the interpolation formula, the approximate parameters of the saddle point:

[0196]

[0197] and the values in Example 2 and In comparison, the values are very close. Compared with the eigenvalue situation in Example 2, the absolute values of the eigenvalues of increment step 8 and increment step 9 are both larger. This shows that the linear interpolation method based on the eigenvalue closest to zero proposed by the present invention can be robustly applied to different increment step lengths, and has the following advantages: a larger step length can be selected without ensuring that the increment step exactly reaches or is very close to the saddle point position (in fact, this is very difficult to achieve unless a very small step length is used, and the reduction of the step length will increase the number of required increment steps and the calculation amount), greatly reducing the calculation amount of the equilibrium path calculation steps.

[0198] Example 4

[0199] Based on the perturbation finite element method proposed by the present invention, since the finite element discretization method is adopted, it can be applied to flexible thin-walled arched structures with non-uniform wall thickness distributions. To illustrate this, the shallow arch loading method and geometric dimensions of this Example 4 are exactly the same as those of Example 1 except that the wall thickness is non-uniformly distributed.

[0200] In Example 1, the wall thickness is uniformly 1 mm, while in this example, the wall thickness is non-uniformly distributed: it varies linearly with the angle within a 10-degree range, with the minimum wall thickness of 0.7 mm at the central end face and the maximum wall thickness of 1 mm at the hinged end. The outer radius of the model remains unchanged at 1 m.

[0201] Comparison Figure 21 and Figure 10 , it can be seen that due to the reduction of the wall thickness, the maximum load that the shallow arch can bear becomes smaller: Figure 21 The maximum load factor in Figure 8 is less than 12, while the load factor in

[0202] Figure 22 shows the changes in the eigenvalue closest to zero and the corresponding parameter at increment steps 7, 8, 9, 10, and 11. Based on the eigenvalues, interpolation can be performed on the results of the normalized eigenvectors of increment step 8 (i.e., ) and the normalized eigenvector of increment step 9 (i.e., ) to obtain the eigenvector at the final saddle point: Figure 24shows the initial shape of the shallow arch, and the shapes of the shallow arch at increment step 8 and increment step 9 after superposing the appropriately scaled eigenvectors (realized by the visualization function of the Abaqus platform in the lightweight system architecture proposed by the present invention), that is, the scaled displacement field is added. (Increment step 8), (Increment step 9), where (for convenient display of deformation). Therefore, the normalized eigenvectors of increment step 8 and increment step 9 are very close, and the displacement field corresponding to the interpolated eigenvector should be between the two.

[0203] Subsequently, based on the obtained by interpolation, the transient jump process of the shallow arch can be obtained. The process is the same as the steps in Embodiment 1 and will not be elaborated here. This embodiment shows that the method proposed by the present invention can adapt to the analysis of shallow arches with non-uniform wall thickness and has good robustness.

[0204] Finally, it should be noted that in order to illustrate that the parameter changes approximately linearly with the eigenvalue, Figure 22 , 23 , Figure 19 , 20 and Figure 16 , 17 The perturbation finite element analysis method of the present invention is used to calculate the parameters of multiple increment steps. , however, in actual implementation, it is only necessary to carry out perturbation finite element analysis and calculation for two increment steps before and after the change of the eigenvalue sign, and the corresponding parameters of the saddle point can be obtained by the interpolation method. Therefore, compared with the case where direct implicit dynamics analysis requires iterative calculations of multiple time increment steps, the calculation amount of the method of the present invention is significantly reduced.

[0205] The saddle point finite element perturbation method for approximately obtaining the dynamic response based on the quasi-static equilibrium path analysis of the shallow arch proposed in the embodiment of the present invention avoids direct dynamic jump analysis and provides an effective analysis tool for the innovative application of flexible shallow arches (see Embodiment 1).

[0206] The embodiment of the present invention combines the perturbation method and the finite element discretization method, adopts the Green strain measure applicable to large deformation, and establishes a saddle point perturbation finite element analysis method for flexible shallow arches. Due to the adaptability of the finite element discretization method in geometric shape, the present invention is applicable to special cases such as shallow arches with non-uniform wall thickness (see Embodiment 3) and provides an analysis tool for the geometric optimization design of flexible shallow arches.

[0207] The interpolation method based on the eigenvalue closest to zero proposed in the embodiment of the present invention effectively solves the problem that the saddle point cannot be accurately reached in the equilibrium path analysis step. Without setting an extremely small incremental step size, the approximate parameters at the saddle point can still be obtained, which significantly reduces the calculation amount of the equilibrium path analysis step (see Example 2 and Example 3).

[0208] Using Lagrange multipliers to add boundary conditions is a conventional method. However, the existence of Lagrange multipliers leads to the existence of multiple non-zero but extremely small eigenvalues ​​in the related stiffness matrix. These eigenvalues ​​are irrelevant to the saddle point and seriously interfere with the extraction of the eigenvalue closest to zero. Therefore, the embodiment of the present invention proposes to multiply the matrix related to the Lagrange multiplier by a large gain coefficient in blocks, thereby avoiding the interference of related eigenvalues ​​(see Example 1) and making eigenvalue-based interpolation possible.

[0209] The present invention proposes a lightweight system framework for shallow arch saddle point perturbation finite element analysis based on Python / C hybrid programming. The Python / C API method is adopted, and the intensive calculation parts such as the stiffness Jacobian matrix and the internal force residual vector of each unit are processed in C language, and compiled into Python modules, which can overcome the problem of low calculation efficiency of Python language. At the same time, the Python API provided by the Abaqus platform is cleverly used to implement the method of using the mature Abaqus platform as a calculation result visualization platform of this system. Finally, a method of solving relevant equations and eigenvalues ​​by using the relevant functions of the scipy library provided by Python is proposed. The system architecture design can ensure the lightweight of the relevant finite element analysis system. At the same time, the relevant modules compiled in C language realize the closed source and source code protection of the system, and can be used for commercial development of relevant systems.

[0210] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modification, equivalent substitution or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A finite element approximate calculation method for a thin-walled arched structure, characterized in that, At least the following steps are included: S1, thin-walled arch geometry model and finite element mesh generation; S2, thin-walled arch equilibrium path tracking and eigenvalue extraction based on Riks method; S3, perturbation finite element calculation of thin-walled arch saddle point and parameter interpolation based on two eigenvalues; S4. Approximate calculation of transient response based on Gaussian hypergeometric function.

2. The finite element approximate calculation method for a thin-walled arched structure according to claim 1, characterized in that, In step S1, a three-dimensional geometric model of a thin-walled arch is established according to the thin-walled arch set parameters, and an eight-node hexahedral solid finite element mesh is divided, and node and unit connection information is output.

3. The finite element approximate calculation method for a thin-walled arched structure according to claim 2, characterized in that, S2 includes the following steps: a. Define the global field variable interpolation function. The specific method is as follows: (1) Using eight-node hexahedral elements, the following Lagrangian interpolation method is adopted; (2)Let the three-dimensional coordinate vector before deformation be , and define in the form of interpolation as ； From the defined and one-to-one correspondence, are all regarded as functions, ; (3) Uniformly number all the nodes , where represents the number of nodes in the finite element model of the thin-walled arch structure. The global field variables within the final thin-walled arch structure are determined by the following interpolation formula: ， Among them, represents the value of the global field variable at the location, is the global field interpolation function, is the value of the field variable at the th node; b. Calculate the equilibrium path under external load; c. Extract the eigenvalue closest to zero for each arc length method increment.

4. The finite element approximate calculation method for a thin-walled arched structure according to claim 3, characterized in that, S3 includes the following steps: a. Calculate the approximate transient response of the saddle point jump of the thin-walled arch structure. b. Based on interpolation of the eigenvalue closest to zero Approximate calculation Obtained at the saddle point by interpolating and approximating based on the positive and negative eigenvalues closest to zero , and the specific implementation method is as follows: After the calculation of the quasi-static equilibrium path, the eigenvalue closest to zero at each increment step is extracted and arranged in sequence. In the sequence, the eigenvalues before and after the sign change are denoted as , ; (2)Respectively take and corresponding increment steps. Assuming that the results of these two increment steps are both in the saddle point state, carry out two perturbation finite element calculations to respectively obtain the corresponding , where and are obtained by using the same normalization method; (3) Use the following linear interpolation formula: ， ， ； Thus, the value of the saddle-point approximation is obtained.

5. The finite element approximate calculation method for a thin-walled arched structure according to claim 4, characterized in that, The specific calculation method in S4 is as follows: Assume the initial conditions , establish and relationship: , The function is a Gaussian hypergeometric function; Taking the limit to obtain the approximate jump time of the perturbed finite element thin-walled arch structure ； Furthermore, from the approximate expression , the transient displacement field and time are related. Finally, from , the relationship between the transient displacement field and the true physical time is established, and the approximate transient response calculation of the saddle point jump of the thin-walled arch structure is realized.

6. A finite element approximate calculation method for a thin-walled arched structure as claimed in claim 5, characterized in that, Obtain the eigenvector field corresponding to the zero eigenvalue , and take the eigenvector corresponding to the eigenvalue closest to zero as the approximation of the eigenvector field in the calculation.

7. A finite element approximate calculation method for a thin-walled arched structure as claimed in claim 6, characterized in that, Normalize the feature vector using maximum normalization, that is, take the maximum value of each component and divide the feature vector by this maximum value to obtain the normalized feature vector.

8. A lightweight system, characterized in that, The finite element approximate calculation method as described in claim 7 is used, based on Python / C hybrid programming, and includes the following modules: input file parsing module, arc length method analysis module, saddle point perturbation finite element analysis module, implicit dynamic analysis module and result visualization module.

9. A lightweight system as claimed in claim 8, characterized in that, The input file parsing module is responsible for inputting the nodes and degrees of freedom directions of load loading, the nodes and degrees of freedom directions of geometric constraints, and the finite element node coordinate table and the node list contained in the unit after the input file is imported into the abaqus platform for modeling and meshing, and is implemented in Python language; the arc length method analysis module is responsible for the Riks analysis step proposed in the present invention. Through the C / Python API method, a C language module with high execution efficiency is compiled and encapsulated into a Python module for Python calling.

10. A lightweight system as claimed in claim 9, characterized in that, The saddle point perturbation finite element analysis module first uses C language to calculate the calculation part of each unit involved in the equation, uses Python's built-in eigsh function to obtain the eigenvector V, and finally uses Python's built-in library function to calculate the Gaussian hypergeometric function; the implicit dynamic analysis module directly calculates the transient dynamic process of the thin-walled arch structure jump, which is mainly used to verify the correctness of the results of the saddle point finite element module; the result visualization module uses the mature Abaqus python development platform to realize the visualization of the results of each analysis module.

Images