Non-linear large deformation calculation method for gradient smooth stable collocation based on super convergence

Abstract

The invention relates to a nonlinear large deformation calculation method for gradient smooth stable collocation based on super convergence. The method comprises the following steps: describing a physical quantity corresponding to a deformed object under a complete Lagrange framework; according to a stable point collocation method, constructing an integral sub-domain by taking each Lagrange particle as a center, and constructing a discrete control equation of an increment iteration format in each integral sub-domain; a smooth function is introduced, the divergence theorem is applied, a gradient smooth sub-domain boundary integral equation of the balance equation is obtained, and the physical boundary integral equation is kept consistent with the stable collocation method; taking the integral sub-domains as gradient smooth domains, and forming two gradient smooth and stable collocation schemes according to whether the integral sub-domains are compatible or not; integral points are arranged according to boundary integral points and Gaussian rules, and weights are selected; and constructing an incremental iteration discrete matrix in a numerical integration form, solving the matrix by adopting a Newton-Rafson iteration method, and performing post-processing. Compared with the prior art, the method can efficiently and accurately solve and calculate the nonlinear large deformation problem.

Description

Technical Field

[0001] The invention relates to the technical field of nonlinear large deformation calculation, and in particular to a nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation. Background Art

[0002] Large deformation is a common phenomenon in fields such as mechanics, materials science, and engineering. When numerically simulating large deformation behavior, the nonlinear strain-displacement relationship under external loads will cause significant displacement and rotation, which poses a major challenge to numerical calculations. Currently, the finite element method (FEM) is the main numerical tool for solving engineering deformation problems. However, when dealing with large deformations, the finite element model (FE) often faces mesh distortion problems due to the rigid connection of the unit to the geometric object, which requires re-meshing. This process involves mapping between the old mesh and the new mesh, which may introduce large errors.

[0003] In recent years, a series of meshless methods have been gradually developed to solve mesh-related problems. In meshless methods, field variables are approximated by the information of multiple particles, and the value of the field variable on each particle can be approximated by the shape function constructed by neighboring particles. In this way, meshless methods do not need to rely on mesh elements, avoiding numerical defects related to elements. Meshless methods can be roughly divided into two categories: Galerkin methods based on weak forms and collocation methods based on strong forms. Galerkin-type meshless methods have attracted widespread attention for the accuracy, stability and convergence of their integral generation, and have been widely used in numerical simulations of large deformation problems. However, these methods may lead to integral errors, high computational costs, and failure to pass piecewise tests in the case of non-uniform particle distribution when using background meshes, rational shape functions, high-order Gaussian integrals, and when the integral domain does not match the support domain.

[0004] In addition, the errors introduced by the insufficient precision of numerical integration may destroy the conservation of the terms on both sides of the discrete large deformation control equation, thereby hindering the convergence of nonlinear iterative calculations. To solve these problems, Chen et al. proposed the stable consistent node integration (SCNI) method that satisfies the first-order integral constraint and applied it to the numerical simulation of nonlinear large deformation. The stable non-compatible node integration (SNNI) method solves the problem of constructing a consistent grid on structurally incomplete fragments. Duan et al. further improved SCNI to satisfy the second-order integral constraint. Subsequently, Chen et al. proposed the variational consistent integration (VCI) method, which can satisfy the integral constraint of any order. Although these improved node integration schemes improve the integral accuracy of the Galerkin-type meshless method, the additional integral terms added to achieve high-order completeness also significantly increase the complexity and computational cost of the method.

[0005] The collocation-based meshless method does not require the special imposition of Dirichlet boundary conditions or integration, and has attracted wide attention due to its simplicity and implementation efficiency. Research shows that this method is more efficient than the Galerkin-based meshless method and the finite element method (FEM). Among them, the radial basis function (RBF) collocation method represented by the direct collocation method (DCM) with global approximation has exponential convergence characteristics and has been used in wave propagation, inverse problems, fluid problems, and even large deformation analysis. However, the traditional RBF-based DCM faces two major problems: one is the instability caused by the single particle representing local characteristics, and the other is the dense and ill-conditioned discrete matrix.

[0006] To solve the instability of DCM, existing research has proposed the least squares collocation method (LSCM) and the improved weighted least squares collocation method (WLSCM). However, although using more collocation points than source points in these methods improves stability, it reduces computational efficiency and fails to solve the problem of matrix ill-conditioning in numerical calculations. To overcome the problem of ill-conditioned matrices, an effective method is to adopt the RK shape function proposed by Liu et al., which satisfies high-order consistency through reproduction conditions. The reconstructed kernel (RK) approximation function with local approximation has promoted the development of the RKCM method due to its high accuracy and stability. To combine the high accuracy and stability of the Galerkin-based meshless method and the efficiency of the collocation method, Wang et al. proposed the SCM with RK approximation. This method performs integration within a regular subdomain centered on particles, which not only retains the high-order consistency of the RK function but also reduces the matrix condition number, significantly improving the accuracy and stability of the algorithm. High-order accuracy can be achieved through low-order Gaussian numerical integration, ensuring an efficient integration process. The initial version of the SCM with RK approximation can pass the patch test under uniform or non-uniform particle distributions and has been applied to linear elasticity problems and heat conduction problems, playing an important role in promoting the development of the collocation-based meshless method.

[0007] However, there are still significant challenges in solving geometrically nonlinear deformation problems in the strong form collocation method, especially when considering the characteristics of hyperelastic materials. This challenge is mainly reflected in two aspects: one is the complex and time-consuming calculation of high-order RK derivatives; the other is the functional relationship between strain and hyperelastic materials involving the deformation gradient, which introduces numerous additional differential terms and a lengthy Newton-Raphson iterative discretization scheme. To address the above problems, Wang and his team have made a series of progress. For example, the gradient reproducing kernel (GRK) approximation is introduced in the framework of the gradient kernel approximation stabilized collocation method (GRKSCM) to avoid direct calculation of high-order RK derivatives; the stabilized Lagrange interpolation collocation method (LISCM) based on Lagrange interpolation provides a simple polynomial shape function; the gradient smoothing collocation method (GSCM) efficiently calculates high-order RK gradients by combining the divergence theorem with gradient smoothing techniques. Wang et al. also proposed a superconvergent gradient smoothing collocation method (SGSMC), which constructs first-order gradient smoothing by performing meshless interpolation on the standard derivatives of the meshless shape function, and the second-order gradient smoothing is represented as two consecutive first-order gradient smoothing operations on the shape function, thus simplifying the numerical implementation process. However, these methods fail to eliminate the additional high-order derivatives in nonlinear deformation and do not achieve a concise discretization formula like the nonlinear Galerkin method. Generally speaking, the collocation method has certain development in numerical solutions, but still cannot overcome its shortcomings in dealing with complex nonlinear problems. Summary of the Invention

[0008] The purpose of the present invention is to overcome the defects of the above-mentioned existing technologies and provide a nonlinear large deformation calculation method based on superconvergent gradient smoothing stabilized collocation, which can efficiently and accurately solve and calculate nonlinear large deformation problems.

[0009] The purpose of the present invention can be achieved through the following technical solutions: A nonlinear large deformation calculation method based on superconvergent gradient smoothing stabilized collocation, comprising the following steps:

[0010] S1. Initialization processing, describing the physical quantities corresponding to the deformed object in the total Lagrangian framework;

[0011] S2. According to the stabilized collocation method, construct an integration subdomain centered on each Lagrangian particle, and construct a discrete control equation in incremental iterative form within each integration subdomain;

[0012] S3. According to the discrete control equation, introduce a smoothing function and apply the divergence theorem to obtain a gradient-smoothed subdomain boundary integral equation of the equilibrium equation, while the physical boundary integral equation remains consistent with the stabilized collocation method;

[0013] S4. Take the integral subdomain as the gradient smoothing domain, and form two gradient smoothing stable collocation schemes according to its compatibility, including the incompatible smoothing domain scheme and the compatible smoothing domain scheme;

[0014] S5. According to the gradient smoothing domain, arrange the integration points and select the weights according to the number of boundary integral points and the Gaussian rule;

[0015] S6. Construct an incremental iterative discrete matrix in the form of numerical integration according to the integration points and weights;

[0016] S7. For the discrete matrix, use the Newton-Raphson iteration method to solve the matrix equation;

[0017] S8. Post-process the solution results of step S7, complete the loading process and large deformation calculation, and obtain the displacements and stresses of the particle points.

[0018] Furthermore, the specific process of step S1 is as follows:

[0019] In a fixed Cartesian coordinate system, let the positions of the particles in the initial configuration and the current configuration be represented as X and x respectively. Assume that the initial occupied region of the deformed body is Ω, and its total boundary is Γ, where Γ = Γ h ∪Γ g , indicating that the total boundary is the union of the Neumann boundary Γ h and the Dirichlet boundary Γ g . Then the large deformation control equation in the total Lagrangian description is:

[0020] (S ij (X)F ki (X)) ,j +b k (X) = 0, in Ω

[0021] S ij (X)F ki (X)N j (X) = h k (X), on Γ h

[0022] u k (X) = g k (X), on Γ g

[0023] where S ij (X) represents the second Piola-Kirchhoff stress, N k (X) is the unit outer normal vector, F ki is the deformation gradient tensor, b k (X) is the body force term, h k (X) is the Neumann boundary external load, gk (X) is the given Dirichlet boundary condition.

[0024] Furthermore, the second Piola-Kirchhoff stress is specifically:

[0025] S ij (X) = D ijkl E kl (X)

[0026]

[0027] where D ijkl is the elastic constant tensor, and E kl (X) represents the Green-Lagrange strain;

[0028] The deformation gradient tensor is specifically:

[0029]

[0030] The incremental large deformation iterative control equation is as follows:

[0031]

[0032] ΔS ij (X) = D ijkl ΔE kl (X)

[0033] ΔS ij,j (X) = D ijkl ΔE kl,j (X)

[0034] Different from the strain formula of small deformation, the Green-Lagrange strain has a non-linear term. Therefore, the solution of the non-linear control equation requires the use of incremental linearization and the Newton-Raphson iteration method, that is, the load step and the iteration step are respectively represented as n and v, and Δ(·) is the increment of the field variable.

[0035] Furthermore, the specific process of step S2 is:

[0036] The particle points arranged on the Neumann boundary Γ h and the Dirichlet boundary Γ g in the problem domain Ω are respectively represented as and where D p , D q and D r represent the number of Lagrangian particle points at the response position. Therefore, the total number of particle points is D = D p + D q + Dr , according to the SCM discrete and sub-domain integration principle, integrate over the sub-domain centered at each particle point to establish the discrete control equation:

[0037]

[0038] where, Ω I , and respectively represent the integration sub-domains on the Neumann boundary Γ h and the Dirichlet boundary in the domain, that is

[0039] Furthermore, in step S3, a piecewise constant function is used as the smoothing function to smooth the gradient of the in-domain integral balance equation, and then integration by parts is performed on the smoothed in-domain integral balance equation and the divergence theorem is applied to obtain the sub-domain boundary integral equation with smoothed gradient.

[0040] Furthermore, the smoothing function in step S3 is specifically:

[0041]

[0042] where, φ I is the smoothing function, and A I represents the area or volume of the sub-domain Ω I .

[0043] Furthermore, the sub-domain boundary integral equation with smoothed gradient in step S3 is specifically:

[0044]

[0045] where, Γ I represents the boundary of the integration sub-domain Ω I , represents the unit outer normal vector on the boundary of the integration sub-domain, represents the m-th segment of the boundary Γ I of the integration sub-domain, and MI is the number of boundary segments of the integration sub-domain.

[0046] Furthermore, the process of arranging integration points in step S5 includes: for and the integration points on and their number and positions follow the traditional SCM. For the integration sub-domain Ω I in the problem domain, in order to meet the requirements of the sub-domain boundary integral of the in-domain balance equation, the integration points are arranged on the boundary of the in-domain integration sub-domain.

[0047] Further, the specific process of step S6 is as follows:

[0048] Using the RK approximation, approximate the displacement field increment Δu k (X) to obtain the approximate displacement field increment as follows: As shown below:

[0049]

[0050] Among them,

[0051] Δa = [Δa 1 , L, Δa D T , Δa I = [Δa 1 I , Δa 2 I , Δa 3 I T

[0052] Φ T = [Φ 1 , L, Φ D T , Φ I = Ψ I I = Ψ I (X)I

[0053] Ψ I (X) represents the RK shape function constructed by particle X within the support domain centered at particle X, matrix I represents the identity matrix, Φ I is the RK shape function matrix generated by all particle X T , and Δa is the increment of the undetermined coefficient;

[0054] Substitute the RK approximation into the equilibrium equation and the boundary condition equation to obtain the following discrete equations: Among them, S

[0055]

[0056] where S a represents the second Piola-Kirchhoff stress matrix, S b is the first derivative matrix of the second Piola-Kirchhoff stress, D a , D b , D c and D d are constitutive matrices in different forms, B a , B b , B c , B d and B​​​e represent different forms of strain-displacement matrices, L a , L b , and L c respectively represent the strain connection matrices formed by the deformation gradient, S c , S d and S c similarly represent different forms of stress matrices, and is the unit outer normal vector of the boundary of the domain integral subdomain, N 1 and N 2 represent different forms of normal vectors on the Neumann boundary, b, h, g, and u respectively represent the body force vector, surface traction vector, boundary constraint displacement vector, and known displacement vector, Δa represents the particle coefficient increment vector, Δa i I represents the coefficient increment of particle I in the i direction. The specific expressions of these matrices are as follows:

[0057]

[0058]

[0059] Under plane stress conditions, D E0 = D E , v 0 = v. For plane strain conditions, D E is the elastic modulus, v is the Poisson's ratio;

[0060] Through numerical integration, the following discrete matrix form is obtained:

[0061]

[0062]

[0063] where, represents the number of boundary integration points, and are the boundary integration points on the integral subdomains centered on particles p I , q I and r I respectively;

[0064] Then it is assembled into the overall matrix form as follows:

[0065]

[0066] where,

[0067]

[0068]

[0069] Among them, is the stiffness matrix of the particle points in the problem domain, is the stiffness matrix of the particle points on the Neumann boundary, is the stiffness matrix of the particle points on the Dirichlet boundary, is the residual matrix of the particle points in the domain, is the residual matrix of the particle points on the Neumann boundary, is the residual matrix of the particle points on the Dirichlet boundary.

[0070] Furthermore, the specific process of step S7 is as follows:

[0071] First, the external load is divided into n incremental steps, and several iterative steps v are required within each incremental step;

[0072] Secondly, preset initial values for the variables, n = 0, v = 0, the undetermined coefficient matrix a = 0, and the particle point displacement matrix u = 0;

[0073] After obtaining the initial values, enter the iterative process. In each iterative step of each load step, according to the undetermined coefficients obtained in the previous step Use the Green-Lagrange strain formula and the constitutive relation to obtain the updated stiffness matrix within the domain Neumann boundary stiffness matrix Dirichlet boundary stiffness matrix Residual matrix within the domain Neumann boundary residual matrix and Dirichlet boundary residual matrix and assemble them into an overall matrix iterative format to calculate the increment of the undetermined coefficients corresponding to the current iterative step

[0074] After that, accumulate and update the undetermined coefficients Update to obtain the current displacement field increment And accumulate and update the displacement field After each iterative step is calculated, perform a convergence judgment. If the norm of the current displacement increment satisfies e is the set error threshold, it means convergence and enter the next load step; otherwise, continue the above iterative steps.

[0075] Compared with the prior art, the present invention has the following advantages:

[0076] In view of the problem of large non - linear deformation, the present invention designs a super - convergent gradient - smoothed stabilized collocation method. In this method, the problem domain is discretized by a series of Lagrangian particles, and the physical fields of the entire deformed body are established in the fully Lagrangian framework. Then, by using the divergence theorem, the traditional sub - domain integration in the stabilized collocation method (SCM) is converted into sub - domain boundary integration, thereby realizing gradient smoothing and simplifying the governing equations. In addition, reducing the derivative order of the equations enables the RK approximation to use low - order basis functions (such as constant basis functions), thus improving the computational efficiency. And in the process of gradient - smoothed numerical integration, compatible and incompatible smoothing domains / integration sub - domains are constructed. Finally, the discrete equations are solved for large - deformation problems through the Newton - Raphson iteration technique, which can efficiently and accurately solve and calculate large - non - linear deformation problems.

[0077] The present invention introduces the divergence theorem to perform gradient smoothing on the strong - form equilibrium equation in the integral form of SCM, converting sub - domain integration into sub - domain boundary integration, avoiding the complex derivative calculations and long non - linear terms in the traditional collocation method. The equations are concise and highly efficient.

[0078] The present invention takes the integration sub - domain as the gradient - smoothing domain, with a simple background grid structure, and can flexibly construct compatible or incompatible sub - domains according to actual needs.

[0079] In the present invention, the integration sub - domain is only related to the particle positions. During the deformation process, the integration sub - domain does not deform, which can avoid the process of re - meshing the integration grid and improve the efficiency of large - deformation analysis.

[0080] The present invention adopts the RK approximation function with low - order basis functions, which can obtain the numerical solutions of high - order strong - form governing equations and further improve the computational efficiency.

[0081] The present invention can achieve accurate integration by using low - order Gaussian integration, improving the accuracy and stability of the collocation method.

[0082] In the present invention, integration sub - domains are constructed for each particle point, which is beneficial to reducing the condition number of the discrete matrix and further improving the stability of the algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

[0083] Figure 1 is a schematic diagram of the method flow of the present invention;

[0084] Figure 2 is a schematic diagram of the application process of the embodiment;

[0085] Figure 3 is a schematic diagram of the SCM discretization scheme;

[0086] Figure 4 is a schematic diagram of the integration sub - domain in the present invention;

[0087] Figure 5Schematic diagram of geometric information and coordinate transformation between an arbitrary integral sub-domain and a reference regular integral sub-domain;

[0088] Figure 6 is a schematic diagram of the discrete scheme of the incompatible gradient smooth stabilized collocation method;

[0089] Figure 7 is a schematic diagram of the discrete scheme of the compatible gradient smooth stabilized collocation method;

[0090] Figure 8a and 8b are the geometric models of the square plate model with a central hole and its quarter-symmetric model in the embodiment;

[0091] Figure 9a and 9b are the sub-domains and particle point distributions of GSSCM-I and GSSCM-II for the square plate with a central hole in the embodiment;

[0092] Figures 10a to 10e are the schematic diagrams of the deformed plate stresses obtained by using FEM, GSSCM-I and GSSCM-II in the embodiment;

[0093] Figure 11a and 11b are the comparisons of the displacement results of DCM, SCM, SGSMC, FEM, GSSCM-I and GSSCM-II in the embodiment;

[0094] Figures 12a to 12d are the comparisons of the errors and convergence of DCM, SCM, SGSMC, GSSCM-I and GSSCM-II in the embodiment;

[0095] Figure 13a and 13b are the comparisons of the condition numbers of DCM, SCM, SGSMC, GSSCM-I and GSSCM-II in the embodiment;

[0096] Figure 14 are the comparisons of the CPU times of DCM, SCM, SGSMC, GSSCM-I and GSSCM-II at the same relative error levels (relative error is about 2.03% when p = 1 and about 1.41% when p = 2) in the embodiment. Specific implementation manners

[0097] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0098] Embodiment

[0099] As Figure 1 shown, a non-linear large deformation calculation method based on superconvergent gradient smooth stabilized collocation includes the following steps:

[0100] S1. Initialization processing, describing the physical quantities corresponding to the deformed object in the total Lagrangian framework;

[0101] S2. According to the stabilized collocation method, construct an integral subdomain centered on each Lagrangian particle, and construct a discrete control equation in incremental iterative form within each integral subdomain;

[0102] S3. According to the discrete control equation, introduce a smoothing function and apply the divergence theorem to obtain a gradient-smoothing subdomain boundary integral equation of the equilibrium equation, while the physical boundary integral equation remains consistent with the stabilized collocation method;

[0103] S4. Take the integral subdomain as the gradient-smoothing domain, and form two gradient-smoothing stabilized collocation schemes according to its compatibility, including an incompatible smoothing domain scheme and a compatible smoothing domain scheme;

[0104] S5. According to the gradient-smoothing domain, arrange the integration points and select the weights according to the number of boundary integration points and the Gaussian rule;

[0105] S6. According to the integration points and weights, construct an incremental iterative discrete matrix in numerical integration form;

[0106] S7. For the discrete matrix, use the Newton-Raphson iteration method to solve the matrix equation;

[0107] S8. Post-process the solution result of step S7, complete the loading process and large deformation calculation, and obtain the displacements and stresses of the particle points.

[0108] This embodiment applies the above scheme, as Figure 2 shown, the main contents are:

[0109] Step 1. Construction of the initial configuration and particle points

[0110] In a fixed Cartesian coordinate system, let the positions of the particles in the initial configuration and the current configuration be represented as X and x respectively. Assume that the object initially occupies the region Ω, and its total boundary is Γ, where Γ = Γ h ∪Γ g , indicating that the total boundary is the union of the Neumann boundary Γ h and the Dirichlet boundary Γ g . Then the large deformation control equation described by the total Lagrangian is:

[0111] (S ij (X)F ki (X)) ,j +b k (X) = 0, in Ω (1)

[0112] S ij (X)F ki (X)Nj ψ(X) = h k ψ(X), on Γ h (2)

[0113] u k φ(X) = g k φ(X), on Γ g (3)

[0114] where S ij (X) represents the second Piola - Kirchhoff stress, N k (X) is the unit outer normal vector, F ki is the deformation gradient tensor, b k (X) is the body force term, h k (X) is the Neumann boundary external load, g k (X) is the given Dirichlet boundary condition. The linear constitutive relation can be expressed as follows:

[0115] S ij (X) = D ijkl E kl (X) (4)

[0116] where D ijkl is the elastic constant tensor, E kl (X) represents the Green - Lagrange strain, and its definition is:

[0117]

[0118] The expression of the deformation gradient tensor is as follows:

[0119]

[0120] Different from the strain formula for small deformations, due to the existence of non - linear terms in the Green - Lagrange strain in Equation (5), the solution of the non - linear control equation requires the use of incremental linearization and the Newton - Raphson iteration method. That is, the load step and the iteration step are represented as n and v respectively. The increment of the field variable is represented by Δ(·). Then the incremental large - deformation iteration control equation is as follows:

[0121]

[0122] ΔS ij (X) = D ijkl ΔE kl (X) (11)

[0123] ΔS ij,j (X) = D ijkl ΔE kl,j (X) (12)

[0124] Step 2. Stabilized collocation method

[0125] In the problem domain Ω, the Neumann boundary Γ h and the Dirichlet boundary Γ g The particle points arranged on are respectively denoted as and where D p , D q and D r represent the number of Lagrangian particle points at the response positions. Therefore, the total number of particle points is D = D p + D q + D r . As Figure 3 shown, according to the SCM discretization and subdomain integration principle, the discrete control equations established in Eqs. (8)–(10) need to be integrated over the subdomains centered at each particle point, and we can obtain:

[0126]

[0127] In the formula, Ω I , and respectively represent the integration subdomains on the Neumann boundary Γ h and the Dirichlet boundary in the domain, that is However, it is not difficult to see from Eqs. (8) and (13) that due to the use of the incremental iteration format and the existence of derivative terms, the equilibrium equations in the large deformation domain constructed by the collocation method are very complex, especially when considering all stress and strain components and the nonlinear terms in the strain matrix. This poses great difficulties for the construction of the final stiffness matrix and load matrix. Therefore, it is very necessary to design this scheme to reduce the order of the derivative of the equation through gradient smoothing to obtain a more concise form.

[0128] Step 3. Construction of the control equations for the smoothed domain and the gradient smoothed stabilized collocation method

[0129] According to the meshless gradient smoothing method, the gradient of the large deformation equilibrium equation of the particles in the subdomain Ω I in the problem domain Ω can be smoothed in the following form through the smoothing function φ I :

[0130]

[0131] This scheme uses a piecewise constant function as the smoothing function, and its definition is as follows:

[0132]

[0133] In the formula, AI Denote the sub - domain Ω I The area or volume. By performing integration by parts on formula (16) and applying the divergence theorem, the following result is obtained:

[0134]

[0135] where Γ I Denote the integration sub - domain Ω I The boundary of, where Denote the unit outer normal vector on the boundary of the integration sub - domain. If the number of boundary segments of the integration sub - domain is MI, then:

[0136]

[0137] where Denote the m - th segment of the boundary Γ I Of the integration sub - domain. The schematic of the integration sub - domain within the problem domain is as shown in Figure 4 As shown. Substituting formula (17) into formula (18), the smooth gradient expression can be obtained:

[0138]

[0139] Combining like terms for the above formula gives:

[0140]

[0141] Retain the SCM boundary integral equations shown in formulas (14) and (15):

[0142]

[0143] To solve the above integral control equation, boundary numerical integration must be performed on the sub - domains and their boundaries within the problem domain. By selecting appropriate Gaussian integration points, accurate boundary numerical integration with p - order accuracy of the RK shape function is achieved. It should be emphasized that Table 1 gives the number of integration points required on the boundary of each regular integration sub - domain. In the one - dimensional case, the integration boundary refers to the endpoints of the one - dimensional line segment; in the two - dimensional case, it corresponds to each side of the polygonal integration sub - domain; and in the three - dimensional case, the integration boundary represents the surface of each integration sub - domain.

[0144] Table 1 Relationship between the number of boundary integration points required for regular integration sub - domains within each domain and the order of the RK basis function

[0145]

[0146] In practice, the shape of the integration sub - domain is arbitrary and irregular, so as to Figure 5Take the two-dimensional quadrilateral integration subdomain in it as an example. The boundary (line segment) of the irregular integration subdomain should be transformed into the integration on the boundary (line segment) of the regular reference subdomain through coordinate transformation:

[0147]

[0148] Among them, ξ and η represent the coordinates of the integration nodes in the regular domain, and the Jacobian determinant is expressed as follows:

[0149]

[0150] The coordinate transformation of linear interpolation is defined as:

[0151]

[0152] Among them, (X m , Y m ) and (X m+1 , Y m+1 ) are the two endpoints of the boundary line segment of the integration subdomain. Similarly, for any three-dimensional hexahedral integration subdomain, its six boundary surfaces must be transformed into the six square boundary planes of the standard cubic integration domain as follows:

[0153]

[0154] Through the sufficient Gauss integration points shown in Table 1 and the exact transformation of the integration subdomain, GSSCM satisfies the following consistency conditions in the integral form in any subdomain within the problem domain:

[0155]

[0156] The consistency conditions for the first and second order gradients of the RK shape function can be expressed as:

[0157]

[0158] To implement the Gauss numerical integration operation, it is necessary to arrange the integration points at appropriate positions on the subdomains corresponding to each particle point. For and the integration points on and their quantities and positions follow the traditional SCM. For the integration subdomain Ω I within the problem domain, in order to meet the requirements of the subdomain boundary integration of the equilibrium equation within the domain, some integration points is arranged on the boundary of the domain integral sub-domain, so the integral sub-domain here is also a gradient smooth domain. The number of boundary integrals is related to the order of the RK basis function, as shown in Table 1. The position requirements of the integration points follow the traditional Gaussian integration requirements. The incompatible smooth domain / integral sub-domain mode of GSSCM-I is as Figure 6a and 6b shown, while the compatible smooth domain / integral sub-domain of GSSCM-II is as Figure 7a and 7b shown. For the sake of illustration, discrete particles, integration nodes, smooth domains, and the compact support domains used to construct the RK approximation function are also shown in the figure.

[0159] In GSSCM-I, the generation of the integral smooth domain is flexible and convenient. Since the smooth domain does not need to strictly conform to the geometry, this method is applicable to cases with complex geometries. For two-dimensional problems, the smooth domain is usually a triangle, square, or other polygon within the domain. However, in the incompatible smooth domain scheme, the integral sub-domains of the particles near the boundary will inevitably partially extend beyond the problem domain. However, the extension of the integral sub-domains of the particles near the problem boundary beyond the problem domain can improve the accuracy of the solution.

[0160] As Figure 7a and 7b shown, the compatible smooth domains in GSSCM-II can be regular or irregular, depending on the requirements of the problem domain and the solution. For two-dimensional compatible smooth domains, continuous triangular and quadrilateral smooth domains can be generated using Delaunay triangulation and Voronoi techniques. The smooth domains are connected and non-overlapping, and each contour forms the integration path of the smooth domain. In addition, the integral sub-domains on the Neumann boundary and the Dirichlet boundary also need to maintain compatibility.

[0161] Step 4. Discrete matrices and iterative formats of GSSCMs

[0162] The RK approximation is an important type of locally supported shape function and is widely used in the meshless method. It has many advantages such as high-order continuity and satisfying the consistency condition. The specific forms of the RK approximation and its derivatives can be referred to the existing technology. Using the RK approximation, the displacement field increment Δu k (X) is approximated to obtain the approximate displacement field increment as follows:

[0163]

[0164] where

[0165] Δa = [Δa 1 , L, Δa D T , Δa I ​= [Δa 1 I , Δa 2 I , Δa 3 I T (34)

[0166] Φ T = [Φ 1 , L, Φ D T , Φ I = Ψ I I = Ψ I (X)I (35)

[0167] Here, Ψ I (X) represents the RK shape function constructed by particle X within the support domain centered on particle X. The matrix I represents the identity matrix, and Φ I is the RK shape function matrix generated by all particle X T . Δa is the increment of the undetermined coefficient. I Substituting the RK approximation into the equilibrium equation and the boundary condition equations (21)-(23), the following discrete equations can be obtained:

[0168]

[0169]

[0170]

[0171] Among them, S a represents the second Piola-Kirchhoff stress matrix, and S b is the first derivative matrix of the second Piola-Kirchhoff stress. D a , D b , D c and D d are constitutive matrices in different forms respectively. B a , B b , B c , B d and B e represent strain-displacement matrices in different forms. L a , L b , and L c represent strain connection matrices formed by the deformation gradient respectively. S c , S d and S c also represent stress matrices in different forms respectively. and are the unit outer normal vectors of the boundary of the domain integration subdomain. N 1 and N​​2 Denote the normal vectors of different forms on the Neumann boundary. b, h, g, and u denote the body force vector, surface traction vector, boundary constraint displacement vector, and known displacement vector, respectively. Δa denotes the particle coefficient increment vector, and Δa i I denotes the coefficient increment of particle I in the i direction. The specific expressions of these matrices are as follows:

[0172]

[0173]

[0174] Under plane stress conditions, D E0 = D E , v 0 = v. For plane strain conditions, D E is the elastic modulus and v is the Poisson's ratio.

[0175] Through numerical integration, equations (36)-(38) are transformed into the following discrete matrix forms:

[0176]

[0177] where, denotes the number of boundary integration points, and are the boundary integration points on the integration subdomains centered at particles p I , q I and r I respectively.

[0178] Equations (53)-(55) are further assembled into the overall matrix form as follows:

[0179]

[0180] where,

[0181]

[0182]

[0183] Equations (57) and (60) show that by converting the subdomain integral into a boundary integral, GSSCMs eliminate the derivative terms of stress and greatly simplify the equilibrium equation.

[0184] Step 5. Superconvergence of GSSCMs

[0185] The local error of the p-th order GSSCMs is:

[0186]

[0187] In the formula,

[0188]

[0189] particle is located within the integral sub-domain, and the particle is located on the boundary of the integral sub-domain. is the operator for smoothing the first-order derivative of :

[0190]

[0191] It can be seen from the above that for GSSCMs using odd-order basis functions in the RK approximation, a (p + 1)-th order superconvergence rate can be obtained.

[0192] Step 6. Newton-Raphson iterative solution process of GSSCMs

[0193] When dealing with large deformation nonlinear problems, the incremental method is needed to linearize the field variable function and the Newton-Raphson method is applied to construct the iterative equation format. First, the external load is divided into n incremental steps, and several iterative steps v are required within each incremental step. Second, preset initial values for some variables, such as n = 0, v = 0, the undetermined coefficient matrix a = 0, and the particle point displacement matrix u = 0. After obtaining the initial values, the iterative process can be entered. Third, in the iterative steps of each load step, according to the undetermined coefficients obtained in the previous step Use the Green-Lagrange strain formula (5) and the constitutive relation formula (4) equations (57)-(62) to obtain the updated in-domain stiffness matrix Neumann boundary stiffness matrix Dirichlet boundary stiffness matrix In-domain residual matrix Neumann boundary residual matrix and Dirichlet boundary residual matrix Third, after assembling into the overall matrix iterative format shown in equation (56), the incremental undetermined coefficients corresponding to the current iterative step can be calculated Third, the undetermined coefficients to be accumulated and updated Third, obtain the current displacement field increment according to equation (33) and accumulate and update the displacement field Third, after each iterative step is calculated, convergence judgment is performed. If the norm of the current displacement increment satisfies (where \(e\) represents the given error threshold), it indicates convergence and the next load step can be entered. If not satisfied, the above iterative steps continue. Finally, when all load steps are calculated, the displacement, stress, and other information of the particle points are output for post-processing.

[0194] To verify the effectiveness of this scheme, in this embodiment, the deformation problem of a flat plate with a central hole is calculated. Figure 8a 、 8b shows a plane strain problem of a square plate with a central hole under a unidirectional tensile load \(T = 2.1\times10\) 10 N / m. Using symmetry, only the upper right quarter of the plate is simulated and modeled. The material properties are Young's modulus \(D\) E \(= 210\times10\) 9 Pa and Poisson's ratio \(v\% = 0.22\), and the aperture radius is set to \(a = 1m\). Since there is no analytical solution for reference, the finite element results with a fine mesh are used for comparison. The number of loading steps is set to 15.

[0195] Figure 9a 、 9b shows the discretization and subdomain configuration of GSSCMs due to irregular geometric features. In this embodiment, the number of particles is \(D = 269,321,436\). According to Table 1, the number of integration points on each boundary of the subdomain is 1 (for \(p = 1\)) and 2 (for \(p = 2\)), respectively. Figures 10a to 10e shows the numerical stress solutions of GSSCM-I and GSSCM-II and compares them with the stress solutions of FEM. The results show that even in the case of a relatively simple mesh and a small number of particles, GSSCM-I and GSSCM-II can obtain a reasonable stress field in this stress concentration problem. Figure 11a 、 11b shows the numerical displacement solutions of the plate problem and the corresponding errors under the same discretization conditions. When \(p = 1\), DCM and SCM cannot converge because they do not meet the consistency conditions required for the second derivative of the RK function. When \(p = 2\), the accuracy of GSSCMs is similar to that of SCM. Generally speaking, under the two basis function conditions, the accuracy of GSSCM-I and GSSCM-II is significantly higher than that of DCM and SGSCM.

[0196] Figures 12a to 12d shows the convergence and error, indicating that the accuracy of GSSCMs is higher than that of DCM and SGSCM, and shows a \((p + 1)\)-th order superconvergence rate at p＝1 , which makes GSSCMs have the same superconvergence as SGSCM. Figure 13a 、 13b shows the condition number. The small condition number of GSSCMs proves its stability, while Figure 14shows the comparison of CPU time under similar error levels. When p = 1, the relative error is approximately 2.03%, and when p = 2, the relative error is approximately 1.41%. These data once again prove that GSSCMs have the highest efficiency and better stability than DCM and SGSCM when obtaining the same accuracy level. And SGSCM is more stable and efficient than DCM, and due to the use of compatible subdomains, the overall performance of GSSCM-II is better than that of GSSCM-I.

[0197] Through the above theoretical analysis and practical examples, it is fully shown that this scheme has high precision, high efficiency and good stability, and can be widely applied to engineering applications with large deformations.

Claims

1. A nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation, characterized in that: The following steps are involved: S1, initialization processing, describing the physical quantities corresponding to the deformed object in the full Lagrangian framework; S2. According to the stable collocation method, an integral subdomain is constructed with each Lagrangian particle as the center, and a discrete control equation in an incremental iterative format is constructed in each integral subdomain; S3. According to the discrete control equation, the smooth function is introduced and the divergence theorem is applied to obtain the gradient smooth subdomain boundary integral equation of the equilibrium equation. The physical boundary integral equation is consistent with the stable collocation method. S4, taking the integral subdomain as a gradient smooth domain, and forming two gradient smooth stable point matching schemes according to whether they are compatible, including an incompatible smooth domain scheme and a compatible smooth domain scheme; S5. According to the gradient smoothing domain, the integration points are arranged according to the number of boundary integration points and the Gaussian rule, and the weights are selected; S6. Constructing an incremental iterative discrete matrix in the form of numerical integration according to the integration points and weights; S7. For discrete matrices, the Newton-Raphson iteration method is used to solve the matrix equation; S8. Post-process the solution of step S7 to complete the loading process and large deformation calculation, and obtain the displacement and stress of the particle point.

2. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 1 is characterized in that: The specific process of step S1 is as follows: In a fixed Cartesian coordinate system, let the positions of the particle in the initial and current configurations be denoted by X and x respectively. Suppose the deformed object initially occupies an area Ω, and its total boundary is Γ, where Γ = Γ h ∪Γ g , indicating that the total boundary is the Neumann boundary Γ h and the Dirichlet boundary Γ g The large deformation governing equations described by the complete Lagrangian are: (S ij (X)F ki (X)) ,j +b k (X)=0,inΩ S ij (X)F ki (X)N j (X)=h k (X),onΓ h u k (X)=g k (X),onΓ g Among them, S ij (X) represents the second type of Piola-Kirchhoff stress, N k (X) is the unit external normal vector, F ki is the deformation gradient tensor, b k (X) is the physical strength item, h k (X) is the Neumann boundary load, g k (X) is the given Dirichlet boundary condition.

3. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 2 is characterized in that: The second type of Piola-Kirchhoff stress is specifically: S ij (X)=D ijkl E kl (X) Among them, D ijkl is the elastic constant tensor, E kl (X) indicates Green-Lagrange strain; The deformation gradient tensor is specifically: The governing equations for the incremental large deformation iteration are as follows: in, ΔS ij (X)=D ijkl ΔE kl (X) ΔS ij,j (X)=D ijkl ΔE kl,j (X) Different from the strain formula for small deformation, the Green-Lagrange strain E kl There are nonlinear terms in , so the solution of the nonlinear control equations requires the use of incremental linearization and Newton-Raphson iteration, that is, the load step and iteration step are represented by n and v, respectively, and Δ(·) is the increment of the field variable.

4. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 3 is characterized in that: The specific process of step S2 is: In the problem domain Ω, Neumann boundary Γ h and the Dirichlet boundary Γ g The particle points arranged on the and Where D p , D q and D r represents the number of Lagrangian particle points at the response position, so the total number of particle points is D = D p +D q +D r According to the principle of discretization and subdomain integration of the stable collocation method (SCM), the discrete control equation is established by integrating on the subdomain centered on each particle point: Among them, Ω I , and They represent the Neumann boundary Γ in the domain respectively. h and the integral subdomain on the Dirichlet boundary, that is 5. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 4 is characterized in that: The step S3 uses a piecewise constant function as a smoothing function, performs smoothing on the gradient of the integral balance equation in the domain, and then performs partial integration on the smoothed integral balance equation in the domain and applies the divergence theorem to obtain a subdomain boundary integral equation with a smooth gradient.

6. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 5, characterized in that: The smoothing function in step S3 is specifically: Among them, φ I is a smooth function, A I Represents subdomain Ω I area or volume.

7. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 6 is characterized in that: The gradient smooth subdomain boundary integral equation in step S3 is specifically: Among them, Γ I Denotes the integral subdomain Ω I The borders of represents the unit external normal vector on the boundary of the integration subdomain, represents the boundary of the integral subdomain Γ I MI is the number of boundary segments of the integral subdomain.

8. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 7, characterized in that: The process of arranging the integral points in step S5 includes: and Points on and Their number The position follows the traditional SCM, for the integral subdomain Ω in the problem domain I In order to meet the boundary integral requirements of the subdomain of the equilibrium equation in the domain, the integration point Arranged on the boundary of the integration subdomain within the domain.

9. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 8, characterized in that: The specific process of step S6 is as follows: Using the RK approximation, the displacement field increment Δu k (X) Approximate operation to obtain approximate displacement field increment As shown below: in, Δa=[Δa 1 ,L,Δa D ] T ,Δa I =[Δa1 I ,Δa2 I ,Δa3 I ] T F T =[Φ1,L,Φ D ] T ,F I =Ψ I I=Ψ I (X)I Ψ I (X) represents particle X I The RK shape function constructed by particle X in the support domain centered at , the matrix I represents the identity matrix, Φ T For all particles X I The generated RK shape function matrix, Δa is the increment of the coefficient to be determined; Substituting the RK approximation into the equilibrium equation and boundary condition equation, we obtain the following discrete equation: Among them, S a represents the second-type Piola-Kirchhoff stress matrix, S b is the first-order derivative matrix of the second-kind Piola-Kirchhoff stress, D a ,D b ,D c and D d They are constitutive matrices of different forms, B a ,B b ,B c ,B d and B e Represents different forms of strain-displacement matrices, L a ,L b , and L c They represent the strain connection matrix formed by the deformation gradient, S c ,S d and S c They also represent different forms of stress matrices, and is the boundary unit external normal vector of the integral subdomain in the domain, N1 and N2 represent the normal vectors of different forms on the Neumann boundary, b, h, g and u represent the body force vector, surface traction vector, boundary constraint displacement vector and known displacement vector respectively, Δa represents the particle coefficient increment vector, Δa i I represents the coefficient increment of particle I in the i direction. The specific expressions of these matrices are as follows: Under plane stress conditions, D E0 =D E ,v0=v. Corresponding to plane strain condition, D E is the elastic modulus, v is the Poisson’s ratio; Through numerical integration, the following discrete matrix form is obtained: Among them, D ￠ represents the number of boundary integration points, and The particles p I ,q I and r I The boundary integration point on the integral subdomain centered at ; Then assemble into the overall matrix form as follows: in, in, is the stiffness matrix of the particle point in the problem domain, is the stiffness matrix of the particle point on the Neumann boundary, is the stiffness matrix of the particle point on the Dirichlet boundary, is the residual matrix at the particle point in the domain, is the residual matrix of the particle points on the Neumann boundary, is the residual matrix of the particle points on the Dirichlet boundary.

10. The nonlinear large deformation calculation method based on superconvergence gradient smooth stable collocation according to claim 9, characterized in that: The specific process of step S7 is as follows: First, the external load is divided into n increments, and a number of iterations v are required in each increment; Secondly, the variables are preset with initial values, n = 0, v = 0, the coefficient matrix to be determined a = 0, and the particle point displacement matrix u = 0; After obtaining the initial value, the iteration process begins. In each iteration step of the load step, the coefficients to be determined are calculated based on the coefficients obtained in the previous step. The updated domain stiffness matrix is ​​obtained using the Green-Lagrange strain formula and the constitutive relation. Neumann boundary stiffness matrix Dirichlet boundary stiffness matrix In-domain residual matrix Neumann boundary residual matrix and the Dirichlet boundary residual matrix And assemble into the overall matrix iteration format to calculate the unknown coefficient increment corresponding to the current iteration step The undetermined coefficients are then cumulatively updated. Get the current displacement field increment And cumulatively update the displacement field After each iteration step is completed, the convergence is judged. If the norm of the current displacement increment satisfies If e is the set error threshold, it means convergence and entering the next load step; otherwise, continue with the above iterative steps.

Images