Simulation of interfacial debonding in hollow particle reinforced composites with VCFEM

The addition of hollow glass microsphere into composites is a method to improve mechanical properties. However, the interfacial debonding of hollow microsphere inevitably causes a decrease in the mechanical properties of the material, which ultimately leads to the failure of the composites. In the numerical simulation of such hollow particle-reinforced composites, the ordinary displacement finite element requires a large number of meshes, which undoubtedly greatly increases the computational cost. In this paper, a new VCFEM is proposed to solve this problem by establishing a two-dimensional Voronoi cell finite element model, deriving the residual energy generalized function of hollow particle-reinforced composites, and calculating the interface debonding. The simulation results are compared with the commercial software MARC, ABAQUS to verify the effectiveness of this VCFEM. The results show that this VCFEM greatly improves the computational efficiency while ensuring the accuracy. Based on this model, this paper also investigates the effect of the generation of interfacial debonding on the overall structure and the effect of different wall thickness of hollow particles on the damage of element debonding.

Hollow particle-reinforced composites have been widely utilized across various industries, such as aviation, military, nuclear power, and construction.Hollow Glass Microspheres (HGM) are characterized by their hollow structure, lightweight, wear resistance, and corrosion resistance.Research indicates that the addition of HGM to composites can enhance their mechanical properties.Swetha, C. 1 conducted a series of quasi-static uniaxial compression experiments to explore the influence of hollow microsphere content and wall thickness on the mechanical characteristics of foams.Scotta, R. 2 investigated the impact of glass microspheres on weight reduction properties and compressive strength of glass microsphere reinforced cementitious composites under quasi-static and high strain rate conditions.Blanco, F. 3 demonstrated that the compressive strength of lightweight cementitious materials surpassed that of conventional perlite lightweight cement upon the incorporation of hollow microspheres.
Tiwari, Vikrant 4 showed that the addition of hollow glass microsphere to the cement matrix significantly improved the acoustic insulation properties of the material.In the investigation of hollow glass microsphere reinforced matrix composites, numerous studies have assessed the mechanical properties by examining the microstructure-property relationship [5][6][7][8] .It has been observed that Eshelby's solution, which is based on the assumption of spatially uniform strain within an elliptical homogeneous envelope, does not accurately predict the elastic properties of the material.Consequently, researchers have incorporated considerations for the nonuniformity of the microstructure to develop more precise mathematical models [9][10][11] .However, challenges such as hollow microsphere aggregation, microsphere cracking, and interface debonding contribute to the deterioration of mechanical properties [12][13][14] , resulting in failure damage of the composites and hindering the model's predictive accuracy.Therefore, in the exploration of hollow glass microsphere reinforced matrix composites, a comprehensive understanding of the influence of damage on the material's mechanical properties is crucial.Given the stress concentration at damaged interfaces and the presence of numerous randomly distributed inclusions of varying sizes and shapes in the composites, conventional displacement-based finite element methods struggle to simulate the complexity of such composites.Consequently, specialized finite element methods have been developed to mitigate computational costs.
The VCFEM (Voronoi Cell Finite Element Method) is a convenient approach that can greatly enhance computational efficiency.In comparison to traditional finite element models, VCFEM is able to reduce the degrees of freedom effectively, leading to a significant improvement in computational productivity.Zhang 15 introduced a hybrid finite element method for analyzing the mechanical response of two-dimensional heterogeneous materials www.nature.com/scientificreports/containing randomly dispersed rigid inclusions of arbitrary sizes.Ghosh and Moorthy proposed the Voronoi cell finite element method [16][17][18] , which has been utilized as a precise and efficient tool for modeling non-uniform microstructures with diverse shapes, sizes, and dispersed heterogeneities.Guo 19 suggested the use of Voronoi cells to simulate creep, thermal, and plastic strains in two-dimensional particle-reinforced composites.Zhang 20 developed a new finite element model, MHSFEM, that simplifies mesh generation and, when combined with VCFEM, enables the analysis of unique microstructures, expanding the application scope of VCFEM.Zhang 21 proposed the SVCFEM to simulate stress concentration at crack tips in homogeneous materials.Ghosh 22 employed VCFEM to study interfacial debonding in fiber-reinforced composites, where the matrix and inclusions are connected through flexible elements known as springs.In contrast, this paper uses shared node pairs to simulate interface debonding.Li and Ghosh 23,24 introduced the XVCFEM to accurately model the propagation of multiple cracks in homogeneous materials.Zhang 21 utilized VCFEM to analyze multi-cracking in homogeneous materials.These methods focus on simulating crack propagation in homogeneous materials.The current study proposes a new VCFEM to simulate interface debonding in hollow particle-reinforced composites.Zhang 25 extended the Voronoi cell finite element method to incorporate arbitrarily distributed inclusion reinforcements with homogeneous phases in composites.Han 26 developed two new elements based on the Voronoi cell finite element model and conducted fracture simulations of materials with large voids under internal pressure.Hao 27 established a twodimensional Voronoi cell finite element model to analyze the fatigue behavior of homogeneous materials with dispersed voids in the matrix.Rao 28 proposed a two-dimensional VCFEM, formulated with plastic, thermal, and creep strain.In conclusion, VCFEM is well-suited for simulating composites with numerous inclusions.Currently, VCFEM is predominantly used for micromechanical damage analysis, such as crack propagation in homogeneous materials, without considering the impact of hollow particle-reinforced composites in crack propagation and interfacial debonding.
Therefore, the objective of this study is to enhance the existing advanced VCFEM to investigate the damage behavior occurring at the interface between the matrix and inclusions in hollow particle-reinforced composites.A novel VCFEM approach is proposed based on a hybrid stress model to characterize the two-phase damage behavior in such composites.A unique cell containing a hollow inclusion is defined based on the principle of complementary residual energy, where stress and displacement are treated as independent fields.The stress field is assumed to satisfy equilibrium equations in advance, represented by a polynomial stress function combined with the inverse stress function within the matrix.A correction function is formulated by introducing Lagrange multipliers into the function.By establishing a relationship between the stress function and the displacement variable, the equations are transformed into a solvable form, and the resulting model is presented following the standard FEM solution format.A new remeshing strategy is implemented to investigate the issue of interfacial debonding in hollow particle-reinforced composites through numerical examples, demonstrating the feasibility of this novel VCFEM approach.

Hollow particle debonding Voronoi cell construction
In the following, the two-dimensional linear elasticity problem is mainly investigated and numerical examples are carried out under plane stress.In order to analyze the effect of debonding between inclusion and matrix on the microevolution of the structure and the macroscopic mechanical properties, we consider a cell model as in Fig. 1.The inclusion and matrix phases of the hollow particle debonding Voronoi cell are denoted as c , m . n e is the direction of the outer normal to the cell boundary ∂� e .n b is the direction of the outer normal to the perfectly bonded matrix-inclusion interface ∂� b .n i is the direction of the outer normal to the hollow inclusion interface ∂� i .n m and n c are the directions of the outer normal at the debonding matrix interface ∂� m and the debonding inclusion interface ∂� c , respectively.
At the bonding interface ∂� b between the inclusion and the matrix, the surface force continuity condition is satisfied.www.nature.com/scientificreports/At the debonding interface ∂� m and ∂� c , the surface force is zero.
At the interface ∂� i between inclusion and hollow part, the surface force is zero.
The modified residual energy generalized function is obtained by applying Lagrange multipliers to apply constraints based on the stress hybrid element generalized function formulation: where u m and u c denote the debonding crack displacements and the cell boundary displacement field is obtained from a linear interpolation function.
Expression of residual energy functional after finite element discretization: In the above equation, the matrix The identification of Gaussian points in the integral computation is realized by dividing multiple quadrangles in the matrix of each phase.P m , P c is the known stress function at the relative point locations of the matrix, inclusion.
By deriving e mc from the stress function β m , β c , we obtain: A weak expression for the relationship between the motions in the element can be obtained: By taking the derivative of e mc with respect to q e , q b , q m , q m , q i , we obtain: (1)

A strategy for matrix-inclusion interface debonding remeshing
In order to realistically simulate the interfacial debonding process of hollow particle reinforced composites, with the change of the debonding region, it is inevitable to change the corresponding mesh structure, so the matrix-inclusion interface remeshing strategy is used to simulate the gradual debonding of the matrix-inclusion interface.Since we mainly study the finite element method, we could have used this very simple and effective strategy for modeling the interface debonding.
Assuming that a particle reinforced composite has produced a section of debonding between the interfaces of matrix and inclusion in an element at the beginning of an incremental step during a certain loading process, according to certain fracture mechanics principles, when the stress at a node reaches or exceeds a critical state, the initial crack tip produces an extension forward at that node extending the interfaces of the matrix and the inclusion.In the first step, the stress state of all the nodes of the interface is determined, and if the node stress reaches or exceeds the critical stress, the interface produces debonding.A second step is connecting the two nodes of the debonding interface to generate a new debonding layer.By the example of A in Fig. 2, when A reaches the critical stress, debonding occurs at the interface.Since the matrix and the inclusion share the same node at point A before the new debonding occurs, this is not a problem in the bonding state.However, with the emergence of the new debonding, the displacement of the part of the matrix at node A and the part of the inclusion at node A are no longer identical.In order to ensure that the crack surface force is zero, we use the node pair A1 and A2 to represent the point A on the inclusion and the point A on the matrix, respectively.The corresponding node variables satisfy the following relationship:

Numerical examples
There are two types of numerical examples in this section: The first type of example is to verify the validity of this VCFEM.The validity is verified by comparing the results with those of the commercial computational software MRAC and ABAQUS.The second type of example is to use this element to simulate a composite material containing a large number of randomly distributed hollow inclusions and to calculate the gradual cracking process at the interface between the matrix and the hollow inclusions under load.Its validity is also verified by comparing certain states of the cracking process with simulation analysis using MARC software.
The same number of polynomials and stress function terms are used in all the calculations.

A matrix-inclusion element
In this example, the VCFEM model consists of a simple Voronoi element, the MARC and ABAQUS models consist of a fine finite element mesh with a large number of elements.The validity of the VCFEM is verified by comparing it with the results of MARC and ABAQUS.The model to be simulated is an 2 × 2mm 2 square matrix containing a hollow inclusion with a radius of 0.2mm and the hollow part of the inclusion with a radius of 0.15mm .The inclusion is located at the center of the square matrix.The inclusion is represented by an octagonal structure in which the right two sides are cracked.For comparative analysis, the VCFEM model, MARC model and ABAQUS model are subjected to exactly the same boundary conditions and displacement loads are shown in Fig. 3a, the vertical displacement on the top boundary is 0, the vertical displacement on the bottom boundary is 0, and the horizontal displacement on the left boundary is 0, and the horizontal displacement on the right boundary is 0.0002 mm.
The MARC model consists of 9856 second-order plane stress triangular elements, with 5107 nodes as shown in Fig. 3c; the ABAQUS model consists of 3337 CPS4R elements and 100 CPS3 elements, with 3509 nodes; the VCFEM model consists of a single Voronoi element as shown in Fig. 3b, with stress functions in both the matrix and Voronoi elements containing 99 terms, composed of 63 polynomial stress functions and 33 interaction stress function terms.The material properties used in the modeling analysis are as follows: Inclusion:E c = 80000 MPa, υ c = 0.25, r c = 0.2 mm, r h = 0.15 mm Matrix:E m = 40000 MPa, υ m = 0.2 The critical normal stress at the interface is 1.05 N and the critical load is 0.0002 m.Comparative analysis of the results in Fig. 4 shows that the three models can be in good agreement with each other.From Figs. 5 and 6 at the centerline position of VCFEM, MARC, and ABAQUS σ x ,σ y comparison, it can be seen that there are different peaks at the matrix-inclusion interface and at the crack, and there are different degrees of stress concentration in the hollow part of the inclusion, and the stress decreases rapidly to the minimum value at the crack.The stress components of the three models are more or less different except at the interface of matrix and inclusion, which is mainly caused by these two reasons: (I) in VCFEM, the matrix part and the inclusion part within each element are assumed to be the stress function, so that describing the stress field with two continuous functions that satisfy the perfect continuity at all the interfaces of the matrix and the inclusion is difficult; (II) in the MARC model, there is only one node at the same location on the interface of matrix and inclusion, and there is only one value of the stress component at the interface node obtained by the software, which does not reflect the stress discontinuity applied at the interface.Whereas, in the VCFEM model and ABAQUS, there are two values of stress at the interface point.

Interfacial debonding of a matrix-inclusion element
In this example, the VCFEM model consists of a simple Voronoi element and the MARC and ABAQUS models consist of a fine finite element mesh with a large number of elements.The validity of the VCFEM is verified by comparing it with the results of MARC and ABAQUS.The model to be simulated is an 2 × 2 mm 2 square matrix containing a hollow inclusion with a radius of 0.2 mm and the hollow part of the inclusion with a radius of 0.15 mm .The inclusion is located in the center of the square substrate and inclusion interface consisted of 16 straight edge segments.For comparative analysis, the VCFEM model, MARC model and ABAQUS model are subjected to exactly the same boundary conditions and displacement loads, the vertical displacement on the top boundary is 0, the vertical displacement on the bottom boundary is 0, and the horizontal displacement on the left boundary is 0, and the horizontal displacement on the right boundary is 0.0002 mm.Initially, each side does not contain debonding, and as the load increases, debonding of different numbers of sides occurs when the stress state at each node of the inclusion reaches or exceeds the critical stress state.
In order to correspond to the debonding process, four MARC models were established for comparative analysis, including 13424 second-order isotropic plane stress triangle elements (containing 6970 nodes); an ABAQUS model consisting of 8316 CPS4R elements and 246 CPS3 elements (containing 8562 nodes); and a VCFEM model consisting of a Voronoi element.Both the matrix and Voronoi elements contain stress functions with 99 terms, composed of 63 polynomial stress function terms and 33 interaction stress function terms.The material properties of the model are the same as above.
The critical normal stress at the interface was 2.31N, and the critical displacement loads for debonding with different numbers of sides were 0.0002mm, 0.00045mm, and 0.00125mm, respectively.Figure 7 gives the stress cloud plots for MARC and VCFEM at complete bonding, 2 × 2 debonded sides, 2 × 4 debonded sides, and 2 × 6 debonded sides.Figure 8 gives the stress distribution of MARC,VCFEM,ABAQUS at the horizontal and vertical

Interfaical debonding of multiple hollow inclusions
In this calculation example, the main simulation is a complex hollow particle-reinforced composite, as shown in Fig. 9a, the target model has 20 randomly distributed elliptical hollow inclusions of arbitrary size and orientation on a square plate of 1 × 1 mm 2 .The 20 randomly distributed inclusions have a thickness of approxi- mately 0.0055 mm, and the range of variation in the semi-major axis length of the hollow elliptical part is about 0.013 mm to 0.022 mm.The material properties of the model are the same as above, fixing the horizontal Initially, all the inclusions interfaces were not debonding and each inclusion interface consisted of 16 straight edge segments.With the increase of load, the normal stresses at some interface nodes were larger than the given critical stresses, resulting in debonding of some interface segments, and the nodes with normal stresses larger www.nature.com/scientificreports/than the given critical stresses were replaced by a node pair after remeshing.The stress functions for both the matrix and the inclusions within the elements of the VCFEM model contain 99 terms The VCFEM model simulates the gradual damage generation of the entire structure under load as follows: a given displacement at the right end of the initial load is 0.0003 mm ; for each subsequent incremental step, the displacement is incremented by 0.00001 mm , and the total number of load steps is 50.Critical normal stress at the interface was 8.05N.
In order to check the validity of the VCFEM, for the damage case obtained from the incremental step VCFEM calculation in step 4 (debonding produced at the ends of elements 5 and 16), a MARC model was built for the same damage case, as shown in Fig. 9b, the MARC model consists of 30,027 elements, which contain 15,974 nodes Fig. 10a,b.Stress σ x cloud plots for MARC and VCFEM, respectively.Figures 11 and 12 show the stress distributions of the VCFEM model and MARC model on the horizontal and vertical centerlines.From these figures, it can be seen that the stress distributions of the two methods are in good agreement.The displacements of the internal nodes are related to the displacements of the external nodes of its own element, and the internal nodes are compressed within the element.Therefore, if no new cracks are created in the element, the calculation of the stiffness of the element is only needed in the first step of the whole calculation process.
The results of the process simulation of the structural damage by the VCFEM model are as follows: the horizontal stresses at step 1 (0.00031 mm), step 6 (0.00036 mm), step 11 (0.00041 mm) and step 16 (0.00046mm) are shown in Fig. 13.For step 1, the interface between the matrix and the hollow inclusions appears to be complete bonding.For step 6, six inclusions debonded at the interface.For step 11, fourteen inclusions experienced debonding at the interface.At step 16, all but one of the nineteen inclusions experienced debonding at the interface.In order to further understand the effect of interface debonding damage on the overall structure, the variation of the right end restraining force with the horizontal displacement of the right end is obtained, as shown in Fig. 14.In this case, the variation of the curve is roughly divided into three stages: (a) the initial linear stage, in which the displacement is less than 0.00033 mm, almost no debonding is generated and the restraining reaction force basically varies linearly with the displacement; (b) the massive damage stage, when the displacement of the right end is increased in the interval from 0.00033 mm to 0.00048 mm, there subtle fluctuations in the increase in confining reaction force due to the production of large amounts of debonding.(c) Linear increase stage, when the displacement of the right end is greater than 0.00048mm, there still exists a small amount of delamination generation.The impact on the overall structure is very small and the restraining reaction force and displacement become linear change relationship again.
Figure 15 shows the debonding of elements with different inclusions wall thickness at the same sixth incremental step, Fig. 15a six inclusions debonding occurred, Fig. 15b sixteen inclusions debonding occurred, Fig. 15c eighteen inclusions debonding occurred.It can be seen that under the same displacement load, the larger the wall thickness of the inclusions in the hollow inclusion element the more prone to debonding damage.

Conclusion
The new VCFEM was employed to model interface debonding of hollow particles reinforced composites.A remeshing strategy was developed to simulate the initiation and propagation of inclusion interface debonding.The efficacy and computational efficiency of this method were validated through comparisons with results obtained from established commercial software packages such as MARC and ABAQUS.Furthermore, this method was utilized to analyze the evolution of interface debonding damage in intricate composite structures.
The investigation also explored the impact of interface damage on horizontal displacement and constraint reaction forces, as well as the influence of hollow microsphere wall thickness on interface debonding.The findings

Figure 1 .
Figure 1.Voronoi cell model with a hollow inclusion, considering interfacial debonding.

Figure 2 .
Figure 2. Mesh remeshing steps for debonding at the matrix-inclusion interface.
https://doi.org/10.1038/s41598-024-57451-xwww.nature.com/scientificreports/centerline positions.As can be seen from the two figures, for each form of debonding, the stress plots computed by the two models are consistent, and the distribution of stress is different only at the interface between the inclusion and the matrix.The stress concentration in the hollow part is more pronounced in MARC and ABAQUS.

Figure 5 .
Figure 5. Distribution of stress (a) σ x , (b) σ y on the horizontal centerline in MARC, VCFEM and ABAQUS models.

Figure 6 .
Figure 6.Distribution of stress (a) σ x , (b) σ y on the vertical centerline in MARC, VCFEM and ABAQUS models.

Figure 8 .
Figure 8. Distribution of horizontal stress on the horizontal centerline and vertical centerline in MARC,VCFEM and ABAQUS models.

Figure 10 .
Figure 10.Comparison of horizontal stress cloud plots between (a) MARC model and (b) VCFEM model.

Figure 11 .
Figure 11.Distribution of stress (a) σ x , (b) σ y on the horizontal centerline in MARC and VCFEM.

Figure 12 .Figure 13 .
Figure 12.Distribution of stress (a) σ x , (b) σ y on the vertical centerline in MARC and VCFEM.

Figure 14 .
Figure 14.The variation of the right end constraint reaction force with the horizontal displacement of the right end during the process of damage generation.

Figure 15 .
Figure 15.Damage profile of (a) 0.0055mm thickness, (b) 0.00825mm thickness and (c) 0.011 mm thickness in the sixth incremental step.
A weak expression for the surface force boundary condition can be obtained: Scientific Reports | (2024) 14:6868 | https://doi.org/10.1038/s41598-024-57451-xwww.nature.com/scientificreports/ When the original crack tip A extends to point B, if the stress state at node B reaches or exceeds the critical state, then the debonding continues at point B, producing a new debonding layer, and the next step continues to judge whether the stress state at the other nodes reaches the critical state.