A VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains

Abstract

The present invention aims to propose a VCFEM numerical simulation method for three-dimensional porous materials that takes into account elastic and thermal strains. This method relates to the technical field of porous material analysis. Utilizing the three-dimensional stress hybrid element method and the variational principle, this application proposes a VCFEM numerical simulation method for three-dimensional porous materials that takes into account elastic and thermal strains. Based on the principle of minimum complementary energy, modified complementary energy functionals are derived for both cases with and without thermal strain. After scaling, each cell is divided into multiple Delaunay tetrahedrons, and integration is performed using the Hammer numerical integration method. A three-dimensional stress function that takes into account an ellipsoidal shape is constructed. An ellipsoidal partitioning method is employed for the pore microstructure, and the transformation relationship between the local and global coordinates of the pores is solved. By assuming and solving a high-order stress field within the cell, the stress calculation results of the VCFEM are compared with those of the MSC MARC method, verifying the accuracy and efficiency of the VCFEM.

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 consisting of interconnected or closed pores, with the pore boundaries or surfaces formed by pillars or plates. Porous materials may contain irregular pore microstructures, which can lead to thermal stress concentrations at the irregular pore boundaries, further affecting the overall performance of the material. Therefore, calculating the overall thermal stress field of porous materials is crucial.

[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. There is an existing technology that 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 based on the hybrid element. The unit grid is divided by generating Voronoi graphics, and there is a particle reinforcement or hole at the center of each unit. The stress hybrid element theory successfully solves the shortcoming that polygonal grids are difficult to establish displacement patterns. Subsequently, the Voronoi unit is combined with the homogenization method to perform multi-scale analysis of 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 shortcomings of the prior art, the present invention provides the following technical solutions:

[0006] A VCFEM numerical simulation method for three-dimensional porous materials considering elastic and thermal strains is proposed. Based on the principle of minimum complementary energy, modified complementary energy functionals are derived when thermal strain is considered and when thermal strain is not considered. After scaling the unit cell, each unit cell 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. Derivation of VCFEM functional without considering thermal strain based on the full-quantity method: The initial complementary energy functional is given based on the assumed stress hybrid finite element model. The complementary energy functional is modified according to the surface force balance between the three-dimensional Voronoi cells containing holes. The modified complementary energy functional is rewritten according to the material phase region and material properties of the three-dimensional hole cell containing only one phase of matrix material in the cell. The stress function of the three-dimensional hole cell matrix material is selected and substituted 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-quantity method: define the thermal strain caused by temperature, substitute the thermal strain into the modified complementary energy functional obtained in S1 without considering the thermal strain, 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, discretizing the three-dimensional space into a geometric space composed of tetrahedrons, and restoring 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 On the surface, the surface force satisfies the equation n + σ=n - The equilibrium condition of σ; the surface force continuity condition is introduced by Lagrange multiplier method to obtain the modified complementary energy functional equation,

[0015]

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

[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 co-energy functional, and we get:

[0019]

[0020] in:

[0021]

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

[0023]

[0024] Rearranging the above equations, we get:

[0025]

[0026] According to the above formula, the stress parameter β is obtained 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, 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] Among them, m represents the number of tetrahedrons after division, n is the Hammer integration point format, 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 is 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 existing technology, 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 modify the functional with an initial thermal strain field, and proposes a three-dimensional Voronoi element with holes that takes into account elasticity and thermal strain. The three-dimensional Voronoi can use fewer elements for numerical simulation. Compared with MSC MARC using 1,776,089 elements to numerically simulate a model with multiple holes, VCFEM can achieve the same accuracy using only 27 elements. The computational efficiency of VCFEM is significantly improved compared to the displacement finite element. In addition, the stress concentration of the Voronoi element 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 Schematic diagram of VCFEM with holes;

[0055] Figure 2 is the Delaunay partitioning diagram of the three-dimensional polyhedral element;

[0056] Figure 3 The natural coordinate diagrams 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 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 for effectiveness verification;

[0061] Figure 8 A schematic diagram of the stress contour slice of the calculation results of MARC and VCFEM for single hole verification under fixed constraints;

[0062] Figure 9 A comparison chart of the stress numerical calculation results of MARC and VCFEM for single hole verification under fixed constraints;

[0063] Figure 10 Schematic diagram of stress contour slices of MARC and VCFEM calculation results for single hole verification under given displacement;

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

[0065] Figure 12 This is a diagram of a multi-element model with holes when verifying a single hole under a given displacement;

[0066] Figure 13 Schematic diagram of stress cloud slices of MARC and VCFEM calculation results for multi-hole verification under given displacement;

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

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

[0069] This application proposes a three-dimensional Voronoi element with holes that accounts for elastic and thermal strains and uses it to simulate the stresses of three-dimensional porous materials under the influence of displacement and temperature fields. By deriving a modified complementary energy functional and utilizing Delaunay partitioning and Hammer numerical integration, the calculation of the three-dimensional Voronoi element with holes is realized. The effectiveness and efficiency of the VCFEM method are verified by comparing the total number of elements and nodes, stress contours, and stress values ​​along selected paths with MSC MARC.

[0070] Example 1:

[0071] 1.1 Derivation of VCFEM functional based on the full-quantity 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 as follows: Figure 1 As shown, the matrix phase material of each unit is Ω m ;n e is the outer normal of the unit 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 Lagrange multiplier method is used to introduce the surface force continuity condition and 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 element contains only one phase of matrix material. The functional is refined according to the characteristics of the hole element material. In the area where the matrix material exists in the unit, Ω m , the displacement boundary of the matrix material Γt e , the force boundary Γt of the entire unit, the interface between the matrix material and the hole Γ i ; The modified complementary energy functional is rewritten as the 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, and u2 is the displacement field on the interface between the matrix and the hole u2; the stress function of the matrix material of the three-dimensional hole unit should be selected as:

[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 takes into account 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 node in the boundary:

[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 condition as follows:

[0100]

[0101] Rearranging the above equations, 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 hole rigid body displacement: 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 corrected single K e It can be expressed as:

[0106]

[0107] 1.1.2 Condensation of degrees of freedom of hole elements:

[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 units' external nodes, and equation (30) that involves only the displacement of each unit's internal nodes and the displacement of its external nodes:

[0115]

[0116] The relationship between internal node displacement, internal node force and external node displacement can be obtained through equation (30) as shown in 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) and substituting them into equation (32), we can obtain:

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

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

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

[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 the partial derivative of the modified co-energy functional equation (39) obtained above, we can obtain the stationary value condition:

[0143]

[0144] Arranging 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 the unit stiffness matrix K e Not strange. The corrected single 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) involving only the displacement of the unit's external nodes and the displacement of other units' external nodes, and equation (57) involving only the displacement of each unit's internal nodes and the displacement of the external nodes:

[0156]

[0157] The relationship between internal node displacement, internal node force and external node displacement can be obtained through equation (57) as shown in 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) and put them into equation (60), we can get:

[0168]

[0169] 1.3 Unit Scaling:

[0170] In order to solve the problem of [H m ] and [H m ] th In order to solve the irreversible problem and avoid large errors in the calculation results, it is necessary to scale the unit to obtain the centroid coordinates during numerical calculations. The method proposed 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 2The paper explains how to use Delaunay triangulation to discretize three-dimensional space into a geometric space composed of tetrahedrons. 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 consists 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 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, as follows 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 Similar to Equation (67), ξ,η are the local coordinates of a point in the two-dimensional problem, and λ is the unit scaling factor. Similarly, the integral of the G matrix in the local coordinate system can be expressed as:

[0190]

[0191] Among them, m represents the number of tetrahedrons after division, n is the Hammer integration point format, 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 for 3D Voronoi element 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) where {σ 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, the number of element nodes, and the rigid body's degrees of freedom. Based on the property that a system of linear equations has a unique solution, the following conclusions can be 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, DF is the degree of freedom of the rigid body displacement of the unit, and the three-dimensional unit is generally taken as 6. In order to meet 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 cell, 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 ellipsoid center point and the ellipsoid system and the actual coordinate system is obtained. The two parts of the deflection angle of the sphere inside the cell in space are used to obtain the generation method of the local inclusion point, which is 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, we need 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 in Figure 1, dividing the surface of the ellipsoidal inclusion requires two division parameters: the number of segments S along the local coordinate x-axis and the number of points 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. Then 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 equations:

[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 x-direction coordinate values ​​of 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 an 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 coordinate values ​​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 ellipsoidal local coordinate system 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 transformation between the hole coordinate system and the matrix coordinate system is actually the transformation 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] Example 2:

[0236] To verify the effectiveness of the 3D Voronoi element, we used the commercial finite element software MSC MARC to construct a model with the same structure and displacement boundary conditions as the Voronoi element. By comparing the calculation results of the two models, we verified the effectiveness and efficiency of VCFEM. We verified the elastic model and the model considering thermal strain respectively, and constructed a 3D single-hole model and a 3D multi-hole model. Due to the characteristics of the Voronoi hybrid element, the 3D single-hole model can only have one unit during meshing. Therefore, the single-hole model and the multi-hole 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 in the center of the model divided into 2 layers with 10 points, 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 MARC and VCFEM meshing 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. The comparison results are shown in Figure 7 .like Figure 7 As shown, the stress trends obtained by VCFEM and MARC are essentially identical, demonstrating the accuracy of VCFEM. Observing the stress values ​​along the Y-axis, VCFEM results show that stress values ​​near the holes are larger than those obtained by MARC. This is due to the use of more polynomial terms in the interpolation function to achieve more accurate results. Furthermore, VCFEM uses far fewer elements and nodes than MARC, shortening the calculation time and significantly improving efficiency. Numerical simulation results for microstructures containing holes demonstrate that VCFEM provides a more accurate simulation of the microstructure.

[0244] 2.2 Thermal strain model verification:

[0245] 2.2.1 Single hole thermal strain calculation 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. 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: Figure 9 The red line shows the stress data comparison results. Figure 9 As shown, the stress trends obtained by VCFEM and MARC are essentially identical, demonstrating the accuracy of VCFEM. Observing the stress values ​​along the Y-axis, VCFEM results show that stress values ​​near the hole are larger than those obtained by MARC. This is due to the use of more polynomial terms in the interpolation function to achieve more accurate results. Furthermore, the number of elements and nodes used by VCFEM is significantly smaller than that used by MARC, shortening the calculation time and significantly improving computational efficiency. In summary, when considering only thermal strain, the stress trends obtained by VCFEM and MARC are identical, with similar values, demonstrating higher computational efficiency, demonstrating 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 are analyzed under given node temperature conditions and given the same displacement boundary conditions as the elastic model, that is, a displacement of 0.002mm is given on the top surface and 0mm displacement is given on the other three surfaces. The specific position working conditions are the same as Figure 5 Considering thermal expansion as 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. Figure 10 It can be shown that the stress distribution and concentration of MARC and VCFEM are in good agreement, which proves the effectiveness of VCFEM. Figure 11 The red line segment 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. Figure 11 .pass Figure 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 multi-hole thermal strain calculation under given displacement: When verifying the multi-element model with holes, the following Figure 12 The porous model with 3×3×3 elements 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 sliced ​​through the points (0, 0, 0), (20, 0, 0), and (40, 0, 0), and the vector (1, 0, 0) is used as the normal to take multi-layer slices. The stress cloud map results of the slices are shown in Figure 13 . Slice the stress cloud map into Figure 14 The red line segment path is the stress path. The stress data is extracted for comparison. The results are shown in Figure 14 The results show that when calculating multi-element models, VCFEM not only retains accuracy, but also significantly improves computational efficiency compared to 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 containing multiple microstructures.

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

[0249] In summary, the three-dimensional Voronoi unit model with regular distribution of pores under thermal strain and displacement field has been verified respectively; in reality, the pores 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 conditions, the following VCFEM is used to construct Figure 15The model shown contains 40 randomly distributed hole units with a size of 40×40×40 mm3. The directions and sizes 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 Figure 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 utilizes the variational principle to modify a functional with an initial thermal strain field, and proposes a three-dimensional Voronoi element with holes that accounts for elastic and thermal strains. Numerical simulation results using the 3D Voronoi element are compared with the finite element software MSC MARC through numerical examples 2.1, 2.2.1, 2.2.2, and 2.2.3, respectively. The effectiveness, high accuracy, and high efficiency of the VCFEM are demonstrated for displacement fields, temperature fields, and simultaneous displacement and temperature fields. In large-scale engineering practice, the 3D Voronoi element can be simulated using fewer elements. For example, in Example 2.2.3, compared to MSC MARC's 1,776,089 elements for a model with multiple holes, the VCFEM achieves the same accuracy with only 27 elements. The VCFEM significantly improves its computational efficiency compared to the displacement finite element method. Furthermore, the VCFEM also demonstrates its effectiveness in analyzing models with randomly distributed holes, as shown in Example 2.3. In addition, the stress concentration of the Voronoi element containing holes at the hole boundary in this paper 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 intended only to help illustrate the present invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the present invention to the specific embodiments described. Obviously, many modifications and variations are possible based on the content of this specification. These embodiments are selected and described in detail in this specification to better explain the principles and practical applications of the present invention, thereby enabling those skilled in the art to better understand and utilize the present invention. 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 by: Based on the minimum complementary energy principle, the modified complementary energy functionals with and without thermal strain are derived. After scaling the cells, each cell is divided into multiple Delaunay tetrahedrons. The Hammer numerical integration method is used for integration to achieve the calculation of three-dimensional Voronoi cells with holes. The specific VCFEM numerical simulation method includes the following steps: S1. Derivation of VCFEM functional without considering thermal strain based on the full-quantity method: The initial complementary energy functional is given based on the assumed stress hybrid finite element model. The complementary energy functional is modified according to the surface force balance between the three-dimensional Voronoi cells containing holes. The modified complementary energy functional is rewritten according to the material phase region and material properties of the three-dimensional hole cell containing only one phase of matrix material in the cell. The stress function of the three-dimensional hole cell matrix material is selected and substituted 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-quantity method: define the thermal strain caused by temperature, substitute the thermal strain into the modified complementary energy functional obtained in S1 without considering the thermal strain, 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, discretizing the three-dimensional space into a geometric space composed of tetrahedrons, and restoring 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.

Images