VCFEM numerical simulation method considering elasticity and thermal strain for three-dimensional porous material

Abstract

The invention aims to provide a VCFEM numerical simulation method considering elasticity and thermal strain for a three-dimensional porous material, relates to the technical field of porous material analysis, and provides the VCFEM numerical simulation method considering elasticity and thermal strain for the three-dimensional porous material by using a three-dimensional stress hybridization element method and a variation principle. The method comprises the following steps of: deducing a modified complementary energy functional when thermal strain is considered and when the thermal strain is not considered based on a minimum complementary energy principle, scaling units, dividing each unit into a plurality of Delaunay tetrahedrons, and performing integration by adopting a Hammer numerical integration method; a three-dimensional stress function considering an ellipsoid shape is constructed. An ellipsoid division method is adopted for the hole microstructure, and the conversion relation between the local coordinates and the overall coordinates of the holes is solved. A high-order stress field in a unit is assumed and solved, and a stress calculation result of the VCFEM is compared with a result of MSC MARC, so that the accuracy and the efficiency of the VCFEM are verified.

Description

Technical Field

[0001] The present invention relates to the technical field related to porous material analysis, and in particular to a VCFEM numerical simulation method for three-dimensional porous materials taking elastic and thermal strains into consideration. Background Art

[0002] Porous materials generally refer to materials with a network structure composed of interconnected or closed pores, with the pore boundaries or surfaces composed of pillars or plates. Porous materials may contain irregular pore microstructures, which leads to thermal stress concentration at the irregular pore boundaries, further affecting the overall performance of the material. Therefore, it is very important to calculate the overall thermal stress field of porous materials.

[0003] Due to the complexity of the hole model structure, ordinary displacement finite element calculations often require a large number of units for numerical simulation. For three-dimensional porous material models, ordinary finite elements require hundreds of thousands of units for numerical simulation to achieve the corresponding accuracy, which greatly reduces the efficiency of the calculation.

[0004] In order to solve the elastic and thermal strain problems of porous materials, it is usually necessary to solve microstructures with complex geometric shapes. The existing technology has created a hypothetical stress hybrid element method. This method constructs units through a hybrid model, and establishes the correct modified complementary energy functional according to the needs of actual conditions through the minimum complementary energy principle. This calculation method is very effective in improving the accuracy of the unit stress field. Later, the Voronoi unit finite element method was established on the basis of the hybrid element. The unit grid is divided by the generation method of Voronoi graphics, and there is a particle reinforcement or hole at the center of each unit. The shortcoming that the polygonal grid is difficult to establish a displacement pattern is successfully solved through the stress hybrid element theory. Subsequently, the Voronoi unit is combined with the homogenization method to perform multi-scale analysis on non-uniform structures with arbitrary microstructures. However, the above studies mostly consider two-dimensional plane problems, while in practice, three-dimensional problems need to be considered. Summary of the invention

[0005] In view of the deficiencies of the prior art, the present invention provides the following technical solutions:

[0006] A VCFEM numerical simulation method for three-dimensional porous materials considering elasticity and thermal strain is proposed. Based on the minimum complementary energy principle, the modified complementary energy functionals when thermal strain is considered and when thermal strain is not considered are derived. After the unit is scaled, each unit is divided into multiple Delaunay tetrahedrons. The Hammer numerical integration method is used for integration to realize the calculation of three-dimensional Voronoi units with holes. The specific VCFEM numerical simulation method includes the following steps:

[0007] S1, derive the VCFEM functional without considering thermal strain based on the full-volume method: give the initial complementary energy functional based on the assumed stress hybrid finite element model, modify the complementary energy functional according to the surface force balance between the three-dimensional Voronoi units with holes, rewrite the modified complementary energy functional according to the material phase region and material properties of the three-dimensional hole unit containing only one phase of matrix material in the unit, select the stress function of the matrix material of the three-dimensional hole unit, and substitute the stress function into the rewritten modified complementary energy functional to obtain the modified complementary energy functional without considering thermal strain;

[0008] S2, derive the VCFEM functional considering thermal strain based on the full energy method: define the thermal strain caused by temperature, substitute the thermal strain into the modified complementary energy functional without considering the thermal strain obtained in S1, and obtain the modified complementary energy functional considering the thermal strain;

[0009] S3, first obtain the centroid coordinates of the unit and scale the unit using the centroid as the origin of the local coordinate system: divide the unit into n tetrahedrons in the Cartesian coordinate system, and then calculate the center coordinates of each tetrahedron through the center of the tetrahedron;

[0010] S4, the Delaunay triangulation method is used to divide the integration region, discretize the three-dimensional space into a geometric space composed of tetrahedrons, and restore the regular Delaunay tetrahedrons, resulting in the irregular complex integration domain being divided into regular geometric regions. In order to realize the integration of complex fields, the divided simple domains are integrated and accumulated during the calculation process;

[0011] S5, integration is performed using the Hammer numerical integration method.

[0012] Furthermore, in S1, the initial coenergy functional is,

[0013]

[0014] The given traction boundary Γt satisfies the equation At a given cell boundary The surface force satisfies the equation n + σ=n - ·σ equilibrium condition; using the Lagrange multiplier method to introduce the surface force continuity condition, we get the modified complementary energy functional equation,

[0015]

[0016] The three-dimensional hole element contains only one phase of matrix material. According to the material phase region and material properties of the element, the modified complementary energy functional is rewritten as an equation,

[0017]

[0018] The stress function of the three-dimensional hole element matrix material is selected as {σ m =[P m ]{β m}, and substitute it into the rewritten modified complementary energy functional, and we get:

[0019]

[0020] in:

[0021]

[0022] For {β m}、u1 e 、u2 e Taking partial derivatives, we can get the stationary value condition as follows:

[0023]

[0024] Rearranging the above equation, we get:

[0025]

[0026] According to the above formula, the stress parameter β m The relationship between nodes,

[0027]

[0028] Furthermore, in S3, it is first assumed that the coordinates of the four vertices of the tetrahedron are:

[0029] p1(x1,y1,z1),p2(x2,y2,z2),p3(x3,y3,z1)

[0030] Then, the center coordinates of the tetrahedron can be calculated as:

[0031]

[0032] If the center coordinates of n tetrahedrons are T1, T2, ...Tn, and their volumes are V1, V2, Vn, then the center coordinates of any three-dimensional polyhedron are:

[0033]

[0034] According to the center of any polyhedral unit obtained, the origin of the coordinate system is moved to the center of the coordinate system and divided by the scale factor to obtain the coordinates in the local coordinate system; the coordinates of any point p (x, y, z) on the element, when scaled, are expressed as

[0035] Furthermore, in S5, when the Hammer integration method is used, the numerical integration of the two-dimensional problem is expressed as:

[0036]

[0037] Where f is the integrand, n is the number of Hammer integration points, l p is the weight coefficient of the Hammer integration point;

[0038] In the local coordinate system,

[0039]

[0040] Among them, ξ p ,η p is the local coordinate of a point in the two-dimensional problem, and λ is the unit scaling factor;

[0041] Similarly, the integral of the G matrix is ​​expressed in the local coordinate system as:

[0042]

[0043] Where m represents the number of tetrahedrons after partitioning, n is the Hammer integration point format, and J i The Jacobian determinant for the volume integral of the i-th tetrahedron is expressed as the equation:

[0044]

[0045] Then the three-dimensional numerical problem can be expressed as:

[0046]

[0047] Furthermore, the stress function of the matrix material is expressed as express:

[0048]

[0049] The interaction stress function of the matrix material is expressed as express:

[0050]

[0051] Near the hole Away from the hole After constructing the matrix stress polynomial, the stress field is derived from the Maxwell stress function. The matrix stress polynomial is substituted into the matrix, interface and inclusion stress fields. The most widely used Maxwell stress function is selected as the stress field introduced in the hybrid stress element.

[0052] Compared with the prior art, the technical solution of this application has the following beneficial effects:

[0053] The present invention is based on the VCFEM method, uses the variational principle to correct the functional with the initial thermal strain field, and proposes a three-dimensional Voronoi unit with holes that takes into account elasticity and thermal strain; the three-dimensional Voronoi can use fewer units for numerical simulation, compared with MSC MARC using 1776089 units to perform numerical simulation on a model with multiple holes, VCFEM can achieve the same accuracy using only 27 units; the computational efficiency of VCFEM is significantly improved compared with the displacement finite element; and the stress concentration of the Voronoi unit with holes at the hole boundary of the present invention is very obvious and clear, which proves that VCFEM is more accurate in analyzing three-dimensional models containing microstructures. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 This is a schematic diagram of VCFEM with holes;

[0055] Figure 2 It is the Delaunay partition diagram of the three-dimensional polyhedral unit;

[0056] Figure 3 The natural coordinate diagram for triangles and tetrahedrons;

[0057] Figure 4 It is the ellipsoid division diagram;

[0058] Figure 5 Structural model diagram for validity verification;

[0059] Figure 6 Schematic diagram of stress cloud map slices of MARC and VCFEM calculation results for validity verification;

[0060] Figure 7 This is a comparison chart of the stress numerical calculation results of MARC and VCFEM during effectiveness verification;

[0061] Figure 8 It is a schematic diagram of the stress cloud map slice of the calculation results of MARC and VCFEM when verifying a single hole under fixed constraints;

[0062] Fig. 9 This is a comparison chart of the stress numerical calculation results of MARC and VCFEM for single hole verification under fixed constraints;

[0063] Fig.10 It is a schematic diagram of the stress cloud map slice of the calculation results of MARC and VCFEM for single hole verification under a given displacement;

[0064] Fig.11A comparison chart of the stress numerical calculation results of MARC and VCFEM for single hole verification under given displacement;

[0065] Fig.12 This is a multi-unit model diagram with holes when verifying a single hole under a given displacement;

[0066] Fig.13 It is a schematic diagram of the stress cloud map slice of the calculation results of MARC and VCFEM for multi-hole verification under given displacement;

[0067] Fig.14 A comparison chart of the stress numerical calculation results of MARC and VCFEM for multi-hole verification under given displacement;

[0068] Fig.15 A model diagram of 40 randomly distributed hole units constructed using VCFEM. DETAILED DESCRIPTION

[0069] This application proposes a three-dimensional Voronoi unit with holes that takes into account elastic strain and thermal strain, and uses it to simulate the stress results of three-dimensional porous materials under the influence of displacement field and temperature field. By deriving the modified complementary energy functional and using the Delaunay partitioning and Hammer numerical integration method, the calculation of the three-dimensional Voronoi unit with holes is realized. And by comparing the total number of units and nodes, stress cloud diagrams and stress values ​​on the selected path with MSC MARC, the effectiveness and efficiency of the VCFEM method calculation are verified.

[0070] Embodiment 1:

[0071] 1.1 Derivation of VCFEM functional based on the full-volume method without considering thermal strain:

[0072] Based on the assumed stress hybrid finite element model, the hole element is constructed as follows: Figure 1 As shown, the initial coenergy functional expression is given:

[0073]

[0074] Where S is the elastic flexibility matrix, is the given boundary condition, T is the surface force of the current boundary, Γt e is the displacement boundary of the unit. Based on the variational principle of the total quantity theory, a typical three-dimensional Voronoi unit with holes is Figure 1 As shown, the matrix phase material of each unit is Ω m ;n e is the outer normal of the cell boundary, n m is the outward normal from the pore to the matrix on the interface between the matrix and the hole. Assuming that the surface forces between the units are balanced, the equation is satisfied at the given surface force boundary Γt:

[0075]

[0076] At a given cell boundary The surface forces satisfy the following equilibrium conditions:

[0077] n + σ=n - ·σ (3)

[0078] The surface force continuity condition is satisfied on the boundary between units. The surface force continuity condition is introduced by using the Lagrange multiplier method to obtain the modified complementary energy functional equation:

[0079]

[0080] The displacement on the boundary can be expressed by the shape function on the boundary as follows:

[0081] u＝Lu e (5)

[0082] The three-dimensional hole unit contains only one phase of matrix material. The functional is refined according to the characteristics of the hole unit material. The area Ω where the matrix material exists in the unit m , displacement boundary Γt of the matrix material e , the force boundary of the whole unit Γt, the interface between the matrix material and the hole Γ i ; The modified complementary energy functional is rewritten as an equation based on the material phase region and material properties of the unit:

[0083]

[0084] Where u is the displacement at any position on the unit boundary, u e is the node displacement on the unit boundary, L is the shape function on the unit boundary; u1 is the displacement field of the matrix material boundary u1, u2 is the displacement field on the interface between the matrix and the hole u2; the selection of the stress function of the matrix material of the three-dimensional hole unit should be:

[0085] {σ m =[P m ]{β m} (7)

[0086] The construction of the stress field can be composed of two parts: the Airy stress function and the interaction function, as shown in the equation:

[0087] σ pol =P pol β pol (8)

[0088] In order to make the stress calculation results more accurate, a high-order stress function is constructed using the Airy stress function of polynomials with different numbers of terms; the interaction force function considers the influence of the hole on the stress field and eliminates the influence of the interaction force at positions far away from the hole:

[0089]

[0090] The boundary displacements of the element can then be expressed as the interpolation of the generalized displacements of each boundary node:

[0091]

[0092] Among them, L1 is the shape function on the boundary of the matrix material, L2 is the shape function on the boundary between the hole and the matrix, and u1 e is the node on the boundary of the matrix material, u2 e are the internal nodes of the unit on the boundary between the matrix material and the hole;

[0093] Substituting equation (7) into the modified complementary energy functional in equation (6), we can obtain

[0094]

[0095] Rearranging equation (11) yields:

[0096]

[0097] in:

[0098]

[0099] For {β m}、u1 e 、u2 e Taking partial derivatives, we can get the stationary value condition as follows:

[0100]

[0101] Rearranging the above equation, we can get:

[0102]

[0103] According to equation (19), the stress parameter β can be obtained m Relationship with nodes:

[0104]

[0105] 1.1.1 Elimination of rigid body displacement of hole: In order to ensure that the displacement of the node on the hole surface is consistent with the displacement of the node outside the unit, an additional constraint φ is introduced by the lagrange multiplier method when considering the rigid body motion, which just eliminates the rigid body displacement of the hole in the unit and ensures the unit stiffness matrix K e Not strange. The modified single rigid K e It can be expressed as:

[0106]

[0107] 1.1.2 Condensation of freedom degree of hole unit:

[0108]

[0109] According to equation (22),

[0110]

[0111] q ei is the external node displacement, q in i is the internal node displacement, F out is the external nodal force, F in is the internal nodal force;

[0112]

[0113] F in =0 (28)

[0114] Combining equations (23), (27), and (28) yields equation (29) that involves only the displacement of the unit's external nodes and the displacement of other unit's external nodes, and equation (30) that involves only the displacement of each unit's internal nodes and the displacement of the external nodes:

[0115]

[0116] Through equation (30), the relationship between internal node displacement, internal node force and external node displacement can be obtained as equation (31):

[0117] q in i =K 22 -1 (F in -K 12 T q ei ) (31)

[0118] Substituting equation (31) into equation (29), we can obtain:

[0119]

[0120] The solution relationship can be expressed as:

[0121]

[0122] in,

[0123] K′ e =(K 11 -K 12 K 22 -1 K 12 T ) (34)

[0124] Equivalent nodal forces:

[0125]

[0126] Then, according to equation (34) and equation (33) into equation (32), we can get:

[0127] ∑ e F=∑ e (F+G T H -1 ) (36).

[0128] 1.2 Derivation of VCFEM functional considering thermal strain based on the full mass method:

[0129] When thermal deformation (expansion or contraction) caused by temperature change is constrained, stress will be generated; this type of stress is usually called thermal stress / temperature stress. Assuming that the temperature change of the porous material is Δt, it can be known that the deformation along a certain direction is:

[0130] Δl=αΔtl(37)

[0131] Where α is the linear expansion coefficient, and the thermal strain caused by temperature is defined as ε th , then the total strain of the porous material during elastic deformation is:

[0132] ε=σS+ε th (38)

[0133] Where S is the flexibility matrix. The selection of stress function is the same as that of unit boundary displacement, and the thermal strain ε caused by temperature is th Substituting into equation (12) we can obtain the modified complementary energy functional considering thermal strain:

[0134]

[0135] in

[0136]

[0137] [H m ] th =[P m ]ε th (41)

[0138]

[0139] ε th =θαI(45)

[0140] In equation (45), θ is the temperature difference, α is the thermal expansion coefficient, and I is a vector matrix that can be expressed as:

[0141] I = [1 0 0 1 0 0 1 0 0] (46)

[0142] Taking partial derivative of the modified complementary energy functional equation (39) obtained above, we can obtain the stationary value condition:

[0143]

[0144] By rearranging equation (47), we can obtain:

[0145]

[0146] According to the equation, we can get:

[0147]

[0148]

[0149] 1.2.1 Elimination of rigid body displacement of holes: In order to ensure that the displacement of the nodes on the hole surface is consistent with the displacement of the nodes outside the unit, an additional constraint φ is introduced by the lagrange multiplier method when considering the rigid body motion, which just eliminates the rigid body displacement of the hole in the unit and ensures that the unit stiffness matrix K e Not strange. The modified single rigid K e It can be expressed as:

[0150]

[0151] 1.2.2 Internal degree of freedom condensation: Equation (50) can be rewritten as:

[0152]

[0153] Where K 11 , K 12 ,(K 12 ) T , K 22 The same as equations (24), (25), and (26), q ei is the external node displacement, qin i is the internal node displacement, F out is the external nodal force, F in is the internal nodal force.

[0154]

[0155] Combining equations (53), (54), and (55) yields equation (56) that involves only the displacement of the unit's external nodes and the displacement of other unit's external nodes, and equation (57) that involves only the displacement of each unit's internal nodes and the displacement of the external nodes:

[0156]

[0157] Through equation (57), the relationship between internal node displacement, internal node force and external node displacement can be obtained as equation (58):

[0158] q in i =K 22 -1 (F in -K 12 T q ei ) (58)

[0159] Substituting equation (58) into equation (56), we can obtain:

[0160]

[0161] The solution relationship can be expressed as:

[0162]

[0163] in,

[0164] K′ e =(K 11 -K 12 K 22 -1 K 12 T ) (61)

[0165] Equivalent nodal forces:

[0166]

[0167] Then, according to equation (62) and equation (61) substituting into equation (60), we can obtain:

[0168]

[0169] 1.3 Unit Scaling:

[0170] In order to solve the problem of [H in equation (40) caused by too large or too small values ​​in the calculation process m ] and [H m ] th In order to avoid large errors in the calculation results, the unit needs to be scaled to obtain the centroid coordinates for numerical calculations. The method given in this application is to divide the unit in the Cartesian coordinate system into n tetrahedrons, and then calculate the center coordinates of each tetrahedron after division. The center of the unit is obtained indirectly through the center of the tetrahedron. First, assume that the coordinates of the four vertices of the tetrahedron are:

[0171] p1(x1,y1,z1),p2(x2,y2,z2),p3(x3,y3,z1) (64)

[0172] Then, the center coordinates of the tetrahedron can be calculated as:

[0173]

[0174] If the center coordinates of n tetrahedrons are T1, T2, ...Tn, and their volumes are V1, V2, Vn, then the center coordinates of any three-dimensional polyhedron are:

[0175]

[0176] According to the center of any polyhedral unit, the origin of the coordinate system is moved to the center of the coordinate system and divided by the scale factor to obtain the coordinates in the local coordinate system. The coordinates of any point p(x,y,z) on the element, when scaled, can be expressed as where ξ,η, They can be expressed as the following equations (67):

[0177]

[0178] Among them, γ is the scaling factor when converting from global coordinates to local coordinates.

[0179] 1.4Hammer numerical integration:

[0180] This application uses the Delaunay triangulation method to divide the integration area. Figure 2It explains how to discretize three-dimensional space into a geometric space composed of tetrahedrons using Delaunay triangulation. The restoration is divided into regular Delaunay tetrahedrons, resulting in the irregular complex integral domain being divided into regular geometric regions; in order to realize the integration of complex fields, the divided simple domains can be integrated and accumulated during the calculation process. In the specific implementation process, after division, each integral domain is composed of tetrahedrons, and the element faces are composed of triangles; it is more convenient to use natural coordinates to describe triangles or tetrahedrons. Taking a triangle as an example, let the Cartesian coordinates of the three vertices be A1, A2, and A3, such as Figure 3 As shown. The coordinates of any point in the triangle can be described as natural coordinates:

[0181] A p =[L2 L2 L3][A1 A2 A3] T (68)

[0182] Among them, L i (i=1,2,3) is the area of ​​the triangle determined by the two points other than point i. In three-dimensional problems, tetrahedral elements can be described in a similar way; assuming that the Cartesian coordinates of the four vertices are A1, A2, A3, and A4, respectively, such as Figure 2 As shown, the natural coordinates of any point in the tetrahedron can be expressed as the equation:

[0183] A p =[H1 H2 H3 H4][A1 A2 A3 A4] T (69)

[0184] Among them, H i (i=1,2,3,4).

[0185] When the Hammer integration method is used, the numerical integration of the two-dimensional problem can be expressed as:

[0186]

[0187] Where f is the integrand, n is the number of Hammer integration points, l p is the weight coefficient of the Hammer integration point. In the local coordinate system, equation (69) can be written as:

[0188]

[0189] Among them, ξ p ,η p That is similar to equation (67), ξ,η is the local coordinate of a point in the two-dimensional problem, and λ is the unit scaling factor. Similarly, the integral of the G matrix can be expressed in the local coordinate system as:

[0190]

[0191] Where m represents the number of tetrahedrons after partitioning, n is the Hammer integration point format, and J i The Jacobian determinant for the volume integral of the i-th tetrahedron is:

[0192]

[0193] Then the three-dimensional numerical problem in this application can be taken as:

[0194]

[0195] Similarly, the integral form of the H matrix in equation (12) in the local coordinate system can be written as follows:

[0196]

[0197] 1.5 Selection of stress function of three-dimensional Voronoi unit with holes:

[0198] In three-dimensional problems, it is necessary to choose a reasonable stress function to accurately reflect the stress field of the unit. The interaction stress function of the matrix is ​​expressed as The construction of polynomial stress function and interaction force function should meet the following requirements: The columns in must be linearly independent. Only when this condition is met can the H matrix be invertible. The selection of the stress function polynomial needs to fully consider the shape of the hole. The interaction stress term approaches zero at places far away from the contact boundary of different phases. The continuity condition of the surface force must be met at the contact boundary of different phases.

[0199]

[0200] This application further simplifies the interface interaction stress function and takes the following expression:

[0201]

[0202] Near the hole Away from the hole After constructing the matrix stress polynomial, the stress field is derived from the Maxwell stress function. Substituting the matrix stress polynomial into the matrix, interface and inclusion stress fields are as follows:

[0203]

[0204] This application uses the most widely used Maxwell stress function as the stress field introduced in the hybrid stress element, taking the third-order complete polynomial as an example:

[0205] Φ(x,y,z)=β1x 3 +β2y 3 +β3z 3 +β4x 2 y+β5x 2 z+β6y 2 x+β7y 2 z+β8z 2 x+β8z 2 y+β9xz 2 +β9yz 2 +β 10 xyz+β 11 x 2 +β 12 y 2 +β 13 z 2 +β 14 xy+β 15 xz+β 16 yz (79)

[0206] Taking partial derivatives of the stress function equation (76), we can obtain the stress components:

[0207]

[0208] Equation (81) can also be expressed in the form of equation (7), that is, {σ m =[P m ]{β m}. Among them, {β m} is the stress parameter. Equation (80) is the stress parameter {β m The number of stress parameters {β m The rationality of the number of} determines the quality of the stress field results. When selecting the number of stress parameters, it is necessary to consider the completeness of the stress polynomial and the relationship between the number of stress parameters and the number of unit nodes and the degree of freedom of the rigid body. According to the property that the linear equation system has a unique solution, the following conclusions are drawn:

[0209] n β ≥n n ·DD F (82)

[0210] Among them, n β That is, the stress parameter {β m}. n n is the number of nodes in the unit, D is the degree of freedom of each node, and DF is the degree of freedom of the rigid body displacement of the unit. The three-dimensional unit is generally taken as 6. In order to satisfy the completeness of the stress function, n is usually β The value obtained by direct calculation cannot be used, but the corresponding fixed value is taken according to the order of the function Φ, see Table 1:

[0211] Table 1 corresponds to the order n of the function Φ β Table of Values

[0212]

[0213]

[0214] 1.6 Generation of internal points of the unit and local coordinate system of the hole:

[0215] To add internal points in a unit, we first need to solve the actual coordinates of all internal points. In this application, a complex local coordinate system is established with the center point of the ellipsoid as the origin and the directions of the three axes of the ellipsoid as coordinates. The conversion relationship between the coordinates of the center point of the ellipsoid and the ellipsoid system and the actual coordinate system is obtained. The two parts of the deflection angle of the sphere inside the unit in space are closely related to the shape of the particle reinforcement phase. Taking the ellipsoidal inclusion as an example, we first need to know the control equation of the inclusion shape. The general form of the ellipsoid equation is:

[0216]

[0217] Where a, b, and c are the lengths of the three axes of the ellipsoid. To solve the actual coordinates of the internal point, it is necessary to find the internal point in the local coordinates of the inclusion based on the governing equations of the ellipsoid and the local coordinates of the inclusion. Figure 4 As shown, two partitioning parameters are required to partition the surface of the ellipsoidal inclusion: the number of segments S to be partitioned along the local coordinate x-axis and the number of points to be partitioned on each circumference.

[0218] Assume that the ellipsoid is divided into S ellipsoidal surfaces along the local coordinate axis X. The number of segments is S+1. The distance between each segment is dS, and dS can be expressed as the following equation:

[0219]

[0220] The local coordinates of each slice in the X direction can be obtained by dividing the number of segments into the following equation:

[0221]

[0222] Where a is the length of the ellipsoid in the x-axis direction. i indicates the number of segments of the ellipsoid currently divided. S is the number of faces divided. After obtaining the local coordinates in the x-direction, continue to divide each ellipsoid face. Substituting the local coordinate values ​​in the x-direction on each face into the shape control equation of the ellipsoid, the elliptical shape equation of each face can be obtained:

[0223]

[0224] The circumference is divided into equal number of regions according to the number of circumferential division points N of each elliptical surface given in advance. The angle dθ in each region is represented by N:

[0225]

[0226] The y and z coordinates of all nodes on the circumference of the ellipse can be expressed in polar coordinates as follows:

[0227]

[0228] Wherein, j is the number of points on the current circumference, j = 1, 2, ..., N; N is the number of points on the circumference used in equation (56). j is the value of the current polar coordinate. Substituting the polar coordinates of y and z into the current ellipse equation (86) can solve the r of each point. j ,

[0229]

[0230] r j Substituting the value of into equation (88) can solve the complete coordinates of an internal point in the local coordinate system. The above division of the ellipsoidal hole surface can be easily realized through the loop algorithm, and the coordinate values ​​of all divided points in the local coordinate system of the ellipsoid can be solved. These points solved by equation (88) are the internal points subsequently added to the three-dimensional hole-containing stress hybrid unit. Before adding the internal points, the coordinates of these points need to be converted from the local coordinate system to the global coordinate system through coordinate transformation. 1.7 Conversion between the hole coordinate system and the matrix coordinate system:

[0231] The so-called conversion between the hole coordinate system and the matrix coordinate system is actually the conversion between the local coordinates and the global coordinates in space. According to equations (54) and (57), we can first assume that the local coordinates of the hole node are (x i ,y j ,z j ), when the origin of the local coordinates coincides with the origin of the global coordinates, the transformation relationship of the spatial coordinate system can be expressed as an equation:

[0232]

[0233]

[0234] In equation (91), α, β, and γ are the angles between the three main axes of the ellipsoid and the global coordinate axes, and x, y, and z are the global coordinates.

[0235] Embodiment 2:

[0236] In order to verify the effectiveness of the three-dimensional Voronoi unit, we used the commercial finite element software MSC MARC to construct a model with the same structure and displacement boundary conditions as the Voronoi unit. By comparing the calculation results of the two models, the effectiveness and efficiency of VCFEM were verified. The elastic model and the model considering thermal strain were verified respectively, and a three-dimensional single-hole model and a three-dimensional porous model were constructed. Due to the characteristics of the Voronoi hybrid element, the three-dimensional single-hole model can only have one unit when the mesh is divided. Therefore, the single-hole model and the porous model can be used to verify the effectiveness of VCFEM. All model material parameters are 2024-T6 aluminum alloy. The specific material parameters are shown in Table 2:

[0237] Table 2 Model material parameters

[0238]

[0239] 2.1 Verification of the validity of the hole model:

[0240] Single hole elasticity example: The structure model is as follows Figure 5 As shown, according to equation (78), there is a hole divided into 2 layers and 10 points in the center of the model, and the overall volume of the model is 20×20×20mm 3 , the top surface of the model is given a 0.002mm displacement boundary condition, and the surfaces perpendicular to the x-axis, y-axis, and z-axis are given a 0mm displacement boundary condition. Thermal strain is not considered in this example; the comparison of meshing between MARC and VCFEM is shown in Table 3:

[0241] Table 3 Comparison of meshing between MARC and VCFEM

[0242]

[0243] The calculation results of MARC and VCFEM are sliced ​​into stress cloud maps as follows Figure 6 As shown in the figure, the normal vector of the slice plane is (1, 0, 0) and passes through the point (0, 0, 0). By observing the stress concentration and distribution in the stress cloud diagram, it can be concluded that the calculation results of VCFEM and MARC are basically consistent, thus verifying the effectiveness of VCFEM; the stress values ​​in the calculation results of MARC and VCFEM are extracted along the Y-axis path and compared specifically. The comparison results are shown in Figure 7 .like Figure 7 As shown in the figure, the stress variation trends of VCFEM and MARC are basically consistent, which proves the accuracy of VCFEM. By observing the stress values ​​along the Y-axis path, it is found that the stress values ​​of VCFEM near the hole are larger than those of MARC. This is because more polynomial terms of the interpolation function are used to make the calculation results more accurate. In addition, the number of units and nodes used by VCFEM is much smaller than that used by MARC, which shortens the calculation time and significantly improves the calculation efficiency. From the numerical simulation results of the unit with microstructure containing holes, VCFEM simulates the microstructure more accurately.

[0244] 2.2 Thermal strain model verification:

[0245] 2.2.1 Single hole thermal strain calculation example verification under fixed constraint: The displacement boundary conditions of the model are as follows: Figure 5 As shown, the overall volume of the model is 20×20×20mm 3 According to equation (81), there is a hole model divided into 2 layers and 6 points in the center of the model. All displacement boundary conditions are 0 mm. This example considers thermal strain, and the given thermal expansion coefficient is shown in Table 2. The calculation results of MARC and VCFEM are sliced ​​into stress cloud maps as shown in Figure 8 As shown in the figure, the normal vector of the slice plane is (1, 0, 0) and passes through the point (0, 0, 0). By observing the stress concentration and distribution in the stress cloud diagram, it can be concluded that the calculation results of VCFEM and MARC are basically consistent, thus verifying the effectiveness of VCFEM; the stress values ​​in the calculation results of MARC and VCFEM are extracted along the Y-axis path and compared specifically. The specific path is as follows Fig. 9 The red line shows the stress data comparison results. Fig. 9 As shown in the figure, the stress change trends of VCFEM and MARC are basically the same, which proves the accuracy of VCFEM. By observing the stress values ​​along the Y-axis path, it is found that the stress values ​​of VCFEM near the hole are larger than those of MARC. This is because more polynomial terms of the interpolation function are used to make the calculation results more accurate. In addition, the number of units and nodes used by VCFEM is much smaller than that used by MARC, which shortens the calculation time and significantly improves the calculation efficiency. In summary, when only thermal strain is considered, the stress change trends of VCFEM and MARC are the same, and the values ​​are relatively close, with higher calculation efficiency, which proves the accuracy and efficiency of VCFEM.

[0246] 2.2.2 Single hole thermal strain calculation example verification under given displacement: In order to verify the correctness of the thermal strain model, the stress distribution and magnitude under given node temperature conditions and given the same displacement boundary conditions as the elastic model, that is, a displacement of 0.002 mm is given on the top surface and 0 mm is given on the other three surfaces. The specific position conditions are the same as Figure 5 Considering that thermal expansion is linear expansion, the given thermal expansion coefficients are shown in Table 2. Since metal materials will produce thermal strain when subjected to temperature changes, but due to the limitations of boundary conditions, metal materials cannot expand freely, so thermal stress will be generated due to the inability to expand freely. Fig.10 It can be shown that the stress distribution and concentration of MARC and VCFEM are well consistent, which proves the effectiveness of VCFEM. Fig.11 The red line in the figure is a specific path. By extracting the stress value on the selected path, a scatter plot of the stress value can be obtained. Fig.11 .pass Fig.11 By comparing the stress value point-line graphs, it can be concluded that when thermal strain and displacement are considered simultaneously, the stress change trends of VCFEM and MARC are the same, and the values ​​are relatively close.

[0247] 2.2.3 Verification of thermal strain under given displacement: When verifying the multi-element model with holes, the following Fig.12 The porous model with 3×3×3 units shown in the figure has the same displacement boundary conditions and node temperatures as the single-hole model in the multi-element case. Figure 5 Similarly, the calculated results are respectively taken through the points (0, 0, 0), (20, 0, 0), and (40, 0, 0), and multi-layer slices are taken with the vector (1, 0, 0) as the normal line. The stress cloud map results of the slices are shown in Fig.13 . Slice the stress cloud map results according to Fig.14 The red line segment path is the stress path. The stress data is extracted for comparison. The results are shown in Fig.14 The results show that when calculating multi-element models, VCFEM not only retains accuracy, but also greatly improves computational efficiency compared with MSC MARC because the number of nodes and elements used is much smaller than that used by MSC MARC. This is one of the advantages of VCFEM when analyzing models with multiple microstructures.

[0248] 2.3 Example of random distribution of multiple holes:

[0249] In summary, the three-dimensional Voronoi unit model considering the regular distribution of holes under the action of thermal strain and displacement field has been verified respectively; in reality, the holes in porous materials are often irregular in shape and randomly distributed in space. In order to simulate the mechanical properties of porous materials in real situations, VCFEM is used to construct the following Fig.15The model shown contains 40 randomly distributed holes. The model size is 40×40×40 mm3. The direction and size of the holes are fixed, and their center positions are randomly distributed in the model. The given displacement boundary conditions are the same as Figure 5 . Numerical simulation of stress Fig.15 It is proved that for the model with randomly distributed multiple holes, the calculation results of VCFEM have good regularity, which proves that VCFEM can obtain reliable stress results and the calculation efficiency of VCFEM is higher than that of displacement finite element.

[0250] Based on the VCFEM method, this paper uses the variational principle to modify the functional with the initial thermal strain field, and proposes a three-dimensional Voronoi element with holes that considers elasticity and thermal strain; the numerical simulation results of the three-dimensional Voronoi element are compared with the finite element software MSC MARC through numerical calculation examples 2.1, 2.2.1, 2.2.2, and 2.2.3, respectively, proving the effectiveness, high precision and high efficiency of VCFEM under the displacement field, temperature field and displacement field and temperature field. In large-scale engineering practice calculations, three-dimensional Voronoi can use fewer units for numerical simulation. For example, in numerical calculation example 2.2.3, compared with MSC MARC using 1776089 units to numerically simulate the model with multiple holes, VCFEM only uses 27 units to achieve the same accuracy; the computational efficiency of VCFEM is significantly improved compared with the displacement finite element; in addition, when analyzing the model with random distribution of multiple holes, as shown in numerical calculation example 2.3, VCFEM also demonstrates its effectiveness. Moreover, the stress concentration of the Voronoi unit containing holes at the hole boundary is very obvious and clear, which proves that VCFEM is more accurate in analyzing three-dimensional models containing microstructures, which is also one of the advantages of the VCFEM method.

[0251] The preferred embodiments of the present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains, characterized in that: Based on the minimum complementary energy principle, the modified complementary energy functionals when thermal strain is considered and when thermal strain is not considered are derived. After the unit is scaled, each unit is divided into multiple Delaunay tetrahedrons. The Hammer numerical integration method is used for integration to realize the calculation of three-dimensional Voronoi units with holes. The specific VCFEM numerical simulation method includes the following steps: S1, derive the VCFEM functional without considering thermal strain based on the full-volume method: give the initial complementary energy functional based on the assumed stress hybrid finite element model, modify the complementary energy functional according to the surface force balance between the three-dimensional Voronoi units with holes, rewrite the modified complementary energy functional according to the material phase region and material properties of the three-dimensional hole unit containing only one phase of matrix material in the unit, select the stress function of the matrix material of the three-dimensional hole unit, and substitute the stress function into the rewritten modified complementary energy functional to obtain the modified complementary energy functional without considering thermal strain; S2, derive the VCFEM functional considering thermal strain based on the full energy method: define the thermal strain caused by temperature, substitute the thermal strain into the modified complementary energy functional without considering the thermal strain obtained in S1, and obtain the modified complementary energy functional considering the thermal strain; S3, first obtain the centroid coordinates of the unit and scale the unit using the centroid as the origin of the local coordinate system: divide the unit into n tetrahedrons in the Cartesian coordinate system, and then calculate the center coordinates of each tetrahedron through the center of the tetrahedron; S4, the Delaunay triangulation method is used to divide the integration region, discretize the three-dimensional space into a geometric space composed of tetrahedrons, and restore the regular Delaunay tetrahedrons, resulting in the irregular complex integration domain being divided into regular geometric regions. In order to realize the integration of complex fields, the divided simple domains are integrated and accumulated during the calculation process; S5, integration is performed using the Hammer numerical integration method.

2. A VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains according to claim 1, characterized in that: In S1, the initial coenergy functional is, The given traction boundary Γt satisfies the equation At a given cell boundary On the surface, the surface force satisfies the equation n + σ=n - · The equilibrium condition of σ; The modified complementary energy functional equation is obtained by introducing the surface force continuity condition using the Lagrange multiplier method. The three-dimensional hole element contains only one phase of matrix material. According to the material phase region and material properties of the element, the modified complementary energy functional is rewritten as an equation, The stress function of the three-dimensional hole element matrix material is selected as {σ m =[P m ]{β m }, and substitute it into the rewritten modified complementary energy functional, and we get: in: For {β m }、u1 e 、u2 e Taking partial derivatives, we can get the stationary value condition as follows: Rearranging the above equation, we can get: According to the above formula, the stress parameter β is obtained m The relationship between nodes, 3. A VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains according to claim 2, characterized in that: In S2, assuming that the temperature change of the porous material is Δt, the deformation along a certain direction is: Δl=αΔtl Where α is the linear expansion coefficient, and the thermal strain caused by temperature is defined as ε th , then the total strain of the porous material during elastic deformation is: ε=σS+ε th Where S is the flexibility matrix, which takes the thermal strain ε caused by temperature as th Substituting into the modified complementary energy functional without considering thermal strain, we get the modified complementary energy functional considering thermal strain: in, [H m ] th =[P m ]e th e th =I will θ is the temperature difference, α is the thermal expansion coefficient, and I is the vector matrix expressed as: I＝[1 0 0 1 0 01 0 0] Taking the partial derivative of the modified complementary energy functional considering thermal strain obtained above, we can obtain the stationary value condition:

4. According to claim 1, a VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains is characterized in that: In S3, firstly, it is assumed that the coordinates of the four vertices of the tetrahedron are: p1(x1,y1,z1),p2(x2,y2,z2),p3(x3,y3,z1) Then, the center coordinates of the tetrahedron can be calculated as: If the center coordinates of n tetrahedrons are T1, T2, ...Tn, and their volumes are V1, V2, Vn, then the center coordinates of any three-dimensional polyhedron are: According to the center of any polyhedral unit obtained, the origin of the coordinate system is moved to the center of the coordinate system and divided by the scale factor to obtain the coordinates in the local coordinate system; the coordinates of any point p (x, y, z) on the element, when scaled, are expressed as where ξ,η, Respectively expressed as: γ is the scaling factor from global coordinates to local coordinates.

5. According to claim 1, a VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains is characterized in that: In S5, when the Hammer integration method is used, the numerical integration of the two-dimensional problem is expressed as: Where f is the integrand, n is the number of Hammer integration points, l p is the weight coefficient of the Hammer integration point; In the local coordinate system, Among them, ξ p ,η p is the local coordinate of a point in the two-dimensional problem, and λ is the unit scaling factor; Similarly, the integral of the G matrix is ​​expressed in the local coordinate system as: Where m represents the number of tetrahedrons after partitioning, n is the Hammer integration point format, and J i The Jacobian determinant for the volume integral of the i-th tetrahedron is expressed as the equation: Then the three-dimensional numerical problem can be expressed as:

6. According to claim 3, a VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains is characterized in that: The stress function of the matrix material is expressed as express: The interaction stress function of the matrix material is expressed as express: Near the hole Away from the hole After constructing the matrix stress polynomial, the stress field is derived from the Maxwell stress function. The matrix stress polynomial is substituted into the matrix, interface and inclusion stress fields. The most widely used Maxwell stress function is selected as the stress field introduced in the hybrid stress element.

Images