Improved SBFEM method based on consistent one-point integral theory

Abstract

An improved SBFEM method based on a consistent one-point integral theory belongs to the field of computational mechanics, and comprises the following steps: firstly, carrying out weighted average on a compatible strain field in a unit domain by adopting a linear strain smoothing technology to obtain a smooth strain field; secondly, representing a linear smooth function in a nodal function derivative form, and establishing a correction equation based on a discrete form of a divergence theorem; thirdly, adopting one-point integration at the geometric center of a sub-triangular domain, introducing a high-order derivative, and correcting the nodal function; fourthly, the corrected first-order partial derivative and the corrected second-order partial derivative of the shape function are substituted into a divergence theorem equation, and a corrected shape function derivative is solved; and finally, substituting the derivative of the nodal function into the stiffness matrix, determining a corrected stiffness matrix of the unit, and expanding the corrected stiffness matrix in a Taylor expansion form to obtain an efficient and stable corrected stiffness matrix. According to the method, the advantages of the consistent integral thought in the meshless method and the boundary discrete characteristic of the SBFEM are synergistically exerted, and elastic and elastic-plastic analysis is supported.

Description

Technical Field

[0001] This invention belongs to the field of computational mechanics and relates to a numerical simulation method based on the scaled boundary finite element method (SBFEM), specifically an SBFEM method improved based on the uniform one-point integral theory. Background Technology

[0002] In the fields of engineering mechanics and numerical simulation, the Scaled Boundary Finite Element Method (SBFEM), as a semi-analytical method combining the advantages of both the finite element method and the boundary element method, significantly improves the computational accuracy and adaptability of complex geometric structures (such as polygonal elements and crack networks) by reducing three-dimensional problems to two-dimensional boundary discretization. Its core advantage lies in combining analytical solutions with numerical discretization, making it particularly suitable for scenarios such as infinite-domain wave analysis and multiphysics coupling. However, as the scale of engineering problems increases, the computational efficiency bottleneck of SBFEM becomes increasingly prominent. Taking the elastoplastic analysis of earth-rock dams as an example, earth-rock dams, as key infrastructure in hydraulic engineering, must withstand complex stress and seepage field coupling during construction, impoundment, and operation. Factors such as the elastoplastic deformation of the soil, nonlinear constitutive relations, and the interaction between the dam body and the foundation all require high-precision numerical simulation to ensure engineering safety.

[0003] When dealing with large and complex structures like earth-rock dams, traditional SBFEM requires decomposing the fundamental domain of an n-sided polygonal element into n sub-triangular regions. Each sub-triangle needs three Gaussian integration points to ensure the numerical stability of the stiffness matrix, resulting in a total of 3n integration points per element. When analyzing earth-rock dam models containing thousands or even tens of thousands of elements, and considering the elastoplastic constitutive relations of the soil under different working conditions, each iteration requires a large number of integration calculations, which drastically increases computational resource consumption and significantly extends the computation cycle. In practical engineering applications, construction parties often need to quickly obtain key data such as stress and deformation of the earth-rock dam to guide construction and operation decisions, and the computational efficiency of traditional SBFEM is clearly insufficient to meet this requirement.

[0004] In existing technologies, the finite element method (FEM) is one of the most widely used methods. For example, Chinese invention patent (application number 202211101083.4) provides a method for soil slope stability analysis based on the precise finite element method. By employing a modified hyperbolic circular arc Mohr-Coulomb yield constitutive model and a high-order explicit Euler stress integral algorithm, it effectively improves the computational accuracy and stress update stability of traditional finite element methods in elastoplastic analysis. However, this method still relies on high-quality structured mesh discretization. When applied to earth-rock dams with complex geometries or uneven material distribution, the complex mesh generation process significantly increases the preprocessing workload. Furthermore, when calculating large deformation problems, the mesh is prone to distortion, which may lead to non-convergence and make it difficult to balance adaptability to complex geometries with computational efficiency. To eliminate the dependence on meshes, meshless methods have been developed. For example, a Chinese invention patent (application number 202111056868.X) provides a solution method and system for complex geometric flow domains with single-layer walls based on a meshless method. Utilizing particle discretization and a smooth, universal wall model, it demonstrates advantages in handling complex geometric boundaries and large deformation flow problems. However, such meshless methods are typically computationally expensive, and the calculation of nodal derivatives and the stability and accuracy control of boundary conditions are quite complex. When solving static problems involving strong seepage-stress coupling, such as earth-rock dams, their computational efficiency and solution accuracy are difficult to meet practical engineering requirements. Furthermore, some hybrid discretization methods (such as the material point method and the coupled Eulerian-Lagrange method) avoid mesh distortion through a strategy of separating the background mesh from the material points. However, the numerical flow of these methods differs significantly from that of traditional finite element methods. There are theoretical and technical barriers to their coupling with the proportional boundary finite element method (SBFEM), which is applicable to regular-irregular mixed structures such as dams. Furthermore, they lack sufficient support for simulating the dynamic coupling process of pore water pressure in saturated porous media, thus limiting their application in the seepage-stress coupling analysis of earth-rock dams throughout their entire life cycle.

[0005] Current optimization research on SBFEM mainly focuses on improving its analytical capabilities (such as mixed-order element extension and scaled boundary perfectly matched layer techniques), while progress in optimizing integration efficiency is limited. Traditional three-point integration, while ensuring numerical stability, results in redundant integration points that cause computational costs to increase linearly with the number of element edges. Single-point integration, although significantly reducing the number of integration points, introduces low-order errors by neglecting higher-order derivative terms, leading to a significant decrease in accuracy. Especially in the large deformation elastoplastic analysis of earth-rock dams, the coupling effect of material nonlinearity and geometric nonlinearity further exacerbates the accumulation of integration errors, making single-point integration schemes unsuitable for meeting engineering accuracy requirements. Therefore, how to overcome the efficiency bottleneck through innovative integration strategies without sacrificing computational stability has become the core challenge for the widespread adoption of SBFEM in the field of high-performance computing. Summary of the Invention

[0006] To address the technical challenge of low integration efficiency in existing SBFEM methods, this invention leverages the advantages of consistent integration in meshless methods and the boundary discretization characteristics of SBFEM. It proposes a stabilized single-point integration framework adapted to complex polygonal elements, establishing a linear consistent one-point integration method for two-dimensional polygonal elements. This is a highly efficient and accurate elastic-elastoplastic analysis method, achieving coupled analysis with the finite element method. The core of this invention lies in reducing the number of integration points for each subtriangle to one through the discrete divergence theorem and Taylor expansion correction techniques, significantly improving computational efficiency while avoiding numerical instability. In other words, this invention solves the problems of low computational efficiency and insufficient numerical stability caused by the large number of integration points in traditional SBFEM, providing a revolutionary numerical tool for high-performance elastoplastic analysis of complex geotechnical engineering projects such as earth-rock dams and underground structures.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0008] An improved SBFEM method based on uniform point integral theory includes the following steps:

[0009] Step S1: The compatible strain field within the element domain is weighted and averaged using linear strain smoothing techniques to obtain a smooth strain field. Specifically:

[0010] Step S11: Construct proportional boundary polygonal units and determine the polygon scaling center O, requiring that all proportional boundaries can be directly seen from O.

[0011] Step S12: Decompose the polygonal unit in step S11 into n sub-triangles. Each sub-triangle shares a common vertex with the scaling center and forms a triangular region with one side of the polygon. Each sub-triangle contains three integration points and one internal integration point O. c (x c ,y c The area of ​​the sub-triangle is A.

[0012] Step S13: In the two-dimensional scaled boundary finite element method (SBFEM), firstly, establish the compatible strain field ε(x) within the element domain for the n sub-triangles obtained in step S12; specifically:

[0013] The formula for calculating the compatible strain field ε(x) is as follows:

[0014] (1)

[0015] in, Represents the strain-displacement matrix; Represents the nodal displacement vector;

[0016] (2)

[0017] in, Let denote the derivative of the nodal shape function, where the subscript i indicates the row number, j indicates the column number, and x indicates the partial derivative with respect to x; The derivative of the nodal function is represented by the subscript i, which indicates the row number, j, which indicates the column number, and y, which indicates the partial derivative with respect to y; m indicates the continuation of the number, i.e., 1, 2, 3, ..., m;

[0018] Step S14: The compatible strain field ε(x) from step S13 is weighted and averaged using a linear smooth function to obtain a smooth strain field. To maintain the computational accuracy of the proportional boundary finite element method;

[0019] The linear smooth function is f(x) = [1 / x, y]. T The smooth strain field As shown in formula (3), the processing procedure is as follows:

[0020] (3)

[0021] in, Represents a smooth strain field; d represents the integration domain; Indicate the area of ​​the boundary infinitesimal element;

[0022] Step S2: Express the linear smooth function in the form of the derivative of a nodal shape function, and establish a corrected equation based on the discrete form of the divergence theorem. Specifically:

[0023] Step S21, the smooth strain field is expressed in the form of the modified derivative of the nodal shape function, as shown in formula (4):

[0024] (4)

[0025] in, Let represent the corrected derivative of the nodal shape function, where the subscript i indicates the row number, j indicates the column number, and x indicates the partial derivative with respect to x;

[0026] Step S22: Based on the discrete form of the divergence theorem, establish the corrected equation, as shown in formula (5):

[0027] (5)

[0028] in, Let represent the corrected derivative of the nodal shape function, where the subscript i indicates the row number, j indicates the column number, and x indicates the partial derivative with respect to x; Indicates the boundary of the integration field; The unit normal vector representing the boundary of the integration domain; the x-direction component; f(x) is a linear smooth function; This represents the partial derivative of a linearly smooth function x;

[0029] Step S3 involves using a point integral at the geometric center of the sub-triangular domain, introducing higher-order derivatives through Taylor expansion, and correcting the nodal shape functions. Specifically:

[0030] Step S31: Fix the integration point of each sub-triangle to the internal integration point O. c (x c ,y c This method configures only one integration point, replacing the traditional three-point integration, thus reducing the number of integration points and the computational load. The processing procedure is shown in formula (6):

[0031] (6)

[0032] Where the coordinates (x1, y1), (x2, y2), and (x3, y3) are the coordinates of the three vertices of the sub-triangle; , For the internal integration point O c The coordinates;

[0033] Step S32, for the nodal shape function At the internal integration point O c (x c ,y c Taylor expansion is performed, and a second derivative term is introduced. The higher-order term is used to compensate for the accuracy loss of single-point integration and maintain the calculation accuracy of the proportional boundary finite element method.

[0034] The specific expansion result of the modified nodal shape function is shown in formula (7):

[0035] (7)

[0036] Where HOT represents the abbreviation for higher-order terms; (x, y) represents the coordinates of any point inside the sub-triangle; Represents nodal shape functions; Indicates the integration point O inside the system. c (x c ,y c The node-shaped function at position () indicates the row number (i) and the column number (j). Indicates the integration point O inside the system. c (x c ,y c The first derivative of the shape function at the node () is given by the subscript i, which indicates the row number, j, which indicates the column number, and x, which indicates the partial derivative with respect to x. Indicates the integration point O inside the system. c (x c ,y c The nodal shape function derivative at point (), where y represents the partial derivative with respect to y; Indicates the integration point O inside the system.c (x c ,y c The derivative of the shape function at the node () is given by the subscript i, which indicates the row number, j, which indicates the column number, and xy, which indicates the first and second partial derivatives with respect to x and y. Indicates the integration point O inside the system. c (x c ,y c The derivative of the shape function at the node () is given by the index i, which indicates the row number, j, which indicates the column number, and yy, which indicates the second-order partial derivative with respect to y. Indicates the integration point O inside the system. c (x c ,y c The derivative of the shape function at the node () is given by the subscript i, which indicates the row number, j, which indicates the column number, and xx, which indicates the second-order partial derivative with respect to x.

[0037] Step S4: Substitute the first and second partial derivatives of the corrected shape function into the divergence theorem equation, and solve for the derivatives of the corrected shape function using Gaussian integrals. Specifically:

[0038] Step S41, substitute the expanded equation shown in formula (7) from step S32 into the corrected equation shown in formula (5) from step S22, and obtain the result as shown in formula (8):

[0039] (8)

[0040] in, Indicates the integration point O inside the system. c (x c ,y c A linear smooth function at (). Indicates the integration point O inside the system. c (x c ,y c The first derivative of the shape function at the node () is given by the subscript i, which indicates the row number, j, which indicates the column number, and x, which indicates the partial derivative with respect to x. Indicates the integration point O inside the system. c (x c ,y c The partial derivative of the linear smooth function x at point (). Indicates the integration point O inside the system. c (x c ,y c The partial derivative of the linear smooth function y at point (). Indicates the integration point O inside the system. c (x c ,y c The derivative of the shape function at the node () is given by the subscript i, which indicates the row number, j, which indicates the column number, and xx, which indicates the second-order partial derivative with respect to x. Let x represent the second-order area moment about x; Let x represent the first-order area moment and y represent the second-order area moment. This represents the second-order area moment about y; Indicates the integration point O inside the system. c (x c ,y c The derivative of the shape function at the modified node at (), where the subscript i indicates the row number, j indicates the column number, and xy indicates the first and second partial derivatives with respect to x and y; A is the area of ​​the subtriangle.

[0041] The process for processing second-order mask moments is as follows:

[0042] (9)

[0043] The corresponding first-order area moment is I x I y The processing procedure is as follows:

[0044] (10)

[0045] Step S42, for f(x) = [1xy] in step S4 T Performing a Taylor expansion, we get:

[0046] (11)

[0047] in, Indicates the integration point O inside the system. c (x c ,y c A linear smooth function at (). Indicates the integration point O inside the system. c (x c ,y c The partial derivative of the linear smooth function x at point (). Indicates the integration point O inside the system. c (x c ,y c The partial derivative of the linear smooth function y at point ().

[0048] Step S43: Integrate the two Gaussian integration points on each side of the sub-triangle to obtain a system of linear equations;

[0049] (12)

[0050] (13)

[0051] (14)

[0052] (15)

[0053] (16)

[0054] in, This represents the first set of data in the new matrix obtained after matrix operations on y, where the index i indicates the row number and j indicates the column number. The second set of data represents the new matrix obtained after matrix operations on y, where the index i indicates the row number and j indicates the column number; The third set of data in the new matrix obtained after matrix operations on y is represented by the index i, which indicates the row number and j, which indicates the column number. Indicates the internal integration point O c The x-coordinate; Indicates the internal integration point O c The y-coordinate; L represents the number of sides in the subtriangle; G represents the number of Gaussian points on each side of the subtriangle; Represents the Gaussian point of the line; Indicates the online Gaussian point The node-shaped function at the position, where the subscript i indicates the row number and j indicates the column number; This represents the x-direction component of the unit normal vector in subtriangle A; This represents the weight of the Gaussian point along the sub-boundary; Represents the x-coordinate of the Gaussian point on the line; This represents the y-coordinate of the Gaussian point on the line; It has no actual meaning; it is just a symbol for equation calculation. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row, j indicating the column, and x indicating the first-order partial derivative with respect to x. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row, j indicating the column, and xx indicating the second-order partial derivative with respect to x. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row and j indicating the column, and xy indicating the first and second partial derivatives with respect to x and y. This represents the first set of data in the new matrix obtained after matrix operations on x, where the index i indicates the row number and j indicates the column number. The second set of data represents the new matrix obtained after matrix operations on x, where the index i indicates the row number and j indicates the column number. The third set of data in the new matrix obtained after matrix operations on x is represented by the index i, which indicates the row number and j, which indicates the column number. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row, j indicating the column, and y indicating the first-order partial derivative with respect to y. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row and j indicating the column, and yx indicating the first and second partial derivatives with respect to y and x. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at point (), where the subscript i indicates the row number, j indicates the column number, and yy indicates the second-order partial derivative with respect to y; This represents the y-direction component of the unit normal vector in the subtriangle;

[0055] Step S44: Solve the linear equation system by Gaussian integration to obtain the modified shape function derivatives, as shown in formulas (17) and (18);

[0056] (17)

[0057] (18)

[0058] in, This represents the first set of data in the new matrix obtained after matrix operations on y, where the index i indicates the row number and j indicates the column number. The second set of data represents the new matrix obtained after matrix operations on y, where the index i indicates the row number and j indicates the column number; The third set of data in the new matrix obtained after matrix operations on y is represented by the index i, which indicates the row number and j, which indicates the column number. This represents the first set of data in the new matrix obtained after matrix operations on x, where the index i indicates the row number and j indicates the column number. The second set of data represents the new matrix obtained after matrix operations on x, where the index i indicates the row number and j indicates the column number. The third set of data in the new matrix obtained after matrix operations on x is represented by the index i, which indicates the row number and j, which indicates the column number.

[0059] Step S5: Substitute the modified nodal shape function derivatives from step S4 into the stiffness matrix to determine the modified stiffness matrix of the element. Expand the modified stiffness matrix using Taylor expansion to obtain an efficient and stable modified stiffness matrix, supporting elastic and elastoplastic analysis. Specifically:

[0060] (19)

[0061] in, Indicates the modified stiffness matrix; Represents the material stiffness matrix; Represents the corrected strain-displacement matrix;

[0062] (20)

[0063] in, Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row, j indicating the column, and x indicating the first-order partial derivative with respect to x. Indicates the integration point O inside the system. c (x c ,y c The modified nodal shape function derivative at position () is given by the subscript i indicating the row, j indicating the column, and y indicating the first-order partial derivative with respect to y.

[0064] (twenty one)

[0065] (twenty two)

[0066] in, This is the corrected strain-displacement matrix; This is the derivative of the corrected strain-displacement matrix in the y-direction; This is the derivative of the corrected strain-displacement matrix in the x-direction.

[0067] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0068] (1) This invention compensates for the loss of single-point integral accuracy by using a weighted average of linear smooth functions to accommodate the strain field, deriving the derivative of the corrected shape function using the divergence theorem, and introducing a second derivative term through Taylor expansion.

[0069] (2) The present invention can improve computational efficiency: each sub-triangle is configured with only 1 integration point, replacing the traditional 3-point integration, reducing the number of integration points and reducing the amount of computation;

[0070] (3) The present invention can construct proportional boundary polygonal units, handle complex geometric structures, and support elastic and elastoplastic analysis;

[0071] (4) This invention improves the SBFEM integration efficiency by more than 40% while ensuring the accuracy of calculation, and does not require additional hardware resources. It is especially suitable for engineering practice of large-scale elastoplastic problems such as earth-rock dams and underground structures.

[0072] In summary, this invention configures only one integration point per sub-triangle, replacing the traditional three-point integration, thus reducing the number of integration points and computational load. By introducing a linear smooth strain field and higher-order correction terms in Taylor expansion, it effectively overcomes the technical bottlenecks of traditional SBFEM, such as low computational efficiency and insufficient numerical stability caused by redundant integration points. It provides an innovative numerical solution for high-performance elastoplastic analysis of complex geotechnical engineering projects such as earth-rock dams and underground structures, achieving coupling with existing finite element software and significantly improving computational efficiency. Attached Figure Description

[0073] Figure 1 This is a schematic diagram of the main process of the method of the present invention;

[0074] Figure 2 This is a one-point integration method in SBFEM;

[0075] Figure 3 A schematic diagram of a one-point integration scheme for a sub-triangle;

[0076] Figure 4 Cook membrane model;

[0077] Figure 5 Schematic diagrams of different mesh structures of Cook membrane; Figure 5 (a) in the diagram represents a 1×1 grid; Figure 5 (b) in the diagram represents a 2×2 grid; Figure 5 (c) in the diagram represents a 4×4 grid. Figure 5 (d) in the diagram represents an 8×8 grid. Figure 5 (e) in the diagram represents a 16×16 grid.

[0078] Figure 6 Schematic diagrams showing displacement results under different conditions; Figure 6 (a) in the figure is a schematic diagram of the displacement results of the 4×4 grid fully integrated scheme (Comparative Example 1); Figure 6 (b) in the figure is a schematic diagram of the displacement results of the 4×4 grid reduced integration scheme (Comparative Example 2); Figure 6 (c) in the figure is a schematic diagram of the displacement results of the one-point integration scheme proposed in this invention with a 4×4 grid. Detailed Implementation

[0079] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, but the scope of protection of the present invention is not limited thereto.

[0080] This embodiment provides an improved SBFEM method based on the uniform one-point integral theory, including the following steps:

[0081] Step S1: The compatible strain field within the element domain is weighted and averaged using linear strain smoothing techniques to obtain a smooth strain field. Specifically:

[0082] Step S11, see Figure 2 Construct proportional boundary polygonal units, determine the polygon scaling center O, and require that all proportional boundaries be directly visible from O.

[0083] Step S12, see Figure 2 In step S11, the polygonal unit is decomposed into n sub-triangles. Each sub-triangle shares a common vertex with the scaling center and forms a triangular region with one side of the polygon. Each sub-triangle contains three integration points and one internal integration point O. c (x c ,y c The area of ​​the sub-triangle is A.

[0084] Step S13: In the two-dimensional scaled boundary finite element method (SBFEM), firstly, a compatible strain field ε(x) is established within the element domain for the n sub-triangles obtained in step S12, as shown in formula (1).

[0085] Step S14: The compatible strain field ε(x) from step S13 is weighted and averaged using a linear smooth function to obtain a smooth strain field. To maintain the computational accuracy of the proportional boundary finite element method; the linear smoothing function is f(x) = [1 / x, y]. T The smooth strain field As shown in formula (3);

[0086] Step S2: Express the linear smooth function in the form of the derivative of a nodal shape function, and establish a corrected equation based on the discrete form of the divergence theorem. Specifically:

[0087] Step S21, the smooth strain field is expressed in the form of the modified derivative of the nodal shape function, as shown in formula (4);

[0088] Step S22: Establish the corrected equation based on the discrete form of the divergence theorem, as shown in formula (5);

[0089] Step S3 involves using a point integral at the geometric center of the sub-triangular domain, introducing higher-order derivatives through Taylor expansion, and correcting the nodal shape functions. Specifically:

[0090] Step S31, see Figure 2 Fix the integration point of each sub-triangle to the internal integration point O. c (xc ,y c This method configures only one integration point, replacing the traditional three-point integration, thus reducing the number of integration points and the amount of computation. The processing procedure is shown in formula (6).

[0091] Step S32, for the nodal shape function At the internal integration point O c (x c ,y c Taylor expansion is performed, and a second derivative term is introduced. The higher-order term is used to compensate for the accuracy loss of single-point integral and maintain the calculation accuracy of the proportional boundary finite element method. The specific expansion result of the modified nodal shape function is shown in formula (7).

[0092] Step S4: Substitute the first and second partial derivatives of the corrected shape function into the divergence theorem equation, and solve for the derivatives of the corrected shape function using Gaussian integrals. Specifically:

[0093] Step S41: Substitute the expanded equation (7) from step S32 into the corrected equation (5) from step S22 to obtain the result (8). The second-order surface moment processing is shown in equation (9); the corresponding first-order surface moment is I. x I y The processing procedure is shown in formula (10).

[0094] Step S42, for f(x) = [1xy] in step S4 T Performing a Taylor expansion, we obtain formula (11);

[0095] Step S43, see Figure 3 Integrating the two Gaussian integration points on each side of the sub-triangle, we obtain a system of linear equations, as shown in formulas (12) to (16).

[0096] Step S44: Solve the linear equation system by Gaussian integration to obtain the modified shape function derivatives, as shown in formulas (17) and (18);

[0097] Step S5: Substitute the modified nodal shape function derivatives from step S4 into the stiffness matrix to determine the modified stiffness matrix of the element. Expand the modified stiffness matrix using Taylor expansion to obtain an efficient and stable modified stiffness matrix that supports elastic and elastoplastic analysis. The processing procedure is shown in formulas (19) to (22).

[0098] To demonstrate the effectiveness of this technical solution, numerical simulations were performed on the Cook membrane problem (Example 1) using this technical solution. Details are as follows:

[0099] Step 1, the geometry and boundary conditions of the Cook membrane are as follows: Figure 4As shown, the beam is tapered, with a length of l2 = 48 mm, a fixed end height of l4 = 16 mm, and a free end height of l3 = 44 mm. The left end of the beam is completely fixed (all degrees of freedom are restricted), and a uniformly distributed shear load of F = 1 N / mm² is applied to the right end. The material behavior is assumed to be linear elastic, with an elastic modulus E = 1 Pa and a Poisson's ratio of ν = 0.33. The problem is solved under plane stress conditions.

[0100] Step 2 involves dividing the grid into different meshes, constructing polygonal units, and determining the polygon scaling center, ensuring that all boundaries are directly visible from the scaling center. Five different mesh configurations with varying numbers of units were used. See [link / reference] Figure 5 The model is discretized into grids of 1×1, 2×2, 4×4, 8×8 and 16×16 quadrilateral elements, respectively.

[0101] Step 3: In the five discrete grids with different numbers of cells, during calculation, the multiple integration points of each cell are converted into a single center point. See [link to relevant documentation] for the processing method. Figure 2 With formula (6);

[0102] Step 4: During the calculation, the strain field after processing with a linear smooth function is calculated at the geometric center point to obtain the processed shape function derivative. The processing method is shown in formulas (7)-(18).

[0103] Step 5: Determine the corrected stiffness matrix of the element using the smoothed nodal derivatives. Expand the corrected stiffness matrix using Taylor expansion, as shown in formulas (19)-(22). Calculate and analyze to obtain the lower right corner of the beam (point A, see...). Figure 4 The vertical displacement of the sample was determined. Results from other established integration schemes in SBFEM were compared to validate the proposed one-point integration method using smoothing functions. (The comparative analysis included: ① a full integration scheme with three integration points; ② a reduced integration scheme using one integration point without strain smoothing techniques. See Table 1 for explanations of the different integration types.)

[0104] Table 1. Detailed information on different integration schemes in SBFEM

[0105]

[0106] Step 6: Using the finite element software GEODYNA 8.0, which includes all the integration schemes of SBFEM listed in Table 1, calculate the above mesh configuration. The displacement calculation results for different integration schemes are as follows: Figure 6As shown, the displacement values ​​of point A are summarized in Table 2, where the reference solution is given by the literature (Y. Long, Y. Xu, Generalized forming triangle membrane element with vertex rigid rotational freedoms. Finite Elem. Anal. Des. 17(4)(1994) 259–271. https: / / doi.org / 10.1016 / 0168-874X(94)90002-7).

[0107] Table 2, Displacement values ​​at point A of the Cook membrane

[0108]

[0109] In the Cook membrane bending problem, using the full integration scheme as the benchmark, as the mesh density increases, both the proposed one-point integration scheme and the full integration scheme converge to the reference value for displacement solutions, while the reduced integration scheme exhibits significant errors. The errors in the reduced integration scheme stem from numerical instability and insufficient displacement field sampling, leading to inaccurate approximations of the integral terms. Our proposed method eliminates these defects by introducing higher-order derivative corrections through linear strain smoothing techniques and Taylor expansion, achieving a maximum error of only 0.16%, comparable to the accuracy of the full integration scheme. This verifies that the proposed one-point integration scheme reliably maintains computational accuracy while significantly reducing the number of integration points (reducing computational cost by 45%-60%), especially demonstrating excellent effectiveness and accuracy under complex stress states such as bending, thus solving the problem of balancing computational efficiency and accuracy in traditional SBFEM.

[0110] The above embodiments are merely illustrative of the implementation methods of the present invention, but should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the protection scope of the present invention.

Claims

1. An improved SBFEM method based on consistent point- integration theory, characterized in that, The SBFEM method comprises the following steps: Step S1, the compatible strain field in the unit domain is weighted and averaged by using linear strain smoothing technique to obtain a smoothed strain field ; Step S2, the linear smooth function is expressed in the form of node shape function derivative, and a modified equation is established based on the discrete form of divergence theorem: Step S3, one-point integration is used at the geometric center of the sub-triangle domain, high-order derivatives are introduced through Taylor expansion, and the node shape function is modified; Step S4, the first-order and second-order partial derivatives of the modified shape function are substituted into the divergence theorem equation, and the Gaussian integral is used to solve the derivative of the modified shape function; Step S5, the derivative of the modified node shape function in step S4 is brought into the stiffness matrix to determine the modified stiffness matrix of the element, and the modified stiffness matrix is expanded in the form of Taylor expansion to obtain a high-efficiency and stable modified stiffness matrix, which supports elastic and elastoplastic analysis.

2. The SBFEM method based on the improved consistent point integration theory according to claim 1, wherein, The step S1 is specifically: Step S11, a proportional boundary polygon element is constructed, a polygon scaling center O is determined, and all proportional boundaries can be directly viewed from O; Step S12, decompose the polygonal element in step S11 into n sub-triangles, each sub-triangle has the scaling center as a common vertex, and forms a triangular area with one side of the polygon, and three integral points are arranged in each sub-triangle, one internal integral point O c (x c ,y c ), and the area of the sub-triangle is A; Step S13, in the two-dimensional proportional boundary finite element method SBFEM, first, a compatible strain field ε(x) in the element domain is established for the n sub-triangles obtained in step S12; Step S14, the compatible strain field ε(x) in step S13 is weighted and averaged by using a linear smoothing function to obtain a smoothed strain field , which maintains the calculation accuracy of the proportional boundary finite element method.

3. The SBFEM method based on the improved consistent point integration theory of claim 2, wherein, In the step S1: In the step S13, the compatible strain field ε(x) is calculated according to the formula: （1）； wherein, denotes the strain-displacement matrix; denotes the nodal displacement vector; （2）； wherein denotes the derivative of the nodal shape function, the index i denotes the row number, j denotes the column number, and x denotes the partial derivative with respect to x; denotes the derivative of the nodal shape function, the index i denotes the row number, j denotes the column number, and y denotes the partial derivative with respect to y; m denotes the continuation of the number; In the step S14, the linear smoothing function is f(x) = [1 x y] T ; the smoothing strain field As shown in equation (3), the process is as follows: （3）； wherein, denotes the smooth strain field; denotes the integration domain; d denotes the boundary surface element area.

4. The SBFEM method based on the improved consistent point integration theory of claim 3, wherein, The step S2 is specifically: Step S21, the smooth strain field is expressed in the form of modified derivative of node shape function; Step S22, a modified equation is established based on the discrete form of divergence theorem.

5. The SBFEM method based on the improved consistent point integration theory according to claim 4, wherein, In the step S2: The step S21 is shown in formula (4): （4）； wherein denotes the modified derivative of the nodal shape function, the index i denotes the row number, j denotes the column number, and x denotes the partial derivative with respect to x; In the step S22, the modified equation is shown in formula (5): （5）； wherein denotes the modified derivative of the nodal shape function, the index i denotes the row number, j denotes the column number, and x denotes the partial derivative with respect to x; denotes the integral domain boundary; denotes the unit normal vector of the integral domain boundary; the x-direction component; f(x) is a linear smoothing function; denotes the partial derivative of the linear smoothing function x.

6. The SBFEM method based on the improved consistent point integration theory of claim 5, wherein, The step S3 is specifically: Step S31, fix the integral point of each sub-triangle at the internal integral point O c (x c ,y c ), only 1 integral point is configured; Step S32, the node shape function At the internal integration point O c (x c , y c ) is Taylor expanded, the second derivative term is introduced, the accuracy loss of single-point integration is compensated by high-order terms, and the calculation accuracy of the proportional boundary finite element method is maintained.

7. The SBFEM method based on the improved consistent point integration theory of claim 6, wherein, In the step S3: The processing process of the step S31 is shown in formula (6): （6）； Wherein, the coordinates of the three vertices of the sub-triangle are (x1, y1), (x2, y2), (x3, y3); , The coordinates of the internal integration point O c ; In the step S32, the specific expansion result of the node shape function is shown in formula (7): （7）； where H.O.T stands for high order term; (x, y) stands for an arbitrary point in the sub-triangle; N stands for nodal shape function; N stands for nodal shape function at internal integration point O c (x c ,y c ); N stands for nodal shape function at internal integration point O c (x c ,y c ); N stands for nodal shape function derivative at internal integration point O c (x c ,y c ); N stands for nodal shape function derivative at internal integration point O c (x c ,y c ); N stands for nodal shape function derivative at internal integration point O c (x c ,y c ); N stands for nodal shape function derivative at internal integration point O c (x c ,y c ).

8. The SBFEM method based on the improved consistent point integration theory of claim 7, wherein, The step S4 is specifically: Step S41, the expansion result shown in formula (7) in step S32 is brought into the modified equation shown in formula (5) in step S22 to obtain a result shown in formula (8): （8）； where represents a linearly smooth function at the interior integration point O c (x c ,y c ); represents a nodal-shaped function first derivative at the interior integration point O c (x c ,y c ); the subscript i represents the row number, j represents the column number, and x represents the partial derivative with respect to x; represents a partial derivative of a linearly smooth function x at the interior integration point O c (x c ,y c ); represents a partial derivative of a linearly smooth function y at the interior integration point O c (x c ,y c ); represents a nodal-shaped function derivative at the interior integration point O c (x c ,y c ); the subscript i represents the row number, j represents the column number, and xx represents the second-order partial derivative with respect to x; represents the second-order area moment with respect to x; represents the first-order area moment with respect to x and the second-order area moment with respect to y; represents the second-order area moment with respect to y; represents a modified nodal-shaped function derivative at the interior integration point O c (x c ,y c ); the subscript i represents the row number, j represents the column number, and xy represents the first- and second-order partial derivatives with respect to x and y; and A is the area of the sub-triangle. The second-order surface mask processing process is as follows: （9）； The corresponding first order area moment is I x , y The process is as follows: （10）； Step S42, for f(x) = [1 x y] in step S4 T Carrying out Taylor expansion, we get: （11）； wherein, denotes a linearly smooth function at the interior integration point O c (x c ,y c ); denotes the partial derivative of the linearly smooth function x at the interior integration point O c (x c ,y c ); denotes the partial derivative of the linearly smooth function y at the interior integration point O c (x c ,y c ); Step S43, two Gaussian integral points on each edge of the sub-triangle are integrated to obtain a linear equation group; （12）； （13）； （14）； （15）； （16）； wherein, represents the first set of data of the new matrix after matrix operation of y; represents the second set of data of the new matrix after matrix operation of y; represents the third set of data of the new matrix after matrix operation of y; represents the x coordinate of the interior integration point O c ; represents the y coordinate of the interior integration point O c ; L represents the number of edges in the sub-triangle; G represents the number of line Gauss points on each edge in the sub-triangle; represents a line Gauss point; represents the node shape function at the line Gauss point ; represents the x direction component of the unit normal vector in the sub-triangle A; represents the weight of the line Gauss point along the sub-edge boundary; represents the x coordinate of the line Gauss point; represents the y coordinate of the line Gauss point; has no actual meaning, and is a symbol for equation calculation; represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); represents the first set of data of the new matrix after matrix operation of x; represents the second set of data of the new matrix after matrix operation of x; represents the third set of data of the new matrix after matrix operation of x; represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); represents the derivative of the modified node shape function at the interior integration point O c (x c ,y c ); y-component of the unit normal vector in the sub-triangle; Step S44, the Gaussian integral is used to solve the linear equation group to obtain the derivative of the modified shape function, which is shown in formula (17) and (18); （17）； （18）； wherein, represents the first set of data of the new matrix resulting from the matrix operation on y; represents the second set of data of the new matrix resulting from the matrix operation on y; represents the third set of data of the new matrix resulting from the matrix operation on y; represents the first set of data of the new matrix resulting from the matrix operation on x; represents the second set of data of the new matrix resulting from the matrix operation on x; represents the third set of data of the new matrix resulting from the matrix operation on x.

9. The SBFEM method based on the improved consistent point integration theory of claim 8, wherein, The step S5 is specifically: （19）； wherein, represents a modified stiffness matrix; represents a material stiffness matrix; represents a modified strain-displacement matrix; （20）； wherein, denotes the modified nodal shape function derivative at the internal integration point O c (x c ,y c ); denotes the modified nodal shape function derivative at the internal integration point O c (x c ,y c ); （21）； （22）； wherein, is the modified strain-displacement matrix; is the derivative of the modified strain-displacement matrix in the y direction; is the derivative of the modified strain-displacement matrix in the x direction.

Images