CN111597649A - Elastic crack problem simulation method based on stable generalized finite element - Google Patents

Abstract

The elastic crack problem simulation method based on the stable generalized finite element improves the precision and condition number of the traditional GFEM/XFAM, so that the numerical simulation of the crack problem is more accurate and stable. In addition, the design of the method is highly automatic and standardized, and the method can be embedded into current structural analysis business software to promote the upgrading of the software and can also be used as a technical collection for the independent development of engineering software. Finally, the technical scheme of the invention is one of the steps of development of key engineering software, and has the advantages of strong algorithm innovation, high simulation precision, clear and reasonable software architecture, strong embeddability and expansibility, high parallelization degree and great development potential and application prospect.

Description

Elastic crack problem simulation method based on stable generalized finite element

Technical Field

The invention relates to the technical field of engineering computational mechanics software application, in particular to an elastic crack problem simulation method based on a stable generalized finite element.

Background

The problem of elastic crack and crack propagation simulation is an important problem in mechanics, and the relevance and importance of the problem are derived from the wide application of numerical fracture mechanics in fatigue life prediction of safety critical parts such as airplane fuselages, pressure vessels, automobile parts and castings. The fatigue failure is usually caused by the generation and propagation of surface or near-surface cracks, for cracks with any shape, only a numerical method can be adopted to simulate the crack propagation problem, and the main link of numerical simulation of the crack problem is numerical solution of an elastic mechanics equation set.

In the prior art, the omega expression is cracked_OAnd a fracture region of the fracture tip O, set as a boundary of Ω,

is the unit of outward normal vector. By essential boundaries_DAnd natural boundary_NComposition, as shown in figure 1. Expressing a function or vector of vector values in space with bold symbols, e.g. u ═ u₁，u₂]^T. The elastic equation set for the crack problem within Ω is as follows:

wherein the content of the first and second substances,

is stress tensor, f is physical strength, g and u₀Respectively a natural boundary condition and an intrinsic boundary condition,

is that_OOutward normal vector. The x and y directions are as shown in fig. 1, (r, θ) are the corresponding polar coordinates.

The strain tensor e (u) is recorded as

The model is specifically represented as:

where the former is the plane strain and the latter is the plane stress, E is the Young's modulus and v is the Poisson's ratio. The main part of the problem (P1) solution u corresponds to the following type I and type II open models u_IAnd u_II[13,10]：

Wherein M is_IAnd M_IIIs a coefficient related to the pressure intensity factor. E, v and kappa are Kolosov constants ([5,13, 10)]). u in the crack_OThe site is discontinuous and singular around the fracture tip O. Therefore, the conventional Finite Element Method (FEM) needs to continuously encrypt and update the grid during the fracture propagation process, and the computation amount is extremely high.

In response to the difficulty, in the last two decades, scholars develop a generalized or extended finite element method (GFEM/XFEM), which adds an enrichment function (enrichment) on the basis of a fixed regular grid FEM to solve the crack problem, and the efficiency is very ideal. GFEM/XFEM has been widely used for a variety of engineering problems such as fracture propagation, material modeling, multiphase flow and fluid-structure interactions etc. [2,7 ]. Some common finite element software has also integrated GFEM/XFAM into the framework, such as Ansys, Abaqus, LS-dyna, etc. Although GFEM/XFAM solves the difficulty of large grid calculation amount and achieves higher calculation precision, the condition number of the rigidity matrix is large due to the linear correlation existing between the original finite element function and the newly added enrichment function, so that a remarkable truncation error can be caused when a linear equation set is solved, and even the failure of a numerical simulation process can be caused. The condition number problem has been one of the major technical difficulties in GFEM/XFEM research [4,7], and the associated research efforts have been numerous [4,11,14,15], but have not been well solved to date. Furthermore, since the program framework of GFEM/XFEM differs from that of general FEM, the parallelization and high performance computation aspects of GFEM/XFEM are also challenging tasks compared to FEM.

Disclosure of Invention

The invention provides an elastic crack problem simulation method based on a stable generalized finite element, aiming at overcoming the technical defects that the existing generalized or extended finite element method has large condition number of a rigidity matrix, and can cause obvious truncation error in solving and even failure in a numerical simulation process.

In order to solve the technical problems, the technical scheme of the invention is as follows:

the elastic crack problem simulation method based on the stable generalized finite element comprises the following steps:

s1: constructing a GFEM/XFEM with an aggregation enrichment strategy according to the crack problem GFEM/XFEM;

s2: correcting an enrichment function in an enrichment strategy by adopting a linear Heaviside function;

s3: changing a finite element PU function of the standard FEM, and constructing a novel SGFEM basic function;

s4: processing the SGFEM basic function based on a local principal component analysis technology, eliminating redundancy caused by a plurality of enrichment functions of a single node, and obtaining an SGFEM model;

s5: and solving the crack problem through the SGFEM model to complete the simulation of the crack problem.

Wherein the step S1 specifically includes:

the region omega is subdivided into quasi-regular triangular or quadrangular grid cells e_sThe dimension parameter of the grid is represented by h, the grid is simple, fixed and connected with the crack_OIndependently, { (x)_i，y_i)：i∈I_hRecording the data as a grid node set; let phi_i，i∈I_hIs a standard finite element function, and

is phi_iThe supporting set of (2); the main idea of GFEM/XFAM is by using the PU function [7,6,10,13 ]]Coupling a finite element function and special functions expressing true solution properties to achieve a high-precision approximation effect, wherein the special functions are called enrichment functions; approximation function u of crack problem GFEM/XFEM_hThe following were used:

wherein I_HNode sets of cells intersecting the fracture, but with the support set removed, nodes containing fracture tips, a_i、b_i、

Denotes the intermediate variable, S_jRepresenting an intermediate function; i is_sThe GFEM/XEM is called GFEM/XEM with geometric enrichment strategy as a node set inside a circle B (O, R) with a radius R and a crack tip as a center; in the formula (2), function

Is a Heaviside function for simulating the discontinuity of the true solution at the crack

Called the crack tip strengthening function, wherein R is the radius of the enrichment range, R, theta are polar coordinates, and the intermediate variables of the mathematical expression; using enrichment functions H and

and completing the construction of GFEM/XFEM.

In the above scheme, the enrichment function H is

The introduction of (a) causes a linear correlation between them and the finite element function, so that the condition ratio of the stiffness matrix is much larger than that of the finite element, thereby possibly causing a significant truncation error when solving the linear equation system and even possibly causing a failure of the numerical simulation process. Therefore, the condition number of the traditional GFEM/XFAM is improved by adopting the idea of a stable generalized finite element method, the approximation error is improved, and the improved means is to modify an enrichment function and changePU function, local principal component analysis.

Wherein the step S2 specifically includes:

in the formula (2), a linear Heaviside function is used

Replaces the Heaviside function H of the conventional GFEM/XFEM, and then

And

the correction is as follows:

wherein

Is a finite element interpolation operator, i.e. for a continuous function F,

wherein the parameters v, phi_jIntermediate variables that are mathematical expressions; and constructing the GFEM/XFEM with the set enrichment strategy based on the corrected enrichment function.

In the above scheme, it can be found from equation (3) that the enrichment functions subtract their finite element interpolation, which is one of the actions to reduce the linear dependence on the finite element functions [1,9,16,17,18 ].

Wherein, the step S3 specifically includes: the standard FEM function phi in the formula (2)_iAs a function of PU, this is also one of the causes of deterioration of the condition number of the stiffness matrix, see [11,17 ]]Discussion of (1). Therefore, the PU function is modified to a high-order polynomial PU function, specifically, let:

Q₀＝(1-ξ)²(1+2ξ)，Q₁＝ξ²(3-2ξ)，ξ∈[0，1]

is a one-dimensional reference unit [0,1 ]]The PU function of (1), from which a two-dimensional reference cell [0,1 ] is obtained]×[0，1]PU function of

The method comprises the following specific steps:

in order to obtain any actual unit e_sThe PU function of (1), let (x, y) be F_s(ξ) is a reference cell to e_sAffine transformation of_sA PU function of

Assembling these unit PU functions according to the nodes obtains the required PU function Q_i，i∈I_h(ii) a Based on the modified enrichment function (3) and the new PU function Q_iObtaining a novel SGFEM basic function expression:

equation (4) is a novel SGFEM basic function expression, and the SGFEM can reach an optimal convergence order mathematically proved according to the ideas in [16,18 ].

Wherein, the step S4 specifically includes: in equation (4), it is found that some nodes are enriched with multiple functions, e.g., I_HEach node in (1) is enriched with three functions, I_SEach node in (1) is enriched with four functions, and I_H∩I_HEach node in (a) is enriched with seven functions. Thus using

Represents each enrichment node (x)_i，y_i) The enrichment function of (2) then has:

obviously corresponds to I ∈ I_H，i∈I_s，i∈I_H∩I_sα equals 3, 4,7 respectively, and the linear correlation between enrichment functions for each node is another source of large condition number of stiffness matrix, for which a local principal component analysis is proposed [ 8]]Local Principal Component Analysis (LPCA) technique to eliminate redundancy caused by multiple enrichment functions at a single node, because the elastic mechanics equation is a vector valued equation for any I ∈ I_H∪I_sDenotes E_iThe basis functions for the x-and y-directions are as follows:

for any I ∈ I_H∪I_sMemory for recording

Is a function of

In the inner product

A covariance matrix of; as will be apparent from the above description,

is a α×α submatrix of the total stiffness matrix A, and carries out PCA (principal component analysis) based on a covariance matrix to obtain a α×α matrix

And

wherein:

each column contains the coefficients of one principal component,

consists of the percentage of each main component; according to the PCA property, the new principal component function is:

are orthogonal under inner product (5) and can be based on vectors

Determining contribution of each principal component, and deleting

λ is a predetermined parameter, where λ is 10^-10This is a very small percentage; the convergence performance is not reduced due to the LPCA, the enrichment function knowledge is recombined, the original approximation space is not changed, only the redundancy with extremely small contribution is deleted, and in addition, the LPCA can be executed according to the rigidity matrix corresponding to the SGFEM basic function (4) without a large amount of extra calculation. In this sense, LPCA can be viewed as a preprocessing algorithm to solve linear systems.

The LPCA is regarded as a preprocessing algorithm for solving a linear system, and the specific implementation process is as follows:

based on enrichment function

Integrating the stiffness matrix A and the load vector b for an arbitrary node (x)_i，y_i)，i∈I_H∪I_sThe subscript of the permutation of multiple enrichment functions is J_i，...，J_i+ α -1, and extracting the submatrix from A

Is composed of

For each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

and updating the rigidity matrix A and the load vector b as follows:

for each (x)_i，y_i)，i∈I_H∪I_sIf, if

Updating the stiffness matrix A and the load vector b as follows:

A(J_i+j-1，：)＝0，A(：，J_i+j-1)＝0，

A(J_i+j-1，J_i+j-1)＝1，b(J_i+j-1，：)＝0

for enrichment function

Executing the same process as the above, updating the rigidity matrix A and the load vector b, and obtaining a new rigidity matrix A and a new load vector b after LPCA, thereby obtainingTo SGFEM model.

In step S5, the SGFEM model is used to solve the crack problem, that is, a linear equation set is solved, specifically:

the system of linear equations is scaled:

let E_ii＝(A_ii)^-1/2Obtaining a diagonal matrix E, and obtaining EAE and Eb after the rigidity matrix A and the load vector b are scaled;

solving a scaled linear system of equations:

EAEu＝Eb (6)

reverse pretreatment:

carrying out inverse scaling on the scaled solution u to obtain u ═ Eu; meanwhile, u is also required to be subjected to inverse pretreatment according to the LPCA pretreatment process:

based on enrichment function

Integrate the solution vector u for any node (x)_i，y_i)，i∈I_H∪I_sSubscript of the arranged multiple enrichment function is J_i，...，J_i+ α -1, and extracting the submatrix from A

Comprises the following steps:

for each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

if it is

Updating

The sum vector u is:

for each (x)_i，y_i)，i∈I_H∪I_sThe updated solution vector u is:

for enrichment function

Executing the same process as the above to obtain a new solution vector u after the inverse LPCA; therefore, the crack problem is solved, and the simulation of the crack problem is completed.

The elastic crack problem simulation method based on the stable generalized finite element further comprises the steps of carrying out parallelization processing on the SGFEM model, solving the crack problem by using the parallelized SGFEM model, and completing the simulation of the crack problem.

The process of the parallel processing of the SGFEM model specifically comprises the following steps: screening out nodes needing to be subjected to LPCA (low power allocation algorithm) of an SGFEM (generalized minimum finite element model) according to existing marks, dispersing the nodes to obtain N units, and dividing the nodes according to the number p of computer processes, wherein N is not required to be divided by the processes p; first p-1 process handling

The p number process processes the remaining N-nel (p-1) units; independently assembling a rigidity matrix A and a load vector b in each process, and finally aggregating to a specific processPerforming the following steps; in LPCA, it is obtained at each node according to A

And updating the matrix A and the vector b, thereby completing the parallelization processing of the SGFEM model.

In the above scheme, the SGFEM model is very different from FEM, so its parallelization algorithm needs a new design. Dispersing the solution area to obtain N units, dividing the area according to the number p of computer processes (the number N of units is not required to be divided by the number p of processes): first p-1 process handling

The p number process processes the remaining N-nel (p-1) cells. The stiffness matrix A and the load vector b are independently assembled in each process and finally aggregated in a specific process (such as process No. 0). Note that no data communication is required between different processes when assembling the stiffness matrix and the load vector. Because the rigidity matrix generated by each sub-area adopts a sparse storage mode, the data communication volume is very small when aggregation operation is finally carried out, and the time spent is far shorter than the numerical calculation.

In the above scheme, only the nodes near the fracture need to be analyzed for LPCA. If all nodes are divided into sessions according to the number of processes, the calculation amount of a part of the processes is small, the calculation amount of a part of the processes is large, and the parallel efficiency is low. Firstly, screening out nodes needing LPCA according to the existing marks, and then dividing the nodes according to the process number by the method. In LPCA, it is available at each node according to A

And updates matrix a and vector b.

The method comprises the steps of solving the SGFEM model after parallelization by adopting a parallelization Krylov subspace method, namely solving a linear system EAEx (equal to Eb) in a parallelization mode to obtain x, solving the crack problem by adopting an inverse LPCA method similar to LPCA after inverse scaling of a solution vector x, and completing simulation of the crack problem.

Compared with the prior art, the technical scheme of the invention has the beneficial effects that:

the elastic crack problem simulation method based on the stable generalized finite element improves the precision and condition number of the traditional GFEM/XFAM, so that the numerical simulation of the crack problem is more accurate and stable. In addition, the design of the method is highly automatic and standardized, and the method can be embedded into current structural analysis business software to promote the upgrading of the software and can also be used as a technical collection for the independent development of engineering software. Finally, the technical scheme of the invention is one of the steps of development of key engineering software, and has the advantages of strong algorithm innovation, high simulation precision, clear and reasonable software architecture, strong embeddability and expansibility, high parallelization degree and great development potential and application prospect.

Drawings

FIG. 1 is a prior art crack problem representation schematic;_Dand_Nrespectively, an intrinsic boundary and a natural boundary. (x, y) and (r, θ) are a rectangular coordinate system and a polar coordinate system, respectively, with an origin O.

FIG. 2 is a schematic flow diagram of the process of the present invention;

FIG. 3 is a flow chart of a SGFEM solution for elastic fracture problem;

FIG. 4 is a schematic representation of a region with a crack in one embodiment; crack (crack)_oNodes ○ are enriched by the Heaviside function, nodes × are enriched by the singular function;

FIG. 5 is a diagram of XFEM/GFEM and SGFEM H in one embodiment¹Error schematic diagram;

FIG. 6 is a SCN schematic of XFEM/GFEM and SGFEM in one embodiment;

FIG. 7 is a diagram of the computational time required for different numbers of processors with a fixed problem-solving scale;

FIG. 8 is a graph showing the time required to solve a problem of different size when the number of processors is fixed.

Detailed Description

The drawings are for illustrative purposes only and are not to be construed as limiting the patent;

for the purpose of better illustrating the embodiments, certain features of the drawings may be omitted, enlarged or reduced, and do not represent the size of an actual product;

it will be understood by those skilled in the art that certain well-known structures in the drawings and descriptions thereof may be omitted.

The technical solution of the present invention is further described below with reference to the accompanying drawings and examples.

Example 1

As shown in FIG. 2, the elastic fracture problem simulation method based on the stable generalized finite element comprises the following steps:

s1: constructing a GFEM/XFEM with an aggregation enrichment strategy according to the crack problem GFEM/XFEM;

s2: correcting an enrichment function in an enrichment strategy by adopting a linear Heaviside function;

s3: changing a finite element PU function of the standard FEM, and constructing a novel SGFEM basic function;

s4: processing the SGFEM basic function based on a local principal component analysis technology, eliminating redundancy caused by a plurality of enrichment functions of a single node, and obtaining an SGFEM model;

s5: and solving the crack problem through the SGFEM model to complete the simulation of the crack problem.

More specifically, the step S1 specifically includes:

the region omega is subdivided into quasi-regular triangular or quadrangular grid cells e_sThe dimension parameter of the grid is represented by h, the grid is simple, fixed and connected with the crack_OIndependently, { (x)_i，y_i)：i∈I_hRecording the data as a grid node set; let phi_i，i∈I_hIs a standard finite element function, and

is phi_iThe supporting set of (2); the main idea of GFEM/XFAM is by using the PU function [7,6,10,13 ]]Rendering finite element functions and expressionsSpecial functions of solution properties are coupled to achieve a high-precision approximation effect, and the special functions are called enrichment functions; approximation function u of crack problem GFEM/XFEM_hThe following were used:

wherein I_HNode sets of cells intersecting the fracture, but with the support set removed, nodes containing fracture tips, a_i、b_i、

Denotes the intermediate variable, S_jRepresenting an intermediate function; i is_sThe GFEM/XEM is called GFEM/XEM with geometric enrichment strategy as a node set inside a circle B (O, R) with a radius R and a crack tip as a center; in the formula (2), function

Is a Heaviside function for simulating the discontinuity of the true solution at the crack

Called the crack tip strengthening function, wherein R is the radius of the enrichment range, R, theta are polar coordinates, and the intermediate variables of the mathematical expression; using enrichment functions H and

and completing the construction of GFEM/XFEM.

In practice, due to enrichment functions H and

leads to a linear correlation between them and finite element functions, making the condition ratio of the stiffness matrix much larger than that of the finite elements, which may lead to significant truncation errors when solving linear equations,and may even result in failure of the numerical simulation process. Therefore, the condition number of the traditional GFEM/XFAM is improved by adopting the idea of a stable generalized finite element method, the approximation error is improved, and the improved means comprises the aspects of correcting an enrichment function, changing a PU function and analyzing a local principal component.

More specifically, the step S2 specifically includes:

in the formula (2), a linear Heaviside function is used

Replaces the Heaviside function H of the conventional GFEM/XFEM, and then

And

the correction is as follows:

wherein

Is a finite element interpolation operator, i.e. for a continuous function F,

wherein the parameters v, phi_jIntermediate variables that are mathematical expressions; and constructing the GFEM/XFEM with the set enrichment strategy based on the corrected enrichment function.

In the implementation, it can be found from equation (3) that the enrichment functions subtract their finite element interpolation, which is one of the actions to reduce the linear dependence on the finite element functions [1,9,16,17,18 ].

More specifically, the step S3 specifically includes: the standard FEM function phi in the formula (2)_iAs a function of PU, this is also one of the reasons why the condition number of the stiffness matrix becomes worseSee [11,17 ]]Discussion of (1). Therefore, the PU function is modified to a high-order polynomial PU function, specifically, let:

Q₀＝(1-ξ)²(1+2ξ)，Q₁＝ξ²(3-2ξ)，ξ∈[0，1]

is a one-dimensional reference unit [0,1 ]]The PU function of (1), from which a two-dimensional reference cell [0,1 ] is obtained]×[0，1]PU function of

The method comprises the following specific steps:

in order to obtain any actual unit e_sThe PU function of (1), let (x, y) be F_s(ξ) is a reference cell to e_sAffine transformation of_sA PU function of

Assembling these unit PU functions according to the nodes obtains the required PU function Q_i，i∈I_h(ii) a Based on the modified enrichment function (3) and the new PU function Q_iObtaining a novel SGFEM basic function expression:

equation (4) is a novel SGFEM basic function expression, and the SGFEM can reach an optimal convergence order mathematically proved according to the ideas in [16,18 ].

More specifically, the step S4 specifically includes: in equation (4), it is found that some nodes are enriched with multiple functions, e.g., I_HEach node in (1) is enriched with three functions, I_SEach node in (1) is enriched with four functions, and I_H∩I_HEach node in (a) is enriched with seven functions. Thus using

Represents each enrichment node (x)_i，y_i) The enrichment function of (2) then has:

obviously corresponds to I ∈ I_H，i∈I_s，i∈I_H∩I_sα equals 3, 4,7 respectively, and the linear correlation between enrichment functions for each node is another source of large condition number of stiffness matrix, for which a local principal component analysis is proposed [ 8]]Local Principal Component Analysis (LPCA) technique to eliminate redundancy caused by multiple enrichment functions at a single node, because the elastic mechanics equation is a vector valued equation for any I ∈ I_H∪I_sDenotes E_iThe basis functions for the x-and y-directions are as follows:

for any I ∈ I_H∪I_sMemory for recording

Is a function of

In the inner product

A covariance matrix of; as will be apparent from the above description,

is a α×α submatrix of the total stiffness matrix A, and carries out PCA (principal component analysis) based on a covariance matrix to obtain a α×α matrix

And

wherein:

each column contains the coefficients of one principal component,

consists of the percentage of each main component; according to the PCA property, the new principal component function is:

are orthogonal under inner product (5) and can be based on vectors

Determining contribution of each principal component, and deleting

λ is a predetermined parameter, where λ is 10^-10This is a very small percentage; the convergence performance is not reduced due to the LPCA, the enrichment function knowledge is recombined, the original approximation space is not changed, only the redundancy with extremely small contribution is deleted, and in addition, the LPCA can be executed according to the rigidity matrix corresponding to the SGFEM basic function (4) without a large amount of extra calculation. In this sense, LPCA can be viewed as a preprocessing algorithm to solve linear systems.

The LPCA is regarded as a preprocessing algorithm for solving a linear system, and the specific implementation process is as follows:

based on enrichment function

Integrating the stiffness matrix A and the load vector b for an arbitrary node (x)_i，y_i)，i∈I_H∪I_sThe subscript of the permutation of multiple enrichment functions is J_i，...，J_i+ α -1, and extracting the submatrix from A

Is composed of

For each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

and updating the rigidity matrix A and the load vector b as follows:

for each (x)_i，y_i)，i∈I_H∪I_sIf, if

Updating the stiffness matrix A and the load vector b as follows:

A(J_i+j-1，：)＝0，A(：，J_i+j-1)＝0，

A(J_i+j-1，J_i+j-1)＝1，b(J_i+j-1，：)＝0

for enrichment function

And executing the same process as the above, updating the rigidity matrix A and the load vector b, and obtaining a new rigidity matrix A and a new load vector b after LPCA, thereby obtaining the SGFEM model.

More specifically, in step S5, the SGFEM model is used to solve the crack problem, that is, a linear equation set is solved, specifically:

the system of linear equations is scaled:

let E_ii＝(A_ii)^-1/2Obtaining a diagonal matrix E, and obtaining EAE and Eb after the rigidity matrix A and the load vector b are scaled;

solving a scaled linear system of equations:

EAEu＝Eb (6)

reverse pretreatment:

carrying out inverse scaling on the scaled solution u to obtain u ═ Eu; meanwhile, u is also required to be subjected to inverse pretreatment according to the LPCA pretreatment process:

based on enrichment function

Integrate the solution vector u for any node (x)_i，y_i)，i∈I_H∪I_sSubscript of the arranged multiple enrichment function is J_i，...，J_i+ α -1, and extracting the submatrix from A

Comprises the following steps:

for each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

if it is

Updating

The sum vector u is:

for each (x)_i，y_i)，i∈I_H∪I_sThe updated solution vector u is:

for enrichment function

Executing the same process as the above to obtain a new solution vector u after the inverse LPCA; therefore, the crack problem is solved, and the simulation of the crack problem is completed.

In the specific implementation process, as shown in fig. 3, the process of solving the elastic crack problem by the SGFEM is intuitively expressed, and the method upgrades the general GFEM/XFEM, so that the advantages of the original GFEM/XFEM, such as simple grids, high approximation precision and the like, are maintained, and the condition number of the stiffness matrix is reduced to the same order as that of the conventional FEM, thereby greatly improving the reliability of numerical simulation of the crack problem. The adopted technology mainly comprises a correction enrichment function, a change unit decomposition (PU) function, local principal component analysis and the like.

In the specific implementation process, on the basis of the technical scheme, an assembly algorithm similar to the FEM framework is developed, and on the basis, an SGFEM parallel algorithm is developed, wherein the parallel processing comprises the processes of rigidity matrix assembly, linear equation solution, preprocessing and the like, so that the parallel efficiency of the SGFEM is greatly improved. The program framework is based on the assembly integration idea similar to the FEM, so that the program framework can be conveniently embedded into the current FEM software.

More specifically, the elastic crack problem simulation method based on the stable generalized finite element further comprises the steps of carrying out parallelization processing on the SGFEM model, solving the crack problem by using the parallelized SGFEM model, and completing the simulation of the crack problem.

More specifically, the process of parallelizing the SGFEM model specifically includes: screening out nodes needing to be subjected to LPCA (low power allocation algorithm) of an SGFEM (generalized minimum finite element model) according to existing marks, dispersing the nodes to obtain N units, and dividing the nodes according to the number p of computer processes, wherein N is not required to be divided by the processes p; first p-1 process handling

The p number process processes the remaining N-nel (p-1) units; independently assembling a rigidity matrix A and a load vector b in each process, and finally aggregating the rigidity matrix A and the load vector b into a certain specific process; in LPCA, it is obtained at each node according to A

And updating the matrix A and the vector b, thereby completing the parallelization processing of the SGFEM model.

In the specific implementation process, the SGFEM model is very different from the FEM model, so that a parallelization algorithm needs to be newly designed. Dispersing the solution area to obtain N units, dividing the area according to the number p of computer processes (the number N of units is not required to be divided by the number p of processes): first p-1 process handling

The p number process processes the remaining N-nel (p-1) cells. The stiffness matrix A and the load vector b are independently assembled in each process and finally aggregated in a specific process (such as process No. 0). Note thatIt is to be appreciated that no data communication is required between different processes when assembling the stiffness matrix and the load vector. Because the rigidity matrix generated by each sub-area adopts a sparse storage mode, the data communication volume is very small when aggregation operation is finally carried out, and the time spent is far shorter than the numerical calculation.

In the implementation, for LPCA, only nodes near the fracture need to be analyzed. If all nodes are divided into sessions according to the number of processes, the calculation amount of a part of the processes is small, the calculation amount of a part of the processes is large, and the parallel efficiency is low. Firstly, screening out nodes needing LPCA according to the existing marks, and then dividing the nodes according to the process number by the method. In LPCA, it is available at each node according to A

And updates matrix a and vector b.

More specifically, a parallelized SGFEM model is solved by adopting a parallelized Krylov subspace method, namely a parallelized linear system EAEx is equal to Eb to obtain x, after the inverse scaling of a solution vector x, the crack problem is solved by adopting an inverse LPCA method, and the simulation of the crack problem is completed.

Example 2

More specifically, on the basis of embodiment 1, the method of the present invention is integrated into an algorithm program for application, and a pseudo code of a main program module is obtained:

algorithm 1 Generation of node coordinates

Inputting: square root of total number of cells nel

And (3) outputting: grid node coordinate g _ coord

Start of

fori＝1：nel+1

Generating an ith node coordinate corresponding to the first column;

forj＝1：nel

generating a jth node coordinate corresponding to the row;

end

end

end up

Algorithm 2 generation of element index finite element part degrees of freedom

Inputting: square root of total number of cells nel

And (3) outputting: element finite element partial degree of freedom g _ num

Start of

fori＝1：nel

Generating the ith unit degree of freedom corresponding to the first row;

end

fori＝1：nel

forj＝2：nel

generating the jth unit degree of freedom in the ith column;

end

end

end up

Algorithm 3 marking corresponding enrichment function type and degree of freedom of node

Inputting: square root of total number of cells nel, Heaviside function type l1, Singular function type l2

And (3) outputting: node and number nf corresponding to enrichment function, total degree of freedom dof, unit enrichment type ef, and intersection point of crack and grid

Coordinates pz1, pz2

Start of

fors＝1：nel

if crack penetration unit

Marking unit nodes as enrichment nodes;

recording the vertical coordinate of the intersection point;

judging the edge of the intersection of the crack and the unit;

recording the intersection points;

end

end

fori ═ 1: nn% total number of nodes

if node enrichment type is Heaviside

Determining a new added degree of freedom;

recording the corresponding number of the node;

end

end

fori ═ 1: nn% total number of nodes

if the node enrichment type is Singular

Determining a new added degree of freedom;

recording the corresponding number of the node;

end

end

updating the total degree of freedom;

end up

Algorithm 4 Assembly stiffness matrix

Start of

Foriel ═ 1: nels% total number of units

if unit is enriched

Generating an integral node and weight of a 12 × 12 Gaussian integral formula

else

Dividing the two smaller blocks into at most 4 triangles;

using a 16 point gaussian integration formula in each triangle;

generating integral nodes and weights;

end

for i ═ 1: NG% number of integration nodes

Calculating a cap function and a cap function derivative of the reference unit;

calculating a reference unit PU function and a PU function derivative;

mapping to an actual cell;

calculating a derivative of the basis function;

assembling a local stiffness matrix;

end

assembling the local stiffness matrix into a total stiffness matrix;

end

end

end up

Algorithm 5LPCA

Start of

fori ═ 1: nn% Total degree of freedom

Judging the node enrichment type as follows: heaviside, singulator, Heaviside and singulator;

extracting a sub-stiffness matrix corresponding to the node;

deriving from sub-matrices

(each column contains coefficients of one principal component) and

(consisting of a percentage of the total variance associated with the principal component);

updating a rigidity matrix A and a load vector b;

judgment of

Updating the diagonal element j of the stiffness matrix to be 1, and the other corresponding row and column elements to be 0;

updating the load vector element j to 0;

end

end up

Algorithm 6 linear equation set preprocessing and solving

Start of

Solving a solution vector u formula (6) according to a diagonal matrix E scaling linear equation set;

updating u according to the diagonal matrix E inverse scaling solution vector;

fori ═ 1: nn% Total degree of freedom

Judging the node enrichment type as follows: heaviside, singulator, Heaviside and singulator;

extracting a sub-stiffness matrix corresponding to the node;

deriving from sub-matrices

(each column contains coefficients of one principal component) and

(consisting of a percentage of the total variance associated with the principal component);

judgment of

Updating

The diagonal element j is 1, and the other corresponding row and column elements are 0;

updating the solution vector element j to 0;

updating the solution vector u according to the 4.2.4 section algorithm;

end

end up

Example 3

More specifically, consider a region Ω ═ 0,1 with a crack]²Crack of_oGiven by the equation ax + by + c ═ 0, where

a is-1-b, c is 1, and the crack apex is at O coordinate (0.5 ), as shown in fig. 4.

Opening the model with type I:

the effect of LPCA-based SGFEM numerical simulations was tested as a true solution to equation (1). Wherein, Young modulus E is 90, Poisson ratio v is 0.28, pressure intensity factor M_I＝1。

Using a grid size of

Uniform n × n quadrilateral mesh discretization area omega, bilinear FE function [ phi ] based on the mesh_i，i∈I_hFinite element parts for XFEM/GFEM and SGFEM. FIG. 4 shows that when n is 17, I_HAnd I_SIn the case of an enriched node, the nodes in sphere B (O, R) having a radius R of 0.2 are enriched by a singular function. In addition to the newly developed LPCA-based SGFEM (4), a conventional XFEM/GFEM (2) will be shown for comparison. Parameter λ 10 in LPCA^-10Relative to H¹The Scale Condition Number (SCN) of the error and stiffness matrix is used as an index for comparison, and is specifically expressed as:

after the SGFEM (4) carries out LPCA, the scale condition number of the rigidity matrix is reduced to be in the same order as that of a finite element method, and the approximation precision is improved compared with that of the original GFEM/XFAM. For any symmetric matrix M, let D be a diagonal matrix, where

D_ii＝(M_ii)^-1/2

Then the Scaled Condition Number (SCN) of M is defined as

Wherein

Is that

The ratio of the largest feature root and the smallest non-zero feature root.

Grid parameters

H of corresponding XFEM/GFEM and SGFEM¹The error and SCN are shown in fig. 5 and 6. It can be seen in fig. 5 that both XFEM/GFEM and SGFEM can achieve the best convergence speed, but the error of the latter is smaller than the former. On the other hand, as is evident in the right diagram of fig. 6, SCN of SGFEM has the same order O (h) as FEM^-2) The SCN of XFAM/GFEM is much larger, and when h is smaller, it is approximately O (h)^-4) This growth order is a much larger SCN growth order than SGFEM.

Numerical experiments show that the SGFEM based on LPCA newly developed in the text achieves the best convergence O (h) and SCNO (h) of the same order as FEM^-2) This is not available with conventional XFAM/GFEM.

Next, the execution efficiency of the SGFEM parallel algorithm was tested. The discrete grid size is N129 × 129. Table 1 shows good parallel acceleration ratios and efficiencies. Note that efficiency will be greater than 1, indicating that the efficiency of the parallel inner serial algorithm has yet to be optimized. FIG. 7 illustrates the computational time required for different numbers of processors with a fixed problem-solving scale. The black solid line is an ideal case, and it can be seen that the rigidity matrix is assembled and the total time is similar to the ideal case. The program is explained to be very extensive, and the problem of fixed size can be solved more quickly when processors are increased. The LPCA acceleration is not ideal when the number of processors is large, mainly because the computation time is small and the inter-process communication time is proportional to the increase. Fig. 8 shows the time required to solve the problem of different sizes when the number of processors is fixed (p-40). As can be seen from the figure, the assembly stiffness matrix, LPCA and total time are all approximately linear, indicating good parallel efficiency.

TABLE 1

It should be understood that the above-described embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. Other variations and modifications will be apparent to persons skilled in the art in light of the above description. And are neither required nor exhaustive of all embodiments. Any modification, equivalent replacement, and improvement made within the spirit and principle of the present invention should be included in the protection scope of the claims of the present invention.

[1]I.Babuska and U.Banerjee,Stable generalized finite element method,Comput.Methods Appl.Mech.Engrg.,(201-204):91-111,2011.

[2]I.Babuska,U.Banerjee,and J.Osborn,Survey of meshless andgeneralized finite element methods:a unified approach.Acta Numerica,12:1-125,2003.

[3]I.Babuska and J.M.Melenk,The partition of unity finite elementmethod,Int.J.Numer.Meth.Engng.,40:727-758,1997.

[4]E.B′echet,H.Minnebo,N.Mo¨es,and B.Burgardt,Improved implementationand robustness study of the X-FEM method for stress analysis around cracks,Int.J.Numer.Meth.Engng.,64:1033-1056,2005.

[5]T.Belytschko and T.Black,Elastic crack growth in finite elementswith minimal remeshing,Int.J.Numer.

Meth.Engng.,45:601-620,1999.

[6]T.P.Fries,A corrected XFEM approximation without problems inblending elements,Int.J.Numer.Meth.

Engng.,75:503-532,2008.

[7]T.P.Fries and T.Belytschko,The extended/generalized finite elementmethod:An overview of the method and its applications,Int.J.Numer.Meth.Engng.,84:253-304,2010.

[8]I.T.Jolli e,Principal Component Analysis,Second Edition,Springer-Verlag,New York,2002.

[9]K.Kergrene,I.Babuska and U.Banerjee,Stable Generalized FiniteElement Method and associated iterative schemes:application to interfaceproblems,Comput.Methods Appl.Mech.Engrg.,305:1-36,2016.

[10]P.Laborde,J.Pommier,Y.Renard,and M.Sala¨un,High order extendedfinite element method for cracked domains,Int.J.Numer.Meth.Engng.,64:354-381,2005.

[11]H.Li,A note on the conditioning of a class of generalized finiteelement methods,Applied Numerical Mathematics,62:754-766,2012.

[12]J.M.Melenk and I.Babuska,The partition of unity finite elementmethod:Theory and application,Comput.MethodsAppl.Mech.Engrg.,139:289-314,1996.

[13]N.Mo¨es,J.Dolbow,and T.Belytschko,A finite element method forcrack without remeshing,Int.J.Numer.Meth.Engng.,46:131-150,1999.

[14]A.G.Sanchez-Rivadeneira,C.A.Duarte,A stable generalized/extendedFEM with discontinuous interpolants for fracture mechanics,Comput.MethodsAppl.Mech.Engrg.,345:876-918,2019.

[15]M.A.Schweitzer,Stable enrichment and local preconditioning in theparticle-partition of unity method,NumerischeMathematik,118:137-170,2011.

[16]Q.Zhang,I.Babuˇska and U.Banerjee,Robustness in stablegeneralized finite element methods(SGFEM)applied to Poisson problems withcrack singularities,Comput.Methods Appl.Mech.Engrg.,311:476-502,2016.

[17]Q.Zhang,U.Banerjee,and I.Babuska,High order stable generalizedfinite element methods,NumerischeMathematik,128:1-29,2014.

[18]Q.Zhang,U.Banerjee,and I.Babuska,Strongly Stable GeneralizedFinite Element Method(SSGFEM)for a non-smooth interface problem,Comput.Methods Appl.Mech.Engrg.,344:538-568,2019.

Claims (9)

1. The elastic crack problem simulation method based on the stable generalized finite element is characterized by comprising the following steps of:

s1: constructing a GFEM/XFEM with an aggregation enrichment strategy according to the crack problem GFEM/XFEM;

s2: correcting an enrichment function in an enrichment strategy by adopting a linear Heaviside function;

s3: changing a finite element PU function of the standard FEM, and constructing a novel SGFEM basic function;

s4: processing the SGFEM basic function based on a local principal component analysis technology, eliminating redundancy caused by a plurality of enrichment functions of a single node, and obtaining an SGFEM model;

s5: and solving the crack problem through the SGFEM model to complete the simulation of the crack problem.

2. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 1, wherein the step S1 is specifically performed by:

the region omega is subdivided into quasi-regular triangular or quadrangular grid cells e_sThe dimension parameter of the grid is represented by h, the grid is simple, fixed and connected with the crack_OIndependently, { (x)_i，y_i)：i∈I_hRecording the data as a grid node set; let phi_i，i∈I_hIs a standard finite element function, and

is phi_iThe supporting set of (2); the main idea of GFEM/XFEM is that a finite element function and special functions expressing true solution properties are coupled by using a PU function to achieve a high-precision approximation effect, and the special functions are called enrichment functions; approximation function u of crack problem GFEM/XFEM_hThe following were used:

wherein I_HNode sets of cells intersecting the fracture, but with the support set removed, nodes containing fracture tips, a_i、b_i、

Denotes the intermediate variable, S_jRepresenting an intermediate function; i is_sThe GFEM/XEM is called GFEM/XEM with geometric enrichment strategy as a node set inside a circle B (O, R) with a radius R and a crack tip as a center; in the formula (2), function

Is a Heaviside function for simulating the discontinuity of the true solution at the crack

Called the crack tip strengthening function, wherein R is the radius of the enrichment range, R, theta are polar coordinates, and the intermediate variables of the mathematical expression; using enrichment functions H and

and completing the construction of GFEM/XFEM.

3. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 2, wherein the step S2 is specifically performed by:

in the formula (2), a linear Heaviside function is used

Replaces the Heaviside function H of the conventional GFEM/XFEM, and then

And

the correction is as follows:

wherein

Is a finite element interpolation operator, i.e. for a continuous function F,

wherein the parameters v, phi_jIntermediate variables that are mathematical expressions; and constructing the GFEM/XFEM with the set enrichment strategy based on the corrected enrichment function.

4. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 3, wherein the step S3 is specifically as follows: modifying the PU function in equation (2) to a high-order polynomial PU function, specifically, letting:

Q₀＝(1-ξ)²(1+2ξ)，Q₁＝ξ²(3-2ξ)，ξ∈[0，1]

is a one-dimensional reference unit [0,1 ]]The PU function of (1), from which a two-dimensional reference cell [0,1 ] is obtained]×[0，1]PU function of

The method comprises the following specific steps:

in order to obtain any actual unit e_sThe PU function of (1), let (x, y) be F_s(ξ) is a reference cell to e_sAffine transformation of_sA PU function of

Assembling these unit PU functions according to the nodes obtains the required PU function Q_i，i∈I_h(ii) a Based on the modified enrichment function (3) and the new PU function Q_iObtaining:

the formula (4) is a novel SGFEM basic function expression.

5. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 4, wherein the step S4 is specifically as follows: in equation (4), it is found that some nodes are enriched with multiple functions, and therefore, are used

Represents each enrichment node (x)_i，y_i) The enrichment function of (2) then has:

obviously corresponds to I ∈ I_H，i∈I_s，i∈I_H∩I_sα equals 3, 4,7, respectively, multiple enriches per nodeLinear correlation between set functions is another source of large condition number of the rigidity matrix, and for this reason, a Local Principal Component Analysis (LPCA) technology is provided to eliminate redundancy caused by a plurality of enrichment functions of a single node, because an elastic mechanical equation is a vector value equation and is used for any I ∈ I_H∪I_sDenotes E_iThe basis functions for the x-and y-directions are as follows:

for any I ∈ I_H∪I_sMemory for recording

Is a function of

In the inner product

A covariance matrix of; as will be apparent from the above description,

is a α×α submatrix of the total stiffness matrix A, and carries out PCA (principal component analysis) based on a covariance matrix to obtain a α×α matrix

And

wherein:

each column contains the coefficients of one principal component,

consists of the percentage of each main component; according to the PCA property, the new principal component function is:

are orthogonal under inner product (5) and can be based on vectors

Determining contribution of each principal component, and deleting

λ is a predetermined parameter, where λ is 10^-10(ii) a The LPCA is regarded as a preprocessing algorithm for solving a linear system, and the specific implementation process is as follows:

based on enrichment function

Integrating the stiffness matrix A and the load vector b for an arbitrary node (x)_i，y_i)，i∈I_H∪I_sThe subscript of the permutation of multiple enrichment functions is J_i，...，J_i+ α -1, and extracting the submatrix from A

Is composed of

For each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

and updating the rigidity matrix A and the load vector b as follows:

for each (x)_i，y_i)，i∈I_H∪I_sIf, if

Updating the stiffness matrix A and the load vector b as follows:

A(J_i+j-1，：)＝0，A(：，J_i+j-1)＝0，

A(J_i+j-1，J_i+j-1)＝1，b(J_i+j-1，：)＝0

for enrichment function

And executing the same process as the above, updating the rigidity matrix A and the load vector b, and obtaining a new rigidity matrix A and a new load vector b after LPCA, thereby obtaining the SGFEM model.

6. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 5, wherein the fracture problem is solved through the SGFEM model in step S5, that is, a linear equation system is solved, specifically:

the system of linear equations is scaled:

let E_ii＝(A_ii)^-1/2To obtain oneScaling a diagonal matrix E, a rigidity matrix A and a load vector b to obtain EAE and Eb;

solving a scaled linear system of equations:

EAEu＝Eb (6)

reverse pretreatment:

carrying out inverse scaling on the scaled solution u to obtain u ═ Eu; meanwhile, u is also required to be subjected to inverse pretreatment according to the LPCA pretreatment process:

based on enrichment function

Integrate the solution vector u for any node (x)_i，y_i)，i∈I_H∪I_sSubscript of the arranged multiple enrichment function is J_i，...，J_i+ α -1, and extracting the submatrix from A

Comprises the following steps:

for each (x)_i，y_i)，i∈I_H∪I_sBased on

Performing PCA acquisition

And

if it is

Updating

The sum vector u is:

u(J_i+j-1，：)＝0；

for each (x)_i，y_i)，i∈I_H∪I_sThe updated solution vector u is:

for enrichment function

Executing the same process as the above to obtain a new solution vector u after the inverse LPCA; therefore, the crack problem is solved, and the simulation of the crack problem is completed.

7. The elastic crack problem simulation method based on stable generalized finite elements according to any one of claims 1-6, further comprising parallelizing the SGFEM model, and solving the crack problem by using the parallelized SGFEM model to complete the simulation of the crack problem.

8. The elastic fracture problem simulation method based on stable generalized finite element according to claim 7, wherein the process of parallelizing the SGFEM model specifically comprises: screening out nodes needing to be subjected to LPCA (low power allocation algorithm) of an SGFEM (generalized minimum finite element model) according to existing marks, dispersing the nodes to obtain N units, and dividing the nodes according to the number p of computer processes, wherein N is not required to be divided by the processes p; first p-1 process handling

The p number process processes the remaining N-nel (p-1) units; at each timeIndependently assembling a rigidity matrix A and a load vector b in each process, and finally gathering the rigidity matrix A and the load vector b into a certain specific process; in LPCA, it is obtained at each node according to A

And updating the matrix A and the vector b, thereby completing the parallelization processing of the SGFEM model.

9. The elastic fracture problem simulation method based on stable generalized finite elements according to claim 7, wherein a parallelized Krylov subspace method is adopted to solve the SGFEM model, that is, a parallelized linear system EAEx ═ Eb obtains x, and after the inverse scaling of the solution vector x, a method similar to LPCA is adopted to solve the fracture problem, thereby completing the simulation of the fracture problem.

Images