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

Abstract

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.

Keywords:
Thermo-elastic properties; representative volume element; finite element analysis; composites

1. Introduction

Ceramic composites are used in specific applications due to certain exceptional properties as high-temperature stability, high hardness and good corrosion resistance¹,². A few examples are magnesium aluminate spinel (MgO/MgAl₂O₄), magnesia-hercynite (MgO/FeAl₂O₄) - both used to replace Chromium in magnesium base refractories - silicon carbides reinforced ceramic matrix composites (SiC/SiC, C/SiC, Al₂O₃/SiC)²-⁴. 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⁵ 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⁶ 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.⁷ studied the influence of cracks caused by CTE mismatch in the Young’s modulus during the manufacturing process of model materials. Joliff et al.⁸ 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.⁹ 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 inclusion 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.

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

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¹⁰-¹², phase debonding¹³,¹⁴ and plasticity¹⁵,¹⁶. In particular, some of them¹²-¹⁴ 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⁹,¹⁷. According to Charalambakis¹⁸, 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¹⁹,²⁰. 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.

2. 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⁶-⁸,²¹-²³. Regarding the temperature-dependent behavior of the Young modulus, Tessier-Doyen et al.⁷ 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₂O₃) inclusions immersed in bore-silicate 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.

Table 1
Material properties (extracted from Tessier-Doyen⁷).

Material	Properties	Temperature (oC)
		25	100	200	300	400	500	600
Alumina	α (10−6 K−1)	4.8	5.3	6.2	6.8	6.8	7.1	7.4
	E (GPa)	340	335	330	325	320	315	310
	Poisson (ν)	0.24
G1	α (10−6 K−1)	4.33
	E (GPa)	68
	Poisson (ν)	0.20
G3	α (10−6 K−1)	5.7	9.7	10.7	10.9	11.3	13.7	-
	E (GPa)	72	71	70	68	64	57	-
	Poisson (ν)	0.23

Figure 2 shows the data from Tessier-Doyen et al.⁷ 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.

Figure 2
Temperature-dependent Young modulus of Alumina/G2. Experimental data obtained by Tessier-Doyen⁷.

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.

Figure 3
Temperature-dependent Young modulus of Alumina/G1 and Alumina/G3. Experimental data obtained by Tessier-Doyen⁷.

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.

3. 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.

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

3.1. 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, Ω_RVE ≡ Ω = {(x, y, z) | x ∈ [0, L], y ∈ [0, L], z ∈ [0, L]}. 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.²⁴, 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²⁵, 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.

Figure 5
Layout of regular inclusion arrangements in the RVE.

For the random microstructure, the Random Sequential Adsorption method is used²⁶. 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:

$d=\sqrt{{\left(x_i^{try}-x_j^{set}\right)}^2+{\left(y_i^{try}-y_j^{set}\right)}^2+{\left(z_i^{try}-z_i^{set}\right)}^2}$ (1)

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:

$$n={\displaystyle \sum }_{i=1}^3[\mathcal{H}\left(x_i+r\right)-\mathcal{H}\left(x_i-r\right)+\mathcal{H}\left(x_i+r-L\right)-\mathcal{H}\left(x_i-r-L\right)]$$ (2)

where $\mathcal{H}\left(x\right)$ 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.

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

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

3.2. 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²⁷. 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.

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

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.

3.3. 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.²⁸, in terms of displacement differences involving nodes on opposite faces of the RVE as:

$$u_i^{+}-u_i^{-}={\displaystyle \sum }_{j=1}^3{\overline{\epsilon }}_{ij}\left(x_j^{+}-x_j^{-}\right)$$ (3)

where u_i denotes the displacement degree of freedom along the direction i, ${\overline{\epsilon }}_{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).

3.4. Effective mechanical properties

3.4.1. 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²⁹, the macroscopic stress and strain can be computed as:

$\overline{\sigma }=\frac{1}{V_{VRE}}{{\displaystyle \int }}_{V_{VRE}}^{}\sigma dV$ (4)

$\overline{\epsilon }=\frac{1}{V_{VRE}}{{\displaystyle \int }}_{V_{VRE}}^{}\epsilon dV$ (5)

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.

3.4.2. 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:

$$\Pi \left(p\right)={\displaystyle \sum }_{i=1}^n{‖\overline{\sigma }-\sigma \left(p\right)‖}^2$$ (6)

where n is the number of numerical analysis performed and $\overline{\sigma }$ is obtained using Equation 4 and σ(p) is written as:

$\sigma \left(p\right)=C\left(p\right):\overline{\epsilon }$ (7)

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.:

$p^{\star }=\arg mi{n}_p\Pi \left(p\right)$ (8)

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.

4. Results

4.1. 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^oC up to 600^oC 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.

Figure 9
Young modulus obtained for Alumina/G1 composites.

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^oC and 500^oC 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.

Figure 10
Young modulus obtained for Alumina/G3 composites.

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^oC to 600^oC for Alumina/G1, and 25^oC to 500^oC for Alumina/G3.

Table 2
Poisson coefficient values for Alumina/G1 and Alumina/G3 using a regular arrangement of inclusions.

Volume fraction	MATERIAL
	Alumina/G1			Alumina/G3
	BCC	FCC	C	BCC	FCC	C
15%	0.161	0.161	0.158	0.183	0.183	0.179
30%	0.165	0.165	0.154	0.184	0.183	0.173
45%	0.169	0.169	0.153	0.184	0.184	0.168

4.2. 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.

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

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

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

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%.

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.

Material	Alumina/G1		Alumina/G3
Volume fraction	15%	30%	15%	30%
E (Young modulus) standard deviation (Mpa)	29	133	44	154

The differences between the numerical values obtained in this research and the experimental values obtained by Tessier-Doyen et al.⁷ 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%.

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

Material	Volume Fraction	Temperature (oC)
		100	200	300	400	500	600
Alumina/G1	15%	-	7.3	2.8	2.6	2.9	4.6
	30%	5.5	2.0	0.2	0.1	0.4	1.8
Alumina/G3	15%	37.4	-0.4	-0.4	-0.9	-1.3	-
	30%	-	3.9	2.8	2.3	0.7	-

5. 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¹⁰-¹⁴. 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.

6. Acknowledgements

The authors would like to thank the financial support from the CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico), process 141041/2019-6, 140250/2020-4 and IFSP - Instituto Federal de São Paulo.

7. References

• ¹
  Ashby MF. Materials selection in mechanical design. Amsterdam: Elsevier; 2005. Chapter 11, Designing hybrid materials. https://doi.org/10.1016/B978-1-85617-663-7.00011-4
  » https://doi.org/10.1016/B978-1-85617-663-7.00011-4
• ²
  Chawla KK. Composite materials. Birmingham: Department of Materials Science and Engineering, University of Alabama; 2012. http://dx.doi.org/10.1007/978-0-387-74365-3
  » http://dx.doi.org/10.1007/978-0-387-74365-3
• ³
  Khlifi I, Pop O, Dupré JC, Doumalin P, Huger M. Investigation of microstructure-property relantionships of magnesia-hercynite refractory composites by a refined digital image correlation technique. J Eur Ceram Soc. 2019;39(13):3893-902. http://dx.doi.org/10.1016/j.jeurceramsoc.2019.05.010
  » http://dx.doi.org/10.1016/j.jeurceramsoc.2019.05.010
• ⁴
  Aksel C, Rand B, Riley FL, Warren PD. Thermal shock behaviour of magnesia-spinel composites. J Eur Ceram Soc. 2004;24(9):2839-45. http://dx.doi.org/10.1016/j.jeurceramsoc.2003.07.017
  » http://dx.doi.org/10.1016/j.jeurceramsoc.2003.07.017
• ⁵
  Selsing J. Internal stresses in ceramics. J Am Ceram Soc. 1961;44(8):419. http://dx.doi.org/10.1111/j.1151-2916.1961.tb15475.x
  » http://dx.doi.org/10.1111/j.1151-2916.1961.tb15475.x
• ⁶
  Davidge RW, Green TD. The strength of two-phase ceramic/glass materials. J Mater Sci. 1968;3(6):629-34. http://dx.doi.org/10.1007/BF00757910
  » http://dx.doi.org/10.1007/BF00757910
• ⁷
  Tessier-Doyen N, Glandus JC, Huger M. Untypical Young’s modulus evolution of model refractories at high temperature. J Eur Ceram Soc. 2006;26(3):289-95. http://dx.doi.org/10.1016/j.jeurceramsoc.2004.10.028
  » http://dx.doi.org/10.1016/j.jeurceramsoc.2004.10.028
• ⁸
  Joliff Y, Absi J, Huger M, Glandus JC. Experimental and numerical study of the elastic modulus vs temperature of debonded model materials. Comput Mater Sci. 2008;44(2):826-31. http://dx.doi.org/10.1016/j.commatsci.2008.04.024
  » http://dx.doi.org/10.1016/j.commatsci.2008.04.024
• ⁹
  Luchini B, Sciuti VF, Angélico RA, Canto RB, Pandolfelli VC. Thermal expansion mismatch interinclusion cracking in ceramic systems. Ceram Int. 2016;42(10):12512-5. http://dx.doi.org/10.1016/j.ceramint.2016.05.013
  » http://dx.doi.org/10.1016/j.ceramint.2016.05.013
• ¹⁰
  He W, Wu YF, Xu Y, Fu TT. A thermodynamically consistent nonlocal damage model for concrete materials with unilateral effects. Comput Methods Appl Mech Eng. 2015;297:371-91. http://dx.doi.org/10.1016/j.cma.2015.09.010
  » http://dx.doi.org/10.1016/j.cma.2015.09.010
• ¹¹
  Pereira L, Weerheijm J, Sluys L. Simulation of dynamic behavior of quasi-brittle materials with new rate dependent damage model. In: 9th International Conference on Fracture Mechanics of Concrete and Concrete Structures (FraMCos-9); 2015; Berkeley, USA. Proceedings. Champs-sur-Marne: RILEM; 2015. p. 14. http://dx.doi.org/10.21012/FC9.036
  » http://dx.doi.org/10.21012/FC9.036
• ¹²
  Bheemreddy V, Chandrashekhara K, Dharani LR, Hilmas G. Computational study of micromechanical damage behavior in continuous fiber-reinforced ceramic composites. J Mater Sci. 2016;51(18):8610-24. http://dx.doi.org/10.1007/s10853-016-0120-4
  » http://dx.doi.org/10.1007/s10853-016-0120-4
• ¹³
  Santos WFd, Pituba JJdC. Yield surfaces of material composed of porous and heterogeneous microstructures considering phase debonding. Lat Am J Solids Struct. 2017;14(8):1387-415. http://dx.doi.org/10.1590/1679-78253776
  » http://dx.doi.org/10.1590/1679-78253776
• ¹⁴
  Azizi R. Micromechanical modeling of damage in periodic composites using strain gradient plasticity. Eng Fract Mech. 2012;92:101-13. http://dx.doi.org/10.1016/j.engfracmech.2012.04.033
  » http://dx.doi.org/10.1016/j.engfracmech.2012.04.033
• ¹⁵
  Grassl P, Jirasek M. Plastic model with non-local damage applied to concrete. Int J Numer Anal Methods Geomech. 2006;30(1):71-90. http://dx.doi.org/10.1002/nag.479
  » http://dx.doi.org/10.1002/nag.479
• ¹⁶
  Xue L, Wierzbicki T. Ductile fracture initiation and propagation modeling using damage plasticity theory. Eng Fract Mech. 2008;75(11):3276-93. http://dx.doi.org/10.1016/j.engfracmech.2007.08.012
  » http://dx.doi.org/10.1016/j.engfracmech.2007.08.012
• ¹⁷
  Kuna M. Finite elements in fracture mechanics: theory - numerics - applications. Dordrecht: Springer, 2013. (Solid Mechanics and Its Applications; 201). http://dx.doi.org/10.1007/978-94-007-6680-8
  » http://dx.doi.org/10.1007/978-94-007-6680-8
• ¹⁸
  Charalambakis N. Homogenization techniques and micromechanics: a survey and perspectives. Appl Mech Rev. 2010;63(3):030803. http://dx.doi.org/10.1115/1.4001911
  » http://dx.doi.org/10.1115/1.4001911
• ¹⁹
  Bensoussan A, Lions JL, Papanicolaou G. Asymptotic analysis for periodic structures. Providence, RI: American Mathematical Society; 2011. (vol. 374). http://dx.doi.org/10.1090/chel/374
  » http://dx.doi.org/10.1090/chel/374
• ²⁰
  Sun CT, Vaidya RS. Prediction of composite properties from a representative volume element. Compos Sci Technol. 1996;56(2):171-9. http://dx.doi.org/10.1016/0266-3538(95)00141-7
  » http://dx.doi.org/10.1016/0266-3538(95)00141-7
• ²¹
  André D, Levraut B, Tessier-Doyen N, Huger M. A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expansion mismatch of two-phase materials. Comput Methods Appl Mech Eng. 2017;318:898-916. http://dx.doi.org/10.1016/j.cma.2017.01.029
  » http://dx.doi.org/10.1016/j.cma.2017.01.029
• ²²
  Bobet JL, Naslain R, Guette A, Ji N, Lebrun JL. Thermal residual stresses in ceramic matrix composites-II. Experimental results for model materials. Acta Metall Mater. 1995;43(6):2255-68. http://dx.doi.org/10.1016/0956-7151(94)00430-7
  » http://dx.doi.org/10.1016/0956-7151(94)00430-7
• ²³
  Stoneham AM, Harding JH. Invited review: mesoscopic modelling: materials at the appropriate scale. Mater Sci Technol. 2009;25(4):460-5. http://dx.doi.org/10.1179/174328408X311062
  » http://dx.doi.org/10.1179/174328408X311062
• ²⁴
  Bargmann S, Klusemann B, Markmann J, Schnabel JE, Schneider K, Soyarslan C, et al. Generation of 3D representative volume elements for heterogeneous materials: a review. Prog Mater Sci. 2018;96:322-84. http://dx.doi.org/10.1016/j.pmatsci.2018.02.003
  » http://dx.doi.org/10.1016/j.pmatsci.2018.02.003
• ²⁵
  Cooper DW. Random-sequential-packing simulations in three dimensions for spheres. Phys Rev A Gen Phys. 1988;38(1):522-4. http://dx.doi.org/10.1103/PhysRevA.38.522
  » http://dx.doi.org/10.1103/PhysRevA.38.522
• ²⁶
  Feder J. Random sequential adsorption. J Theor Biol. 1980;87(2):237-54. http://dx.doi.org/10.1016/0022-5193(80)90358-6
  » http://dx.doi.org/10.1016/0022-5193(80)90358-6
• ²⁷
  Geuzaine C, Remacle JF. Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities. Int J Numer Methods Eng. 2009;79(11):1309-31. http://dx.doi.org/10.1002/nme.2579
  » http://dx.doi.org/10.1002/nme.2579
• ²⁸
  Tian W, Qi L, Chao X, Liang J, Fu M. Periodic boundary condition and its numerical implementation algorithm for the evaluation of effective mechanical properties of the composites with complicated micro-structures. Compos, Part B Eng. 2019;162:1-10. http://dx.doi.org/10.1016/j.compositesb.2018.10.053
  » http://dx.doi.org/10.1016/j.compositesb.2018.10.053
• ²⁹
  Nemat-Nasser S, Hori M. Micromechanics: overall properties of heterogeneous materials. Amsterdam: Elsevier; 1993. http://dx.doi.org/10.1016/C2009-0-09128-4
  » http://dx.doi.org/10.1016/C2009-0-09128-4

Publication Dates

• Publication in this collection
  01 May 2023
• Date of issue
  2023

History

• Received
  07 Nov 2022
• Reviewed
  11 Feb 2023
• Accepted
  06 Mar 2023