Three-dimensional Voronoi unit simulation method for solving elasticity problem of interface phase-containing particle reinforced composite material

Abstract

The invention designs a three-dimensional Voronoi unit simulation method aiming at the elasticity problem of a particle reinforced composite material containing an interface phase, and a three-dimensional correction complementary energy functional containing inclusions, a matrix and an interface is deduced, so that the method can be used for calculating a three-dimensional three-phase composite material; besides, stress functions are constructed in the three integral domains respectively, Maxwell complete polynomials are adopted and comprise a stress function polynomial and an interaction force function polynomial, and good convergence and continuity are achieved for three-dimensional interface units with different node numbers. Compared with commercial software, the method is faster and more efficient when the Voronoi model of the three-dimensional interface is calculated, and results of hundreds of thousands of units in the commercial software can be obtained with fewer units when the same model is calculated. Compared with a traditional displacement finite element method, the method is more targeted when solving the stress problem of the particle composite reinforced material.

Description

Technical Field

[0001] The invention belongs to the technical field of material simulation analysis, and in particular relates to a three-dimensional Voronoi unit simulation method for the elasticity problem of a composite material reinforced with interface phase particles. Background Art

[0002] Particle-reinforced composites have attracted much attention due to their wide applications in automotive engineering, aerospace, military, electronics, and nuclear engineering. As the particle size of materials continues to decrease with technological advances, the influence of the interface layer - a thin layer between the matrix and the inclusions - has become increasingly significant. The formation mechanisms of the interface are diverse, including chemical reactions between the reinforcing particles and the matrix, and surface coatings of the particles, all of which have a profound impact on the mechanical properties of the composite. Therefore, in-depth research on the interface characteristics is crucial to accurately predict the overall performance of the composite.

[0003] In the study of particle reinforced composites, analytical methods and finite element method (FEM) are two commonly used methods. The homogenization theory proposed by K. Ding et al. established an elastic-plastic framework for composites with ductile interface characteristics, revealing the significant influence of interface ductility on the performance of composites. Carrère et al. used FEM to analyze the influence of the interface on the deflection of matrix cracks in ceramic matrix composites, verifying the key role of the interface in the crack propagation mechanism. Gu, JH et al. explored the damping performance of particle reinforced metal matrix composites through the FEM model and found that the elastic modulus of the interface phase has a significant effect on damping;

[0004] However, traditional FEM has limitations in dealing with constraint problems, singularity problems and mesh distortion, resulting in high computational cost and limited accuracy. To solve these problems, researchers have developed a variety of improved methods. For example, the Trefftz computational particle (TCGS) proposed by LTDong and Atluri improves the accuracy and efficiency of the calculation by satisfying the displacement field of the control differential equation. Fan, JP used FEM to study the influence of interface performance on Young's modulus and established a three-dimensional cell model to optimize the interface design. Hassanzadeh-Aghdam MK explored the influence of interface layer thickness and stiffness on the mechanical properties of composites and found that the interface effect weakened with the increase of filler aspect ratio. In addition, the Voronoi cell finite element method (VCFEM) has attracted attention in the study of particle reinforced composites due to its advantages such as simple unit division and fast calculation speed. This method can accurately reflect the random shape, size and spatial distribution of inclusion particles, and is suitable for analyzing models with a large number of inclusions. R. Guo and R. Zhang proposed a modified complementary energy functional based on VCFEM to study the microscopic and macroscopic mechanical behaviors of dissimilar materials. WYHao et al. used VCFEM to analyze the plasticity, creep and thermal strain problems of homogeneous materials and verified the effectiveness of VCFEM.

[0005] Although these methods have made significant progress in the study of particle-reinforced composites, challenges remain. In particular, for three-dimensional Voronoi cells with complex interfaces, there is currently a lack of effective and accurate numerical methods. As a result, the computational efficiency for large-scale complex microstructures is low, which provides an effective solution to larger-scale microstructural problems. Summary of the invention

[0006] In order to solve the above technical problems, the present invention designs a three-dimensional Voronoi unit simulation method for the elastic problem of composite materials reinforced with interface phase particles, which regards stress and displacement as independent fields. Each Voronoi unit is divided into multiple Delaunay tetrahedrons, and the Hammer integral is used for integral solution. By interpolating the surface of each tetrahedron, the internal stress and strain of the tetrahedron can be obtained. Subsequently, the displacement, stress and strain of the internal points of all units can be obtained. Compared with the finite element software ABAQUS, VCFEM can achieve the calculation accuracy of hundreds of thousands of units of ABAQUS with fewer units, and has stronger pertinence to the problem of three-dimensional composite reinforced materials containing interface particles; it effectively improves the calculation efficiency of complex microstructures and provides an effective solution for larger-scale microstructure problems.

[0007] In order to achieve the above technical effects, on the one hand, the present invention discloses a three-dimensional Voronoi unit simulation method for the elastic problem of a composite material reinforced with particles containing an interface phase, comprising the following steps:

[0008] S1: Establish a model based on the actual problem and discretize the model into interconnected Voronoi cells;

[0009] S2: Divide each unit into multiple Delaunay tetrahedrons;

[0010] S3: Construct units based on the minimum complementary energy principle to modify the complementary energy functional;

[0011] S4: Use high-order perfect polynomials as unit stress fields and construct suitable trial functions;

[0012] S5: Obtain the element stiffness matrix by solving the G matrix and the H matrix;

[0013] S6: According to the relationship between the external node displacement and the stiffness matrix, the node stress and internal node displacement are obtained;

[0014] S7: Calculate β based on internal node displacements;

[0015] S8: By taking partial derivatives of the stress function, we obtain the coefficient matrix ρ and the final result of stress, (σ=ρβ).

[0016] Furthermore, the process of correcting the complementary energy functional in S3 comprises the following steps:

[0017] S3.1: Construct a three-dimensional phase interface Voronoi unit based on the existing assumed stress hybrid element method equilibrium model I and equilibrium model II; the unit consists of three phases, including matrix, interphase and inclusions, on a given traction boundary:

[0018]

[0019] At the boundaries of inclusions and at interfaces within cells:

[0020] n c ·σ i -n c ·σ c =0#(2)

[0021] On the internal interfaces of the element and on the boundaries of the matrix:

[0022] n i ·σ m -n i ·σ i =0#(3)

[0023] S3.2: The modified complementary energy functional is obtained by introducing the Lagrange multiplier: It can be specifically expressed as:

[0024]

[0025] Γ is the surface force boundary; Ω m is the integration area of ​​the matrix, Ω i is the integral area of ​​the interface phase, Ω c is the integral region of inclusion;

[0026] S3.3: Introducing self-balanced stress field in the element in the finite element:

[0027] σ=pβ#(5)

[0028] σ is a strong vector with 6 stress components, β is a column vector with m unknown stress coefficients, P is a 3×m matrix; the boundary displacement u is the interpolation of the generalized displacement d of the node:

[0029] u=Ld#(6)

[0030] Substituting (5) and (6) into (4), we obtain:

[0031]

[0032] After discretization, the weak form of the system complementary energy can be obtained, and according to the modified complementary energy stationary value condition:

[0033]

[0034] in,

[0035] H m =∫ Ω P mT S m P m dΩ H i =∫ Ω P iT s m P i dΩ H c =∫ Ω P cT S m P c dΩ#(10)

[0036] G mm =∫ Γ P mT n mT LdG ii =∫ Γ p iT n iT LdG cc =∫ Γ P cT n cT LdΓ#(11

[0037] G im =∫Γ P iT n mT LdG ci =∫ Γ P cT n iT LdΓ#(12)

[0038] The above formula can be written as the displacement expression of stress parameter:

[0039] β=H -1 Gd#(13)

[0040] The total energy of the system Π e The first-order variation of the node displacement d can give the corresponding weak expression of the traction boundary condition:

[0041]

[0042] Substitute (14) into (13),

[0043]

[0044] The equations for solving the generalized displacement are obtained, where the element stiffness matrix is:

[0045]

[0046] Furthermore, in S4, the construction process of the stress function is:

[0047] S4.1: The matrix stress function is used The inclusion stress function is expressed as The interface stress function is expressed as The interaction stress function of the matrix is ​​expressed as The interface interaction force function and stress function are expressed as express;

[0048] S4.2: The interfacial interaction stress function obtained by Ghosh is further simplified to obtain the following expression:

[0049]

[0050] In the formula, β m ,β i and β c are column vectors of stress function coefficients of matrix, interface and inclusion respectively; α is the ellipsoid equation, which is a known quantity derived from the coordinates of the model nodes; α 1 is the description boundary Γ c The parameters of the ellipsoid, and α 2 is the description boundary Γ i Parameters of the ellipsoid;

[0051] S4.3: Near the interface Away from the interface In the interface After constructing the stress polynomials of the matrix and inclusions, the stress field is derived from the Maxwell stress function. Substituting the matrix stress polynomial into the matrix, interface and inclusion stress fields is as follows:

[0052]

[0053]

[0054]

[0055] in are scaled (x,y,z) coordinates, and when 1 / α i =1, β pqri The number of terms is equal to the number of terms in the polynomial stress function;

[0056] S4.4: Select the most widely used Maxwell stress function. When assuming that the stress is a complete polynomial, the stiffness matrix must be reversible. According to the relationship between the number of terms in the stress function and the number of degrees of freedom of the rigid body, select d>bc, where d represents the total number of terms in the stress function, b represents the total number of degrees of freedom of the rigid body, and c represents the degrees of freedom of the rigid body. The basic number of terms in the Maxwell stress function used to calculate the three-dimensional arbitrary polyhedral unit is at least 31. Introduce the Maxwell stress function that satisfies the self-balance in the unit:

[0057]

[0058] According to this stress function, the six Maxwell stress components can be obtained:

[0059]

[0060] get

[0061] [σ x σ y σ z τ xy τ yz τ zx ]=Pβ#(222)

[0062] P is a coefficient matrix with 6 rows and m columns, which is only related to the coordinates of the points. The β matrix is ​​an m-dimensional column vector, where m represents the number of terms selected by the polynomial.

[0063] Furthermore, the division process of the Delaunay tetrahedron in S2 is:

[0064] S2.1: Input unit node information (node ​​coordinates, node number); input surface information (surface node information, surface quantity);

[0065] S2.2: Determine which phase the nodes and faces belong to based on the node information, and construct Delaunay triangles in space based on the node numbers;

[0066] S2.3: Connect Delaunay triangles to construct tetrahedrons;

[0067] S2.4: Determine whether the minimum circumscribed sphere of the tetrahedron generated by the four points satisfies the empty circle property. If so, the Delaunay tetrahedron partitioning is completed. If not, repeat S2.3 and subsequent steps until it is satisfied.

[0068] Furthermore, in S5, the stiffness matrix of the unit obtained by solving is:

[0069]

[0070] You can get H * =H(λ * ) 3 andG * =G(λ * ) 2 , then there will be an additional unit scaling factor when calculating the stiffness matrix, namely:

[0071]

[0072] The beneficial effects of the present invention are:

[0073] The present invention designs a three-dimensional Voronoi unit simulation method for the elastic problem of particle-reinforced composite materials containing interface phases. Based on the hybrid finite element method and the variational principle, a three-dimensional modified complementary energy functional including inclusions, matrix and interface is derived, so that the method can be used to calculate three-dimensional three-phase composite materials;

[0074] In addition, the present invention constructs stress functions in three integral domains respectively, and adopts Maxwell complete polynomials, including stress function polynomials and interaction force function polynomials, which have good convergence and continuity for three-dimensional interface units with different numbers of nodes. Compared with commercial software, this method is faster and more efficient in calculating three-dimensional interface Voronoi models. When calculating the same model, the results of hundreds of thousands of units in commercial software can be obtained with fewer units. Compared with the traditional displacement finite element method, this method is more targeted in solving stress problems of particle composite reinforced materials; and the calculation efficiency of large-scale complex microstructures is higher than that of traditional finite element analysis, and the calculation accuracy of hundreds of thousands of units can be achieved with fewer units, which is more targeted for three-dimensional composite reinforced material problems containing interface particles. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for describing the embodiments are briefly introduced below.

[0076] Figure 1 is a Voronoi model diagram with an interface of the present invention;

[0077] Figure 2 Schematic diagram of the hollow circle characteristic of the present invention, in which a satisfies the hollow circle characteristic diagram, and b satisfies the non-hollow circle characteristic diagram;

[0078] Figure 3 is an example diagram of tetrahedron partitioning of the present invention;

[0079] Figure 4 It is a schematic diagram of boundary conditions of the present invention;

[0080] Figure 5 Schematic diagram of the polygonal approximate ellipsoid of the present invention, wherein (a) is a polygonal approximate ellipsoid; (b) single layer 6 nodes (c) single layer 12 nodes (d) divided into 24 layers;

[0081] Figure 6 1 is a mesh division diagram of the present invention, in which (a) is ABAQUS spherical inclusion (361105 units; (b) VCFEM 6-node polygonal inclusion (1 unit); (c) VCFEM 12-node polygonal inclusion (1 unit); (d) VCFEM 6-node no interface (1 unit);

[0082] Figure 7It is a comparison diagram of the calculation results of ABAQUS and VCFEM of the present invention, in which (a) is ABAQUS spherical inclusion (361105 units); (b) VCFEM 6-node polygonal inclusion (1 unit); (c) VCFEM 12-node polygonal inclusion (1 unit); (d) VCFEM 6-node no interface (1 unit); (e) VCFEM 6-node 0.65mm interface (1 unit);

[0083] Figure 8 1 is a stress curve diagram of three models and extraction positions of the present invention, wherein (a) is for effectiveness verification and (b) is for interface influence verification;

[0084] Fig. 9 is a boundary condition diagram of the present invention;

[0085] Fig.10 It is the multi-element model establishment and boundary conditions of the present invention, wherein (a) is ABAQUS spherical inclusion; (b) VCFEM 6-node polygonal inclusion; (c) VCFEM 12-node polygonal inclusion;

[0086] Fig.11 : is the mesh division diagram of the present invention, wherein (a) is ABAQUS spherical inclusion (1708174 units); (b) VCFEM 6-node polygonal inclusion (27 units); (c) VCFEM 12-node polygonal inclusion (27 units);

[0087] Fig.12 It is the multi-element stress cloud diagram of the present invention, where (a) is the ABAQUS spherical inclusion result; (b) is the VCFEM 6-node polygonal inclusion result; (c) VCFEM 12-node polygonal inclusion;

[0088] Fig.13 are stress curves of three models and pull-out positions of the multi-unit of the present invention;

[0089] Fig.14 It is the model establishment and boundary condition diagram when the unit of the present invention is expanded to 4 through array, wherein (a) is ABAQUS, and (b) is VCFEM;

[0090] Fig.15 The unit of the present invention is meshed when the array is expanded to 4, where (a) is ABAQUS spherical inclusion (512423); (b) is VCFEM 6-node polygonal inclusion (4 elements);

[0091] Fig.16 are stress cloud diagrams of two calculation results of the present invention;

[0092] Fig.17 is a diagram of stress conditions on the path of the present invention;

[0093] Fig.18 is a random distribution model and boundary condition diagram of the present invention;

[0094] Fig.19 It is the stress cloud diagram of the random distribution model of the present invention. DETAILED DESCRIPTION

[0095] The present invention discloses a three-dimensional Voronoi unit simulation method for elastic problems of composite materials reinforced with interface phase particles.

[0096] According to the hypothetical stress hybrid element method created by Mr. Bian Xueqi, the equilibrium model in the hybrid element method can be divided into equilibrium model I and equilibrium model II. The equilibrium model I hybrid element is a finite element model based on the minimum complementary energy principle and using the stress function as the field variable. The expression of the complementary energy functional is:

[0097]

[0098] Where S is the elastic flexibility tensor, and Γ is the given displacement boundary. The stress of the unit is expressed by the stress function, with the node stress function as the basic unknown quantity. The stress function is guaranteed to be continuous on the boundary between units, so the equilibrium condition can be met everywhere. The disadvantage of equilibrium model I is that the distribution of stress and strain is obtained by differentiating the stress function, and the displacement solution can only be obtained by integrating the geometric equation, which depends on the selection of the integral path. Therefore, the displacement solution obtained in this way is not unique. In addition, the stress function is not convenient for intuitively explaining the physical meaning of the boundary conditions.

[0099] Equilibrium model II was proposed by Fraeijs in 1964 and is used to analyze many plate bending problems. This model is based on the traditional complementary energy functional and uses the Lagrange multiplier method to introduce the equilibrium condition of the inter-element surface force continuity into the functional.

[0100] For the original problem, the traction boundary conditions are given:

[0101] n is the normal of the boundary, and the equilibrium condition of the surface forces between the elements is as follows:

[0102] n + ·σ+n - ·σ=0#(2-3)

[0103] n + is the outer normal of the boundary, n - is the internal normal of the boundary. By introducing the above two constraints through Lagrange multipliers, the modified functional can be obtained. And the boundary surface force can be written as T = n·σ. The given displacement also satisfies the unit displacement interpolation, then This function is usually written as:

[0104]

[0105] It can be seen that the stress hybridization model introduces two unknown variables inside the unit and at the unit boundary, namely the stress field variable σ and the boundary displacement variable u. The unit node displacement interpolation can satisfy all edge coordination conditions of the entire unit. The advantage of the stress hybrid element method for calculating boundary displacements is that it is not limited by the number of unit edges. In the study of two-dimensional Voronoi units, Guo proposed a new modified complementary energy functional for interface, delamination and plasticity problems, which provided a basis for the interface problems of three-dimensional Voronoi units. Based on this principle, a three-dimensional phase interface Voronoi unit is constructed. The Voronoi model with interface used in the present invention is as follows Figure 1 shown.

[0106] The unit cell consists of three phases, including matrix, interphase and inclusions. On a given traction boundary:

[0107]

[0108] At the boundaries of inclusions and at interfaces within cells:

[0109] n c ·σ i -n c ·σ c =0#(2-6)

[0110] On the internal interface of the unit and the boundary of the matrix: n i ·σ m -n i ·σ i =0#(2-7

[0111] The modified complementary energy functional is obtained by introducing the Lagrange multiplier:

[0112]

[0113] In (2-8)Γ is the surface force boundary. Ω m is the integration area of ​​the matrix, Ω i is the integral area of ​​the interface phase, Ω c is the inclusion integral region.

[0114] Introducing a self-balanced stress field within the element into the finite element:

[0115] σ=pβ#(2-9)

[0116] σ is a stress vector with 6 stress components, β is a column vector with m unknown stress coefficients, and P is a 3×m matrix.

[0117] The boundary displacement u is the interpolation of the node generalized displacement d:

[0118] u=Ld#(2-10)

[0119] Substituting (2-9) and (2-10) into (2-8), we obtain:

[0120]

[0121]

[0122] After discretization, the weak form of the system complementary energy can be obtained. According to the modified complementary energy stationary value condition:

[0123]

[0124] Weak expressions for the kinematic relationships are obtained in each element.

[0125] in,

[0126] H m =∫ Ω P mT S m P m dΩ H i =∫ Ω P iT S m P i dΩ H c =∫ Ω P cT S m P c dΩ#(2-14)

[0127] G mm =∫ Γ P mT n mT LdG ii =∫ Γ P iT n iT LdG cc =∫ Γ P cT n cT LdΓ#(2-15

[0128] G im =∫ Γ P iT n mT LdG ci =∫ Γ PcT n iT LdΓ#(2-16)

[0129] The above formula can be written as the displacement expression of stress parameter:

[0130] β=H -1 Gd#(2-17)

[0131] The total energy of the system Π e The first-order variation of the node displacement d can give the corresponding weak expression of the traction boundary condition:

[0132]

[0133] Bringing (2-17) into (2-18)

[0134]

[0135] The system of equations for solving the generalized displacement is obtained, where the element stiffness matrix is

[0136]

[0137] Equation (2-8) shows that the unit displacement variable appears on all boundaries of the unit, so the displacement mode of the unit can be interpolated to obtain the displacement of any point. In order to continue the derivation of the equation, it is necessary to obtain the expression of the stress function, the normal vector of the unit surface and the cosine formula. Next, the calculation and solution methods of the relevant variables are introduced.

[0138] In three-dimensional problems, it is necessary to choose a reasonable stress function to accurately reflect the stress field of the unit. The inclusion stress function is expressed as The interface stress function is expressed as The interaction stress function of the matrix is ​​expressed as The interface interaction force function and stress function are expressed as The construction of polynomial stress function and interaction force function should meet the following requirements: and 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 inclusion. The interaction stress term approaches zero at places far away from the contact boundary of different phases, and the continuity condition of the surface force must be met at the contact boundary of different phases.

[0139]

[0140] β m ,β i and β care column vectors of stress function coefficients of matrix, interface and inclusion respectively. In this paper, the interface interaction stress function obtained by Ghosh is further simplified and expressed as follows:

[0141]

[0142] α is the equation of the ellipsoid, which is a known quantity derived from the coordinates of the model nodes. 1 is the description boundary Γ c The parameters of the ellipsoid, and α 2 is the description boundary Γ i Parameters of the ellipsoid. Figure 1 As shown. Near the interface Away from the interface In the interface After constructing the stress polynomials of the matrix and inclusions, the stress field is derived from the Maxwell stress function. Substituting the matrix stress polynomial into the matrix, interface and inclusion stress fields is as follows:

[0143]

[0144]

[0145]

[0146] Equations 2-23 and 2-24 have the same form. In the model, they correspond to the stress fields of the matrix and interface phases, respectively. Due to the different parameters α of the two ellipsoids and the different choices of the stress function polynomials, the column vector β is unknown. m and β i The final stress field is different due to the different number of terms. are scaled (x,y,z) coordinates, and when 1 / α i =1, β pqri The number of terms is equal to the number of terms of the polynomial stress function. This paper selects the most widely used Maxwell stress function. When assuming a complete stress polynomial, the stiffness matrix must be reversible. According to the relationship between the number of terms in the stress function and the number of degrees of freedom of the rigid body, d>bc can be selected, where d represents the total number of terms in the stress function, b represents the total number of degrees of freedom of the rigid body, and c represents the degrees of freedom of the rigid body. The basic number of terms of the Maxwell stress function used to calculate arbitrary three-dimensional polyhedral units is at least 31.

[0147] Introduce the Maxwell stress function that satisfies the self-equilibrium in the unit:

[0148]

[0149] This function is the sum of complete polynomials of all orders. We know that there is 1 term in the 0th order, 3 terms in the 1st order, and 6 terms in the 2nd order. Therefore, the number of terms in the mth (m>0)th order complete polynomial is Then the total number of terms in the stress function Φ is Excluding the 0th and 1st order terms, there are 31 terms in total in equation (2-21). Based on this stress function, the six Maxwell stress components can be obtained:

[0150]

[0151] Then you can get

[0152] [σ x σ y σ z τ xy τ yz τ zx ]=Pβ#(2-28)

[0153] P is a coefficient matrix with 6 rows and m columns, which is only related to the coordinates of the points. The β matrix is ​​an m-dimensional column vector, called the stress coefficient. m represents the number of terms in the polynomial selection.

[0154] like Figure 2 As shown in a, in the circumscribed circle of each triangle, there are only three points of the triangle and no other points. If there are other points in the circumscribed circle, the empty circle property is not satisfied, such as Figure 2 b. In three-dimensional problems, the circumscribed sphere of each tetrahedron contains only four points of the tetrahedron and no other points. For example, Figure 3 The hexahedron in is divided into tetrahedrons. After constructing Delaunay triangles on each face, find the fourth point that satisfies the empty circle feature. The tetrahedron is split successfully, and each tetrahedron satisfies the empty circle feature.

[0155] Next, we calculate the displacement interpolation of the tetrahedron surface, considering Figure 2 A space triangle is numbered as i, m, j. In the rectangular coordinate system, the coordinates of each node are (x i ,y i ,z i ),(x m ,y m ,z m ),(x j ,y j ,z j ) As shown in Equation 2-29, the boundary displacement of the unit surface can be obtained by interpolation of the node displacement.

[0156]

[0157] In the formula Δ is the area of ​​the triangle, which can be calculated based on the coordinates of the three points. i is the area of ​​the small triangle formed by the two corner points on the sides and the internal point p. Rewrite equation (2-29) into matrix form: u = Ld#(2-30)

[0158] in

[0159] d=[uvw] T =[u i u j u m v i v j v m w i w j w m ] T #(2-31)

[0160]

[0161] Any point P(x,y,z) inside the triangle and the vertices of the triangle line divide the large triangle into three smaller triangles, each with an area of ​​Δ i ,Δ j and Δ m , the area of ​​the large triangle and the area of ​​the small triangle Δ can both be determined by the node coordinates.

[0162]

[0163] You can get:

[0164]

[0165] From the properties of the isoparametric element, we can also obtain the coordinates of any point on the triangular surface:

[0166]

[0167] When solving the G matrix, you also need to solve the cosine of the normal vector of each face of all tetrahedrons separately to form the normal vector matrix of the triangle. The normal vector formula for each face of the triangle can be obtained here:

[0168]

[0169] Then the direction cosine matrix of the triangle is:

[0170]

[0171] where n i , n j n m It is expressed as:

[0172]

[0173] In order to solve the problem of irreversibility of the matrix caused by too large or too small values ​​during the calculation process and avoid excessive errors in the calculation results, it is necessary to scale the unit to obtain local coordinates for numerical calculations. Therefore, the coordinates of the center of mass of the unit should be obtained, and the unit should be scaled with the center of mass as the origin of the local coordinate system. However, the shape of the polyhedral unit is complex and the number of faces is uncertain, so its center of mass cannot be directly determined. The method given in this paper is to divide the unit into Delaunay tetrahedrons in the Cartesian coordinate system, 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.

[0174] According to the tetrahedron center coordinate formula, the center coordinates of each small tetrahedron can be obtained. Here, it is assumed that the coordinates of the four vertices of the tetrahedron are

[0175] p 1 (x 1 ,y 1 ,z 1 ),p 2 (x 2 ,y 2 ,z 2 ),p 3 (x 3 ,y 3 ,z 3 ),p 4 (x 4 ,y 4 ,z 4 )#(2-39)

[0176] Assume that an arbitrary polyhedron is divided into n tetrahedrons, and their center coordinates are:

[0177]

[0178] The volumes are:

[0179]

[0180] Then we can get the center coordinates of any three-dimensional polyhedron

[0181]

[0182] Once the center coordinates of the polyhedron unit are obtained, scaling can be performed. The center coordinates of the polyhedron unit are used as the origin of the local coordinates, also called the scaling center. The scaled local coordinates of any point (x, y, z) on the polyhedron can be obtained.

[0183]

[0184] Since the stress function is closely related to the element coordinates, the stress function can be written in terms of scaled local coordinates as Related complete polynomial functions.

[0185]

[0186] According to the functional of the three-dimensional Voronoi unit with interface, it can be known that calculating the stiffness matrix of the unit is actually calculating the H matrix and the G matrix. The calculation formula is:

[0187]

[0188] is the integration area Ω e The boundary of this method mainly uses Hammer integral. The Hammer integral formula is:

[0189]

[0190] Where ω k is the weight coefficient of the Hammer integration point. The integral of the H matrix is ​​integrated in the tetrahedron domain, so the integral accumulation of each small tetrahedron can get the result of the entire unit, that is:

[0191]

[0192] The local coordinates after scaling are:

[0193]

[0194] Where n tet Indicates the number of tetrahedrons obtained by partitioning, n int represents the number of tetrahedral Hammer integration points used, ω j Represents the weight coefficient of the jth integration point. i represents the Jacobian of the i-th tetrahedron. Specifically:

[0195]

[0196] P(x,y,z) is a coefficient matrix related only to the coordinates of the integration point. Similarly, the G matrix needs to be accumulated when it is integrated on the triangular surface. That is:

[0197]

[0198] It is expressed in scaled local coordinates as:

[0199]

[0200] You can get H * =H(λ * ) 3 and G * =G(λ * ) 2 , then there will be an additional unit scaling factor when calculating the stiffness matrix, namely:

[0201]

[0202] Example 2

[0203] The simulation method constructed in Example 1 is further verified by comparing VCFEM with a single unit of commercial finite element analysis software; the size of the model and the radius of the internal structure (ellipsoidal inclusions and interfaces) are determined. Then, according to the stratification of inclusions and interfaces and the number of points required for each layer, the information of all points in VCFEM can be obtained, and the continuity conditions of surface forces are met at the interfaces of different materials. Both models use displacement constraints, and in the same coordinate system, the two models apply the same boundary conditions to the nodes on the same surface. The stress distribution of the two models in tension is calculated using the static structural analysis method.

[0204] ABAQUS model establishment: Generate a cube entity with a side length of 20mm, and then create two spherical entities with radii of 2mm and 2.35mm respectively. Move the entities so that their centers of mass coincide. Through Boolean operations, you can obtain entities that maintain their contours, such as Figure 1 The material parameters specified in Table 1 are applied to the corresponding areas. Multiple materials coexist in a single entity, which is consistent with the characteristics of composite materials.

[0205] Table 1 Material parameters

[0206]

[0207] Therefore, there is no need to establish additional contact or connection conditions between different material phases. The entire solid is then meshed and displacement boundary conditions are applied, such as Figure 4 The stress results can then be obtained through static structural analysis.

[0208] VCFEM modeling: 3D Voronoi cells with interfaces such as Figure 4 The unit consists of three parts, namely matrix phase, inclusion phase and interface phase. The unit side length is 20 mm. The inclusion and interface are approximated by the inscribed polygon of the ellipsoid, such as Figure 5As shown, the radius of the inclusion ellipsoid is 2mm, the interface thickness is 0.35mm, and the corresponding material properties are shown in Table 1. In order to study the influence of the interface layer on the composite material, we set up an unbounded model and a model with an interface layer thickness of 0.65mm. The other model parameters are consistent with the VCFEM six-node model. The boundary conditions are as follows: the displacement of the X, Y, and Z surfaces outside the matrix in the corresponding X, Y, and Z directions is 0mm, and the displacement of the Y surface in the Y direction is 0.02mm. The specific model is as follows Figure 4 As shown:

[0209] Table 3-1 Material parameters

[0210]

[0211] One of the advantages of the method proposed in this paper is that the results of tens of thousands of units of the traditional displacement finite element method can be obtained with a very small number of units. Figure 5 As shown:

[0212] Computational models such as Figure 6 As shown in Figure 2, the VCFEM model has only one Voronoi unit with an interface, and the ABAQUS model consists of 361,105 4-node tetrahedral units. The VCFEM6-node model and the VCFEM12-node model have 8 external nodes, 48 ​​internal nodes (inclusion nodes and interface nodes), and 96 internal nodes (inclusion nodes and interface nodes), respectively. The number of stress function items is shown in Table 2:

[0213] Table 2 Stress function item number table

[0214]

[0215] The calculation results are compared and verified, and the result cloud diagram is as follows Figure 7 As shown. This paper compares the calculation results of ABAQUS and VCFEM, and verifies the convergence, accuracy and feasibility of polygonal modeling of the calculation results of this method. The stress data on the vertical middle path are extracted from the cloud map. The stress curves of the three models and the pull-out position are shown in Figure 8 As shown;

[0216] In the study of three-dimensional composite materials with interfacial phases, it is pointed out that the exact solution of uniform cubic deformation can be expressed as:

[0217]

[0218] u 1 ,u 2 , and u 3 is the displacement of the node in the X, Y and Z directions, p is the total stress at the node, v is Poisson's ratio, E is the elastic modulus of the material, x 1 ,x 2 , and x3 are the coordinates of the nodes before the model is deformed.

[0219] The three material parameters of the VCFEM 6-node polygonal inclusion model (1 element) were modified to 72000 MPa Young's modulus and 0.33 Poisson's ratio. The boundary conditions are as follows Figure 3-6 As shown, nodes 2, 4, 6 and 8 apply a displacement boundary condition of 0.01 mm in the Y direction, and nodes 1, 3, 5 and 7 constrain the Y direction. Substitute the stress in the calculation result into formula 3-1 to verify whether the calculation of stress and node displacement in this method is correct. 100 sampling points are selected in the plane of nodes 2, 4, 6, and 8. The verification results are shown in Table 3.

[0220] Table 3 Node displacement and error

[0221]

[0222]

[0223] Among the 100 sampling points, the error is between 0% and 1.25%, and the average node displacement is 0.010068503mm, which is very close to the 0.01mm given by us. This shows that the method proposed in this paper is correct in the calculation of composite materials containing interface phases.

[0224] Figure 7 and Figure 8 The normal stress in the Y direction for the ABAQUS spherical model and the VCFEM polygonal model under uniaxial tension in the Y direction is shown. The results obtained by the stress program and the stress path curve show a high degree of consistency, indicating that the proposed method is able to obtain a computational accuracy comparable to that of the traditional displacement finite element method using hundreds of thousands of units using a single unit. Although this method uses fewer points to approximate the ellipsoid, it still treats the polygon as an ellipsoid when processing the stress function. The computational results further show that the polygon can achieve the expected results of the ellipsoid model. By increasing the number of points within the polygon, results closer to the true solution can be obtained. The analysis of the model without an interface shows that the presence of the interface significantly affects the stress of the entire model. In particular, its absence leads to larger boundary stress fluctuations, which significantly affects the material damage. In actual composite materials, the interface phase is usually a coating with a thickness of 1 / 6 of the inclusion radius, and occasionally there is an artificially added interface phase with a thicker thickness. The stress results of the model with an interface layer thickness of 0.65 mm are compared with the stress results of the model without an interface layer. The maximum stress differs by 15% and the minimum stress differs by 5%. Under the same boundary conditions, as the thickness of the interface phase in the model increases, the overall stress value of the model also increases.

[0225] Example 3

[0226] In this embodiment, the VCFEM established by this method is compared and verified with multiple units of commercial finite element analysis software. The specific verification process is as follows:

[0227] After verifying the stress results of a single Voronoi unit, this paper will verify multiple units. Using an array, the unit is expanded to 27 units, the unit size remains unchanged, and the material properties are the same as in Table 1. The boundary conditions are: the displacements of the X, Y, and Z surfaces on the outside of the matrix are fixed to 0 mm in the corresponding X, Y, and Z directions, and the displacement of the Y surface in the Y direction is given to be 0.02 mm. The specific model is as follows Fig.10 As shown: The grid division of the calculation model is as follows Fig.11 As shown in Table 4, the VCFEM model has only 27 Voronoi units with interfaces, and the ABAQUS model consists of 1708174 4-node tetrahedral units. The VCFEM6-node model and the VCFEM12-node model have 64 external nodes, 1296 internal nodes (inclusion nodes and interface nodes) and 2592 internal nodes, respectively. The number of stress function items is shown in Table 4, where pol: stress function, rec: interaction force function:

[0228] Table 4. Number of terms in stress function

[0229]

[0230] Comparison and verification results cloud diagram as follows Fig.12 As shown in Figure 1, the stress data on the vertical centerline path are extracted from the program. The stress data on the vertical middle path are extracted from the cloud map. The stress curves under the three models and the pull-out position are shown in Figure 1. Fig.13 As shown, Fig.12 and Fig.13 As shown in the figure, the stress distribution in the Y direction of the ABAQUS spherical model and the VCFEM polygonal model under uniaxial tension along the Y axis is given. The calculation results of the stress program show a high degree of consistency with the stress path curve, indicating that the method proposed in this paper can obtain a calculation accuracy comparable to that of the traditional displacement finite element method (more than 1 million units) using only dozens of units.

[0231] Example 4

[0232] In this embodiment, complete polynomials with different numbers of terms are selected to verify the influence of the number of stress function terms on the results;

[0233] In this embodiment, the unit is expanded to 4 through the array, and the unit size and material properties are the same as in Table 1. The boundary conditions are as follows: the displacement of the three fixed X, Y, and Z edges outside the matrix in the corresponding X, Y, and Z directions is 0 mm, and the displacement of the given Y direction on the Y plane is 0.02 mm. The specific model and mesh division are as follows Fig.14 , Fig.15As shown;

[0234] Comparison and verification of calculation results Table 5 lists the number of stress items. Fig.16 and Fig.17 The calculation results are shown;

[0235] Table 5 Stress item number table

[0236]

[0237] In the table, pol: stress function, rec: interaction force function.

[0238] Select the stress on a path on the cloud map for analysis. The path is as follows: Fig.16 As shown in Figure 2, the stress on the path is as follows: Fig.17 As shown; From the stress cloud diagram and stress path curve diagram, it can be seen that the results of the five examples are relatively close, and very similar solutions can be obtained using complete polynomials with different numbers of terms. It can be concluded that when the number of stress terms reaches a certain number, the calculation results tend to be stable, and this method can produce correct results within a certain range of stress terms; the relationship between the number of terms in the stress function and the number of degrees of freedom of the rigid body is selected as d>bc.

[0239] Example 5

[0240] In actual particulate composite reinforced materials, in most cases, the distribution of inclusions is random. In order to simulate the actual situation, a random distribution model is used in this embodiment. A 6-node polygon model is adopted, in which 20 inclusions are randomly distributed in a regular hexahedron with a size of 20mm*20mm*20mm, and an interface phase with uniform thickness is distributed between the matrix and the inclusions. The number of Voronoi units containing the interface is 20, the number of external nodes is 105, the number of internal nodes is 960, and the boundary conditions are: the displacement of the fixed X, Y, and Z edges in the corresponding X, Y, and Z directions is 0mm outside the matrix, and the displacement of the Y edge in the given Y direction is 0.02mm. The boundary conditions and mesh division are as follows Fig.18 The number of stress items is shown in Table 6, and the unit node information is shown in Table 7 and Table 8:

[0241] Table 6 Stress item numbers for random distribution verification

[0242]

[0243] Table 7 Unit node information

[0244]

[0245] Table 8 Unit node information

[0246]

[0247]

[0248] The stress cloud diagram of the calculation results is as follows Fig.19 As shown in the figure, the stress cloud diagram shows that even under random distribution conditions, the number of nodes, faces and boundary conditions of each unit is different. However, the preset complete polynomial can still effectively capture the complex details of the stress distribution within the unit. In this model, the unit has a maximum of 48 internal nodes and 24 external nodes. The total number of degrees of freedom of these units is calculated as (48+24)×3=216. Considering the elimination of 12 rigid body displacements, the final cumulative total number of degrees of freedom is 228. The number of terms in the stress function is set to 600. After multiple evaluations on the selection of stress function terms, we selected the most ideal calculation results for display. When selecting the complete polynomial, the number of terms must exceed the sum of the degrees of freedom of all nodes of each unit plus the degrees of freedom corresponding to twelve rigid bodies. When a certain number is reached, the same stress function shows good convergence on different units, and the calculation results tend to be stable; it can be effectively shown that this method provides a reference for dealing with models that cannot be calculated using commercial software or verified by experimental means, and provides a new basis for studying composite materials with random distribution structures.

[0249] The preferred embodiments of the present invention disclosed above are only used to help explain the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to only the specific implementation methods described.

Claims

1. A three-dimensional Voronoi element simulation method for the elastic problem of composite materials reinforced with interfacial phase particles, characterized in that: The following steps are involved: S1: Establish a model based on the actual problem and discretize the model into interconnected Voronoi cells; S2: Divide each unit into multiple Delaunay tetrahedrons; S3: Construct units based on the minimum complementary energy principle to modify the complementary energy functional; S4: Use high-order perfect polynomials as unit stress fields and construct suitable trial functions; S5: Obtain the element stiffness matrix by solving the G matrix and the H matrix; S6: According to the relationship between the external node displacement and the stiffness matrix, the node stress and internal node displacement are obtained; S7: Calculate β based on internal node displacements; S8: By taking partial derivatives of the stress function, we obtain the coefficient matrix ρ and the final result of stress, (σ=ρβ).

2. A three-dimensional Voronoi unit simulation method for elastic problems of composite materials with interfacial particles according to claim 1, characterized in that: The process of correcting the complementary energy functional in S3 comprises the following steps: S3.1: Construct a three-dimensional phase interface Voronoi unit based on the existing assumed stress hybrid element method equilibrium model I and equilibrium model II; the unit consists of three phases, including matrix, interphase and inclusions, on a given traction boundary: At the boundaries of inclusions and at interfaces within cells: n c ·s i -n c ·s c =0#(2) On the internal interfaces of the element and on the boundaries of the matrix: n i ·s m -n i ·s i =0#(3) S3.2: The modified complementary energy functional is obtained by introducing the Lagrange multiplier: It can be specifically expressed as: Γ is the surface force boundary; Ω m is the integration area of ​​the matrix, Ω i is the integral area of ​​the interface phase, Ω c is the integral region of inclusion; S3.3: Introducing self-balanced stress field in the element in the finite element: σ=pβ#(5) σ is a strong vector with 6 stress components, β is a column vector with m unknown stress coefficients, P is a 3×m matrix; the boundary displacement u is the interpolation of the generalized displacement d of the node: u=Ld#(6) Substituting (5) and (6) into (4), we obtain: After discretization, the weak form of the system complementary energy can be obtained, and according to the modified complementary energy stationary value condition: in, H m =∫ Ω P mT S m P m dΩH i =∫ Ω P iT S m P i dΩH c =∫ Ω P cT S m P c dΩ#(10) G mm =∫ Γ P mT n mT LdΓG ii =∫ Γ P iT n iT LdΓG cc =∫ Γ P cT n cT LdΓ#(11) G im =∫ Γ P iT n m TLdΓG ci =∫ Γ P cT n iT LdΓ#(12) The above formula can be written as the displacement expression of stress parameter: β=H -1 Gd#(13) The total energy of the system Π e The first-order variation of the node displacement d can give the corresponding weak expression of the traction boundary condition: Substitute (14) into (13), The system of equations for solving the generalized displacement is obtained, where the element stiffness matrix is:

3. A three-dimensional Voronoi unit simulation method for elastic problems of composite materials reinforced with interfacial phases according to claim 1, characterized in that: In S4, the construction process of the stress function is: S4.1: The matrix stress function is used The inclusion stress function is expressed as The interface stress function is expressed as The interaction stress function of the matrix is ​​expressed as The interface interaction force function and stress function are expressed as express; S4.2: The interfacial interaction stress function obtained by Ghosh is further simplified to obtain the following expression: In the formula, β m ,β i and β c are column vectors of stress function coefficients of matrix, interface and inclusion respectively; α is the ellipsoid equation, which is a known quantity derived from the coordinates of the model nodes; α1 is the boundary Γ c The parameters of the ellipsoid, and α2 is the boundary Γ i Parameters of the ellipsoid; S4.3: Near the interface Away from the interface In the interface After constructing the stress polynomials of the matrix and inclusions, the stress field is derived from the Maxwell stress function. Substituting the matrix stress polynomial into the matrix, interface and inclusion stress fields is as follows: in are scaled (x,y,z) coordinates, and when 1 / α i =1, β pqri The number of terms is equal to the number of terms in the polynomial stress function; S4.4: Select the most widely used Maxwell stress function. When assuming that the stress is a complete polynomial, the stiffness matrix must be reversible. According to the relationship between the number of terms in the stress function and the number of degrees of freedom of the rigid body, select d>bc, where d represents the total number of terms in the stress function, b represents the total number of degrees of freedom of the rigid body, and c represents the degrees of freedom of the rigid body. The basic number of terms in the Maxwell stress function used to calculate the three-dimensional arbitrary polyhedral unit is at least 31. Introduce the Maxwell stress function that satisfies the self-balance in the unit: Φ(x,y,z)=a0+a1x+a2y+a3z+β1x 2 +β2y 2 +β3z 2 +β4xy+β5yz+β6zx+β7x 3 +β8y 3 +β9z 3 +b 10 x 2 y+β 11 xyz 2 +b 12 y 2 z+β 13 yz 2 +b 14 z 2 x+b 15 zx 2 +b 16 xyz+b 17 x 4 +b 18 y 4 #(2-26) +β 19 z 4 +β 20 x 3 y+β 21 x 2 y 2 +β 22 xy 3 +β 26 y 3 z+β 27 y 2 z 2 +β 28 yz 3 +β 23 z 3 x 3 +β 24 z 2 +β 25 zx 3 +β 29 x 2 yz+β 30 xy 2 z+β 31 xyz 2 +… According to this stress function, the six Maxwell stress components can be obtained: get [s x s y s z t xy t yz t zx ]=Pβ#(22) P is a coefficient matrix with 6 rows and m columns, which is only related to the coordinates of the points. The β matrix is ​​an m-dimensional column vector, where m represents the number of terms selected by the polynomial.

4. A three-dimensional Voronoi unit simulation method for elastic problems of composite materials reinforced with interfacial phases according to claim 1, characterized in that: The division process of the Delaunay tetrahedron in S2 is: S2.1: Input unit node information (node ​​coordinates, node number); input surface information (surface node information, surface quantity); S2.2: Determine which phase the nodes and faces belong to based on the node information, and construct Delaunay triangles in space based on the node numbers; S2.3: Connect Delaunay triangles to construct tetrahedrons; S2.4: Determine whether the minimum circumscribed sphere of the tetrahedron generated by the four points satisfies the empty circle property. If so, the Delaunay tetrahedron partitioning is completed. If not, repeat S2.3 and subsequent steps until it is satisfied.

5. A three-dimensional Voronoi unit simulation method for elastic problems of composite materials reinforced with interfacial phases according to claim 1, characterized in that: Furthermore, in S5, the stiffness matrix of the unit obtained by solving is: You can get H * =H(λ * ) 3 andG * =G(λ * ) 2 , then there will be an additional unit scaling factor when calculating the stiffness matrix, namely:

Images