Novel finite element optimization solution algorithm

Abstract

The invention discloses a novel finite element optimization solution algorithm, which comprises the following steps of: preprocessing a real image scanned by CT (Computed Tomography) to obtain an image gray value; reading an image gray value to obtain an image finite element model; establishing a mapping relationship between the unit node number and the overall degree-of-freedom number; constructing a Jacobian matrix and a geometric matrix according to the unit node numbers; calculating an element stiffness matrix and a stress function matrix corresponding to each material attribute; calculating an equivalent node force according to the element stiffness matrix, and assembling a node force vector; solving node displacement under a given boundary condition; according to the displacement field, stress components corresponding to all units under the given boundary condition are calculated and averaged, and finally the homogenized constitutive tensor is obtained. The invention provides a novel finite element optimization solution algorithm, a traditional finite element preprocessing process is simplified, the calculation efficiency is remarkably improved, and calculation resources are greatly saved.

Description

Technical Field

[0001] The present invention relates to the field of numerical calculations of heterogeneous materials, and particularly to a new finite element optimization solution algorithm. Background Art

[0002] Finite element analysis is one of the most commonly used numerical simulation methods. By discretizing the original problem into a corresponding algebraic system, an approximate solution of the partial differential equation can be obtained. However, in large-scale computing applications, traditional finite element analysis relies on auxiliary data structures such as storing the global stiffness matrix and node connectivity. Especially when dealing with complex-shaped domains and unstructured grids, the memory requirements increase significantly. For materials with complex microstructures, traditional finite element analysis is difficult to capture complex micro details and cannot accurately describe the behavior of the materials. In order to obtain reliable effective properties, it is often necessary to continuously refine the grid to improve the accuracy of the solution, which in turn leads to an exponential increase in the stiffness matrix, consuming not only excessive computing time but also a sharp rise in memory requirements. Computational resources have become one of the important factors restricting the rapid development of traditional finite element solution techniques. Summary of the Invention

[0003] In view of the above-mentioned defects of the prior art, the technical problem to be solved by the present invention is that traditional finite element analysis has high storage requirements and relies on the support of computing resources. The present invention provides a new finite element optimization solution algorithm, which simplifies the traditional finite element preprocessing process for two-dimensional and three-dimensional image data, significantly improves the computing efficiency, and greatly saves computing resources.

[0004] To achieve the above object, the present invention provides a new finite element optimization solution algorithm, including the following steps:

[0005] Step 1: Preprocess the real image scanned by CT, represent it with a series of gray values, and obtain the image gray values;

[0006] Step 2: Read the image gray values and assign corresponding material numbers to each pixel unit in the model to obtain an image finite element model;

[0007] Step 3: Under periodic boundary conditions, establish the mapping relationship between the pixel unit node numbers and the overall degree-of-freedom numbers;

[0008] Step 4: Calculate the four-node rectangular element shape functions and their partial derivatives according to the pixel unit node numbers obtained in Step 3, and construct the Jacobian matrix J and the geometric matrix B;

[0009] Step 5: Calculate the element stiffness matrix K corresponding to each material property according to the material parameters e and the stress function matrix S;

[0010] Step 6: Apply boundary displacement loads, according to the element stiffness matrix K eCalculate the equivalent nodal forces and assemble the right - hand side vector f;

[0011] Step 7: Solve the displacement field. Transform the calculation of the linear equation system matrix - vector product into the convolution process of a convolutional neural network to solve the nodal displacements under the given boundary conditions;

[0012] Step 8: Calculate the stress components corresponding to each element under the given boundary conditions according to the displacement field;

[0013] Step 9: Average the stress components of all elements to obtain the homogenized constitutive tensor.

[0014] Furthermore, pre - process the real image scanned by CT, which is represented by a series of gray values to obtain the image gray values. Specifically, it includes pre - processing operations such as segmentation and denoising on the real image scanned by CT, dividing different material regions into discrete material categories and storing them in an 8 - bit array, and representing them with a series of gray values.

[0015] Furthermore, assign gray values to the heterogeneous materials, where the gray value of the model matrix filler is 0 and the gray value of the cylindrical fiber is 255.

[0016] Furthermore, in Step 2: Read the image gray values and assign corresponding material numbers to each pixel unit in the model to obtain the image finite - element model. Specifically, it includes reading the image gray values and converting them into the material numbers corresponding to each unit, and constructing the mapping relationship between the unit and the corresponding material properties for each pixel unit to obtain the image finite - element model.

[0017] Furthermore, in Step 3: Under the periodic boundary conditions, establish the mapping relationship between the pixel unit node numbers and the global degree - of - freedom numbers. Specifically, according to the pixel unit grid structure, calculate the degree - of - freedom index corresponding to each node under the periodic boundary conditions according to the spatial position relationship between the nodes, and establish the mapping relationship from the pixel unit node numbers to the global degree - of - freedom numbers.

[0018] Furthermore, in Step 4: Calculate the four - node rectangular element shape functions and their partial derivatives according to the pixel unit node numbers obtained in Step 3, and construct the Jacobian matrix J and the geometric matrix B. Specifically, use the natural coordinate system (ξ, η) to map the actual element to the standard element, and calculate the shape function N in the natural coordinate system according to Equation (1) i , and construct the Jacobian matrix J and the geometric matrix B through the shape function and its partial derivatives, as shown in Equations (2) and (3):

[0019]

[0020] where ξ and η represent the coordinates of the natural coordinate system.

[0021]

[0022] where x and y are coordinates in the physical coordinate system; ξ and η represent coordinates in the natural coordinate system; N i is the element shape function; J -1 is the inverse matrix of the Jacobian matrix.

[0023] Furthermore, in Step Five, calculate the element stiffness matrix K e and the stress function matrix S corresponding to each material property according to the material parameters; specifically as follows:

[0024] Calculate the element stiffness matrix K of each material at the Gauss integration points according to the element node numbers e :

[0025]

[0026] where B is the geometric matrix; B T is the transpose of the geometric matrix; |J| is the Jacobian determinant; t represents the element thickness, and in two-dimensional problems, t = 1.

[0027] Furthermore, calculate the stress function matrix S according to the element node numbers and Equation (3):

[0028] S = CB (5)

[0029] where C is the constitutive matrix.

[0030] Furthermore, in Step Six, apply the boundary displacement load and calculate the equivalent nodal forces according to the element stiffness matrix K e to assemble the right-hand side vector f; specifically including applying the displacement load at the boundary. In two-dimensional problems, consider the cases in the x-direction, y-direction, and xy-direction. By traversing the specific boundary elements and using the element stiffness matrix K e to convert it into equivalent nodal forces. Each node contains two degrees of freedom in the x and y directions. Through the degree-of-freedom mapping, the results calculated from the node numbers are assembled into the global right-hand vector f.

[0031] Furthermore, in Step Seven, solve the displacement field and transform the calculation of the linear equation matrix-vector product into the convolution process of a convolutional neural network to solve the nodal displacements under the given boundary conditions; specifically including using the node-based Node-by-Node matrix-free method to directly calculate the stiffness contributions at the node level, assembling the element stiffness matrix according to the assembly principle, setting the local stiffness matrix K batch consisting of m nodes as a minibatch, calculating the linear equation matrix-vector product, and finally transforming it into the convolution process of a convolutional neural network to solve the nodal displacement components under the given boundary conditions.

[0032] Technical effects

[0033] A new finite element optimization solution algorithm provided by the present invention directly generates a pixel-based grid using a high-resolution image, eliminating the cumbersome geometric modeling and mesh generation steps, and transforming the calculation of the matrix-vector product of a linear equation system into a convolution process of a convolutional neural network for solution. The solution process is based on parallel computing accelerated by computing hardware adapted to the Pytorch package, assembling nodes batch by batch without the need to assemble the global stiffness matrix, and repeatedly calling according to rules using the characteristics of the pixel finite element structured grid, effectively saving computing resources and improving computing efficiency.

[0034] The concept, specific structure and technical effects of the present invention will be further described below in conjunction with the accompanying drawings to fully understand the purpose, features and effects of the present invention. Description of the Drawings

[0035] Figure 1 is a schematic flow chart of a new finite element optimization solution algorithm according to a preferred embodiment of the present invention;

[0036] Figure 2 is a schematic diagram of the distribution of the element stiffness matrix of a new finite element optimization solution algorithm according to a preferred embodiment of the present invention (K in step five e );

[0037] Figure 3 is a schematic diagram of the node numbering of a new finite element optimization solution algorithm according to a preferred embodiment of the present invention (step three);

[0038] Figure 4 is a schematic diagram of collecting the element contributions to the global degrees of freedom by the Node-by-Node method of a new finite element optimization solution algorithm according to a preferred embodiment of the present invention (step seven). Detailed Embodiment

[0039] In order to make the technical problems, technical solutions and beneficial effects to be solved by the present invention clearer, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0040] In the following description, specific details such as specific internal programs and technologies are proposed for the purpose of illustration rather than limitation, so as to thoroughly understand the embodiments of the present invention. However, those skilled in the art should clearly understand that the present invention can also be implemented in other embodiments without these specific details. In other cases, the detailed descriptions of well-known systems, devices, circuits and methods are omitted to avoid unnecessary details from hindering the description of the present invention.

[0041] As Figure 1As shown in the figure, the present invention provides a new finite element optimization solution algorithm, which is a matrix-free finite element solution method for heterogeneous materials convolutional neural network based on images, and includes the following steps:

[0042] Step 1: Preprocess the real image scanned by CT, represent it with a series of gray values, and obtain the image gray values; specifically include: perform segmentation and denoising preprocessing operations on the real image scanned by CT, divide different material regions into discrete material categories and store them in an 8-bit array, and represent them with a series of gray values.

[0043] Assign gray values to heterogeneous materials, where the gray value of the model matrix filler is 0, and the gray value of the cylindrical fiber is 255. Use the Matplotlib library of Python to process the image, and finally store the preprocessed image in the RAW file format.

[0044] Step 2: Read the image gray values and assign corresponding material numbers to each pixel unit in the model to obtain an image finite element model, specifically including reading the image gray values and converting them into material numbers corresponding to each unit, and constructing a mapping relationship between the unit and the corresponding material properties for each pixel unit according to the material numbers to obtain an image finite element model. Among them, the number of pixel grids in the x direction is 100, the number of pixel grids in the y direction is 100, the total number of grids is 10000, the number of grid nodes is 10201, the degree of freedom value is 20000, the pixel gray values are 0 and 255, the number of material types is 2, and the material properties are shown in Table 1.

[0045] Table 1 Material parameters

[0046]

[0047] Step 3: Under periodic boundary conditions, establish a mapping relationship between the node numbers of pixel units and the global degree of freedom numbers, specifically including calculating the corresponding degree of freedom index of each node under periodic boundary conditions according to the spatial position relationship between nodes based on the pixel unit grid structure, and establishing a mapping relationship from the unit node numbers to the global degree of freedom numbers.

[0048] As Figure 3 shown, number the nodes in the pixel unit in the order from top to bottom and from left to right. Each node has degrees of freedom in both the x and y directions. For an x×y two-dimensional rectangular grid, there are (x + 1)×(y + 1) nodes. Considering the position of the nodes in the grid, convert the consecutive node numbers into the corresponding unique consecutive degree of freedom numbers.

[0049] Step 4: Calculate the shape functions and their partial derivatives of the four-node rectangular element based on the pixel unit node numbers obtained in Step 3, and construct the Jacobian matrix J and the geometric matrix B. Specifically, map the actual element to the standard element using the natural coordinate system (ξ, η), and calculate the shape functions N in the natural coordinate system according to Equation (1). i Construct the Jacobian matrix J and the geometric matrix B through the shape functions and their partial derivatives, as shown in Equations (2) and (3).

[0050]

[0051] where ξ and η represent the coordinates of the natural coordinate system; N i is the element shape function.

[0052]

[0053] where x and y are the coordinates of the physical coordinate system; ξ and η represent the coordinates of the natural coordinate system; N i is the element shape function; J is the Jacobian matrix; J -1 is the inverse of the Jacobian matrix.

[0054] Step 5: Calculate the element stiffness matrix K corresponding to each material property according to the material parameters e and the stress function matrix S; as Figure 2 shown, use the Jacobian matrix to perform coordinate transformation according to Equation (2), adopt 2×2 Gauss integration points on the standard element for numerical integration, and accumulate and calculate at the Gauss points to obtain the element stiffness matrix K e and the stress function matrix S, as shown in Equations (4) and (5).

[0055]

[0056] S = CB (5)

[0057] where B is the geometric matrix; B T is the transpose of the geometric matrix; |J| is the Jacobian determinant; t represents the element thickness, and in two-dimensional problems, t = 1; C is the constitutive matrix.

[0058] Step 6: Apply the boundary displacement load, calculate the equivalent nodal forces according to the element stiffness matrix K e and assemble the right-hand side vector f; apply the displacement load Δ at the boundary. In two-dimensional problems, consider the cases in the x direction, y direction, and xy direction. By traversing the specific boundary elements and converting them into equivalent nodal forces through the element stiffness matrix K e as shown in Equation (6). Each node contains two degrees of freedom in the x and y directions. Assemble the results calculated from the node numbers into the global right-hand vector f through the degree-of-freedom mapping i in.

[0059] K ij Δ = f i (6)

[0060] Among them, K ij is the form of the degree - of - freedom component of the stiffness matrix corresponding to the boundary element; Δ is the displacement load vector; f i is the form of the degree - of - freedom component of the nodal force vector.

[0061] Step Seven: Solve the displacement field, transform the calculation of the matrix - vector product of the linear equations into the convolution process of the convolutional neural network to solve the nodal displacements under the given boundary conditions. As Figure 4 shown, in the structured periodic grid, each node is surrounded by 4 elements. Utilizing this property, a node - by - node matrix - free method based on nodes is used to directly calculate the contribution of the stiffness matrix at the node level, avoiding the construction and storage of the global stiffness matrix.

[0062] First, traverse m nodes in the order of node numbers, set them as a batch, collect the element numbers around the nodes and call the corresponding element stiffness matrices, assemble the element stiffness matrices according to the assembly principle, and set the local stiffness matrix K batch of the m nodes as a minibatch:

[0063]

[0064] Among them, K batch is the stiffness matrix of a minibatch size, is the element stiffness matrix corresponding to the material type, Θ is the assembly operator, and E is the number of pixel elements under each minibatch.

[0065] Construct the local Jacobi preconditioner M according to the size of the minibatch to preprocess the stiffness matrix:

[0066] M = diag{K 11 , K 22 … K nn} (8)

[0067] Among them, diag{} is the diagonal matrix.

[0068] Convert the original solution equation (9) to (10):

[0069] K batch u i = f i (9)

[0070] M -1 K batch u i = M-1 f i (10)

[0071] Among them, M -1 is the inverse matrix of the local Jacobi preconditioner, u i is the nodal displacement vector under each minibatch, and f i is the nodal force vector under each minibatch.

[0072] Secondly, set the displacement field to be solved as a 1-row and n-column random initial convolution kernel u 0 , set the local stiffness matrix K bat composed of m nodes as the input tensor, set the network structure as a 1-layer convolutional layer, set the target loss value to 10 -7 , use the convolution kernel to convolve K bat Keep iterating until the target loss value is reached to stop training, and then proceed to the assembly of the local stiffness matrix K batch for the next batch:

[0073]

[0074] Among them, represents the convolution process, is the force vector in the iterative form corresponding to K batc , and u i is the displacement vector in the iterative form of the convolutional neural network.

[0075] Step Eight: Under the boundary displacement loading condition, calculate the stress components σ x , σ y , τ xy , including the stress generated by the boundary displacement and the internal response. Combine the stress function matrix S of all elements and the actual displacement field to calculate the complete stress response:

[0076] σ = Su (12)

[0077] Step Nine: Traverse all elements to accumulate the stress σ x , σ y , τ xy generated by the boundary displacement and the internal response and average them for all elements to obtain the corresponding column of the homogenized constitutive tensor. The test results are shown in Table 2:

[0078]

[0079] Among them, is the stress accumulation value under the corresponding boundary displacement condition, and E is the total number of pixel units.

[0080] Table 2 Solving Tests with Different Pixel Sizes

[0081]

[0082] Table 2 selects two different pixel sizes, 20×20 and 100×100, for testing. The gray value of the matrix filler of the model is 0, and the gray value of the cylindrical fiber is 255, to verify the accuracy comparison between the convolutional neural network matrix-free finite element solution method and the direct method under different pixel sizes.

[0083] The preferred specific embodiments of the present invention have been described in detail above. It should be understood that those of ordinary skill in the art can make many modifications and variations based on the concept of the present invention without creative work. Therefore, all technical solutions that can be obtained by those skilled in the art in this technical field based on the concept of the present invention through logical analysis, reasoning, or limited experiments on the basis of the prior art should fall within the protection scope determined by the claims.

Claims

1. A new finite element optimization solution algorithm, characterized in that, Including the following steps: Step 1: Preprocess the real image scanned by CT, represent it with a series of gray values, and obtain the image gray values; Step 2: Read the image gray values and assign corresponding material numbers to each pixel unit in the model to obtain the image finite element model; Step 3: Under periodic boundary conditions, establish the mapping relationship between the pixel unit node numbers and the global degree-of-freedom numbers; Step 4: Calculate the four-node rectangular element shape functions and their partial derivatives according to the pixel unit node numbers obtained in Step 3, and construct the Jacobian matrix J and the geometric matrix B; Step 5. Calculate the element stiffness matrix K e corresponding to each material property and the stress function matrix S according to the material parameters; Step 6: Apply the boundary displacement load, and calculate the equivalent nodal forces according to the element stiffness matrix K e to assemble the right-hand side vector f; Step 7: Solve the displacement field, and transform the matrix-vector product calculation of the linear equations into the convolution process of the convolutional neural network to solve the node displacements under the given boundary conditions; Step 8: Calculate the stress components corresponding to each element under the given boundary conditions according to the displacement field; Step 9: Average the stress components of all elements to obtain the homogenized constitutive tensor.

2. The new finite element optimization solution algorithm according to claim 1, characterized in that, Preprocess the real image scanned by CT, represent it with a series of gray values, and obtain the image gray values. Specifically, it includes performing segmentation and denoising preprocessing operations on the real image scanned by CT, dividing different material regions into discrete material categories and storing them in an 8-bit array, and representing them with a series of gray values.

3. The new finite element optimization solution algorithm according to claim 2, characterized in that, Assign gray values to the heterogeneous materials, where the gray value of the model matrix filler is 0 and the gray value of the cylindrical fiber is 255.

4. A new finite element optimization solution algorithm as claimed in claim 1, characterized in that, Step 2: Read the image gray values and assign corresponding material numbers to each unit in the model to obtain the image finite element model. Specifically, it includes reading the image gray values and converting them into the material numbers corresponding to each unit, and constructing the mapping relationship between the unit and the corresponding material properties for each pixel unit to obtain the image finite element model.

5. The novel finite element optimization solution algorithm according to claim 1, wherein Step 3: Under periodic boundary conditions, establish the mapping relationship between the pixel unit node numbers and the global degree-of-freedom numbers. Specifically, it includes calculating the degree-of-freedom index corresponding to each node under the periodic boundary conditions according to the pixel unit grid structure and the spatial position relationship between the nodes, and establishing the mapping relationship from the pixel unit node numbers to the global degree-of-freedom numbers.

6. The novel finite element optimization solution algorithm according to claim 1, characterized in that, Step 4: Calculate the shape functions and their partial derivatives of the four-node rectangular element based on the pixel unit node numbers obtained in Step 3, and construct the Jacobian matrix J and the geometric matrix B. Specifically, map the actual element to the standard element using the natural coordinate system (ξ, η), and calculate the shape function N in the natural coordinate system according to Equation (1). i , and construct the Jacobian matrix J and the geometric matrix B through the shape functions and their partial derivatives, as shown in Equations (2) and (3): Among them, ξ and η represent the coordinates of the natural coordinate system. where x and y are coordinates in the physical coordinate system; ξ and η represent coordinates in the natural coordinate system; N i is the element shape function; J -1 is the inverse of the Jacobian matrix.

7. The novel finite element optimization solution algorithm according to claim 6, characterized in that, Step 5. Calculate the element stiffness matrix K corresponding to each material property according to the material parameters e and the stress function matrix S; specifically as follows: Calculate the element stiffness matrix K of each material at the Gauss integration points according to the element node numbers e : where B is a geometric matrix; B T is the transpose of the geometric matrix; |J| is the Jacobian determinant; t represents the element thickness, and in two-dimensional problems, t = 1.

8. A new finite element optimization solution algorithm as claimed in claim 7, characterized in that, Calculate the stress function matrix S according to the unit node numbers and Equation (3): S = CB (5) Among them, C is the constitutive matrix.

9. The novel finite element optimization solution algorithm according to claim 8, characterized in that, Step 6. Apply boundary displacement loads, and calculate the equivalent nodal forces according to the element stiffness matrix K e to assemble the right-hand side vector f; specifically, displacement loads are applied at the boundary. For two-dimensional problems, the cases in the x-direction, y-direction, and xy-direction are considered. By traversing specific boundary elements and using the element stiffness matrix K e to convert them into equivalent nodal forces. Each node has degrees of freedom in both the x and y directions. The results calculated by node numbers are assembled into the global right-hand vector f through degree-of-freedom mapping.

10. A new finite element optimization solution algorithm as claimed in claim 8, characterized in that, Step 7. Solve the displacement field, transform the linear equation system matrix-vector product calculation into a convolutional neural network convolution process to solve the nodal displacements under given boundary conditions; specifically, use the Node-by-Node matrix-free method based on nodes to directly calculate the stiffness contributions at the nodal level, assemble the element stiffness matrices according to the assembly principle, and set the local stiffness matrix K batch formed by m nodes as a minibatch, calculate the linear equation system matrix-vector product, and finally transform it into a convolutional neural network convolution process to solve the nodal displacement components under given boundary conditions.

Images