Three-dimensional thermo-elastic-plastic voronoi cell modeling of composites with interfacial phase particles

Abstract

The present application relates to the technical field of interfacial phase particle reinforced composite material, and particularly relates to a three-dimensional thermo-elastic-plastic Voronoi unit simulation method of interfacial phase particle reinforced composite material, a three-dimensional Voronoi unit model generated by a polygon approximation method of an embedded sphere can accurately reproduce the real geometric configuration of the matrix, inclusions and interfacial phase inside the particle reinforced composite material. By introducing a modified complementary energy functional of Lagrange multipliers, the strictness of force balance and energy conservation in numerical calculation is ensured, thereby the numerical stability of simulation is enhanced. In the simulation of thermo-elastic-plastic behavior, the present application combines the temperature increment and the thermal expansion coefficient of the material, adopts the Mises yield condition, and adopts the equivalent plastic strain expansion to specifically describe, thereby the reaction of the material under different thermal mechanical conditions is comprehensively modeled. The present application significantly improves the calculation efficiency, and exhibits the advantage of VCFEM in processing complex material models.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of interfacial phase particle reinforced composites, and particularly relates to a three-dimensional thermo-elastic-plastic Voronoi element simulation method of interfacial phase particle reinforced composites. BACKGROUND

[0002] Particle reinforced composites are widely used. The formation of the interface layer is due to the chemical reaction between the reinforcing particles and the matrix and the coating on the surface of the particles. Existing research results show that the interface has an important influence on the overall mechanical properties of the composite material. K. Ding et al. found that even if the interface phase is very thin, it will have a very significant influence on the overall mechanical properties of the composite material; Gu, JH et al. research shows that the damping capacity of the elastic interface composite material increases with the increase of the elastic modulus. For the plastic interface phase, weak interface phase or hard interface phase is beneficial to improve the overall damping of the composite material; HassanzadehAghdamM K et al. research results show that the influence of the interface phase on the transverse thermal performance of the fiber composite material is more significant than the influence on the longitudinal thermal performance of the fiber composite material.

[0003] In the research of thermal problems, Buryachenko VA's research results show the relationship between thermal expansion, stored energy and average thermal elastic strain in the assembly. Wang H et al. proposed a hybrid finite element scheme with basic solutions as kernel functions, considering a representative volume element with inclusions, and studied the influence of inclusions on the global temperature distribution and the interaction between inclusions; Kaur et al. studied the photo-thermo-elastic interaction of a solid cylinder in a strong magnetic field, and illustrated the influence of the thermo-elastic theory on the Hall current displacement, temperature and thermal stress through expressions in the physical domain.

[0004] With the development of two-dimensional particle reinforced composites, Guo proposed a new two-dimensional modified complementary energy function considering thermal strain and plastic strain, and proved that the modified complementary energy function is correct and effective in studying the mechanical properties of composite materials. Zhang Rui established a two-dimensional Voronoi element finite element model based on linear elastic theory to study the mesoscopic and macroscopic mechanical behavior of heterogeneous materials. The results of VCFEM are compared with the analytical solution and numerical results of standard finite element analysis to confirm its effectiveness. Hao et al. numerically simulated the fatigue behavior of porous materials under different conditions by using a modified complementary energy function composed of plastic strain and thermal strain; Mao Chao et al. proposed a new extended VCFEM method to simulate the generation of particle damage, which embodies the accuracy of determining stress concentration and coarse grid discretization. Rao et al. studied the plastic strain, thermal strain and creep strain of particle composites by VCFEM.

[0005] In the study of three-dimensional particle-reinforced composites, Dong and Atluri proposed the concept of Treffutz computational particles (TCGs). Each particle is geometrically a polyhedron composed of three phases: inclusions, coating surroundings (or intervals), and the matrix. They used radial basis functions and the Treffutz formulation for multifunctional composites. The newly developed TCG involves only boundary integrals, thus ensuring the accuracy and efficiency of micromechanical calculations. Wang Lihui et al. proposed the HF-Fem and HTS-Fem methods based on three-dimensional polyhedral octrees to calculate steady-state thermal conduction and thermal stresses in particle-reinforced composites and verified their effectiveness. Rui Xu et al. established a multiscale analytical model for fiber-reinforced composites using the multiscale finite element method. The properties of SMA fibers are closely related to internal stress. When the matrix stiffness is high, the internal stress during the unloading phase may eliminate some of the residual strain in the SMA fibers, resulting in weakening of the SME.

[0006] Currently, no researchers have used the stress hybrid element method to analyze the nonlinearity of particle-reinforced composites with interfaces. Compared with the traditional displacement finite element method, which requires a large number of elements to obtain accurate results, this method is superior in three-dimensional nonlinear research. Summary of the Invention

[0007] The purpose of the present invention is to provide a three-dimensional thermo-elastic-plastic Voronoi unit simulation method for particle-reinforced composite materials containing an interface phase. The three-dimensional Voronoi unit model generated by the polygonal approximation method of embedded spheres can accurately reproduce the true geometric configuration of the matrix, inclusions and interface phase inside the particle-reinforced composite material.

[0008] In order to achieve the above technical objectives and the above technical effects, the present invention is implemented through the following technical solutions:

[0009] A three-dimensional thermo-elastic-plastic Voronoi element simulation method for a composite material reinforced with particles containing an interface phase includes the following steps:

[0010] S1: Establish a calculation model based on dimensions, boundary conditions, and material parameters;

[0011] S2: Discretize the model into several Voronoi cells containing one inclusion;

[0012] S3: Divide each Voronoi cell into several Delaunay tetrahedrons;

[0013] S4: Derivation of modified complementary energy functional based on the minimum complementary energy principle and thermoplastic constitutive equation;

[0014] S5: Select high-order complete polynomials as trial functions to construct stress field functions of three-phase materials considering ellipsoidal interfaces;

[0015] S6: Obtain the element stiffness matrix by varying and solving the G matrix and H matrix;

[0016] S7: Determine the iterative expressions for stress parameter β and displacement d;

[0017] S8: Calculate the node stress based on the external node displacement and the element stiffness matrix to obtain the internal node displacement;

[0018] S9: Calculate stress parameter dβ = -H -1 G th +H -1 Gdq;

[0019] S10: Take the partial derivative of the stress function to obtain the matrix P and calculate the full-field stress:

[0020] [Δσ x Δσ y Δσ z Δτ xy Δτ yz Δτ zx ]=PΔβ

[0021] Complete the simulation of three-dimensional thermo-elastic-plastic Voronoi elements.

[0022] Beneficial effects of the present invention:

[0023] The three-dimensional Voronoi unit model generated by the polygonal approximation method of embedded spheres in the present invention can accurately reproduce the real geometric configuration of the matrix, inclusions and interface phases inside the particle-reinforced composite material. This approximation method allows the microstructural characteristics of the interface layer, including chemical reaction products and physical coatings on the particle surface, to be taken into account in the model. Since the presence of the interface phase significantly affects the mechanical properties of the composite material, even if its thickness is very small, accurate geometric modeling is crucial to the accuracy of the overall simulation. In this multiphase interface structure, by combining the assumed stress hybridization unit method with the Maxwell stress function, the complex stress field changes between different material phases can be described, and the interface stress concentration phenomenon and force transfer behavior can be accurately captured.

[0024] The present invention ensures the strictness of force balance and energy conservation in numerical calculations by introducing the modified complementary energy functional of the Lagrange multiplier, thereby enhancing the numerical stability of the simulation. This method improves the calculation deficiencies of traditional finite element analysis under nonlinear conditions and avoids possible numerical divergence problems. During the discretization process, the model is meshed using Delaunay triangulation technology, and numerical integration is performed within each tetrahedral unit to effectively construct a global stiffness matrix. This method can not only reduce the accumulation of errors in the calculation process, but also reduce the resources and time required for calculation. Compared with traditional finite element analysis, the calculation efficiency is greatly improved and it is suitable for processing larger-scale three-dimensional structure simulations.

[0025] In the simulation of thermo-elastic-plastic behavior, the present invention combines temperature increment with the material's thermal expansion coefficient, employs the Mises yield condition, and develops a detailed description using equivalent plastic strain, comprehensively modeling the material's response under different thermo-mechanical conditions. By selecting appropriate displacement interpolation functions to ensure the accuracy of boundary conditions, the model accurately reflects the deformation and strain distribution of the material under the combined action of heat and force. This comprehensive thermo-elastic-plastic analysis capability enables the model to more realistically reflect the complex mechanical behavior of materials in actual engineering, especially when the temperature-induced stress and deformation are highly correlated, providing scientific guidance for material design and optimization. The obtained results were compared with those calculated using the commercial finite element software ABAQUS and MARC. The accuracy and effectiveness of the method were verified by comparing stress contours and stress path diagrams of the thermo-elastic and plastic results. Calculations implemented in Fortran demonstrate the method's outstanding performance under complex models, particularly random distribution models, further verifying its universality and practicality in practical engineering applications. The entire method significantly improves computational efficiency while ensuring accuracy, demonstrating the advantages of VCFEM in handling complex material models.

[0026] Of course, any product implementing the present invention does not necessarily need to achieve all of the advantages described above at the same time. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0028] Figure 1 Schematic diagram of Voronoi cell with interface;

[0029] Figure 2 Schematic diagram of Delaunay tetrahedron partitioning;

[0030] Figure 3 It is a schematic diagram of a space triangle;

[0031] Figure 4 It is a schematic diagram of a polygonal approximation ellipsoid;

[0032] Figure 5 Schematic diagram of boundary conditions;

[0033] Figure 6 Schematic diagram of mesh division;

[0034] Figure 7 is the stress-Y diagram of the thermal model;

[0035] Figure 8 Schematic diagram of stress-Y (model does not yield) of the plastic model;

[0036] Figure 9 Schematic diagram of stress-Y (matrix yield) of the plasticity model;

[0037] Figure 10 Schematic diagram of stress-Y (matrix and interface yield) of the plasticity model;

[0038] Figure 11 Schematic diagram of stress path location;

[0039] Figure 12 is a schematic diagram of the stress path;

[0040] Figure 13 Schematic diagram of boundary conditions;

[0041] Figure 14 Schematic diagram of model mesh division;

[0042] Figure 15 is the stress-Y diagram of the thermal model;

[0043] Figure 16 Schematic diagram of stress-Y (model does not yield) of the plastic model;

[0044] Figure 17 Schematic diagram of stress-Y (matrix yield) of the plasticity model;

[0045] Figure 18 Schematic diagram of stress-Y (matrix and interface yield) of the plasticity model;

[0046] Figure 19 Schematic diagram of stress path location;

[0047] Figure 20 is a schematic diagram of the stress path;

[0048] Figure 21Provide a schematic diagram of the model establishment and boundary conditions;

[0049] Figure 22 Schematic diagram of stress cloud diagram. DETAILED DESCRIPTION

[0050] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts shall fall within the scope of protection of the present invention.

[0051] Example 1

[0052] The three-dimensional thermo-elastic-plastic Voronoi element simulation method of the composite material containing interfacial particles described in this embodiment includes the following steps:

[0053] S1: Establish a calculation model based on dimensions, boundary conditions, and material parameters;

[0054] S2: Discretize the model into several Voronoi cells containing one inclusion;

[0055] S3: Divide each Voronoi cell into several Delaunay tetrahedrons;

[0056] S4: Derivation of modified complementary energy functional based on the minimum complementary energy principle and thermoplastic constitutive equation;

[0057] S5: Select high-order complete polynomials as trial functions to construct stress field functions of three-phase materials considering ellipsoidal interfaces;

[0058] S6: Obtain the element stiffness matrix by varying and solving the G matrix and H matrix;

[0059] S7: Determine the iterative expressions for stress parameter β and displacement d;

[0060] S8: Calculate the node stress based on the external node displacement and the element stiffness matrix to obtain the internal node displacement;

[0061] S9: Calculate stress parameter dβ = -H -1 G th +H -1 Gdq;

[0062] S10: Take the partial derivative of the stress function to obtain the matrix P and calculate the full-field stress:

[0063] [Δσ x Δσ y Δσ zΔτ xy Δτ yz Δτ zx ]=PΔβ

[0064] Complete the simulation of three-dimensional thermo-elastic-plastic Voronoi elements.

[0065] Example 2

[0066] Derivation of the modified complementary energy functional in incremental form

[0067] According to the assumed stress hybrid element method, 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 takes the stress function as the field variable. The expression of the energy compensation function is

[0068]

[0069] Where S is the elastic flexibility tensor, Γ is the given displacement boundary, T is the surface force, Ω is the unit integration domain, is the displacement. The unit stress is expressed by the stress function, and the node stress function is used as the basic unknown. The stress function is guaranteed to be continuous on the boundary between units, so the equilibrium condition can be satisfied everywhere. The disadvantage of equilibrium model I is that the distribution of stress and strain is obtained by differentiating the stress function, while the displacement solution is obtained by integrating the geometric equation, which depends on the choice of the integration path. Therefore, the displacement solution obtained by this method is not unique. In addition, the stress function is not convenient for intuitively explaining the physical meaning of the boundary conditions. The Voronoi unit with interface used in this embodiment is as follows Figure 1 shown.

[0070] Based on the two-dimensional modified complementary energy function, the unit consists of three stages, including the matrix stage, the interface stage and the inclusion stage. At a given traction boundary:

[0071]

[0072] n m is the outer normal of the cell boundary, is the given surface force.

[0073] On the boundary between the inclusion phase and the interface phase inside the unit:

[0074] n c ·(σ i +Δσ i )-n c ·(σ c +Δσ c )=0 (3)

[0075] At the boundary between the interfacial phase and the matrix phase:

[0076] n i ·(σ m +Δσ m )-n i ·(σ i +Δσ i )=0 (4)

[0077] n c is the external normal of the inclusion phase, n i is the external normal of the interface phase.

[0078] By introducing the Lagrange multiplier, the modified coenergy functional can be obtained:

[0079]

[0080]

[0081] The residual energy increment ΔB satisfies:

[0082]

[0083] The self-equilibrium stress field in the element is introduced in the finite element method:

[0084] Δσ=pΔβ(7)

[0085] Δσ is a column vector containing six stress components, and Δβ is a column vector containing m unknown stress coefficients.

[0086] Δu=LΔq(8)

[0087] The boundary displacement Δu is an interpolation of the generalized displacement Δq of the node.

[0088]

[0089] Formula 9 can be simplified to

[0090]

[0091]

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

[0093]

[0094] Can get

[0095]

[0096] Where H and G are obtained by integrating the tetrahedron, N Tis the transposed matrix of the external normal.

[0097] In the process of solving equations 12, 13 and 14, the stress coefficient increment of step i can be expressed as the stress coefficient increment of step i-1 plus the stress coefficient change of step i.

[0098]

[0099] Similarly, the displacement increment can be expressed as

[0100]

[0101] A weak expression of the kinematic relationship can be obtained in each unit:

[0102]

[0103] In formulas 12 to 14

[0104] 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Ω (18)

[0105] G mm =∫ Γ P mT N mT LdΓ. G ii =∫ Γ P iT N iT LdΓ. G cc =∫ Γ P cT N cT LdΓ (19)

[0106] G mi =∫ Γ P mT N iT LdΓ. G ic =∫ Γ P iT N cT LdΓ (20)

[0107]

[0108] G th is the G matrix related to heat, and Δε can be calculated from the constitutive equation.

[0109] The displacement expression of stress parameter can be written from the above formula:

[0110]

[0111] According to the first-order variation of the total energy of the system with respect to the node displacement dq, Πe can obtain the weak expression of the corresponding surface force boundary condition:

[0112]

[0113] Substituting into formula 22, we get

[0114]

[0115] This results in a system of equations for solving the generalized displacement, where the element stiffness matrix is

[0116]

[0117] dF is the node stress increment, which is calculated by decomposing the element stiffness matrix

[0118]

[0119]

[0120] is the nodal stress increment of the interface phase in the i-th incremental step, is the nodal stress increment of the inclusion phase in the i-th increment, K 22 -1 is the inverse matrix of the stiffness matrix of the interface and the inclusion, K 12 T is the inverse matrix of the stiffness matrix of the inclusion and interface influence on the matrix displacement, and the node equivalent stress increment is expressed as

[0121]

[0122] Thus, in Formula 29, the calculation only involves the displacement of the external nodes of the element, which greatly reduces the calculation amount of the entire stiffness matrix.

[0123] Constitutive Relationship between Thermal Strain and Plasticity

[0124] Elastic-plastic constitutive relation

[0125] This embodiment adopts the strain increment theory or the flow theory of plasticity. In this case, it can be assumed that the principal axis of the strain increment coincides with the principal axis of the stress deviation at a certain moment, so the component of the stress deviation is the same as the component of the strain increment. The principal axis of the stress deviation coincides with the principal axis of the stress deviation at a certain moment, so the component of the stress deviation is proportional to the component of the strain increment at that moment. The Mises yield condition can be expressed as follows

[0126]

[0127] σ so is the material yield stress, s ij =σ ij -σ m δ ij is the stress deviator tensor, σ m =(σ 11 +σ 22 +σ 33 ) / 3 is the hydrostatic pressure, s ij is the Kronecker triangular tensor.

[0128] The flow rule specifies the relationship between the components of the plastic strain increment and the components of the stress and stress increment. Using the Mises yield condition, it can be written as

[0129]

[0130] Where dλ is an undetermined finite positive quantity called the plastic multiplier. Combined with the definition of equivalent stress, we can get

[0131]

[0132]

[0133] The Mises yield condition is introduced as the equivalent plastic strain increment

[0134]

[0135] The isotropic hardening rule states that after a material enters plastic deformation, the loading surface expands uniformly in all directions, while its center of mass and its orientation in stress space remain unchanged. If the Mises buckling condition is used, the yield function after isotropic hardening can be expressed as

[0136]

[0137] where σ ij Current elastic-plastic stress, which is a function of the equivalent plastic strain function, It can be obtained from the σ-ε curve of the uniaxial tensile test of the material.

[0138]

[0139] H ′ is the plastic modulus of the material, also known as the hardening coefficient. Its relationship with the elastic modulus E and the tangent modulus is

[0140]

[0141] If the elastic strain is of the same order of magnitude as the plastic strain, ignoring the elastic strain will often lead to large errors. Therefore, it is recommended to consider the elastic strain in the plastic zone, that is, the total strain increment consists of two parts, which are expressed as

[0142]

[0143] dσ ij and Satisfying the generalized Hooke's law, in elastic stress and strain

[0144]

[0145] According to the hardening criteria

[0146]

[0147] The relationship between the stress increment matrix and the strain increment matrix is

[0148] [dε]=D[dσ] (44)

[0149] D=D e +D p (45)

[0150] in

[0151]

[0152] thermal strain

[0153] In the classical thermoelastic-plastic theory, the temperature strain increment can be expressed as

[0154] Δε th ij =αΔTδ ij (48)

[0155] α is the coefficient of thermal expansion, ΔT is the temperature increment;

[0156] Tectonic stress function

[0157] In three-dimensional problems, it is necessary to consider how to reasonably select stress functions to construct the stress field of the unit. For stress functions, the matrix stress function is expressed as represents the matrix stress function The interface stress function is Regarding the interaction force function, the interaction force stress function of the matrix is ​​expressed as The interaction force stress function of the interface is expressed as The construction of polynomial stress function and interaction force function should meet the following requirements: and The columns in are linearly independent, and only if this condition is met can the matrix H be invertible. The shape of the enclosing volume should be considered when selecting the stress function polynomial. The interaction stress term should approach zero far from the contact boundary between the two phases, and the continuity condition of the traction should be met at the contact boundary between the two phases.

[0158]

[0159] The interfacial interaction stress function obtained by Ghosh is further simplified as follows:

[0160]

[0161] Close to the interface, Stay away from the interface, Within the interface, after constructing the stress polynomials for the matrix and inclusions, the stress field is derived from the Maxwell stress function. Substituting the matrix stress polynomial into the matrix stress polynomial, the stress fields for the matrix, interface layer, and inclusions are obtained as follows:

[0162]

[0163]

[0164] in is the scaled (x, y, z) coordinate. This embodiment uses the Maxwell stress function. Assuming that the stress function is a perfect polynomial, the stiffness matrix should be reversible. According to the relationship between the number of stress function terms and the number of degrees of freedom of the rigid body, d>bc can be selected, where d is the total number of stress function terms, b is the total number of degrees of freedom of the rigid body, and c is the degree of freedom. According to previous calculations and studies, the basic number of Maxwell stress functions used for three-dimensional element calculations is at least 31.

[0165] Φ(x,y,z)=a0+a1x+a2y+a3z+β1x 2 +β2y 2 +β3z 2 +β4xy+β5yz+β6zx+β7x 3 +β8y3 +β9z 3 +β 10 x 2 y+β 11 xy 2 +β 12 y 2 z+β 13 yz 2 +β 14 z 2 x+β 15 zx 2 +β 16 xyz+β 17 x 4 +β 16 y 4 +β 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 +…(54)

[0166] This function takes the form of the sum of perfect polynomials of all orders, and the total number of terms in the stress function is Φ, which is 6+10+15+...+(m+1)(m+2) / 2. Excluding the 0th- and 1st-order terms, the number of terms in Equation 54 is 31. Based on this stress function, the six Maxwell stress components can be obtained:

[0167]

[0168] By taking partial derivatives of the Maxwell stress function Φ, the stress components in the element can be obtained.

[0169] [Δσ x Δσ yΔσ z Δτ xy Δτ yz Δτ zx ]=PΔβ (56)

[0170] Numerical integration of H and G matrices

[0171] According to equations 18 and 19, the integral of the H matrix is ​​integrated into the tetrahedral domain, so the integrals of each small tetrahedron are accumulated to obtain the H matrix of the entire unit:

[0172]

[0173] The local coordinates are expressed as

[0174]

[0175] In the formula, n tet Represents the number of tetrahedrons obtained after mesh division, n int represents the number of Hammer integration points of the tetrahedron used, and ωj represents the weight coefficient of the integration point numbered j. Ji represents the Jacobian coefficient of the tetrahedron numbered i;

[0176] P(x, y, z) is a coefficient matrix that depends only on the coordinates of the integration points. Similarly, the triangular surface integral of the G matrix on the tetrahedron also needs to be accumulated:

[0177]

[0178] The local coordinates are expressed as

[0179]

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

[0181]

[0182] Displacement interpolation function

[0183] This method divides the three-dimensional Voronoi cell containing the interface layer into multiple Delaunay tetrahedrons, such as Figure 2 As shown, each tetrahedron has four triangular surfaces. When calculating the displacement interpolation of the tetrahedron surface, consider the spatial triangles, such as Figure 3As shown, they are numbered i, m, j in the rectangular coordinate system. The coordinates of each node are (xi, yi, zi), (xm, ym, zm) and (xj, yj, zj). By interpolating the node displacement, the boundary displacement of the element surface can be obtained.

[0184]

[0185] in Δ is the area of ​​the triangle, which can be calculated based on the coordinates of the three points, and Δ i is the area of ​​the small triangle formed by the two corner points on the edge. u, v, w are the displacements of the node in the x, y, and z directions, and L is the shape function. Rewrite Equation 62 in matrix form:

[0186] u=Ld (63)

[0187] Example 3

[0188] Based on VCFEM, incremental modified plasticity and thermal residual energy functionals of three-dimensional Voronoi elements containing interfaces are constructed to solve nonlinear problems of particle reinforced composite materials containing interfaces. In this embodiment, the Voronoi elements containing interfaces are programmed using Fortran language, and the results are compared with those of commercial finite element analysis software ABAQUS and MARC. With the rapid development of modern science and technology, large-scale engineering problems are gradually increasing, and more and more people are using commercial software to model and analyze actual engineering problems. ABAQUS and MARC have become the most advanced software in the world in the engineering field due to their excellent simulation performance and huge solution functions. Similarly, ABAQUS and MARC are also widely used in research institutions in various countries. Therefore, it is very convincing to choose ABAQUS and MARC to compare the calculation results of the examples.

[0189] Verifying the Thermoelastic-Plastic Results of a Single Element

[0190] Modeling

[0191] First, calculate a single cell and place the center point of the inclusion at the center point of the cell. Figure 4 As shown; Figure 4 (a): This unit consists of three parts, namely the matrix phase, the inclusion phase and the interface phase. The unit is a regular hexahedron with a side length of 20 mm. The inclusion phase and the interstitial phase are approximately represented by the inscribed polygon of the sphere, as shown in Figure 4 (b) shown. Figure 4The radius of the inclusion is 2 mm, the thickness of the mesophase is 0.35 mm, and the corresponding material properties are shown in Table 1. The boundary conditions are: when considering plasticity, the displacement in the Y direction is 0.002 mm, and each increment is 10 steps; when considering thermal strain, the temperature increment of the entire unit is 40 degrees Celsius. The specific model boundary conditions are as follows Figure 5 .

[0192] Table 1 Material properties

[0193]

[0194]

[0195] In this embodiment, the results of tens of thousands of units of the traditional displacement finite element method can be obtained with only a few units, which is also the advantage of VCFEM. Figure 6 As shown in Figure 2, the VCFEM model has only one Voronoi element with an interface, while the ABAQUS model and the MARC model consist of 361,105 4-node tetrahedral elements. The number of stress function terms in the VCFEM model is shown in Table 2:

[0196] Table 2 Number of stress function terms

[0197]

[0198] In Table 2, pol: stress function, rec: interaction force function.

[0199] Comparison and verification of calculation results

[0200] The stress cloud diagram of the calculation results is as follows Figures 7 to 10 This example compares the calculation results of ABAQUS, MARC and VCFEM respectively to verify the reliability of the calculation method proposed in the present invention.

[0201] Extract stress data on the vertical centerline path from the stress cloud diagram. The path position is as follows Figure 11 The stress path diagrams of the above four groups of cloud diagrams are as follows: Figure 12 .

[0202] Figure 10 and Figure 12 The calculation results of VCFEM and commercial analysis software are shown. The comparison results show that the calculation results of this method are consistent with those of the commercial software. This proves that the method proposed in this paper is correct for stress calculation of three-dimensional Voronoi elements with interfaces, has high computational efficiency, and does not require solving tens of thousands of elements. Given the same number of cores and threads, VCFEM requires less time. This demonstrates the advantages of VCFEM in terms of simplicity and speed.

[0203] Verifying Multi-element Thermoelastic-Plastic Results

[0204] Modeling

[0205] This embodiment will calculate multiple units based on a single unit with the center point of the inclusion located at the center point of the unit to verify whether the method is still applicable. The multi-unit calculation model uses the unit model in this embodiment for array arrangement, and the unit size and material parameters are the same as the model modeled in this embodiment. The boundary conditions are: when considering plasticity, the given displacement in the Y direction is 0.006 mm, and each increment is 10 steps; when considering thermal strain, the temperature increment of the entire element is 120 degrees Celsius. The specific model boundary conditions are as follows Figure 13 .

[0206] The grid division of the computational model is as follows Figure 4 The VCFEM model has only 27 Voronoi elements with interfaces, while the ABAQUS and MARC models consist of 906,012 four-node tetrahedral elements. The VCFEM model has 64 external nodes and 1,296 internal nodes (inclusion nodes and interface nodes). The number of stress function terms in the VCFEM model is shown in Table 3.

[0207] Table 3 Number of stress function terms

[0208]

[0209] In Table 3, pol: stress function, rec: interaction force function.

[0210] Comparison and verification of calculation results

[0211] The cloud diagram of the calculation results is as follows Figures 15 to 18 This example compares the calculation result distributions of ABAQUS, MARC, and VCFEM to verify the reliability of the calculation method proposed in this example in the multi-element model.

[0212] Extract stress data on the vertical centerline path from the cloud map. The path position is as follows Figure 19 The stress path diagrams of the above four cloud diagrams are as follows: Figure 10 shown.

[0213] like Figures 15 to 20 As shown, the VCFEM results are consistent with those of the traditional displacement finite element method. A comparison of the elements also demonstrates that the method performs well when calculating multi-element models. Continuity conditions are satisfied both within and between elements. The multi-element model demonstrates the VCFEM's superior computational efficiency compared to the traditional displacement finite element method.

[0214] Calculation of random distribution models

[0215] Modeling

[0216] In order to simulate the random distribution of inclusions in actual composite reinforcement materials with interface particles, this embodiment adopts a random distribution model. The model is still modeled with 6 polygons. 20 inclusions are randomly distributed in a regular hexahedron with a side length of 20 mm, and the interface phase with uniform thickness is distributed between the matrix and the inclusions. The mesh after model meshing is as follows Figure 11 As shown. In this model, the number of Voronoi elements with interfaces is 20, the number of external nodes is 109, and the number of internal nodes is 960. The boundary conditions are the same as those in the thermoelastic-plastic result modeling section for verifying a single element in Example 3. The boundary conditions in the surface normal direction apply to all nodes on the corresponding surface. The displacement on the given displacement boundary in each incremental layout is 0.001 mm. There are 15 steps in total. The temperature increment is 40 degrees Celsius. The number of stress terms is shown in Table 6. The unit node information is shown in Tables 4 and 5.

[0217] Table 4 Unit node information

[0218]

[0219] Table 5 Unit node information

[0220]

[0221] Table 6 Number of stress function terms

[0222]

[0223] In Table 6, pol: stress function, rec: interaction force function.

[0224] Calculation results display

[0225] Stress cloud diagram Figure 22 As shown, it can be seen from the stress cloud diagram of the calculation results that even in the case of random distribution, when the number of nodes and the number of faces of each unit are not the same, each unit uses the same preset high-order stress field. Although the complete polynomials solved are the same, the stress distribution details of each unit can still be captured. In this model, the number of internal points of all units is 48, the maximum number of external points is 24, and the minimum is 10. The maximum degree of freedom of a unit is 216, and the minimum degree of freedom is 174. Even if the determination conditions of the high-order stress field required for the unit are quite different, according to the selection conditions of the stress function terms in Example 2, the selection of the stress function takes into account the unit with the largest degree of freedom, and the stress distribution results of all units in the model also conform to actual rules. Practice has proved that the method of the present invention has universality in the selection of stress function and can be used to solve the analysis problem of composite reinforced materials containing interface particles in actual engineering.

[0226] The preferred embodiments of the application disclosed above are only to facilitate the elucidation of the application. The preferred embodiments do not describe all the details of the application and limit the application to the specific embodiments described. Obviously, many modifications and variations can be made in light of the teachings above. The description is chosen and described in order to best explain the principles of the application and its practical application to thereby enable others skilled in the art to best utilize the application and get the best results from the application. The application is only limited by the claims and their full scope and equivalents.

Claims

1. A three-dimensional thermo-elastic-plastic Voronoi element simulation method for composite materials containing interfacial phase particles, characterized in that: The following steps are involved: S1: Establish a calculation model based on dimensions, boundary conditions, and material parameters; S2: Discretize the model into several Voronoi cells containing one inclusion; S3: Divide each Voronoi cell into several Delaunay tetrahedrons; S4: Derivation of modified complementary energy functional based on the minimum complementary energy principle and thermoplastic constitutive equation; S5: Select high-order complete polynomials as trial functions to construct stress field functions of three-phase materials considering ellipsoidal interfaces; S6: Obtain the element stiffness matrix by varying and solving the G matrix and H matrix; S7: Determine stress parameters and displacement The iterative expression of S8: Calculate the node stress based on the external node displacement and the element stiffness matrix to obtain the internal node displacement; S9: Calculate stress parameters ; in is the stress parameter, is the G matrix integral, is the node displacement; S10: Take the partial derivative of the stress function to obtain the matrix P and calculate the full-field stress: ； in is the coordinate related P matrix, is the stress parameter increment; is the normal stress increment in the x-direction, is the normal stress increment in the y direction, is the z-direction normal stress increment, 、 、 is the shear stress increment; Complete the simulation of three-dimensional thermo-elastic-plastic Voronoi elements.

2. The three-dimensional thermo-elastic-plastic Voronoi element simulation method of a particle-reinforced composite material containing an interface phase according to claim 1, characterized in that: The steps S1-S3 specifically include: According to the assumed stress hybrid element method, the equilibrium model in the hybrid element method is 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 takes the stress function as the field variable; the expression of the energy compensation function is （1）； in is the elastic compliance tensor, is a given displacement boundary, For surface force, is the unit integration domain, , the element stress is expressed as a stress function, and the nodal stress function is used as the basic unknown; the stress function is guaranteed to be continuous on the boundaries between elements, so the equilibrium condition can be satisfied everywhere; Based on the two-dimensional modified complementary energy function, the unit consists of three stages, including the matrix stage, the interface stage and the inclusion stage; at a given traction boundary: （2）； is the outer normal of the cell boundary, is the given surface force; On the boundary between the inclusion phase and the interface phase inside the unit: （3）； is the external normal of the inclusion phase, is the external normal of the interface phase; At the boundary between the interfacial phase and the matrix phase: （4）； By introducing the Lagrange multiplier, the modified coenergy functional can be obtained: ； ； ； （5）； Residual energy increment satisfy: （6）； The self-equilibrium stress field in the element is introduced in the finite element method: （7）； is a column vector containing six stress components, is a column vector containing m unknown stress coefficients; （8）； Boundary displacement is the generalized displacement of the node Interpolation: ； ； ； ； ； （9）； Formula 9 can be simplified to ； ； ； ； ； ； After discretization, the weak form of the system complementary energy can be obtained; according to the modified complementary energy stationary condition ； Can get ； ； ； Where H and G are obtained by integrating the tetrahedron, N T is the transposed matrix of the external normal; In solving equations 12, 13, and 14, the stress coefficient increment at step i can be expressed as the stress coefficient increment at step i-1 plus the stress coefficient change at step i; ； Similarly, the displacement increment can be expressed as ； A weak expression of the kinematic relationship can be obtained in each unit: ； In formula 12-14 ； ； ； ； is the G matrix related to heat, It can be calculated by the constitutive equation; The displacement expression of stress parameter can be written from the above formula: ； According to the first-order variation of the total energy of the system with respect to the node displacement dq, Π𝑒 can obtain the weak expression of the corresponding surface force boundary condition ； Substituting into formula 22, we get ； This results in a system of equations for solving the generalized displacement, where the element stiffness matrix is ; ； ； is the node stress increment, by decomposing the element stiffness matrix ； ； ； is the nodal stress increment of the interface phase in the i-th incremental step, is the nodal stress increment of the inclusion phase in the i-th increment, is the inverse matrix of the stiffness matrix of the interface and the inclusion, is the inverse matrix of the stiffness matrix of the inclusion and interface influence on the matrix displacement, and the node equivalent stress increment is expressed as 。 3. The three-dimensional thermo-elastic-plastic Voronoi element simulation method of a composite material containing an interfacial phase particle reinforcement according to claim 1, characterized in that: The step S4 specifically includes: using the strain increment theory or the flow theory of plasticity; setting the principal axis of the strain increment to coincide with the principal axis of the stress deviation at a certain moment, so that the component of the stress deviation is the same as the component of the strain increment; the strain component increment coincides with the principal axis of the stress deviation at a certain moment, so the component of the stress deviation is proportional to the component of the strain increment at that moment; the Mises yield condition can be expressed as follows ； is the material yield stress, is the stress deviator tensor, is the hydrostatic pressure, is the Kronecker triangular tensor; The flow rule specifies the relationship between the components of the plastic strain increment and the components of the stress and stress increment; using the Mises yield condition, it can be written as ； Where d𝜆 is an undetermined finite positive quantity called the plastic multiplier; combined with the definition of equivalent stress, we can get ； ； ； The Mises yield condition is introduced as the equivalent plastic strain increment ； The isotropic hardening rule states that after the material enters plastic deformation, the loading surface expands uniformly in all directions, while its shape center and its direction in stress space remain unchanged; if the Mises buckling condition is used, the yield function after isotropic hardening can be expressed as ； in Current elastic-plastic stress, which is a function of the equivalent plastic strain function, c can be obtained from the uniaxial tensile test of the material Obtained from the curve; ； 𝐻 ′ is the plastic modulus of the material, also known as the hardening coefficient; its relationship with the elastic modulus E and the tangent modulus is : ； If the elastic strain is of the same order of magnitude as the plastic strain, ignoring the elastic strain will often lead to large errors. Therefore, it is recommended to consider the elastic strain in the plastic zone, that is, the total strain increment consists of two parts, which are expressed as ； and Satisfying the generalized Hooke's law, in elastic stress and strain ； ； According to the hardening criteria ； The relationship between the stress increment matrix and the strain increment matrix is ； ； in ； ； In the classical thermoelastic-plastic theory, the temperature strain increment is expressed as ； is the coefficient of thermal expansion, is the temperature increment.

4. The three-dimensional thermo-elastic-plastic Voronoi element simulation method of a composite material containing interfacial particles according to claim 1, characterized in that: The step S5 specifically includes: For the stress function, the matrix stress function is expressed as represents the matrix stress function The interface stress function is ; Regarding the interaction force function, the interaction force stress function of the matrix is ​​expressed as The interaction force stress function of the interface is expressed as ; The construction of polynomial stress function and interaction force function meets the following requirements: , and The columns in are linearly independent. Only when this condition is met can the matrix H be reversible. When selecting the stress function polynomial, the shape of the containing body should be considered. The interaction stress term should approach zero far away from the contact boundary of the two phases. The continuity condition of the surface force should be met at the contact boundary of the two phases. ； The interfacial interaction stress function obtained by Ghosh is further simplified as follows: ； Close to the interface, →0 is far away from the interface, >1In 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 stress field of the matrix, interface layer and inclusions, the stress fields are as follows: ; ； ； Where (𝜉, 𝜂, 𝜍) are scaled (x, y, z) coordinates; the Maxwell stress function is selected; setting the stress function to a perfect polynomial should satisfy the condition of the stiffness matrix being invertible; based on the relationship between the number of stress function terms and the number of degrees of freedom of the rigid body, d>bc can be selected, where d is the total number of stress function terms, b is the total number of degrees of freedom of the rigid body, and c is the number of degrees of freedom; based on previous calculations and research, the basic number of Maxwell stress functions used for three-dimensional element calculations is at least 31; ； ； ； ； . ； This function takes the form of the sum of perfect polynomials of all orders. The total number of terms in the stress function is Φ6+10+15+...+(m+1)(m+2) / 2. Excluding the 0th- and 1st-order terms, the number of terms in Equation 54 is 31. Based on this stress function, the six Maxwell stress components can be obtained: ； By using the Maxwell stress function By taking partial derivatives, we can get the stress components in the element: 。 5. The three-dimensional thermo-elastic-plastic Voronoi element simulation method of a composite material containing interfacial particles as claimed in claim 1, characterized in that: The step S6 specifically includes: According to formulas 18 and 19, the integral of the H matrix is ​​integrated into the tetrahedron domain, so the integrals of each small tetrahedron are accumulated to obtain the H matrix of the entire unit: ； The local coordinates are expressed as ； In the formula 𝑛 𝑡𝑒𝑡 represents the number of tetrahedrons obtained after mesh division, 𝑛 𝑖𝑛𝑡 represents the number of Hammer integration points of the tetrahedron used, and 𝜔𝑗 represents the weight coefficient of the integration point numbered j; 𝐽𝑖 represents the Jacobian coefficient of the tetrahedron numbered i; P(x, y, z) is a coefficient matrix that depends only on the coordinates of the integration points; similarly, the triangular surface integral of the G matrix on the tetrahedron also needs to be accumulated: ； The local coordinates are expressed as ； Can get , , then there will be an additional element scaling factor when computing the stiffness matrix: 。 6. The three-dimensional thermo-elastic-plastic Voronoi element simulation method of a particle-reinforced composite material containing an interfacial phase according to claim 1, characterized in that: Steps S7 and S8 specifically include: The three-dimensional Voronoi element containing the interface layer is divided into multiple Delaunay tetrahedrons, each of which has four triangular surfaces. When calculating the displacement interpolation of the tetrahedron surface, the spatial triangles are considered and numbered as i, m, and j in the rectangular coordinate system. The coordinates of each node are (𝑥𝑖,𝑦𝑖,𝑧𝑖), (𝑥𝑚,𝑦𝑚,𝑧𝑚) and (𝑥𝑗,𝑦𝑗,𝑧𝑗). By interpolating the node displacements, the boundary displacement of the element surface can be obtained. ； in , Δ is the area of ​​the triangle, which can be calculated based on the coordinates of the three points, and is the area of ​​the small triangle formed by the two corner points on the side; rewrite formula 62 into matrix form: 。

Images