Prediction of Elastic Parameters of Particle Reinforced Composites Using Finite Element Simulations

The macroscopic properties of composite materials depend on the microscopic properties of the constituents and the geometric arrangement of their phases. Therefore, it is essential to predict heterogeneous materials’ mechanical properties by simulating microstructural finite element models. The present article aims to analyze particle reinforced composites composed of spherical alumina inclusions surrounded by a glass matrix using a tridimensional representative volume element. Herein, microstructures are artificially created considering a regular or random arrangement of the inclusions. Two materials systems previously studied in the literature were analyzed. The discretization of the models was performed to have periodic mesh, thus enabling the use of periodic boundary conditions. A finite element model is created using Abaqus software. Numerical results show that the macroscopic properties can be estimated with high accuracy for the temperature where linear matrix behavior stands. The predictions were compared to experimental data from the literature. The models with a regular arrangement of inclusions show a difference inferior to 10%, while random arrangements show a difference inferior to 3.9%. The developed numerical algorithms can be modified to include new features, such as other dispersed phase arrangements or nonlinear material behavior.


Introduction
Ceramic composites are used in specific applications due to certain exceptional properties as high-temperature stability, high hardness and good corrosion resistance 1,2 .A few examples are magnesium aluminate spinel (MgO/MgAl 2 O 4 ), magnesia-hercynite (MgO/FeAl 2 O 4 ) -both used to replace Chromium in magnesium base refractories -silicon carbides reinforced ceramic matrix composites (SiC/SiC, C/SiC, Al 2 O 3 /SiC) [2][3][4] .These materials can be designed to obtain an increased elastic modulus and improved mechanical properties, as strength and toughness, and during this phase it is of the most importance to consider the acting thermal loads all along the material processing and its application.This concern occurs because high temperature gradients or temperature variations can lead to the appearance of microcracks, a phenomenon associated with the different thermal expansion coefficients (CTE) between the dispersed and continuous phases.Therefore, understanding how the phases interacts with each other is essential to develop new materials with better elastic properties.
The influence of temperature variations in the mechanical strength of ceramic composites has been a topic of study for many researchers.Selsing 5 wrote an expression to evaluate the internal stress in biphasic composites with spherical inclusions.His model considers a single spherical inclusion surrounded by an infinite matrix and a perfect interface between the phases, so it is applicable for composites with a low volume fraction.Based on Selsing's model, Davidge and Green 6 studied the strength of two-phase ceramic materials and proposed an analytical model that established a critical inclusion size for the crack initiation process depending on the thermomechanical properties and the thermal loading.To better understand this phenomenon, Tessier-Doyen et al. 7 studied the influence of cracks caused by CTE mismatch in the Young's modulus during the manufacturing process of model materials.Joliff et al. 8 studied a system with two inclusions in order to understand the stress distribution between inclusions, concluding that the closer the inclusions, the higher is the interfacial stress.Luchini et al. 9 did a computational investigation over the matrix inter inclusion cracks in composite systems, observing that for high volume fractions (ϕ > 35%), the mid-distance between two inclusions is the most stressed point.
According to the works previously mentioned, the coefficient of thermal expansion mismatch between the *e-mail: fernando.schiavon@ifsp.edu.brinclusion and matrix phase leads essentially to two kinds of cracks, viz.radial and orthoradial.The occurrence of one or another depends on the relation ∆α ∆T, where ∆α is the difference between the inclusion and matrix CTE and ∆T is the temperature variation.For ∆α ∆T < 0, the tension would appear in a radial direction and compression in an orthoradial direction resulting in debonding or orthoradial cracks.For ∆α ∆T > 0, compression would occur in a radial direction and tension in an orthoradial direction, thus creating cracks in a radial direction.This last scenario is considered the most critical case, because while circumferential cracks end when the phases are debonded, radial cracks can growth as long as the system provides energy for its propagation.The two types of cracks can be seen in Figure 1.
As can be seen, the design of ceramic composites is complex and demands a great deal of analysis and studies.This process can be assisted by numerical simulations, which reduce the experimental working and, when associated with experimental results, enables a better comprehension of deformation mechanisms.However, phenomena such as debonding and crack propagation are very non-linear, making their modeling more complex.Therefore, many researchers dedicate themselves to simulating the damage [10][11][12] , phase debonding 13,14 and plasticity 15,16 .In particular, some of them [12][13][14] take into account the nonlinear phenomena in the homogenization process.The present study is limited to the linear analysis of RVEs considering different geometric arrangements in a two-phase system.
Predictions via numerical methods have become a powerful tool to analyze the material properties of heterogeneous materials 9,17 .According to Charalambakis 18 , through these simulations it is possible to select and identify the appropriate materials to achieve the desired thermomechanical properties.The numerical analysis of these composites requires the creation of a Representative Volume Element (RVE), i.e. a small portion of the material providing results that are in agreement with the entire media's behavior.By averaging the stress and strains in the RVE volume, the macroscopic stress and strains can be calculated.These procedures to obtain effective macroscopic properties based on the results of a RVE are called homogenization methods 19,20 .Although computational power has increased a lot in the last decades, it is convenient to simplify the RVE geometry when working with numerical simulations.
In the context of the prediction of the effective mechanical properties using numerical simulations, the present paper aims to estimate elastic properties of biphasic model composites using regular and random periodic RVEs.The analysis are conducted at several temperatures and the simulations results are compared with experimental data.The RVE creation and the application of Periodic Boundary Conditions (PBC) are described and the results are presented.

Materials
The deformation mechanisms of ceramic composites are complex because of the volume fraction, the number and shape of phases, and the interfacial properties between the dispersed phase(s) and continuous one.Some authors developed model materials to highlight phenomena of interest and better comprehend deformation mechanisms [6][7][8][21][22][23] . Regardng the temperature-dependent behavior of the Young modulus, Tessier-Doyen et al. 7 developed model materials to emphasize the influence of the CTE mismatch.For that, the authors developed a ceramic composite material composed of spherical Alumina (Al 2 O 3 ) inclusions immersed in boresilicate glasses matrix with volume fractions of 15%, 30%, and 45%.Three different types of glass materials were mixed with the dispersed alumina phase to generate three typical microstructure configurations, Alumina/G1 (∆α > 0), Alumina/G2 (∆α ≈ 0) and Alumina / G3 (∆α < 0).The same materials used for Tessier-Doyen will be simulated in the present paper.Material properties are shown in Table 1.
The material systems will be nominated as GgAlϕ where g  denotes the glass matrix, which can be 1, 2, or 3, and ϕ is the volume fraction of inclusions.
Figure 2 shows the data from Tessier-Doyen et al. 7 for Alumina/G2 composites.The mechanical properties were obtained on a cylindrical sample of the material pressed and subsequently sintered at a high temperature.The sample is then subjected to heat treatment, while ultrasonic methods measure the Young modulus.The blue curve shows the behavior of the Young modulus for the pure G2 glass matrix during heating and cooling.The green curve represents the G2Al15.The heating cycle starts at A, where the Young modulus is almost constant.The temperature rises to 800°C.The Young modulus starts to fall because the glass transition temperature has been reached (region B).During the cooling cycle, the material Young modulus increases because of the matrix solidification and remains almost constant up to room temperature.G2Al30 and G2Al45, representing the materials with 30 and 45% of alumina volume fraction, respectively, have almost the same behavior, with the only difference in the mean value of the resistance of the material been higher for higher volume fractions.This behavior is expected once the CTE difference between alumina and G2 is very small, producing little to no damage in the Alumina/G2 composites by thermal loads.
The curves for the composites Alumina/G1 and Alumina/ G3 are shown in Figure 3, which presents an untypical Young modulus behavior when compared to the pure glassy material.Pure materials G1 and G3 have the same behavior as material G2 with the difference in the glass transition temperature.In the green curve G1Al15, six regions are highlighted (A, B, C, D, E, and F) and they refer to different material stages.In the beginning of the heating (region A) no significant change in the Young modulus occurs.In region B, the thermal dilatation starts to close orthoradial cracks and an increment of the modulus of elasticity can be seen, followed by a decreasing in region C, where the glass transition temperature is attained.At the beginning of the cooling cycle (point D), the resistance grows as the matrix starts to solidify, and its values are higher than those in region C because the interface is re-established.During the cooling process the Young modulus keeps increasing up to a maximal value at region E after which cracks and debonding reappear, causing it to decrease.This same behavior can be seen in materials G1Al30 and G1Al45 with higher maximal modulus values for higher volume fractions.Another phenomenon worth of mention is the rise in the difference between the Young modulus values during heating and cooling with the increasing of volume fraction.This can be explained by a larger amount of Alumina in higher volume fractions, which leads to higher internal stress, making the material more susceptible to cracking.
Material Alumina/G3 macroscopic behavior has some similarities with Alumina/G1 system but has different deformation mechanisms because of the different CTE mismatch.Opposite of what was seen for Alumina/G2, the Young modulus increases during the heating phase (region B) during a small temperature range.This occurs because the cracks are primarily radial, creating a new interface in the composite and therefore restoring the adhesion between the cracked surfaces.Then, the Young modulus stays in a plateau (region C) until the end of the heating cycle.The cooling cycle starts in region D with an increasing in Young modulus up to region E, where there is an abrupt decrease (region F).This decrease in the Young modulus is due to the reappearance of the cracks, and it is more accentuated than in material Alumina/G1 because in this case the cracks are radial.Also, the Young modulus difference between the heating and cooling cycles is higher than in Alumina/G1.It is interesting to highlight that the pure matrix G3 has a higher Young modulus at low temperatures than their composites.
In the present study, only the combinations that present a mismatch of the CTE are investigated, so Alumina/G1 and Alumina/G3 materials systems.Besides, the material is only analyzed in the cooling stage.The composite system was considered to be undamaged and with no residual stresses at the glass transition temperature.As the composite system cools, the possibility of failure increases due to rising stress levels.At a certain level of temperature variation, the damaging process due to radial cracks and debonding between inclusion and matrix leads to a reduction of the macroscopic Young modulus.The only region of the experimental curve where the material can be considered to have an almost elastic behavior is from the glass transition temperature up to the decrease of the Young modulus.

Numerical model
An overview of the steps concerning the creation and analysis of the numerical model is shown in Figure 4.The process of obtaining the effective properties consists on the generation of artificial microstructures, their analysis via computer simulations, and evaluation of the macroscopic values of mechanical properties from the stress and strain fields.The first step comprises the generation of a periodic microstructure with regular or random dispersed phase distribution.Then, meshes with same geometry are generated on opposite sides of the cube's domain.The final step to constructing the RVE is applying constraints on boundary opposite nodes to enforce a periodicity in the stress and strain fields.Further details on each step of the developing process are presented in the following sections.

RVE geometry
The RVE consists of a cube with spherical inclusions of radii r.A cube with edge L positioned at the coordinate system origin represents the boundary of the RVE domain, that is, The matrix complements the spherical inclusions.The ratio r/L defines the number of inclusions necessary to fulfill the volume fraction of the inclusion phase.According to Bargmann et al. 24 , for computational efficiency purposes, it is sought the smallest RVE volume that allows computation of effective properties such as elastic stiffness without reverting to a modeling attempt of the whole material domain.
Herein, two geometric types of inclusion distribution are considered: regular and random.The regular distribution assumes that the inclusion positions are predefined in three possible arrangements: C -simple cubic structure, BCC -body centered cubic structure, and FCC -face centered cubic structure, as shown in Figure 5. RVEs with regular geometric  arrangements were modeled with volume fractions of 15, 30, and 45%.For the RVE with random distribution, only fractions of 15 and 30% could be modeled.The 45% volume fraction could not be achieved due to a limitation in the algorithm employed.As shown by Cooper 25 , the maximum fraction enabled by the random adsortion method with particles of the same radius is 0.385 ± 0.010.In his study, Cooper performed thousands of runs with ratios L/r going up to 40.In the random arrangements, a ratio L/r = 10 was adopted.
For the random microstructure, the Random Sequential Adsorption method is used 26 .It inserts the inclusions into a random position in the domain without to collide with the existing particles.The collision detection check is made between the trial inclusion and the ones already settled in the RVE by measuring the distance between them: where try stands for the trial position, and set for the inclusions already inserted in the RVE.In case none of them collide, that is, if for all i particles d > 2r, the inclusions belonging to be part of the RVE.After that, a check is performed to verify how many faces the inserted inclusion intersects.Four cases are possible: I -three faces, II -two faces, III -one face and IV -no face intersection (inner).The number n of intersections is defined as: where ( ) x  is null for x ≤ 0 and the unity otherwise.Depending upon the n value, one, two, four or eight particles will be stored on a temporary variable as shown in Figure 6 according to the drawn particle position (not shown for inner particle).It is important to mention that the sum of all the particles volumes that are inside the cube and are stored in the auxiliary variable is precisely equal to one sphere volume.Figure 7 depicts the whole process of generating a periodic RVE of a random microstructure.A MATLAB script was developed to define the random position of the inclusions.

RVE finite element meshing
Once the inclusions position have been established, the next step is to discretize the domain into finite elements.Gmsh was used to generate the mesh to be posterly imported in Abaqus to perform the finite element analysis.Gmsh enables meshing the domain verifying the periodicity, and also it has a Python API (Application Programming Interface) that allows a relatively simple connection with Abaqus.
Gmsh is an open-source 3D finite element mesh generator with two CAD engine kernels (Built-in and OpenCascade) that enables the creation of complex mesh geometries with several types of elements 27 .The Built-in kernel has a bottom-up approach, while the second offers a constructive solid geometry approach.Therefore, considering the desired models, the OpenCascade (OCC) kernel is the most adequate to create the geometries.A particle reinforced composite can be generated using simple boolean operations involving a cube and spheres (representing the dispersed phase) as can be seen in Figure 8.
Regarding mesh generation, Gmsh offers a feature that allows creating periodic meshes by establishing an affine transformation matrix.This matrix creates a geometric transformation between two entities, and it preserves points, straight lines, and planes.Particularly, for parallelepipedal RVEs with edges parallel to the coordinate system, the translation transform is enough to apply the periodicity along with the three directions.Figure 8 illustrates the periodicity of the mesh between opposite sides of a cubic RVE.The meshing process does not create an interface region between inclusion and matrix.The model considers that nodes in the inclusion boundaries are shared with matrix and inclusion elements.

Periodic boundary constraints
After the RVE creation and meshing, the mesh is imported in Abaqus finite element software, where the analyses are performed.Since the model is imported as an orphan mesh, the attribution of materials, sections, and boundary conditions are performed using mesh entities.The application of periodic boundary conditions is split into three parts: (i) creation of sets, (ii) identification of node pairs, and (iii) creation of nodal constraints.The imposition of nodal constraints is made, as suggested by Tian et al. 28 , in terms of displacement differences involving nodes on opposite faces of the RVE as: ( ) where u i denotes the displacement degree of freedom along the direction i, ij ε is the macroscopic strain prescribed, and the subscripts plus and minus denote the RVE face with negative and positive normal, respectively.In practice, the components of the macroscopic strain tensor are applied using Reference Points (RPs) in Abaqus software.The identification and constraint creation were automated using Python routines and Abaqus Scripting Interface (ASI).

Homogenization method
In the last part, the macroscopic strains are applied on the RVE, and the effective material parameters are identified after the integration of stress and strain tensors over the RVE volume.According to the homogenization method 29 , the macroscopic stress and strain can be computed as: where ε and σ are the strain and stress tensors, respectively.
The integration was performed using the stress and strain components evaluated at each element integration point.From the macroscopic components obtained using Equations 4 and 5,  the effective elastic parameters can be identified submitting the RVE to different loading cases.

Identification of the elastic constants
The constitutive tensor (C) is identified by a minimization process between the obtained stress state and the predicted one with a set of values p.The process consists of minimizing the functional Π defined as: where n is the number of numerical analysis performed and σ is obtained using Equation 4and σ(p) is written as: ( ) ( ): The tensor C is written in terms of the vector of parameters p.The number of unknowns in p depends on the type of materials system considered.For instance, herein the material is considered to be isotropic, thus the constitutive tensor is defined using two parameters.The vector of parameters that minimizes Π, p ⋆ is the one used to determine the constitutive tensor, i.e.: ( ) In order to obtain the elastic properties of the RVE, six mechanical loading cases are applied on the RVE, with three under pure axial deformation and other three under pure shear deformation: (1) ε xx = 0.001, ε yy = 0, ε zz = 0, (2) ε xx = 0, ε yy = 0.001, ε zz = 0, (3) ε xx = 0, ε yy = 0, ε zz = 0.001, (4) ε xy = 0.001, ε xz = 0, ε yz = 0, (5) ε xy = 0, ε xz = 0.001, ε yz = 0, (6) ε xy = 0, ε xz = 0, ε yz = 0.001.The numerical analysis performed at each temperature consists of linear elastic computation, so it requires only a single step to obtain the fields of interest.
The generated RVEs were analysed using the finite element software Abaqus for the cooling stage of Tessier's experiment.All samples were analysed for several temperatures between the glass transition, viz.600 ºC for Alumina/G1 and 500 ºC for Alumina/G3, and room temperature (25 ºC).Below the glass transition temperature regarding the cooling phase of the heat treatment, the matrix behavior is governed by a linear elastic constitutive equation.Thus, the glass transition temperature is the maximum temperature simulated.Also, it is important to remark that the interface between inclusion and the surrounding matrix is perfect.There is no contact model in the interface, and it is modeled in such a way that neighbor elements share the nodes that belong to the interface.Inclusion and matrix materials are considered elastic.The analyses were carried out using the Abaqus Standard solver.

RVEs with a regular arrangement of inclusions
The results for the case with a regular arrangement for the material Alumina/G1 are depicted in Figure 9.It is possible to see that BCC (Figure 9a) and FCC (Figure 9b) have very similar values, while C structure presented higher values.The RVE with 30% of volume fraction is the one that presented values nearer experimental values.Regarding the RVE structure type, the C structure has shown high discrepancy.This divergence can be associated with the distance between the particles, since C structure has the smallest distance with 0.049L for 45% while for BCC is 0.111L and 0.108L for FCC.The model predicts the values quite well in the range of temperature 300 o C up to 600 o C with error below 10% for the fractions of 15 and 30%.The model does not predict values with accuracy for higher fractions as 45%, with values around 10% of error for structures type BCC and FCC and around 15% for the structure type C.
The results for the case with material Alumina/G3 are depicted in Figure 10.For this material, the lower the volume fraction, the lower is the error.Again, the structure type C with fraction of 45% is the one that presented the worst predictions for the Young modulus.Between the temperatures of 200 o C and 500 o C and for the fraction of 15% the higher error is about 7% for the structures BCC and FCC.C structure presents a maximum error of 10%.For the fraction of 45% the errors arise for 15% for structures BCC and FCC and up to 20% for the C structure type.
The composite material elastic properties can be defined by three variables, the Poisson coefficient, the Young modulus and the shear modulus.The obtained Young modulus for the regular cases is plotted in Figures 9 and 10.The obtained Poisson coefficient is described in Table 2.The Poisson ratio does not change significantly with the temperature, because of that, only the average values are shown in Table 2.As the variations for the Young modulus with temperature are below 0.7% (for instance, 0.001 variation in a value of 0.161 for the case 15%, BCC structure and Alumina/G1 material), the Young modulus can be considered constant in the temperature range from 25 o C to 600 o C for Alumina/G1, and 25 o C to 500 o C for Alumina/G3.

RVEs with a random arrangement of inclusions
The RVE with a random geometry particles distribution differs from the regular geometry by the fact that the particles are not in the same position and also that the RVE contains much more particles.The properties are calculated on an average taken on five different samples.As periodic boundary conditions are imposed over the boundaries of the domain, there is continuity of displacement and stresses between the left/right, rear/frontal and upper/bottom sides.Figure 11 illustrates this with nine RVEs put side by side on horizontal and vertical directions (z direction not shown) and with the original RVE in the center.Images of the RVE under applied normal deformation in the directions xx, yy and zz are depicted in Figure 12.Pure shear stresses fields due to shear deformations xy, xz and yz are depicted in Figure 13.
The results obtained for the random distribution can be seen in Figure 14 that also shows the results obtained by experimental means.For the Alumina/G1 material, the model predicts the values quite well above the temperature of 300°C up to 600°C.For the Alumina/G3 the temperature range for a good prediction starts from 200°C up to 500°C.Below 200°C for the Alumina/G1 and 300°C for the Alumina/ G3 material, the CTE mismatch can induce the failure of the matrix/inclusion interface or cracks throughout the matrix, reducing the Young modulus considerably.This phenomenon becomes evident for the lower temperatures, corresponding to more significant temperature variations.As the developed numerical model has no matrix failure model implemented, there was a significant difference in observed values for lower temperatures.It is shown in Table 3 the standard deviation of the Young modulus obtained by virtual tests measured over the five samples taken from each volume fraction batch.As can be seen, the more significant value is 154 MPa, that compared with measured values of E, gives a variation of about 0.15%.The differences between the numerical values obtained in this research and the experimental values obtained by Tessier-Doyen et al. 7 are described in Table 4.The lowest differences considering regions with no expected cracks were obtained for Alumina/G1 30% and Alumina/G3 15% materials.Higher temperatures have increasing differences because matrix material starts to achieve the glass transition temperature.For the Alumina/G1 composite material in the range of 200 up to 600°C, the maximum difference is about 7.3%.For the Alumina/G3 in the range of 200 up to 500°C, the maximum difference is 3.9%.

Conclusions
The present paper evaluated the temperature-dependent macroscopic elastic properties of particle-reinforced composites using computer simulations.The simulations were conducted using the finite element method in periodic artificial two types of microstructures, viz.regular and random inclusion arrangements.The predictions have shown a good agreement with the experimental measurements, especially for the random inclusion distribution.For volume fractions of 15% and 30%, body centered (BCC) and face centered (FCC) arrangements have shown a difference inferior to 10%.The simple cube arrangement (C) presented the worst results.For the same volume fractions, the random inclusion arrangement presented a difference inferior to 3.9%.
The best-fitted region corresponds to the region where no damage occurs in the material, i.e., before the Young modulus decreases.A nonlinear interface and matrix behavior must be incorporated into the numerical model to capture this damaging process, such as proposed in the literature [10][11][12][13][14] .The identification of model parameters of further models is complex and, generally, involves indirect measurements.The simulation of the temperature dependent Young modulus in the entire cycle remains a defiance since it must combine damaging and healing models.Finally, the developed numerical algorithms can be modified to include new features as other disperse phase arrangements or nonlinear material behavior.

Figure 1 .
Figure 1.Influence of the CTE mismatch and temperature variation in the expected failure modes for a biphasic composite.The subindexes i and m are referent to the inclusion and matrix respectively.Adapted from Tessier-Doyen 7 .

Figure 4 .
Figure 4. Overview of the computational tool from the RVE creation up to the estimative of elastic properties.

Figure 5 .
Figure 5. Layout of regular inclusion arrangements in the RVE.

Figure 6 .
Figure 6.Replication of the inclusions depending on the number of intercepted faces.

Figure 7 .
Figure 7. Flow chart describing the process to insert the inclusions into the cubic domain.

Figure 8 .
Figure 8. Periodic meshes between opposite sides of a cubic RVE.

Figure 11 .
Figure 11.Periodic condition with continuity of stress field on left/ right and bottom/top sides.

Figure 12 .
Figure 12.Principal stress fields on RVE under different tension loading states in the deformed configuration.

Figure 13 .
Figure 13.Principal stress fields on RVE under different shearing loading states in the deformed configuration.

Figure 14 .
Figure 14.Young modulus obtained through computer simulation with tridimensional random RVE for the materials Alumina/G1 and Alumina/G3.

Table 3 .
Standard deviation for the calculated Young modulus evaluated over five samples.

Table 4 .
Percentual difference between the predicted Young modulus with respect to the experimental values obtained by Tessier-Doyen 7 .