A finite element approximation calculation method for a thin-walled arched structure and a lightweight system

Abstract

The present invention relates to the field of simulation analysis of the dynamic jumping behavior of thin-walled arched structures, and particularly relates to a finite element approximation calculation method and a lightweight system for thin-walled arched structures. Specifically, the present invention proposes a finite element approximation calculation method and a lightweight system for thin-walled arched structures. First, a one-time quasi-static equilibrium path analysis is carried out to obtain the eigenvalue closest to zero at each increment step, identify the eigenvalues before and after the sign change, and carry out two perturbation finite element analyses to obtain approximate transient dynamic jumping parameters and eigenvectors. Based on the eigenvalue closest to zero, approximate jumping parameters and eigenvectors at the saddle point are interpolated. Finally, an approximate transient jumping response is obtained by the Gaussian hypergeometric function, and a Python / C hybrid programming method is proposed to realize the lightweight of the finite element analysis system, which can effectively reduce the calculation cost and analysis difficulty of the transient jumping analysis of the flexible shallow arch saddle point.

Description

Technical Field

[0001] The present invention relates to the technical field of microelectromechanical intelligent applications and finite element simulations of flexible structures, 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 relatively low material stiffness and prone to elastic deformation. Such structures are prone to deformation, and this deformation does not represent structural failure but is what is desired artificially. It 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, which is the load-deformation equilibrium path diagram of a 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 occurs at the center of the thin-walled arch, and the system remains in an equilibrium state. The deformation 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 flexible thin-walled arch jump makes it convenient 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 through 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, 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 a flexible thin-walled arch saddle-point type jump and a related lightweight finite element system development method, which can obtain an approximate transient dynamic jump response without conventional dynamic analysis, only through 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 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:

[0006] S1. Generation of a geometric model and finite element mesh of the thin-walled arch;

[0007] S2. Tracking of the equilibrium path and extraction of eigenvalues of the thin-walled arch based on the Riks method;

[0008] S3. Finite element calculation of saddle point perturbation of the thin-walled arch and parameter interpolation based on two eigenvalues;

[0009] S4. Approximate calculation of the transient response based on the Gaussian hypergeometric function.

[0010] 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.

[0011] Preferably, S2 includes the following steps:

[0012] a. Define the interpolation function of the global field variable, and the specific method is as follows:

[0013] (1) Use the eight-node hexahedron element and adopt the following Lagrangian interpolation method:

[0014] ,

[0015] ,

[0016] ,

[0017] ,

[0018] ,

[0019] 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;

[0020] (2) Let the three-dimensional coordinate vector before deformation be , and define The interpolation form is

[0021] ;

[0022] The one-to-one correspondence defined by Equation [2] and is regarded as a function of ; ;

[0023] (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 in the final thin-walled arch structure are determined by the following interpolation formula:

[0024] ,

[0025] where represents the value of the global field variable at , is the global field interpolation function, is the field variable value at the -th node;

[0026] b. Calculate the equilibrium path under the external load. The specific calculation method is as follows:

[0027] Adopt the Green-Lagrange strain , where is the deformation gradient tensor, represents the transpose of the tensor, represents the unit tensor;

[0028] Based on the principle of virtual work, establish the variational equation:

[0029] ,

[0030] where represents the degrees of freedom in the X, Y, Z directions, represents the node number, is the region occupied by the thin-walled arch structure before deformation, represents the variational operation, is the isotropic linear elastic stiffness fourth-order tensor, and ":" represents the contraction operation of the tensor, 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, The numerical value of the degree of freedom defined by the boundary conditions represents the load acting force value applied in the direction of the -th degree of freedom of the -th node, is the load factor characterizing 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 to which the load acting force is applied;

[0031] Based on the variational principle, a discrete equilibrium equation is established:

[0032] ,

[0033] where, represents the tensor product operation. When , when , when , , when, take ; defines the constraint conditions; taking the load factor and the displacement degree of freedom as unknowns simultaneously, the arc-length method is used to solve the discrete equilibrium equation to obtain the equilibrium path under the external load ;

[0034] c. Extract the eigenvalue closest to zero for each arc-length method increment step. The specific method is as follows:

[0035] Based on the and linearization, a tensor equation in the following form is obtained (for all repeated indices, the Einstein summation convention is adopted, is the Kronecker symbol ( ), represents the differential):

[0036] When :

[0037] ,

[0038] ;

[0039] When :

[0040] ,

[0041] Through the flattening and reshaping operations of the tensor, it is transformed into matrix form

[0042]

[0043] Among them, represents the flattened 1D column vector, and let the dimension be , represents the remaining flattened 1D column vector, and let the dimension be , represents the flattened 1D column vector, and the dimension is 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 applied force;

[0044] Set a large gain coefficient , and for each increment step, solve the eigenvalues of the modified matrix , eliminate the interference of the extremely small non-zero eigenvalues caused by the Lagrange multiplier term, and extract the eigenvalue closest to zero at this increment step.

[0045] Preferably, S3 includes the following steps:

[0046] a. Calculate the approximate transient response of the saddle point jump of the thin-walled arch structure, and the specific method is as follows:

[0047] Based on D'Alembert's principle, establish a variational equation:

[0048] ,

[0049] Among them, 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, represents the determinant of the deformation gradient tensor;

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

[0051] ,

[0052] Assume that the Lagrange multipliers satisfy the perturbation expansion:

[0053] ,

[0054] where, is the flattened vector of all and is denoted as , is the flattened vector of the Lagrange multipliers at the saddle point, is the coefficient vector of each term in the perturbation expansion, is the physical real time and the perturbation parameter multiplied to obtain the scaled time parameter, is the transient displacement vector field, is the displacement vector field at the saddle point, is the coefficient of the perturbation expansion;

[0055] Based on the order expansion variational equation, the corresponding equations of each order are obtained:

[0056] (1) Equation of order:

[0057] and ;

[0058] (2) Equation of order:

[0059] ,

[0060] where, , are variational components, represents the unit vectors in the X, Y, Z directions; Fs is the deformation gradient tensor at the saddle point;

[0061] Convert the equation order equation to the eigenvalue form, and for all repeated indices, use the Einstein summation convention to sum:

[0062] , when , ;

[0063] Based on The equation of the order, after node numbering, tensor reshaping, and flattening, is transformed into an eigenvalue problem of a matrix to obtain the eigenvector field corresponding to the zero eigenvalue , that is

[0064] ,

[0065] where, is the eigenvector of the eigenvalue problem, is the eigenvector field; after normalizing the eigenvector , the normalized eigenvector field is denoted as ;

[0066] (3) The equation of the order is:

[0067]

[0068] , ,

[0069] where, is the determinant of the deformation gradient tensor at the saddle point;

[0070] Decompose into ; based on the solvability condition of the equation of the order, obtain the ordinary differential equation about :

[0071] , where,

[0072] ,

[0073] ;

[0074] b. Based on the interpolation of the eigenvalue closest to zero approximate calculation

[0075] 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:

[0076] (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 , ;

[0077] (2)Respectively take and For the corresponding increment steps, assuming that the results of these two increment steps are both in the saddle point state, two perturbation finite element calculations are carried out to obtain the corresponding , where and are obtained by using the same normalization method;

[0078] (3) The following linear interpolation formula is used:

[0079] , , ;

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

[0081] Preferably, the specific calculation method in S4 is as follows:

[0082] Assume the initial condition , establish the relationship between and :

[0083] , where the function is the Gaussian hypergeometric function;

[0084] Take the limit , and obtain the approximate jump time

[0085] of the perturbed finite element for the thin-walled arched structure;

[0086] Furthermore, from the approximate expression , establish the relationship between the transient displacement field and time , and finally establish the relationship between the transient displacement field and the real physical time

[0087] Preferably, obtain the eigenvector field corresponding to the zero eigenvalue, and in the calculation, take the eigenvector corresponding to the eigenvalue closest to zero as the approximation of the eigenvector field .

[0088] Preferably, for the normalization of the eigenvector , 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.

[0089] To achieve the above objectives, the present invention proposes a lightweight system that uses the above-mentioned finite element approximate calculation method. Based on Python / C hybrid programming, the system 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.

[0090] 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.

[0091] 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 module of 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.

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

[0093] 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 the Lagrange multiplier term. The closest eigenvalue to zero is obtained for each incremental step, and the eigenvalues before and after the sign change are identified. Two perturbation finite element analyses are carried out to obtain approximate transient dynamic jump parameters and eigenvectors. Based on the eigenvalues closest to zero, the approximate jump parameters and eigenvectors at the saddle point are obtained by interpolation. Finally, the approximate transient jump response is obtained by the Gaussian hypergeometric function, which can effectively reduce the computational cost and analysis difficulty of the transient jump analysis of the flexible thin-walled arch saddle point. A Python / C hybrid programming method is proposed, which can effectively improve the computational efficiency of computationally intensive parts such as the unit stiffness matrix. Matrix calculations are implemented through the Python scientific computing library, and the system results are visualized using the Abaqus platform, providing a certain reference value for the lightweighting of related finite element analysis systems.

[0094] 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.

[0095] 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

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

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

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

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

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

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

[0102] Figure 7 . Node diagram of constraints in the symmetric boundary condition.

[0103] Figure 8 . Node diagram of constraints in the hinged boundary condition.

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

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

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

[0107] Figure 12 . Gain coefficient K = 10 4 Eigenvalue distribution at that time.

[0108] Figure 13 . Comparison of the calculation results of the perturbation finite element and the implicit finite element.

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

[0110] 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.

[0111] Figure 16 . The c0 values for different increment steps.

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

[0113] Figure 18 . Relationship diagram between 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).

[0114] Figure 19 . The eigenvalue and the parameter c0 at increment steps 7, 8, 9, and 10.

[0115] Figure 20 . The eigenvalue and the parameter c1 at increment steps 7, 8, 9, and 10.

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

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

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

[0119] Figure 24 . Initial shape of the shallow arch, and the shapes of the shallow arch at increment step 8 and increment step 9 after superposing the scaled eigenvectors. Specific implementation mode

[0120] To make the purpose, technical solutions, and advantages of the present invention clearer, the present invention will be further described in detail below through 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 description of well-known structures and technologies is omitted to avoid unnecessarily confusing the concept of the present invention.

[0121] The 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 approximation calculation method for flexible thin-walled arched structures is adopted, and the thin-walled arched structure can also be called a shallow arch structure.

[0122] Refer toFigure 2 , the finite element approximate calculation method for the flexible thin-walled arch structure at least includes the following steps:

[0123] 【1】 Generation of the geometric model and finite element mesh of the flexible shallow arch;

[0124] 【2】 Tracking of the equilibrium path of the flexible shallow arch and extraction of the eigenvalue closest to zero based on the Riks method;

[0125] 【3】 Finite element calculation of the saddle point perturbation of the flexible shallow arch and parameter interpolation based on two eigenvalues;

[0126] 【4】 Approximate calculation of the transient response based on the Gaussian hypergeometric function.

[0127] Step 【1】

[0128] 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 included in the element, the loading force information, and the constraint node information into the perturbation finite element analysis system of the present invention.

[0129] Step 【2】 specifically includes Step 201, Step 202 and Step 203.

[0130] Step 201. Field variable interpolation based on the eight-node hexahedron element

[0131] Adopt the eight-node hexahedron element, such as Figure 4 shows a schematic diagram of the standard element of a regular hexahedron with a side length of 2, 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.

[0132] Adopt the following Lagrangian interpolation formula

[0133] 【1】

[0134] ,

[0135] ,

[0136] ,

[0137] ,

[0138] Among them, is the Lagrangian coordinate of the standard element, is the field variable at the coordinate , are the field variable values at 8 nodes, and are the interpolation functions of each node inside the element.

[0139] For non - standard elements, let the three - dimensional coordinate vector before deformation be , The interpolation form is the same as that of the field variable:

[0140]

[0141] Inside each element, Equation 【2】defines a and one - to - one correspondence. Thus, the field variable and the interpolation function inside the element can both be regarded as functions, that is: .

[0142] By uniformly numbering all the nodes , where, represents the number of nodes in the model. The global field variable inside the shallow arch can be obtained by the following interpolation formula

[0143] , 【3】

[0144] For each element, each node has internal numbers 1, 2, 3, 4, 5, 6, 7, 8 inside this element.

[0145] However, the finite - element model consists of multiple elements, and all the nodes have global numbers .

[0146] For each node , find all the elements that contain this node (denote the element numbers as ).

[0147] Inside each element , the internal number corresponding to the node in the element is a number among 1, 2, 3, 4, 5, 6, 7, 8. Therefore, the internal number of this node in the element is a definite function of and , which may be denoted as . Obviously, .

[0148] Then inside the element , take (This function is determined by the aforementioned element internal interpolation function). Note that the domain of this function is the element internally, and its domain can be extended to the entire finite element model by zero-value extension, that is, when is within the element internally, define , and when is not within the element internally, ; in this way is a function defined on the entire finite element model;

[0149] Define the global interpolation function of node , that is, traverse all elements (denoted by the number ) that contain node and accumulate the corresponding global function . Obviously is a function related to the global node number defined on the entire finite element model. For example, assume the global node , and the numbers of the elements containing this node are (that is, element 2 and element 3 contain this node), then define ;

[0150] where 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.

[0151] Step 202. Calculation of the equilibrium path under external loads

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

[0153] Adopt Green-Lagrange strain , where is the deformation gradient tensor, represents the transpose of the tensor, represents the unit tensor. From the principle of virtual work, we get:

[0154] , 【4】

[0155] where is the region occupied by the shallow arch before deformation, represents a variational operation, is the isotropic linear elastic stiffness tensor (uniquely determined by the elastic modulus and Poisson's ratio), and ":" represents the contraction operation of tensors, represents the constraint on the - th node's - th degree of freedom's Lagrange multiplier, represents the - th node's - th displacement degree of freedom, is the numerical value of the displacement degree of freedom defined by the boundary condition, represents the load acting force applied in the direction of the - th node's - th degree of freedom (unchanged during the equilibrium path), represents the load factor of the concentrated force (changing during the equilibrium path, characterizing the change in the magnitude of the concentrated force), represents the set of node numbers and degree - of - freedom pairs corresponding to the constraints in the boundary condition, represents the set of node numbers and degree - of - freedom pairs where the load acting force is applied. The calculation of the equilibrium path requires determining the deformation of the shallow arch under static equilibrium under different .

[0156] From Equation [4] and the variational principle, we get:

[0157] ,

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

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

[0160] From the variational principle, we obtain the equation:

[0161] , 【5】

[0162] where, represents the tensor - product operation. When , , when , , similarly, when , , when, take ;

[0163] In addition, due to the arbitrariness, the equations corresponding to the boundary conditions are obtained, that is, when ,

[0164] . 【6】

[0165] Taking the load factor and the displacement field as unknowns, the arc-length method (Riks) is used to obtain the solutions of Equation 【5】 and Equation 【6】, and the equilibrium path under the external load is obtained.

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

[0167] Based on the saddle point being a singular point, the global 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.

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

[0169] When ,

[0170] 【7】

[0171] ; 【8】

[0172] When : 【9】

[0173] The above equations are written in the following matrix form:

[0174]

[10]

[0175] Among them, represents the 1D column vector obtained by flattening the constrained (assuming the dimension is ), represents the 1D column vector obtained by flattening the remaining (assuming the dimension is , represents The flattened 1D column vector (the dimension is equal to ), is the identity matrix, is a square matrix ( order), is order matrix, is order matrix, is order matrix, is a zero matrix block, is the flattened one-dimensional column vector related to the loading force.

[0176] The finite element model consists of nodes, and each node has three degrees of freedom. For the convenience of expression, therefore, the degrees of freedom involved in the above equations , that is, there are degrees of freedom, 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, all non-homogeneous linear equations can be written in matrix form );

[0177] For the purpose of solving, the equations involved need to be transformed into matrix form, that is, in the form of an equation like (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);

[0178] 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 independent of , can also be flattened into the form of a column vector . The following explains the relevant steps.

[0179] The specific steps are as follows:

[0180] Number the elements of set B as follows:

[0181] Among them, #B represents the number of elements in the set;

[0182] Define the set of all degrees of freedom as , take , and "-" represents the set difference operation;

[0183] For in, the element numbers are as follows, , where #D represents the number of elements in set D.

[0184] Define , represents the transpose, is a column vector of dimension #D;

[0185] Similarly, , is a column vector of dimension #B; similarly, write the Lagrange multipliers as a column vector, , so, is a column vector of dimension #B;

[0186] (4) is the column vector formed by arranging all the elements of in sequence, that is .

[0187] Take , , (including #B zero elements), , that is is the column vector obtained by multiplying the column vector composed of by .

[0188] (6) The elements of correspond to the equations in Equation Set 【7】(see the specific implementation part) in sequence, the elements of correspond to Equation Set 【9】, the elements of correspond to the equations in Equation Set 【8】; after arranging the above equations 【7】, 【8】, 【9】 in the corresponding order, extract the corresponding coefficients of the elements regarding in each equation (in the order of elements) to form the matrix , so, each equation corresponds to a row of the matrix , and each variable corresponds to a column of .

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

[10] , and Equation

[10] is the so-called equation 。

[0190] 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 with extremely small values 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 ), thus, the present invention proposes a method for tracking the eigenvalue closest to zero, which is as follows:

[0191] 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.

[0192]

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

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

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

[0196] ,

[11]

[0197] 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 's product, is the material density, is the determinant of the deformation gradient tensor.

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

[0199]

[12]

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

[0201]

[13]

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

[0203] Substitute the perturbation expansion (Equations

[12] and

[13] ) into the variational equation

[11] , and based on order, we get:

[0204] (1) For order,

[0205] [14 - 1]

[0206] and ; [14 - 2]

[0207] (2) For order,

[0208] [15 - 1]

[0209] [15 - 2]

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

[0211] Using the arbitrariness, write Equation [15 - 1] in the form of an eigenvalue problem (for all repeated indices, sum according to the Einstein summation convention):

[0212] [15 - 3]

[0213] , where, when , when at that time, .

[0214] Equations [15 - 3] and [15 - 2] define a zero eigenvalue problem, which is transformed into an eigenvalue problem of a matrix through node numbering and flattening, and the eigenvector field corresponding to the zero eigenvalue is obtained, 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 the eigenvector corresponding to the eigenvalue with the smallest absolute value is used for replacement.

[0215] The eigenvector After normalization (as an option, maximum value normalization is adopted, that is, the maximum value of each component is taken, and the eigenvector is divided by this maximum value to obtain the normalized eigenvector), the normalized eigenvector field is denoted as .

[0216] (3) For order,

[0217] 【16-1】

[0218] , , 【16-2】

[0219] where, is the determinant of the deformation gradient tensor at the saddle point, which characterizes the volume ratio of the saddle point state to the state before deformation.

[0220] satisfies Equation 【15-3】 and can be decomposed into ;

[0221] Based on the solvability conditions of Equations 【16-1】 and 【16-2】, an ordinary differential equation about is obtained:

[0222]

[0223] where, 【17-1】

[0224] 【17-2】

[0225] 【17-3】

[0226] Step 302. Based on the interpolation of the eigenvalue closest to zero, and approximately calculate

[0227] As can be seen from the process of Step 301, it is necessary to extract the parameters and , and the approximate position of the saddle point is calculated from the equilibrium path in 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 computational cost.

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

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

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

[0231] (3) Use the following interpolation formula:

[0232] ,

[0233] ,

[0234] .

[0235] Thus, the approximate value of the saddle point is obtained.

[0236] Step [4]

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

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

[0239] , [17 - 4]

[0240] , where the function is the Gaussian hypergeometric function. Take the limit , and obtain the theoretical shallow arch jump time

[0241] 。 [17-5]

[0242] Thus, 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 true physical time is established to calculate the approximate transient response of the shallow arch saddle point jump.

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

[0244] Module (1) is responsible for inputting the nodes and degrees of freedom directions of load application, 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 element. This module is completely implemented in the Python language;

[0245] ​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 internal force residual vector of each unit (corresponding to equation [5]). Python combines the Jacobian matrix and internal force residual vector of the unit into a global Jacobian matrix and a global internal force residual vector, and converts the matrix into a sparse matrix form. The matrix equation is solved by the spsolve function of the scipy library function provided by Python, and the eigsh function of the sparse module of the scipy library provided by Python is used to calculate the eigenvalue 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. The C language is first used to calculate the equation [15-3] involved The calculation part of each unit is carried out, and the eigenvector V is obtained by using the eigsh function provided by Python, and the part involving each unit in equations [17-2] and [17-3] is calculated by C language. Finally, the Gaussian hypergeometric function is calculated using the scipy.special.hyp2f1 library function provided by Python; Module (4) is a module proposed by the present invention for directly calculating the transient dynamic process of shallow arch jump, which is mainly used to verify the correctness of the results of the saddle point finite element module proposed by the present invention. It is divided into an initial acceleration calculation step and an implicit dynamic analysis step based on the GN22 algorithm; Module (5) is used for visualization of the results of each analysis module. In order to avoid repeated development, the mature Abaqus python development platform is used to realize the visualization of the results of the perturbation finite element analysis system involved in the present invention.

[0246] Taking the visualization of displacement or deformation in the calculation results as an example, the specific implementation method of module (5) is explained:

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

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

[0249] (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;

[0250] (4) After saving the odb file, the deformation nephogram and other information of the shallow arch can be viewed using the visualization function of abaqus.

[0251] Example 1

[0252] Take Figure 5 the circular shallow arch with hinged ends shown 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 connection line between the end face and the center of the circle forms 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 .

[0253] 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 verification between the shallow arch jump transient response given by the perturbation finite element method and the direct finite element dynamic analysis.

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

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

[0256] In the Abaqus platform, due to the structural symmetry, establish half of the finite element mesh model as shown in Figure 5 using 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:

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

[0258] 1, 0.173474535,0.983822942, 0.00999999978

[0259] 2, 0.17261605, 0.98397547, 0.00999999978

[0260] The unit connection relationship list includes the unit node numbers and the numbers of the 8 included nodes. For example:

[0261] 1, 1006, 1007, 1208, 1207, 1, 2, 203, 202

[0262] 2, 1007, 1008, 1209, 1208, 2, 3, 204, 203

[0263] 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 node positions 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; The 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, such as Figure 9 .

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

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

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

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

[0268] ,

[0269] where is the area occupied by the shallow arch before deformation, that is, the volume area 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:

[0270] , ":" represents the contraction operation of the tensor, 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 degree of freedom defined by the boundary conditions. represents the reference concentrated force (which does not change along the equilibrium path) applied in the direction of the -th degree of freedom of the -th node, that is ( , see Figure 9 ). represents the load factor of the concentrated force (which changes along the equilibrium path and characterizes the change in the magnitude of the concentrated force). represents the set of node numbers and degree-of-freedom pairs corresponding to the constraints in the boundary conditions, determined by the nodes and related degrees of freedom shown in Figure 7 and Figure 8 . represents the set of node numbers and degree-of-freedom pairs where the concentrated force is applied, that is .

[0271] Derived from the above equations:

[0272] ,

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

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

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

[0276]

[0277] , where represents the tensor product operation. When , , when , , and similarly, when , . At this time, take .

[0278] In addition, due to being arbitrary, the corresponding equation of the boundary condition is obtained, that is, when At this time, .

[0279] The arc-length method (Riks) is used 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 clear physical meaning obtained by 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 for each step is fixed at 6.

[0280] 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.

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

[0282] According to the method of tracking the eigenvalue closest to 0 proposed by the present invention, the specific implementation is as follows:

[0283] Figure 11 shows Figure 10 the distribution of the eigenvalues of the global matrix in each increment step (taking 50 load steps, corresponding to the data points in Figure 10 ), including the absolute value of the negative eigenvalue and the positive eigenvalue.

[0284] 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 .

[0285] From Figure 12 it can be known 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 indicates some of the smallest eigenvalues of the gain steps obtained therefrom. This embodiment shows that by setting a larger gain , the eigenvalue closest to zero can be effectively obtained, and the approximate position of the saddle point can be obtained.

[0286] The load factor value at the saddle point of this embodiment is approximately (corresponding to Figure 10 14th incremental step in the process), that is, when the Y-axis direction of the center end face is applied When the saddle point is reached, the shallow arch is in a critical equilibrium state. At this time, the displacement three-dimensional vector of the central node of the central end face is calculated as ,Right now Figure 9 Node 2413 in the original undeformed configuration is as follows Figure 14 As shown in configuration 1 in , the deformed configuration at the saddle point is as follows Figure 14 As shown in configuration 2.

[0287] (4) Saddle point jump parameter extraction based on perturbation finite element

[0288] When the load at the saddle point is , further increased to When, is the additional load magnitude coefficient relative to the saddle point, is a fixed quantity related to the additional load direction ( represents the vector of the additional load), in this case ,Right now ( ,See Figure 9 ) ,and ( ). Obviously, The size of the additional load is determined, here we take .

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

[0290] ,

[0291] Lagrange multiplier perturbation expansion:

[0292] ,

[0293] in, For all The flattened vector, and remember , is the vector flattened by the Lagrange multiplier at the saddle point, is the coefficient vector of the perturbation expansion, Physical real time and perturbation parameters The scaled time parameter obtained by multiplying, is the displacement vector field, is the displacement vector field at the saddle point, For perturbation expansion and time Related coefficients.

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

[0295] 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;

[0296] The magnitude of the eigenvector is arbitrary. Therefore, normalization is required. In this embodiment, 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 Y-direction displacement degree of freedom 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.

[0297] Finally, carry out the perturbation finite element analysis step, and obtain the parameters from Equation [17-1], Equation [17-2], and Equation [17-3] , . To clarify The practicality of the calculation, carry out a direct implicit dynamics finite element analysis, and adopt the GN22 algorithm, that is, the relationships of displacement, velocity, and acceleration satisfy Equation [18-1] and Equation [18-2].

[0298] [18-1]

[0299] [18-2]

[0300] Among them, take , Represents 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.

[0301] Thus, only the unknown quantity needs to be calculated to obtain Displacement and velocity at a moment.

[0302] Due to the additional load applied at the saddle point which results 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. The specific method is as follows:

[0303] (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:

[0304]

[19]

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

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

[0307] (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 thus the displacement field and velocity field at each moment can be obtained.

[0308] Figure 13 respectively show the comparison diagrams of the radial displacement at the center of the shallow arch in the perturbation finite element calculation and the direct implicit dynamics calculation of the present invention ( Figure 9 the displacement increment in the -Y direction of node 2413 relative to the saddle point state in ). Among them, , , is the load at the saddle point ( ). Thus, the load during the jump increases by compared with the saddle point load.

[0309] Figure 14shows the shape of the global deformation obtained from the direct implicit dynamics finite element calculation, corresponding to Figure 13 configurations 2, 3, 4, 5, and 6 marked in

[0310] From Figure 13 the implicit dynamics finite element calculation in it can be seen that at time 0 s, it corresponds to the saddle point; as time increases, a dynamic jump occurs 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 significant. During this time interval, the perturbation finite element method proposed in 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

[0311] when the displacement first reaches the maximum match well. The present invention is applicable to different Figure 15 without the 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 the displacement-time path at

[0312] 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 and the corresponding time when the displacement first reaches the maximum are close, indicating that the perturbation finite element method of the present invention can approximately obtain the dynamic response. It should be noted that a series of results can be obtained by the perturbation finite element calculation without the need to carry out perturbation finite element analysis multiple times for

[0313] each one. For direct dynamic finite element analysis, calculations need to be carried out for each

[0314] Example 2

[0315] 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 perturbation finite element analysis calculations. However, in actual calculations, 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. In this case, a large number of increment steps and a high computational cost are required.

[0316] To further improve the accuracy of the saddle point, the present invention proposes a method of linear interpolation of eigenvalues, which is as follows:

[0317] 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:

[0318] ,

[0319] ,

[0320] .

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

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

[0323] 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

[0324]

[0325]

[0326] Due to 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.

[0327] Embodiment 3

[0328] To further illustrate the effectiveness of the interpolation method based on the eigenvalue closest to zero proposed by the present invention, the same parameters and model as in Embodiment 1 are used, but the arc length of each equilibrium path calculation step is enlarged (the increment step size is larger), making it difficult to accurately obtain the saddle point artificially. 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.

[0329] Figure 18 Fig. 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 taking a larger arc length increment, therefore, compared with Embodiment 2, the eigenvalue is larger, and using any calculation point as an approximate saddle point will inevitably bring a larger error.

[0330] 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 of the . 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 .

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

[0332] In increment step 8, the eigenvalue is , , ;

[0333] In increment step 9, the eigenvalue is , , 。

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

[0335]

[0336]

[0337] and the values in Example 2 and are compared. The values are very close. Compared with the eigenvalue situation in Example 2, the absolute values of the eigenvalues at 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 in the present invention can be robustly applicable to different increment step sizes and has the following advantages: a larger step size can be selected without ensuring that the increment step exactly reaches the saddle point position or is very close to the saddle point (in fact, this is very difficult to achieve unless a very small step size is used, and the reduction of the step size will increase the number of required increment steps and the amount of calculation), greatly reducing the calculation amount of the equilibrium path calculation steps.

[0338] Example 4

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

[0340] In Example 1, the wall thickness is uniformly 1 mm, while in this example, the wall thickness is non-uniformly distributed: it linearly varies 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.

[0341] Comparison Figure 21 and Figure 10 shows that since the wall thickness decreases, 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

[0342] Figure 22 gives 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 carried out on the results of the normalized eigenvectors at increment step 8 (i.e., ) and the normalized eigenvector at increment step 9 (i.e., ) to obtain the eigenvector at the final saddle point:Figure 24 shows 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 is taken (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.

[0343] Subsequently, based on the interpolated 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.

[0344] Finally, it should be noted that, in order to illustrate that the parameter varies 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 in . However, in actual implementation, only the perturbation finite element analysis needs to be carried out 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.

[0345] 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).

[0346] The embodiment of the present invention combines the perturbation method and the finite element discretization method, adopts the Green strain measure applicable to large deformations, 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 terms of 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.

[0347] 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 generally 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 computational complexity of the equilibrium path analysis step (see Examples 2 and 3).

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

[0349] 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 stiffness Jacobian matrix, internal force residual vector and other intensive calculation parts of each unit are processed in C language. The modules are compiled into Python modules to overcome the problem of low computational efficiency of Python language. At the same time, the Python API provided by the Abaqus platform is cleverly utilized to implement the mature Abaqus platform as a calculation result visualization platform for this system. Finally, a method for solving relevant equations and eigenvalues 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.

[0350] 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 include the following steps: S1. Generation of thin-walled arch geometric model and finite element mesh; S2. Tracking of the equilibrium path of the thin-walled arch based on the Riks method and extraction of eigenvalues; S3. Finite element calculation of thin-walled arch saddle point perturbation and parameter interpolation based on two eigenvalues; S4. Approximate calculation of transient response based on the Gaussian hypergeometric function; Among them, S2 includes the following steps: a. Define the interpolation function of the global field variable; b. Calculate the equilibrium path under the external load; c. Set the gain coefficient , eliminate the interference of the minimum non-zero eigenvalue caused by the Lagrange multiplier term, and extract the eigenvalue closest to zero for each arc-length method increment step; Among them, 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 For approximate calculation, perform perturbation finite element analysis for each increment step to obtain the parameter , where the parameter is the eigenvector field; approximately obtain the at the saddle point based on interpolation approximation of the positive and negative eigenvalues closest to zero. 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 change of positive and negative signs are respectively 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, and 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; The specific calculation method in S4 is as follows: Assume initial conditions , establish and relationship: , The function is a Gaussian hypergeometric function; Take the limit , and obtain the jump time of the approximate thin-walled arch structure of the perturbed finite element ； 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 to calculate the approximate transient response of the saddle-point jump of the thin-walled arched structure. Among them, is the scaled time parameter obtained by multiplying the physical true time and the perturbation parameter . is the transient displacement vector field, is the displacement vector field at the saddle point.

2. A finite element approximation calculation method for a thin-walled arched structure as described in claim 1, characterized in that In step S1, according to the thin-walled arch set parameters, 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.

3. A finite element approximate calculation method for a thin-walled arch structure according to claim 2, characterized in that Define the global field variable interpolation function, and the specific method is as follows: (1) Use eight-node hexahedron elements and the Lagrangian interpolation method to calculate the field variables and interpolation functions; (2) Let the three-dimensional coordinate vector before deformation be , and define in the interpolation form as ; is the Lagrangian coordinate of the standard element, are the interpolation functions of the respective nodes inside the element; From the defined and one-to-one correspondence, are all regarded as functions, ; are the numerical values of the field variables at the eight nodes; (3) Uniformly number all the nodes , where represents the number of nodes in the finite element model of the thin-walled arched structure. The global field variables in the final thin-walled arched structure are determined by the following interpolation formula: ， Among them, 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.

4. A finite element approximation calculation method for a thin-walled arched structure according to claim 1, 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.

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

6. A lightweight system, characterized in that, The finite element approximate calculation method described in claim 5 is used, and based on Python / C mixed programming, it 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.

7. A lightweight system according to claim 6, characterized in that, The input file parsing module is responsible for inputting the nodes and degrees of freedom directions of load application, the nodes and degrees of freedom directions of geometric constraints, the finite element node coordinate table and the node list included in the elements after the input file is imported into the abaqus platform for modeling and mesh division, and is implemented in the python language; the arc length method analysis module is responsible for the Riks analysis step, 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 to call.

8. A lightweight system according to claim 7, characterized in that, The saddle point perturbation finite element analysis module first calculates the calculation part of each element involved in the equation in C language, obtains the eigenvector V using the eigsh function in Python, and finally calculates the Gaussian hypergeometric function using the library function in Python; the implicit dynamic analysis module directly calculates the module for the transient dynamic process of the thin-walled arch structure jump, mainly used to verify the correctness of the results of the saddle point finite element module; the result visualization module realizes the result visualization of each analysis module with the help of the mature Abaqus python development platform.

Images