A three-dimensional voronoi cell simulation method for elastic problems of particle reinforced composites with interfacial phase

Abstract

The application designs a three-dimensional Voronoi unit simulation method for elastic problems of particle reinforced composite materials containing interface phase, deduces a three-dimensional modified complementary energy functional containing inclusions, matrix and interface, so that the method can be used for calculating three-dimensional three-phase composite materials; in addition, stress functions are constructed in three integral domains respectively, Maxwell complete polynomials are adopted, including stress function polynomials and interaction force function polynomials, and the three-dimensional interface units with different node numbers have good convergence and continuity. Compared with the commercial software, the method is more efficient when calculating the three-dimensional interface Voronoi model, and when calculating the same model, the results of hundreds of thousands of units in the commercial software can be obtained by using fewer units. Compared with the traditional displacement finite element method, the method is more targeted when solving the stress problems of particle composite reinforced materials.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of material simulation analysis, and particularly relates to a three-dimensional Voronoi cell simulation method for elastic problems of particle reinforced composite materials containing interface phases. BACKGROUND

[0002] Particle reinforced composites have attracted much attention due to their wide applications in automotive engineering, aerospace, military, electronics and nuclear engineering, etc. With the progress of technology, the particle size of materials is continuously reduced, and the influence of the interface layer, which is a thin layer between the matrix and the inclusion, becomes more and more significant. The formation mechanism of the interface is diverse, including chemical reaction of the reinforced particles and the matrix, surface coating of the particles, etc., which have a profound impact on the mechanical properties of the composite materials. Therefore, in-depth study of the interface characteristics is crucial for accurately predicting the overall performance of the composite materials.

[0003] In the study of particle reinforced composites, analytical method and finite element method (FEM) are two commonly used methods. K. Ding et al. proposed a homogenization theory, which established an elastic-plastic framework for composite materials with ductile interface characteristics, revealing the significant influence of interface ductility on the performance of the composite materials. Carrere et al. analyzed the influence of the interface phase on the crack deflection of the matrix in ceramic matrix composites using FEM, 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 FEM models, finding that the elastic modulus of the interface phase has a significant impact on the damping;

[0004] However, traditional FEM has limitations in handling constraint problems, singularity problems and mesh distortion, resulting in high computational cost and limited accuracy. To solve these problems, researchers have developed various improved methods. For example, TCGS proposed by L.T. Dong and Atluri improves the accuracy and efficiency of the calculation by satisfying the displacement field that controls the differential equation. Fan, J.P. 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 M K explored the influence of interface layer thickness and stiffness on the mechanical properties of composite materials, and found that the interface effect decreases with the increase of the aspect ratio of the filler. In addition, Voronoi cell finite element method (VCFEM) has attracted attention in the study of particle reinforced composites due to its simple cell division and fast calculation speed. This method can accurately reflect the random shape, size and spatial distribution of inclusions, and is suitable for analyzing models with a large number of inclusions. R. Guo and R. Zhang based on VCFEM proposed a modified complementary energy functional to study the micro and macro mechanical behavior of dissimilar materials. W.Y. Hao et al. used VCFEM to analyze the plasticity, creep and thermal strain of homogeneous materials, and verified the effectiveness of VCFEM.

[0005] Although these methods have made significant progress in the study of particle reinforced composites, there are still challenges. In particular, for three-dimensional Voronoi cells with complex interfaces, there is currently a lack of effective and accurate numerical methods. Therefore, the computational efficiency of large-scale complex microstructures is low, and an effective solution is provided for larger-scale microstructure problems. SUMMARY

[0006] To solve the above technical problems, the present application designs a three-dimensional Voronoi cell simulation method for elastic problems of particle reinforced composite materials containing interface phases, which regards stress and displacement as independent fields. Each Voronoi cell is divided into multiple Delaunay tetrahedrons, and Hammer integration 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 all internal points of the cell can be obtained. Compared with the finite element software ABAQUS, VCFEM can achieve the calculation accuracy of ABAQUS with hundreds of thousands of elements with fewer elements, and has stronger pertinence for three-dimensional composite reinforced material problems containing interface particles; effectively improves the calculation efficiency of complex microstructures, and provides an effective solution for larger-scale microstructure problems.

[0007] To achieve the above technical effects, on the one hand, the present application discloses a three-dimensional Voronoi cell simulation method for elastic problems of particle reinforced composite materials containing interface phases, comprising the following steps:

[0008] S1: Establishing a model according to the actual problem, and discretizing the model into Voronoi cells connected with each other;

[0009] S2: Dividing each cell into a plurality of Delaunay tetrahedrons;

[0010] S3: Constructing a cell modified complementary energy functional according to the minimum complementary energy principle;

[0011] S4: Adopting a high-order complete polynomial as a cell stress field, and constructing a suitable trial function;

[0012] S5: Obtaining a cell stiffness matrix by solving G matrix and H matrix;

[0013] S6: Obtaining node stress and internal node displacement according to the relationship between external node displacement and the stiffness matrix;

[0014] S7: Calculating β according to the internal node displacement;

[0015] S8: Obtaining the final result of stress by calculating the coefficient matrix ρ through the partial derivative of the stress function, (σ = ρβ).

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

[0017] S3.1: Constructing a three-dimensional phase interface Voronoi cell according to the existing assumed stress hybrid element method balance model I and balance model II; the cell is composed of three phases, including a matrix, an interphase and an inclusion, and on a given traction boundary:

[0018]

[0019] On the boundary of the inclusion and the interface inside the cell:

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

[0021] On the internal interface of the cell and the boundary of the matrix:

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

[0023] S3.2: Obtaining a modified complementary energy functional by introducing a Lagrange multiplier, which can be specifically expressed as:

[0024]

[0025] Γ is the boundary of Ω m Ω is the integral region of the matrix i Ω is the integral region of the interface phase c Ω is the integral region of the inclusion

[0026] S3.3: Introducing self-balanced stress field within element in 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 3xm matrix; the boundary displacement u is the interpolation of the nodal generalized displacement d:

[0029] u = Ld#(6)

[0030] Substituting (5) and (6) into (4) gives:

[0031]

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

[0033]

[0034] where,

[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 LdΓ G ii = ∫ Γ p iT n iT LdΓ G cc = ∫ Γ P cT n cT LdΓ#(11

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

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

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

[0040] The first variation of the total energy of the system Π e with respect to the nodal displacement d, the weak expression of the boundary condition of the surface force can be obtained:

[0041]

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

[0043]

[0044] The equation group of the generalized displacement is obtained, in which the element stiffness matrix is:

[0045]

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

[0047] S4.1: the stress function of the matrix part is represented by , the stress function of the inclusion is represented by , the stress function of the interface is represented by , the stress function of the interaction force of the matrix part is represented by , and the stress function of the interaction force of the interface is represented by ;

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

[0049]

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

[0051] S4.3: at the vicinity of the interface at the far field of the interface at the interface After the stress polynomial of the matrix and the inclusion is constructed, the stress field is derived from the Maxwell stress function, and the stress field of the matrix, the interface and the inclusion is obtained by substituting the stress polynomial of the matrix into the Maxwell stress function as follows:

[0052]

[0053]

[0054]

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

[0056] S4.4: the most widely used Maxwell stress function is selected, and when the stress complete polynomial is assumed, the condition of the reversible stiffness matrix is met, d > b-c is selected according to the relationship between the number of terms of the stress function and the number of degrees of freedom of the rigid body, wherein d represents the total number of terms of 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 for calculating a three-dimensional arbitrary polyhedral element is at least 31 terms; the Maxwell stress function satisfying the self-balancing of the element is introduced:

[0057]

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

[0059]

[0060]

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

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

[0063] Further, the partition process of the Delaunay tetrahedron in S2 is as follows:

[0064] ​S2.1: input unit node information (node coordinates, node number); input face information (node information on the face, the number of faces);

[0065] S2.2: according to the node information, determine which phase the node and the face belong to, and construct a Delaunay triangle in the space according to the node number;

[0066] S2.3: connect the Delaunay triangle to build a tetrahedron;

[0067] S2.4: determine whether the minimum circumscribed sphere of the tetrahedron generated by the four points satisfies the empty circle property, if it satisfies, the Delaunay tetrahedron partition is completed, if it does not satisfy, repeat S2.3 and the subsequent steps until it is satisfied.

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

[0069]

[0070] H can be obtained * = H (lambda * ) 3 and G * = G (lambda * ) 2 When calculating the stiffness matrix, an additional unit scaling coefficient will be added, that is:

[0071]

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

[0073] The present application designs a kind of three-dimensional Voronoi unit simulation method for the elastic problem of interface phase particle reinforced composite material based on hybrid finite element method and variation principle, deduces three-dimensional modified complementary energy functional including inclusion, matrix and interface, so that the method can be used to calculate three-dimensional three-phase composite material;

[0074] In addition, the stress function is constructed in three integral domains respectively, the Maxwell complete polynomial is adopted, the stress function polynomial and the interaction force function polynomial are included, and the three-dimensional interface element with different node numbers has good convergence and continuity. Compared with the commercial software, the method is more efficient when calculating the three-dimensional interface Voronoi model, and when calculating the same model, a few elements can obtain the result of hundreds of thousands of elements in the commercial software. Compared with the traditional displacement finite element method, the method is more targeted when solving the stress problem of the particle composite reinforced material; and the calculation efficiency of the large-scale complex microstructure is higher than that of the traditional finite element analysis, and a few elements can achieve the calculation precision of hundreds of thousands of elements, and the three-dimensional composite reinforced material problem containing the interface particle has stronger pertinence. BRIEF DESCRIPTION OF DRAWINGS

[0075] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings needed to be used for the embodiment description will be briefly introduced as follows.

[0076] Figure 1 is a Voronoi model with interface of the present application;

[0077] Figure 2 is an empty circle characteristic diagram of the present application, wherein a satisfies the empty circle characteristic diagram, and b satisfies the non-empty circle characteristic diagram;

[0078] Figure 3 is a tetrahedron partition example diagram of the present application;

[0079] Figure 4 is a boundary condition diagram of the present application;

[0080] Figure 5 is a polygon approximate ellipsoid diagram of the present application, wherein (a) is a polygon approximate ellipsoid; (b) is a single layer 6 node; (c) is a single layer 12 node; and (d) is divided into 24 layers;

[0081] Figure 6 is a mesh partition diagram of the present application, wherein (a) is an ABAQUS spherical inclusion (361105 elements); (b) is a VCFEM 6 node polygon inclusion (1 element); (c) is a VCFEM 12 node polygon inclusion (1 element); and (d) is a VCFEM 6 node without interface (1 element);

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

[0083] Figure 8 Figure is the stress curve chart of three models and pull-out position of the application, in which (a) is validity verification, (b) is interface influence verification;

[0084] Figure 9 Figure is the boundary condition chart of the application;

[0085] Figure 10 Figure is the multi-unit model establishment and boundary condition of the application, in which (a) is ABAQUS spherical inclusion; (b) is VCFEM 6-node polygonal inclusion; (c) is VCFEM 12-node polygonal inclusion;

[0086] Figure 11 Figure is the meshing chart of the application, in which (a) is ABAQUS spherical inclusion (1708174 units); (b) is VCFEM 6-node polygonal inclusion (27 units); (c) is VCFEM 12-node polygonal inclusion (27 units);

[0087] Figure 12 Figure is the multi-unit stress cloud chart of the application, in which (a) is ABAQUS spherical inclusion result; (b) is VCFEM 6-node polygonal inclusion result; (c) is VCFEM 12-node polygonal inclusion;

[0088] Figure 13 Figure is the stress curve of three models and pull-out position of the multi-unit of the application;

[0089] Figure 14 Figure is the model establishment and boundary condition chart of the unit of the application when the unit is expanded to 4 through array, in which (a) is ABAQUS, (b) is VCFEM;

[0090] Figure 15 Figure is the meshing of the unit of the application when the unit is expanded to 4 through array, in which (a) is ABAQUS spherical inclusion (512423); (b) is VCFEM 6-node polygonal inclusion (4 units);

[0091] Figure 16 Figure is the stress cloud chart of two calculation results of the application;

[0092] Figure 17 Figure is the stress condition chart on the path of the application;

[0093] Figure 18 is the random distribution model and boundary condition diagram of the present application;

[0094] Figure 19 is the stress nephogram of the random distribution model of the present application. DETAILED DESCRIPTION

[0095] The application discloses a three-dimensional Voronoi cell simulation method for elastic problems of particle reinforced composite materials containing interface phases,

[0096] According to the assumed stress hybrid element method created by Mr. Bian Xuecheng, 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 taking stress function as a field variable according to the minimum complementary energy principle. The expression of the complementary energy functional is as follows:

[0097]

[0098] In the formula, S is an elastic compliance tensor, and Γ is a given displacement boundary. The stress of the element is expressed by using a stress function, and the stress function is taken as a basic unknown quantity at a node; the stress function is ensured to be continuous on the boundary between elements, and therefore the equilibrium condition can be met everywhere. The equilibrium model I has the following defects: the distribution of stress and strain is obtained by derivation of the stress function, and the displacement solution is further obtained by integration of the geometric equation, which depends on the selection of the integration path, therefore, the displacement solution obtained in this way is not unique, and in addition, the stress function is not convenient for intuitively explaining the physical meaning of the boundary condition.

[0099] The equilibrium model II is proposed by Fraeijs in 1964 and is used for analyzing many plate bending problems. The model is obtained by introducing the equilibrium condition of the surface force continuity between elements into the functional based on the traditional complementary energy functional by using the Lagrange multiplier method.

[0100] For the original problem, there is a given surface force boundary condition:

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

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

[0103] n + is the outer normal of the boundary, and n - is the inner normal of the boundary. The two constraint conditions can be introduced by using the Lagrange multiplier to obtain a modified functional. And the boundary surface force can be written as T = n·σ. Given displacement also satisfies the element displacement interpolation, and therefore The function is usually written as:

[0104]

[0105] It can be seen that the stress hybrid model introduces two unknown variables, namely the stress field variable σ and the boundary displacement variable u, in the unit interior and unit boundary respectively. The unit node displacement interpolation can satisfy all edge coordination conditions of the entire unit. The stress hybrid element method has the advantage of calculating the boundary displacement without being 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 provides a basis for three-dimensional Voronoi unit interface problems. Based on this principle, a three-dimensional phase interface Voronoi unit is constructed. The Voronoi model with an interface used in the present application is as shown in Figure 1 .

[0106] The unit is composed of three phases, including a matrix, an interphase and an inclusion. On the given traction boundary:

[0107]

[0108] On the boundary of the inclusion and the interface in the unit interior:

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

[0110] On the internal interface of the unit and the boundary of the matrix: 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 integral region of the matrix, Ω i is the integral region of the interface phase, and Ω c is the integral region of the inclusion.

[0114] The self-balanced stress field in the unit is introduced in the finite element:

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

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

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

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

[0119] Substituting (2-9) and (2-10) into (2-8) gives:

[0120]

[0121]

[0122] The weak form of the system residual energy can be obtained after discretization. According to the modified residual energy stationary condition:

[0123]

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

[0125] where,

[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 LdΓ G ii =∫ Γ P iT n iT LdΓ G cc =∫ Γ P cT n cT LdΓ#(2-15

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

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

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

[0131] The total energy of the system Π e The first variation of the node displacement d gives the weak expression of the boundary condition of surface force:

[0132]

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

[0134]

[0135] The equation group of generalized displacement is obtained, where the element stiffness matrix is

[0136]

[0137] Equation (2-8) shows that the element displacement variable appears on all boundaries of the element, so the displacement mode of the element can be obtained by interpolation to get the displacement of any point. In order to continue the derivation of the equation, the expression of the stress function, the normal vector of the element surface and the cosine formula are needed. The calculation and solving method of related variables are introduced next.

[0138] In three-dimensional problems, the reasonable selection of stress function and the accurate representation of the stress field of the element are problems that need to be considered. The matrix of the base stress function is represented by , the matrix of the inclusion stress function is represented by , the matrix of the interface stress function is represented by , the matrix of the base interaction stress function is represented by , and the matrix of the interface interaction stress function is represented by The construction of the polynomial stress function and the interaction function should meet the following requirements: and The columns in the matrix H should meet the condition of linear independence, and only when this condition is met, the matrix H is invertible. The selection of the stress function polynomial needs to fully consider the shape of the inclusion, and the interaction stress term should approach zero far away from the contact boundary of different phases, and the surface force continuity condition should be met at the contact boundary of different phases.

[0139]

[0140] β m ,β i and β care column vectors of the stress function coefficients of the 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 ellipsoid equation, which is a known quantity derived from the coordinates of the model nodes. α1 is the equation describing the boundary Γ c The parameters of the ellipsoid, and α2 is the boundary Γ i Parameters of the ellipsoid. Figure 1 As shown. Near the interface Away from the interface In the interface After constructing the matrix and inclusion stress polynomials, 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 phase, respectively. Due to the different parameters α of the two ellipsoids, the choice of stress function polynomial is different, resulting in the unknown column vector β m and β i The final stress field is different depending on the number of terms in . is the scaled (x,y,z) coordinate, and when 1 / α i =1, β pqri The number of terms in the polynomial stress function is equal to the number of terms in the stress function. This paper uses the most widely used Maxwell stress function. When assuming a complete stress polynomial, the stiffness matrix must be reversible. Based on 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. Here, 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 number of basic terms in the Maxwell stress function used to calculate arbitrary three-dimensional polyhedral elements must be 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 one term in the 0th order, three terms in the 1st order, and six terms in the 2nd order. Therefore, the number of terms in the mth order (m>0) complete polynomial is The total number of terms of the stress function Φ is Removing the 0th and 1st order terms, there are 31 terms in total shown in equation (2-21). According to this stress function, the 6 Maxwell stress components can be obtained:

[0150]

[0151] Then we can get

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

[0153] P is a 6-row m-column coefficient matrix, which is only related to the coordinates of the points, and β is an m-dimensional column vector, called the stress coefficient. m represents the number of terms selected by the polynomial.

[0154] As shown in a in Figure 2 , in the circumcircle of each triangle, only the three points of the triangle are present, and no other points are present. If other points are present in the circumcircle, the empty circle property is not satisfied, as shown in b in Figure 2 . In a three-dimensional problem, only the four points of a tetrahedron are present in the circumsphere of each tetrahedron, and no other points are present. For example, Figure 3 , a hexahedron in is divided into tetrahedrons. After constructing the Delaunay triangle on each face, the fourth point that satisfies the empty circle feature is found. The tetrahedron partitioning is successful, and each tetrahedron satisfies the empty circle property.

[0155] Figure 2 Next, the displacement interpolation on the surface of the tetrahedron is calculated, considering a spatial triangle as shown in , which is numbered i, m, j. The coordinates of the respective nodes in the rectangular coordinate system are (x i ,y i ,z i ), (x m ,y m ,z m ), and (x j ,y j ,z j ), respectively. The boundary displacement on the surface of the element can be obtained by interpolating the node displacement, as shown in equation 2-29.

[0156]

[0157] In the equation, Δ is the area of the triangle, which can be calculated according to the coordinates of the three points, and Δ iThe area of the small triangle formed by the two corner points on the edge and the interior point p. Rewrite equation (2-29) in matrix form: u = Ld# (2-30)

[0158] where

[0159] d = [u v w] 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 will 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 be determined by the node coordinates.

[0162]

[0163] We 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, we also need to solve the normal vector cosine of each face of all tetrahedrons respectively, to form the normal vector matrix of the triangle. The formula of the normal vector of 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 is represented as:

[0172]

[0173] In order to solve the problem of too large or too small numerical value leading to the irreversibility of the matrix in the calculation process, and to avoid too large error of the calculation result, it is necessary to scale the unit to obtain the local coordinates for numerical calculation. Therefore, the centroid coordinates of the unit should be obtained, and the unit is scaled with the centroid 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 it is impossible to directly determine the centroid thereof. 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 indirectly obtained through the centers of the tetrahedrons.

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

[0175] p1(x1,y1,z1),p2(x2,y2,z2),p3(x3,y3,z3),p4(x4,y4,z4)#(2-39)

[0176] It is assumed that an arbitrary polyhedron is divided into n tetrahedrons, and the center coordinates of the tetrahedrons are respectively

[0177]

[0178] The volumes are respectively

[0179]

[0180] Therefore, the center coordinates of the three-dimensional arbitrary polyhedron can be obtained as

[0181]

[0182] After obtaining the center coordinates of the polyhedral unit, scaling can be performed. With the center coordinates of the polyhedral unit as the origin of the local coordinates, also called the scaling center, the local coordinates of any point (x, y, z) on the polyhedron after scaling can be obtained

[0183]

[0184] Since the stress function is closely related to the unit coordinates, the stress function can be written as a complete polynomial function related to ξ, η, and ζ according to the scaled local coordinates.

[0185]

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

[0187]

[0188] is the boundary of the integral region Ω e . The main method is Hammer integration. Hammer integration formula is:

[0189]

[0190] In the formula, ω k is the weight coefficient of Hammer integration point. The integration of H matrix is carried out in the tetrahedron domain, so the result of the whole unit can be obtained by accumulating the integration of each small tetrahedron, that is:

[0191]

[0192] The scaled local coordinates are:

[0193]

[0194] In the formula, n tet represents the number of tetrahedrons obtained by division, n int represents the number of tetrahedron Hammer integration points adopted, ω j represents the weight coefficient of the jth integration point. J i represents the Jacobian of the ith 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 also needs to be accumulated in the triangular surface integration, that is:

[0197]

[0198] The scaled local coordinates are:

[0199]

[0200] H * = H(λ * ) 3 and G * = G(λ * ) 2 , so there is an additional unit scaling coefficient when calculating the stiffness matrix, that is:

[0201]

[0202] Example 2

[0203] To further validate the simulation method constructed in Example 1, the VCFEM was compared with a commercial finite element analysis software ANSYS® single element; the size of the model and the radius of the internal structure (ellipsoidal inclusions and interface) were determined. Then, according to the layered situation of the inclusion and interface and the required number of points for each layer, the information of all points in the VCFEM could be obtained, and the continuity condition of the surface force at the interface of different materials was satisfied. Both models used displacement constraints, and the same boundary conditions were applied to the nodes on the same face in the same coordinate system. The stress distribution of the two models under tensile state was calculated using the static structural analysis method.

[0204] ABAQUS model establishment: a cube entity with a side length of 20 mm was generated, and then two spherical entities with radii of 2 mm and 2.35 mm were created. The entities were moved so that their centroids coincided. By Boolean operation, an entity that retains its outline can be obtained, as shown in Figure 1 The material parameters specified in Table 1 were applied to the corresponding regions. Multiple materials coexist in a single entity, consistent with the characteristics of composite materials.

[0205] Table 1 Material parameters

[0206]

[0207] Therefore, it is not necessary to establish additional contact or connection conditions between different material phases. Then mesh the entire entity and apply displacement boundary conditions, as shown in Figure 4 The stress results can be obtained by static structural analysis.

[0208] VCFEM modeling: a three-dimensional interface-containing Voronoi cell is as shown in Figure 4 The cell is composed of three parts, namely the matrix phase, the inclusion phase and the interface phase, and the cell edge length is 20 mm. The inclusion and interface are approximated by the inscribed polygon of an ellipsoid, as shown in Figure 5 The radius of the inclusion ellipsoid is 2 mm, and the interface thickness is 0.35 mm. The corresponding material properties are shown in Table 1. To study the effect of the interface layer on the composite material, we set up a model without an interface layer and a model with an interface layer thickness of 0.65 mm, and other model parameters are consistent with the VCFEM six-node model. The boundary conditions are as follows: the displacement of the X, Y, Z faces on the outer side of the matrix in the corresponding X, Y, Z directions is 0 mm, and the displacement of the Y face in the Y direction is 0.02 mm. The specific model is as shown in Figure 4

[0209] Table 3-1 Material parameters

[0210]

[0211] ​One of the advantages of the method presented in this paper is that it can get the result of traditional displacement finite element method with several ten thousand elements by using very few elements. The model meshing is shown in Figure 5

[0212] The calculation model is shown in Figure 6 The VCFEM model has only one Voronoi cell with interface, and the ABAQUS model is composed of 361105 four-node tetrahedral elements. The VCFEM 6-node model and the VCFEM 12-node model have 8 external nodes, and the internal nodes (inclusion nodes and interface nodes) are 48 and 96 respectively. The number of stress function terms is shown in Table 2:

[0213] Table 2 Number of stress function terms

[0214]

[0215] The result comparison and verification is shown in the result cloud map Figure 7 This paper compares the calculation results of ABAQUS and VCFEM, verifies the convergence, accuracy and polygon modeling feasibility of the calculation results of the method. The stress data on the vertical middle path is extracted from the cloud map. The stress curves of the three models and the pull-out position are shown in Figure 8

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

[0217]

[0218] u1, u2, and u3 are the displacements of the nodes in the X, Y and Z directions, p is the total stress of the node, v is the Poisson's ratio, E is the elastic modulus of the material, x1, x2, and x3 are the coordinates of the node before the model deformation.

[0219] The three material parameters of the model VCFEM 6-node polygon inclusion (1 cell) are modified to set the Young's modulus to 72000 MPa and the Poisson's ratio to 0.33. The boundary conditions are shown in Figures 3-6 The nodes 2, 4, 6 and 8 are applied with a displacement boundary condition of 0.01 mm in the Y direction, and the nodes 1, 3, 5 and 7 are constrained in the Y direction. The stress in the calculation result is brought into formula 3-1 to verify whether the calculation of stress and node displacement in the method is correct. 100 sampling points are selected in the plane of nodes 2, 4, 6 and 8, and the verification results are shown in Table 3.

[0220] Table 3 Node displacement and error

[0221]

[0222] ​​

[0223] In 100 sampling points, the error is between 0% and 1.25%, and the average node displacement is 0.010068503 mm, which is very close to the 0.01 mm 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 of the ABAQUS spherical model and the VCFEM polygonal model under uniaxial tension in the Y direction is shown. The results obtained from the stress program and the stress path curve show a high degree of consistency, indicating that the method proposed in this paper can obtain a calculation precision comparable to that of the traditional displacement finite element method using hundreds of thousands of elements. Although this method uses fewer points to approximate the ellipsoid, it still treats the polygon as an ellipsoid when dealing with the stress function. The calculation results further show that the polygon can achieve the expected results of the ellipsoid model. By increasing the number of points inside the polygon, a result 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, and in particular, its absence leads to greater boundary stress fluctuations, thereby significantly affecting material damage. In actual composites, the interface phase is usually a coating, and its thickness is 1 / 6 of the inclusion radius, and occasionally a man-made interface phase is added, which has 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, and the maximum stress difference is 15%, and the minimum stress difference is 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 the method is compared and verified with multiple elements of a commercial finite element analysis software, and the specific verification process is as follows:

[0227] After verifying the stress results of a single Voronoi cell, this paper will verify multiple cells. In an array manner, the cells are expanded to 27, and the cell size remains unchanged, and the material properties are the same as in Table 1. The boundary conditions are: fixing the X, Y, and Z sides of the matrix outside in the corresponding X, Y, and Z directions to 0 mm, and giving the Y direction displacement on the Y surface to 0.02 mm, and the specific model is as shown in Figure 10 The grid partitioning of the calculation model is as shown in Figure 11The VCFEM model has only 27 Voronoi cells with interfaces, while the ABAQUS model is composed of 1708174 4-node tetrahedral elements. The VCFEM6 and VCFEM12 models have 64 external nodes, 1296 and 2592 internal nodes (inclusion nodes and interface nodes), respectively. The number of stress function terms is shown in Table 4, where pol: stress function, rec: interaction force function.

[0228] Table 4 Stress function term table

[0229]

[0230] The comparative verification result cloud map is shown in Figure 12 , and the stress data on the vertical center line path is extracted from the program. The stress data on the vertical center path is extracted from the cloud map. The stress curves of the three models and the pull-out position are shown in Figure 13 , as shown in Figure 12 and Figure 13 , the stress distribution 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 and the stress path curve show a high consistency, indicating that the method proposed in this paper can obtain the calculation accuracy comparable to the traditional displacement finite element method (more than 1 million elements) using only a few tens of elements.

[0231] Example 4

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

[0233] In this embodiment, the element is expanded to 4 through an array, and the element 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, Z edges on the outer side of the matrix in the corresponding X, Y, Z directions is 0 mm, and the given Y direction displacement in the Y plane is 0.02 mm. The specific model and meshing are shown in Figure 14 , Figure 15 .

[0234] The comparison of the calculation results is shown in Table 5, which lists the number of stress terms. Figure 16 and Figure 17 show the calculation results;

[0235] Table 5 Stress term table

[0236]

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

[0238] A path on the stress cloud is selected for analysis, as shown in Figure 16 , and the stress on the path is shown in Figure 17 . From the stress cloud and stress path curve, it can be seen that the results of the five examples are relatively close, and the use of different numbers of complete polynomials can also obtain very similar solutions. Thus, when the stress term number reaches a certain number, the calculation result tends to be stable, and the stress term number within a certain range can obtain correct results. The relationship between the stress function term number and the number of rigid body degrees of freedom is selected as d > b-c.

[0239] Example 5

[0240] In actual particle composite reinforced materials, the distribution of inclusions is random in most cases. In order to simulate the actual situation, a random distribution model is used in this embodiment. A 6-node polygon is used for modeling, in which 20 inclusions are randomly distributed in a regular hexahedron with a size of 20mm*20mm*20mm, and a uniform thickness interface phase is distributed between the matrix and the inclusion. The number of Voronoi cells containing the interface is 20, the number of external nodes is 105, and the number of internal nodes is 960. The boundary conditions are: the displacement of the fixed X, Y and Z edges in the corresponding X, Y and Z directions is 0mm on the outside of the matrix, and the displacement of the Y edge in the given Y direction is 0.02mm. The boundary conditions and meshing are shown in Figure 18 , the stress term number is shown in Table 6, and the element node information is shown in Tables 7 and 8:

[0241] Table 6 Stress term number for random distribution verification

[0242]

[0243] Table 7 Element node information

[0244]

[0245] Table 8 Element node information

[0246]

[0247]

[0248] The stress cloud of the calculation result is shown in Figure 19As shown, the stress contour shows that even in the case of random distribution, the number of nodes, faces and boundary conditions of each element is different, however, the preset complete polynomial can still effectively capture the complex details of the stress distribution within the element. In this model, the element has a maximum of 48 internal nodes and 24 external nodes. The total number of degrees of freedom of these elements is calculated as (48+24) x 3 = 216. Considering the elimination of 12 rigid body displacements, the final total number of accumulated degrees of freedom is 228. The number of terms of the stress function is set to 600. After several evaluations on the selection of the terms of the stress function, 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 each node of the element plus the degrees of freedom corresponding to the twelve rigid bodies. When a certain number is reached, the same stress function exhibits better convergence on different elements, and the calculation results tend to be stable; it can effectively show that the present method provides a reference when dealing with models that cannot be calculated using commercial software or verified through experimental means, and provides a new basis for studying composite materials with randomly distributed structures.

[0249] The preferred embodiments of the application disclosed above are only used to help illustrate the application, the preferred embodiments do not describe all the details, nor limit the application to the specific embodiments described.

Claims

1. A three-dimensional Voronoi element simulation method for the elasticity problem of composite materials reinforced with interfacial 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 element stress fields and construct appropriate trial functions; S5: Obtain the element stiffness matrix by solving the G matrix and the H matrix; S6: Based on the relationship between the external node displacement and the stiffness matrix, the node stress and internal node displacement are obtained; S7: Calculated based on internal node displacements ; S8: Obtain the coefficient matrix by taking partial derivatives of the stress function , get the final result of stress, , is the node displacement, is the coefficient matrix.

2. A three-dimensional Voronoi element simulation method for elasticity problems of composite materials containing interfacial particles according to claim 1, characterized in that: The process of correcting the complementary energy functional in S3 includes the following steps: S3.1: Construct a three-dimensional phase interface Voronoi cell based on the existing hypothetical stress hybrid element method equilibrium model I and equilibrium model II. The cell consists of three phases, including matrix, interphase, and inclusions. On a given traction boundary: ； At the boundaries of inclusions and interfaces within cells: ； On the internal interfaces of the element and on the boundaries of the matrix: ； S3.2: The modified complementary energy functional is obtained by introducing the Lagrange multiplier: Specifically, it can be expressed as: ； ； ； is the surface force boundary; is the integration region of the matrix, is the integrated area of ​​the interface phase, is the inclusion integral region; S3.3: Introducing self-equilibrium stress fields within the element in the finite element: ； is a strong vector with 6 stress components, is a column vector with m unknown stress coefficients, P is a 3×m matrix, and the boundary displacement u is an interpolation of the generalized displacement d of the node: ； Substituting (5) and (6) into (4) yields: ； ； ； After discretization, the weak form of the system complementary energy can be obtained, and according to the modified complementary energy stationary value condition: ； ； in, ； ； ； The above formula can be written as the displacement expression of stress parameter: ； The total energy of the system For the first-order variation of the node displacement d, the corresponding weak expression of the traction boundary condition can be obtained: ； Substitute (14) into (13), ； The equations for solving the generalized displacement are obtained, where the element stiffness matrix is: 。 3. The three-dimensional Voronoi element simulation method for elasticity problems of composite materials reinforced with interfacial particles according to claim 1, characterized in that: In S4, the construction process of the stress function is: S4.1: Matrix stress function 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: ； Where, , and are column vectors of stress function coefficients of matrix, interface and inclusion respectively; is the equation of the ellipsoid, which is a known quantity derived from the coordinates of the model nodes; Describes the boundary The parameters of the ellipsoid, and Describes the boundary Parameters of the ellipsoid; S4.3: Near the interface , away from the interface , in the interface After constructing the matrix and inclusion stress polynomials, 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 hour, The number of terms in is equal to the number of terms in the polynomial stress function; S4.4: Use the most widely used Maxwell stress function. When assuming a complete stress polynomial, the stiffness matrix must be reversible. Based on 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. For calculating arbitrary three-dimensional polyhedral elements, the basic number of terms in the Maxwell stress function must be at least 31. Introduce a Maxwell stress function that satisfies self-equilibrium within the element: ； ； ； ； ； According to this stress function, the six Maxwell stress components can be obtained: ； get ； Is a 6-row m-column coefficient matrix, which is only related to the coordinates of the point. A matrix is ​​an m-dimensional column vector, where m represents the number of terms in the polynomial selection.

4. The three-dimensional Voronoi element simulation method for elasticity problems of composite materials containing interfacial particles according to claim 1, characterized in that: The partitioning process of the Delaunay tetrahedron in S2 is: S2.1: Input unit node information, including node coordinates and node numbers; input surface information, including surface node information and 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. The three-dimensional Voronoi element simulation method for elasticity problems of composite materials containing interfacial particles according to claim 1, characterized in that: Furthermore, in S5, the stiffness matrix of the element obtained by solving is: ； Can get and , then there will be an additional unit scaling factor when calculating the stiffness matrix, namely: 。

Images