Simulation method of three-dimensional thermal-elastic-plastic Voronoi unit containing interface phase particle reinforced composite material

Abstract

The invention relates to the technical field of a particle reinforced composite material containing an interface phase, in particular to a simulation method of a three-dimensional thermal-elastic-plastic Voronoi unit containing the particle reinforced composite material containing the interface phase, which comprises the following steps of: generating a three-dimensional Voronoi unit model through a polygon approximation method of an embedded sphere; and the real geometric configuration of the internal matrix, inclusions and interface phase of the particle reinforced composite material can be accurately reproduced. By introducing the corrected complementary energy functional of the Lagrange multiplier, the strictness of force balance and energy conservation in numerical calculation is ensured, so that the numerical stability of simulation is enhanced. In the simulation of thermal-elastic-plastic behaviors, the method combines the temperature increment and the thermal expansion coefficient of the material, adopts the Mises yield condition and the equivalent plastic strain expansion specific description, and performs comprehensive modeling on the reaction of the material under different thermal mechanical conditions. According to the method, the calculation efficiency is remarkably improved, and the advantages of the VCFEM in processing a complex material model are shown.

Description

Technical Field

[0001] The invention relates to the technical field of composite materials containing interface phase particles, and in particular to a three-dimensional thermo-elastic-plastic Voronoi unit simulation method of composite materials containing interface phase particles. Background Art

[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 effect on the overall mechanical properties of the composite material; Gu, JH et al.'s research shows that the damping capacity of elastic interface composite materials increases with the increase of elastic modulus. For plastic interface phases, weak interface phases or hard interface phases are beneficial to improving the overall damping of composite materials; Hassanzadeh Aghdam M K et al.'s research results show that the effect of the interface on the transverse thermal properties of fiber composites is more significant than the effect on the longitudinal thermal properties of fiber composites.

[0003] In the study of thermal problems, the results of Buryachenko VA showed the relationship between thermal expansion, stored energy and average thermoelastic strain in components. Wang H et al. proposed a hybrid finite element scheme with the basic solution as the kernel function, considered a representative volume unit with inclusions, and studied the influence of inclusions on the global temperature distribution and the interaction between inclusions; Kaur et al. studied the photothermoelastic interaction of solid cylinders in strong magnetic fields, and explained the influence of thermoelastic theory on 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. Based on the linear elastic theory, Zhang Rui established a two-dimensional Voronoi unit finite element model to study the mesoscopic and macroscopic mechanical behavior of heterogeneous materials. The results of VCFEM were compared with the analytical solutions and numerical results of standard finite element analysis to confirm its effectiveness. Hao et al. used a modified complementary energy function composed of plastic strain and thermal strain to numerically simulate the fatigue behavior of porous materials under different conditions; Mao Chao et al. proposed a new extended VCFEM method to simulate the generation of particle damage, which reflects the accuracy of determining stress concentration and coarse mesh discretization. Rao et al. used VCFEM to study the plastic strain, thermal strain and creep strain of particle composites.

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

[0006] At present, no scholars have used the stress hybrid element method to analyze the nonlinear analysis of particle reinforced composites with interfaces. Compared with the traditional displacement finite element method, the latter requires a large number of elements to obtain accurate results, which shows the superiority of this method 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 a particle-reinforced composite material containing an interface phase. The three-dimensional Voronoi unit model generated by the polygonal approximation method of embedded spheres can accurately reproduce the real 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 containing an interface phase particle reinforcement 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 unit 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 according to 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 3D 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 surface of the particles, 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 hybrid 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 calculation by introducing the modified complementary energy functional of Lagrange multipliers, 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. In the discretization process, the model is meshed using Delaunay triangulation technology, and numerical integration is performed in each tetrahedral unit to effectively construct a global stiffness matrix. This method can not only reduce the error accumulation 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 the temperature increment with the thermal expansion coefficient of the material, adopts the Mises yield condition, and the equivalent plastic strain to develop a specific description, and comprehensively models the response of the material under different thermo-mechanical conditions. By selecting a suitable displacement interpolation function to ensure the accuracy of the boundary conditions, the model can correctly reflect 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 the design and optimization of materials. The obtained results are compared with the calculation results of commercial finite element software ABAQUS and MARC. The accuracy and effectiveness of the method are verified by comparing the stress cloud map and stress path map of the thermoelastic and plastic results. The calculations implemented by Fortran show the excellent performance of this method under complex models, especially in random distribution models, which further verifies its universality and practicality in actual engineering applications. While ensuring accuracy, the entire method significantly improves the calculation efficiency, showing the advantages of VCFEM in processing 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 accompanying drawings required for describing the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying creative work.

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

[0029] Figure 2 This is a schematic diagram of Delaunay tetrahedron partitioning;

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

[0031] Figure 4 This is a schematic diagram of a polygonal approximation to an ellipsoid;

[0032] Figure 5 is a schematic diagram of boundary conditions;

[0033] Figure 6 This is a schematic diagram of grid division;

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

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

[0036] Fig. 9 It is the schematic diagram of stress-Y (matrix yield) of the plasticity model;

[0037] Fig.10 It is a schematic diagram of stress-Y (matrix and interface yield) of the plasticity model;

[0038] Fig.11 is a schematic diagram of the stress path location;

[0039] Fig.12 is a schematic diagram of stress path;

[0040] Fig.13 is a schematic diagram of boundary conditions;

[0041] Fig.14 Schematic diagram of model mesh division;

[0042] Fig.15 It is the stress-Y schematic diagram of the thermal model;

[0043] Fig.16 It is the stress-Y diagram of the plastic model (the model does not yield);

[0044] Fig.17 It is the schematic diagram of stress-Y (matrix yield) of the plasticity model;

[0045] Fig.18 It is a schematic diagram of stress-Y (matrix and interface yield) of the plasticity model;

[0046] Fig.19 is a schematic diagram of the stress path location;

[0047] Fig. 20 is a schematic diagram of stress path;

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

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

[0050] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0051] Example 1

[0052] The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing interface phase 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 unit 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 according to 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 3D 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 with 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 matrix stage, interface stage and inclusion stage. At a given traction boundary:

[0071]

[0072] n m is the outer normal of the cell boundary, is a 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 interface 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 complementary energy 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 nodes.

[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 at the i-th step can be expressed as the stress coefficient increment at the i-1-th step plus the stress coefficient change at the i-th step.

[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 ith incremental step, is the node stress increment of the inclusion phase in the ith increment step, 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 effects 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, thereby greatly reducing the amount of calculation 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 plastic mechanics. 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 rules specify 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 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

[0136]

[0137] where σ ij Current elastic-plastic stress, which is a function of the equivalent plastic strain The function of 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. 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, and the continuity condition of the surface force should be met at the contact boundary of 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. In the interface, after constructing the stress polynomials of the matrix and inclusions, the stress field is derived from the Maxwell stress function. Substituting the stress polynomial of the matrix into the stress field of the matrix, interface layer and inclusions, the following is obtained:

[0162]

[0163]

[0164] in is the scaled (x, y, z) coordinate. The present 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 terms in the stress function and the number of degrees of freedom of the rigid body, d>bc can be selected, where d is the total number of terms in the stress function, 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 terms in the Maxwell stress function used for three-dimensional unit 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 adding complete polynomials of all orders, and 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 formula 54 is 31. Based on this stress function, 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 meshing, 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, the spatial triangle is considered, 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 unit 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, z directions, and L is the shape function. Rewrite formula 62 into matrix form:

[0186] u=Ld (63)

[0187] Example 3

[0188] On the basis of VCFEM, the incremental modified plasticity and thermal residual energy functionals of the three-dimensional Voronoi unit containing the interface are constructed to solve the nonlinear problem of the particle reinforced composite material containing the interface. The present embodiment uses Fortran language to program the Voronoi unit containing the interface, and compares its results with the results 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 use commercial software to model and analyze practical engineering problems. ABAQUS and MARC have become the most advanced software in the world in the engineering field because of their excellent simulation performance and huge solution function. Similarly, ABAQUS and MARC are also widely used in research institutions of various countries. Therefore, it is very convincing to choose ABAQUS and MARC to compare the calculation results of the examples.

[0189] Verifying Thermoelastic-Plastic Results for a Single Element

[0190] Modeling

[0191] First, a single cell is calculated, with the center point of the inclusion placed at the center point of the cell. Figure 4 As shown; Figure 4 (a): This unit consists of three parts, namely, matrix phase, inclusion phase and interface phase. This unit is a regular hexahedron with a side length of 20 mm. The inclusion phase and interstitial phase are approximately represented by the inscribed polygon of a sphere, as shown in Figure 4 (b) as shown. Figure 4The radius of the inclusion is 2 mm, the thickness of the intermediate phase 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 in the traditional displacement finite element method can be obtained with very few units, which is also the advantage of VCFEM. Figure 6 As shown in Figure 2, the VCFEM model has only one Voronoi unit with an interface, while the ABAQUS model and the MARC model consist of 361,105 4-node tetrahedral units. The number of stress function items of 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 embodiment 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 path perpendicular to the center line from the stress cloud diagram. The path location is as follows: Fig.11 The stress path diagrams of the above four groups of cloud diagrams are as follows: Fig.12 .

[0202] Fig.10 and Fig.12 The calculation results of VCFEM and commercial analysis software are shown. From the analysis of the comparison results, the calculation results of this method are consistent with those of the commercial software. This proves that the method proposed in the present invention is correct for the stress calculation of three-dimensional Voronoi cells with interfaces, has high calculation efficiency, and does not require solving tens of thousands of cells. Under the same number of cores and threads, VCFEM takes less time. This reflects the advantages of VCFEM in simplicity and speed.

[0203] Verifying Multi-element Thermoelastic-Plastic Results

[0204] Modeling

[0205] This example 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 example for array arrangement, and the unit size and material parameters are the same as the model modeled in this example. 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 Fig.13 .

[0206] The meshing of the computational model is as follows: Figure 4 As shown in Table 3, the VCFEM model has only 27 Voronoi elements with interfaces, while the ABAQUS model and the MARC model consist of 906,012 4-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 items of 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-unit model.

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

[0213] like Figures 15 to 20 As shown in the figure, the results of VCFEM are consistent with those of the traditional displacement finite element method. Combined with the comparison of the elements, this method also has excellent performance when calculating multi-element models. The continuity conditions are met both within and between elements. Compared with the traditional displacement finite element method, the multi-element calculation model shows the advantages of the VCFEM method in terms of computational efficiency.

[0214] Calculation of random distribution models

[0215] Modeling

[0216] In order to simulate the random distribution of inclusions in the actual composite reinforcement material 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 Fig.11 As shown. In this model, the number of Voronoi cells 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 plasticity result modeling section for verifying a single unit in Example 3. The boundary conditions in the surface normal direction apply to all nodes on the corresponding surface. The displacement on a 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 Fig. 22 As shown, it can be seen from the stress cloud diagram of the calculation result that even in the case of random distribution, when the number of nodes and the number of faces of each unit are different, 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 the unit is 216, and the minimum degree of freedom is 174. Even if the determination conditions of the high-order stress field required by 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 the actual law. 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 present invention disclosed above are only used to help illustrate the present invention. The preferred embodiments do not describe all the details in detail, nor do they limit the invention to the specific implementation methods described. Obviously, many modifications and changes can be made according to the content of this specification. This specification selects and specifically describes these embodiments in order to better explain the principles and practical applications of the present invention, so that those skilled in the art can understand and use the present invention well. The present invention is limited only by the claims and their full scope and equivalents.

Claims

1. A three-dimensional thermo-elastic-plastic Voronoi unit simulation method for composite materials with 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 unit 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 the iterative expressions for stress parameter β and displacement d; S8: Calculate the node stress according to the external node displacement and the element stiffness matrix to obtain the internal node displacement; S9: Calculate stress parameter dβ = -H -1 G th +H -1 Gdq; S10: Take the partial derivative of the stress function to obtain the matrix P and calculate the full-field stress: [D.O. x Board y Board z Dt xy Dt yz Dt zx ]=PΔβ Complete the simulation of 3D thermo-elastic-plastic Voronoi elements.

2. The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing the interface phase particles as claimed in 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 Where S is the elastic flexibility tensor, Γ is the given displacement boundary, T is the surface force, Ω is the unit integration domain, The displacement is the unit stress, 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; Based on the two-dimensional modified complementary energy function, the unit consists of three stages, including matrix stage, interface stage and inclusion stage; at a given traction boundary: n m 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: n c ·(s i +Ds i )-n c ·(s c +Ds c )=0 (3) n c is the external normal of the inclusion phase, n i is the external normal of the interface phase; At the boundary between the interface phase and the matrix phase: n i ·(s m +Ds m )-n i ·(s i +Ds i )=0 (4) By introducing the Lagrange multiplier, the modified complementary energy functional can be obtained: The residual energy increment ΔB satisfies: The self-equilibrium stress field in the element is introduced in the finite element method: Δσ=pΔβ(7) Δσ is a column vector containing six stress components, and Δβ is a column vector containing m unknown stress coefficients; Δu=LΔq(8) The boundary displacement Δu is an interpolation of the generalized displacement Δd of the node: 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 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) G mm =∫ Γ P mT N mT LdΓ.G ii =∫ Γ P iT N iT LdΓ.G cc =∫ Γ P cT N cT LdΓ (19) G mi =∫ Γ P mT N iT LdΓ.G ic =∫ Γ P iT N cT LdΓ (20) G th is the G matrix related to heat, Δε can be calculated from 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, Πe 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 dF is the node stress increment, which is calculated by decomposing the element stiffness matrix is the nodal stress increment of the interface phase in the ith incremental step, is the node stress increment of the inclusion phase in the ith increment step, 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 effects on the matrix displacement, and the node equivalent stress increment is expressed as 3. The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing the interface phase particles as claimed in claim 1, characterized in that: The step S4 specifically includes: using strain increment theory or flow theory of plastic mechanics; setting the principal axis of strain increment to coincide with the principal axis of stress deviation at a certain moment, so that the component of stress deviation is the same as the component of strain increment; strain component increment coincides with the principal axis of the stress deviation at a certain moment, so the component of the stress deviation is in the same proportion as the component of the strain increment at that moment; the Mises yield condition can be expressed as follows σ 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; The flow rules specify the relationship between the components of the plastic strain increment and the components of the stress and stress increment; using the Mises yield condition, this 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 outward 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 where σ ij Current elastic-plastic stress, which is a function of the equivalent plastic strain The function of It can be obtained from the σ-ε curve of the uniaxial tensile test of the material; 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 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 dσ ij 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 [dε]=D[dσ] (44) D=D e +D p (45) in In the classical thermoelastic-plastic theory, the temperature strain increment is expressed as No th ij =αΔTδ ij (48) α is the coefficient of thermal expansion and ΔT is the temperature increment.

4. The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing the interface phase particles as claimed in 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, Stay away from the interface. In the interface, after constructing the stress polynomials of the matrix and inclusions, the stress field is derived from the Maxwell stress function; substituting the stress polynomial of the matrix into the stress field of the matrix, interface layer and inclusions, the stress field is obtained as follows: in is the scaled (x, y, z) coordinate; select the Maxwell stress function; set the stress function to be a complete polynomial, then 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 research, the basic number of Maxwell stress functions used for three-dimensional unit calculations is at least 31; Φ(x,y,z)=a0+a1x+a2y+a3z+β1x 2 +β2u 2 +β3z 2 +β4xy+β5yz+β6zx+β7x 3 +β8y 3 +β9z 3 +b 10 x 2 y+β 11 xyz 2 +b 12 y 2 z+β 13 yz 2 +b 14 z 2 x+b 15 zx 2 +b 16 xyz+b 17 x 4 +b 16 y 4 +b 19 z 4 +b 20 x 3 y+β 21 x 2 y 2 +b 22 xyz 3 +b 26 y 3 z+β 27 y 2 z 2 +b 28 yz 3 +b 23 z 3 x 3 +b 24 z 2 +b 25 zx 3 +b 29 x 2 yz+β 30 xyz 2 z+β 31 x.y.z 2 +....(54) This function takes the form of adding complete 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 formula 54 is 31. Based on this stress function, six Maxwell stress components can be obtained: By taking partial derivatives of the Maxwell stress function Φ, the stress components in the element can be obtained: [Δσ x Δσ y Δσ z Δτ xy Δτ yz Δτ zx ]=PΔβ (56).

5. The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing interfacial phase particles as claimed in claim 1, characterized in that: The step S6 specifically includes: 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: The local coordinates are expressed as In the formula n tet Represents the number of tetrahedrons obtained after meshing, 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; 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 over the tetrahedron also needs to be accumulated: The local coordinates are expressed as You can get H * =H(λ * ) 3 , G * =G(λ * ) 2 , then there will be an additional element scaling factor when computing the stiffness matrix:

6. The three-dimensional thermo-elastic-plastic Voronoi unit simulation method of the composite material containing the interface phase particles as claimed in claim 1, characterized in that: Steps S7 and S8 specifically include: The three-dimensional Voronoi unit 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, j in the rectangular coordinate system. The coordinates of each node are (xi, yi, zi), (xm, ym, zm) and (xj, yj, zj), respectively. By interpolating the node displacement, the boundary displacement of the unit surface can be obtained. 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 side; rewrite formula 62 into matrix form: u=Ld(63).

Images