Thermomechanical Analysis of Ceramic Composites Using Object Oriented Finite Element Analysis

Abstract

This chapter discussed the object oriented finite element (OOF2)-based studies for ceramic composites. OOF2 is an effective method that uses an actual microstructure image of the material/composites for simulation. The effect of filler inclusions on the thermomechanical properties (coefficient of thermal expansion, thermal conductivity, Young’s modulus, stress and strain) is discussed. For this purpose, various ceramics composites (thermal barrier coating and ferroelectric based) are considered at homogeneous and heterogeneous temperature/stress conditions. The maximum stress is found at the interface of the filler/matrix due to their mismatch of thermal expansion coefficient. Further, residual and localized interface stress distributions are evaluated to analyze the composite’s failure behavior. The possible integration of OOF2 with other simulation techniques is also explored.

1. Introduction

Ceramics composites are a significant field of research for the industrialist and researcher [1, 2, 3, 4, 5, 6]. This is because of their wide variety of applications and better mechanical properties such as higher strength, toughness or fracture, etc. Apart from the superior mechanical properties, ceramic composites are extensively used for the thermal barrier coating (TBC) [7, 8, 9] and nuclear fuel cell [6] applications due to high thermal expansion coefficient (α) [10, 11, 12], thermal shock resistance [3, 6] and thermal conductivity (λ) [1, 3, 5, 6, 11]. Many researchers have studied ceramics composites’ effective thermal and mechanical properties with changes in compositions using different modeling approaches [4, 5, 13, 14]. However, the prediction of precise thermal and mechanical properties by various numerical or simulation techniques has limited success so far [4]. Several modeling methods are used to obtain ceramics composites’ thermal and mechanical properties [4, 5, 13, 14]. These methods are categorized in two ways (i) macroscopic and (ii) microscopic. The macroscopic approach is easy to implement and predicts the average or global response of composites. Thus, it considers the volumetric effect of individual phases in the composites. However, the effect of size, shape, orientation, and arrangement of individual phases related to microstructure is wholly ignored. This parameter plays a vital role in the effective thermal and mechanical properties of composites. Therefore, many finite element methods (FEM) are widely applied to predict/analyze and improve the thermomechanical properties of composite materials [4, 5, 13, 14]. FEM can easily integrate the microstructural details of composites and are also computationally intensive [4].

In this direction, the Object-oriented finite element method (OOFEM) can be an effective method in the FEM analysis. This method takes care of complexities of microstructure such as the effect of size, shape, orientation, and arrangement of individual phase, particularly in multiple component microstructure [2]. OOFEM is considered the actual microstructure of the materials compared to conventional FEM, where a “unit cell” model is used to predict material properties [2]. Moreover, the microstructure boundary conditions can be easy to implement. In this context, various researchers have used open-source software OOF2 2D version [1, 3, 4, 7, 15, 16]. The OOF2 software is developed by the National Institute of Standard and Technology (NIST) USA, an open-source tool [http://www.ctems.nist.govt./oof/oof2/index:html#download]. The OOF2 modeling uses microstructural images as input and considered individual phase grain size, shape, local orientation and distribution, etc., with their mechanical and physical properties for analysis [4, 5, 10, 15, 16]. It is used to predict composite thermal and mechanical properties such as λ, α, Young’s modulus (Y) and thermal/mechanical stress–strain contour on microstructure images. Similarly, the residual thermal stresses (σ_r), the effect of thickness, λ and cracking in TBC due to stress relaxation are also investigated [3, 8, 9, 11, 17]. Additionally, these (σ_r, Y, λ and α) properties of composite with respect to filler content and operating conditions are analyzed. In this direction, numerous composites, i.e., 0.94Na_1/2Bi_1/2TiO₃–0.06BaTiO₃/ZnO (NBT-6BT-ZnO), Ni-Al₂O₃, Al/B₄C, Al-TiB₂ Al-MgO, WC-Al₂O₃, AlN-TiN, and Al-TiO₂ are investigated and results are supported by various other techniques [2, 4, 10, 12, 18, 19, 20, 21].

S. Patel et al. [2, 19, 20, 21] predicted the mechanical and thermal properties of NBT-6BT-ZnO, Al-MgO, WC-Al₂O₃ and AlN-TiN composites by the OOF2 method. The simulation results are validated with different analytical models. They show that an increase in the filler content increases the local stress concentration, which can be the start point of failure in the material. S. Patel et al. [2, 19, 20, 21] also studied filler orientation, gradient and uniform temperature environment effect on the thermomechanical properties. Neeraj et al. [4, 10] obtained Y of Ni-Al₂O₃ composite by OOF2 and compared it with the ultrasonic measurement with the possibility of crack initiation due to stress distribution. Andrew et al. [16, 22] studied the complex microstructure using the OOF2 FE tool and reported that the quality of the results depends on the quality of the generated mesh. The set of generics is modified, which results in improve the quality of the 2D mesh. Elomari et al. [23] have analyzed thermal expansion behavior between matrix and reinforcement in terms of silica layer formed during oxidation, size, and thermal stress. Chawla et al. [13, 14] performed a microstructure-based simulation to predict the thermomechanical behavior of composite. Plasma-sprayed Y₂O₃-stabilized ZrO₂ analysis was carried out to measure Y of the actual coating and compared with experimentally observed values [24]. The magnitudes of σ_r in textured and untextured Al₂O₃ are obtained concerning grain orientations with the help of OOF2 [17]. Thermal shock resistance and λ of yttria-stabilized zirconia (YSZ)-Al₂O₃ and 3Al₂O₃·2SiO₂ are computed and compared by analytical methods [3]. Further, λ is analyzed for 4-phase and 3-phase composite of Y₂O₃ stabilized ZrO₂-Al₂O₃-MgAl₂O₄-LaPO₄ and CeO₂-MgAl₂O₄- CeMgAl₁₁O₁₉, respectively to use them with nuclear fuel [6]. Moreover, porous W/CuCrZr composites micrographs with OOF2 are simulated for tensile deformation and thermal conduction behavior [7].

Recently, OOF2 was used to obtain the stress distribution and Young’s modulus in the thermally treated Mg-9 wt.%Li-7 wt.%Al-1 wt.%Sn alloy to enhance the wear resistance [24, 25, 26]. Further, continuously reinforced concrete pavement α correlation with spalling, performance and mechanical behavior is also studied with OOF2 and results are compared with commercial FE software [27]. Furthermore, the λ and thermal effusivity are simulated for the different regions of TBCs as synthesized by thermal spray processes and suspension plasma spray for gas turbine applications [28]. It was found that modeling provides good results [28]. Moreover, Si₃N₄ welding with 316L stainless steel is designed using Mo/Ag composite as interlayer and the residual tensile stress in the interlayer is performed [29]. It can be said that the OOF2 provides a better transition between the mechanical/thermal behavior of a heterogeneous material at the macro-scale and the mechanical/thermal response of its constituent phases. It provides adequate information about microstructures/phases (volume fraction, distribution, orientation) to ensure a realistic assessment of thermo-mechanical properties [30, 31, 32]. More recently, a study was performed on ferritic-pearlitic based steel and found that the predicted results was in good agreement with the experimental results [31]. They found that the microstructural morphology plays a vital role in strain partitioning, strain localization and formability of the ferritic-pearlitic steels [31].

Recently, OOF2 is also used for meshing and FE solutions are obtained by integrating meshing with ABAQUS or MATLAB. In this direction, many researchers have used OOF2-ABAQUS to get the shear stress distribution [26, 33], elastic modulus, crack distribution in the materials/composites [32, 34], residual stresses [29], and thermal properties [28]. Devi lal et al. [24, 26] fabricated 7 wt.% Y₂O₃ stabilized ZrO₂ beam to study the presence of microcracks/pores during bending by OOF2 with ABAQUS. Similarly, crack analysis of duplex stainless steel (ferrite + austenite phases) strength and hydrogen diffusion characteristics are investigated at the microstructure scale [30, 32]. Most recently, the space charge distribution between two-grain boundaries or interface in bulk is also studied [35]. Moreover, to generate input data sets for machine learning, an OOF2 based simulation is performed with a variety and radii of pores for brittle porous materials failure analysis [36]. As discussed above, OOF2 is used by various researchers for numerous applications and integration with other software. Hence, this work focused on the ceramic composites’ thermal/mechanical properties prediction.

This chapter discussed the detailed OOF2 analysis procedure with boundary conditions and assumptions for ceramics composites. The thermal (λ and α) and mechanical (Y) properties are predicted and compared with other analytical methods. In FE analysis, thermal stress–strain contour and heat flux are used to indicate the different temperature conditions on microstructure images. The local stress and σ_r distribution interpenetrating phase and particle-reinforced structure are also studied.

2. Method

2.1 Microstructure image and mapping

Figure 1(a) shows the workflow procedure of OOF2 analysis. The complete OOF2 simulation and analysis process consists of six steps. In the first step, a microstructure image is needed as input; for this purpose, materials scanning electron microscopy (SEM) or back-scattered electron (BSE) images can be used. A BSE microstructure image of NBT-6BT/ZnO composites contains 4% ZnO by volume is used for analysis, as shown in Figure 1(b). The image consists of a length of 24 μm and size of 645 × 484 pixels. In the BSE microstructure, contrast is used to differentiate NBT-6BT (matrix) and ZnO (inclusion or filler phase). The single-phase material pixels are selected by pixel selection and grouped to give each phase thermal and mechanical properties. Figure 1(c) illustrates that the red color pixel belongs to NBT-6BT and the gray color pixel represents the ZnO. It is to be noted that the composite is not fully dense; thus, pores (∼4%) are also obtained, as shown in Figure 1(c) in white color. Then selected single phase pixel materials are assigned with their thermal and mechanical properties, i.e., Poisson ratio, Y (GPa), λ (W/m-K) and α (K⁻¹).

Figure 1.

(a) Outline of OOF2 analysis procedure (b) microstructural image of composites having NBT-6BT-4% ZnO, (c) single-phase material pixel selection, assign group and properties. (Reprint with permission from Ref. [2]. Copyright© 2020 Elsevier BV).

Similar kinds of pixel selection and properties assignment are used by many researchers where two-phase materials or composites such as Ni-Al [37], Al-MgO [19], WC-Al₂O₃ [20] and AlN-TiN [21] for various applications. Moreover, a varying ZnO composition from 4–10% is also studied by Manish et al. [2] for NBT-6BT-ZnO composites. In the next step, microstructural finite element meshing is created by a skeleton, as shown in Figure 2. The geometry of the mesh explained the geometry of the skeleton. It does not cover the information about the finite element shape functions, equations, fields, etc. The skeleton represents the FE discretization of the microstructural. It is composed of quadrilateral and triangular elements, nodes, and segments. In the OOF2, the skeleton is refined by various tools, as discussed in Section 2.2.

Figure 2.

(a) and (b) Microstructural meshing with node movement, (Reprint with permission from Ref. [2] Copyright © 2014 Wiley Publishing Ltd) (c) energy functional approach used during skeleton mesh refinement process [38], (d) final mesh with boundary conditions of NBT-6BT-4% ZnO (Reprint with permission from Ref. [21] Copyright© 2020 Elsevier BV).

2.2 Skeleton modification and refinement

Figure 2(a) and (b) show the nodes of a skeleton element where nodes are on the element’s corner and edges. However, the other FEM tools generally do not consider the nodes along their edges or interiors. In the refinement process, nodes of the edge move along the edge only. Moreover, the interior nodes are free to move anywhere. These node’s movement can also be restricted by the pin node tool. In the Figure 2(a), an element (A, B, C) is shown, which can create an illegal element via movement, as shown in Figure 2(b). The node ‘B’ motion (see Figure 2(a)) generates two illegal elements (Figure 2b). The illegal element refers to an element that breaks the ordering; three collinear nodes also create illegal elements (nonconvex quadrilaterals). The highlighted quadrilateral element represents the illegal element because its nodes (A, B, C) are not convex (Figure 2b).

This type of node movement and element creation depends on the mesh energy parameters. Further, the FE skeleton elements refinement process increases the mesh quality and reduces functional energy parameter (Ё) [36]. Two types of effective energy function maintain the mesh quality as shape energy (Ё_S) and homogeneity energy (Ё_H). The Ё_H is given as [15]:

$${\mathrm{Ё}}_H=1-{\mathrm{Ё}}_{\mathit{HM}}$$E1

where Ё_HM is microstructural components homogeneity energy. Ё_H can be reduced when all components of the microstructure are encompassing a single-phase material. In the case of triangular element Ё_S is defined as [13]:

$${\mathrm{Ё}}_S=1-4\sqrt{3}{A}_T/l^2$$E2

where A_T and l are triangle element area and side length, respectively. For the equilateral triangular component, ${Ё}_S$ is zero and equal to one when all the vertices are collinear (zero aspect ratio). It can also be estimated for quadrilateral elements using a quality factor ${\gamma }_i$ for each corner i as [16]:

$${\gamma }_i=\frac{A_{\mathit{par}}}{l_0^2+l_1^2}$$E3

where A_par is, the area of the parallelogram formed by l₀ and l₁ edge length element converges at node i. ${\gamma }_i<1$ at a corner where the two converging edges have different lengths or meet at an acute/obtuse angle. ${\gamma }_i=0$ when edges are collinear or the length of one edge is zero (degenerate cases). Shape energy for a quad element is estimated by weighted average from the ${\gamma }_S$ at the corner with the minimum γ (${\gamma }_M$) and the corner opposite to it (${\gamma }_0$). The weighted average method is utilized to make sure that each node movement affects the total energy and quadrilateral element Ё_S can be defined as [13, 15]:

$${\mathrm{Ё}}_S=1-\left[\left(1-\eta \right)\bullet {\gamma }_m+\eta {\gamma }_0\right],$$E4

where η = 10⁻⁵ is an arbitrary small parameter. It is needed to avoid pathologies that arise due to the shape energy has no dependence on the position of one of the nodes. The energy functional parameter ($Ё$) is expressed by [13, 16]:

$$\mathrm{Ё}=\mathrm{\delta }{\mathrm{Ё}}_H+\left(1-\mathrm{\delta }\right){\mathrm{Ё}}_S,$$E5

For δ = 0, Ё_S will note benefaction in Ё of the elements. There will be a trade-off between Ё_S and Ё_H when 0 < δ > 1. Some skeleton refinement processes are used to achieve the optimum homogeneity and shape of the element. These energy functions elements are depicted in Figure 2(c).

The homogeneity energy (Ё_H) of the elements is contributing the total Ё when δ = 1 and skeleton refining δ is several according to homogeneity index. For the mesh refinement, firstly, we focus on reduce Ё_H of an element by assigning a higher value of δ. Figure 2(c) displays different types of Ё_S and Ё_H components. In the equilateral triangular element 1, as shown in Figure 2(c), the contribution of Ё_S is zero. Further, element 1 also contain single material; hence the Ё_H value of the element is also zero. Ё_S rise for elements 2 and 3 as the triangular shape deviates from that of an equilateral triangle. Similarly, element 2 also has Ё_H as zero because of homogeneous phase content. Therefore, Ё_S and Ё_H of elements 2 and 3 contribute to the total energy function parameter (Ё).

Subsequently, the skeleton is generated and modified by various refinement tools, i.e., snap refine, anneal, snap node, merge triangle, refine, smooth, etc. In the snap anneal process, each node tries to move and turns to a pixel boundary. Snap refine skelton modifier tool combines the most likable refine and snap nodes features in a single method. It is subdivided into components along pixel boundaries. Snap refine introduces additional edges and new nodes that can be made to follow pixel category boundaries. So, it can also be more flexible than snap nodes in fitting a skeleton, as shown in Figure 2(c). In a smooth tool, nodes move to the average positions of their nearby element. The node satisfied the acceptance criteria and the move is accepted. Smoothing of the skeleton requires several iterations. The merge triangle modifier tool is used to merge neighboring homogeneous triangles to form a quadrilateral element. Rationalize is used to fix an irregular-shaped component of a microstructure skeleton by removing their immediate neighbors. Relax tool moves nodes and improves element shape and homogeneity. Thus, each refinement tool has its refinement process. The final skeleton explains only the mesh’s geometry; it does not cover any information about the finite element shape functions, field, equations, etc. The skeleton creates FE mesh, and the final mesh generation is shown in Figure 2(d). The accuracy of FE solutions be dependent on the homogeneity index and mesh quality. The FE mesh has a 0.98 homogeneity index with 6039 nodes and 8373 elements [38].

2.3 Micro-scale fields and equations

In the case of temperature and displacement applied on the image, the force balance equation and temperature field should be defined in-plane. Fields are defined on two and more graphically overlapping subproblems on a FE mesh. For each subproblem, the field’s value is the same. These values are only stored on a FE mesh. OOF2 used 3D material but solved 2D problems. OOF2 has a generalized plane strain and plane stress approach, implemented by an in-plane elastic modulus. In this software, two types of equations are used plane-flux and divergence. It is activated and deactivated on subproblems. The static divergence equation is expressed as [39]:

$$\nabla ·\mathit{\text{flux}}+\mathit{\text{force}}=0$$E6

where force is generalized force, i.e., heat equation ($\nabla$·J = −∂U/∂T, where -∂U/∂T is a change of energy density) and force balance equation ($\nabla$·σ = f, where f is actual force). The nonstatic equation is a general form of divergence equation (time-dependent version) and given by [39]:

$$M\frac{{\partial }^2}{\partial t^2}\mathit{\text{field}}+C\frac{\partial }{\partial t}\mathit{\text{field}}+\nabla ·\mathit{\text{flux}}+\mathit{\text{force}}=0,$$E7

Heat equation for the first-order problem, M term will be zero. C is the second problem that may be zero. Plane flux equation is a generalization of plane stress and the out of plane component flux is zero and can be written as [39]:

$${\left(\mathit{\text{Flux}}\right)}_z=0$$E8

It can solve nonlinear Eqs. (6)-(8) if the nonlinearity is limited to force and flux terms.

2.4 Boundary conditions

OOF2 FE tool provides five boundary conditions: Dirichlet, floating, periodic, Navman, and generalized force. This work uses Dirichlet boundary conditions on the microstructure image, as shown in Figure 2(d). The Dirichlet boundary condition has identified the value of one component of a field along with a boundary. In both cases, the temperature and displacement application force balance equation and temperature field are defined in-plane. However, when plane strain is used, then the displacement is defined, active, and in-plane. Similarly, for thermal stress force balance equation is used when the temperature is defined, active, and in-plane. The periodic boundary conditions are applied on the left and right vertical edges. The microstructure edge is considered at a constant temperature and/or displacement. In the present work, boundary conditions are depicted in the Figure 2(d), where the left side edge is fixed or at zero displacements (i.e., U|_left = U₀ = 0). When a temperature gradient is used, the same edge is considered under a low temperature (T|_left = T_L). However, the right edge is kept as defined displacement (U|_right = U₁). When a temperature gradient is used, the same edge is considered for higher temperature (T|_right = T_R)_. Moreover, the governing equations are solved iteratively and the convergence criterion is fixed as 10⁻¹³. Finally, thermal or elastic stress/strain contours are plotted for the analysis (discussed in Section 4). It is to be mentioned that these steps are previously discussed in detail by several researchers [2, 13, 16, 19, 20, 40].

2.5 Assumption

The matrix and individual filler properties are considered isotropic, linear elastic and homogeneous. It is considered that matrix and filler interfaces have good bonding. The homogeneity index of mesh is obtained as 98% and the same is used for simulation. In the selected temperature range, material properties are assumed to be constant. In the simulation 2D microstructure is used; thus, the effect of thickness is not considered. The OOF2 software maintained vector (displacement), scalar (temperature) and tensor fields on two and more graphically overlapping sub-problems of a FE mesh. In the case of uniform temperature application, all the edges temperature fixed at the same temperature. The basic iterative matrix equation solver is applied to obtain the solution. The symmetric and asymmetric matrixes are solved using Conjugate Gradient and generalized minimum residual method, respectively. In both cases, a Bi-conjugate gradient with the incomplete lower upper (ILU) preconditioner is used for FE simulation. Once field equations achieve the solution, thermal and mechanical properties are obtained with the help of averaging the applied field.

3. Comparison of analytical and simulation analysis

3.1 Young’s modulus

Young’s modulus (Y) is predicated using OOF2 analysis and results are compared with the analytical method. It is well known that the NBT-6BT ceramics is brittle in nature and have lower fracture toughness. Therefore, ZnO addition can improve the $Y$ and elevate other properties of the composite. The left edge (${\left.d\right|}_{\mathit{\text{left}}}$) of microstructure is considered at zero displacement and right edge (${\left.d\right|}_{\mathit{\text{right}}}$) is given a uniform displacement as:

$${\left.d\right|}_{\mathit{\text{left}}}=d_0\left(d_x=0\right)$$E9

$${\left.d\right|}_{\mathit{\text{right}}}=d_1\left(d_x=0.001\right)$$E10

A similar boundary condition is depicted in the Figure 2(d). Y of composites is calculated by:

$$Y_{\mathit{avg}}=\frac{\mathit{\text{stress}}}{\mathit{\text{strain}}}=\frac{{\sigma }_{\mathit{avg}}}{{\mathrm{Ɛ}}_{\mathit{avg}}}$$E11

Where ɛ_avg is average strain and σ_avg is average stress. The small displacement is applied to the microstructure; hence, it can be assumed that only elastic deformation occurs. The calculated Y is given in the Table 1.

Methods	Young’s modulus (GPa)	Methods	The effective α (×10^−6°C^−1)
OOF2	132.80	OOF2	6.41
Voigt	136.08	Rule of Mixture	6.89
Reuss	134.98	Kerner’s Models	6.86
Hashin and Shtrikman lower limit	135.53	Turner’ Models	6.81
Hashin and Shtrikman upper limit	135.41	Scharpery’ Models	6.86

Table 1.

Comparison of Young’s modulus (GPa) and the effective α (×10⁻⁶°C⁻¹) obtained using OOF2 analysis and various theoretical methods [2].

The predicated Y is also compared with various theoretical models, i.e., Reuss rule of mixture [41], Voigt rule mixture [42] and Hashin and Shtrikman [43] are also given in Table 1. Many researchers considered these methods to estimate Y of composites and compared them with OOF2 results [2, 13, 19, 44, 45]. The estimated value by OOF2 is in close agreement with that of theoretical methods values. Similar studies are also considered with various compositions by variation of filler content [2, 13, 19, 44, 45]. As the filler volume increases, the effective Y of composites also increases in most cases. However, the composites Y variations depend on the difference between matrix and filler Y.

3.2 Thermal expansion coefficient and thermal conductivity analysis

In order to predict the thermal properties of NBT-6BT- ZnO composites, the α is estimated by OOF2 analysis. The α is predicted when the microstructure edges are at a 10°C temperature difference. The thermal boundary conditions are depicted in the Figure 2(d), where the left side edge is kept at low temperature (T|_left = T_L = 30°C). However, the right edge is kept at a higher temperature (T|_right = T_R = 40°C). Moreover, the bottom and top edges are kept as adiabatic. The α is obtained by:

$$\alpha =\frac{{\mathrm{Ɛ}}_{\mathit{avg}}}{∆t}$$E12

where ΔT is the temperature difference. The α predicted with the help of OOF2 analysis is given in Table 1. The α is also estimated with the help of different theoretical methods, i.e., Kerner’s model [46], Rule of a mixture, Turner’s model [47], and Schapey’s model [48] and given in Table 1. It can be said that the OOF2 predicted value of α is in good agreement with different theoretical methods. A similar analysis is performed with different filler (ZnO) compositions based composited and found that as the filler content increases, α decreases [2]. This is because the filler’s α is lower than the matrix α. A similar study is also performed to analyze the behavior of α in composites with different filler content by many researchers [1, 2, 12, 13, 49].

Figure 3(a) and (b) show the boundary conditions for WC-Al₂O₃ and AlN-TiN based composites to obtain the λ. In this case, the top and bottom edges are considered insulating to ensure heat transfer in one direction only.

Figure 3.

(a) Boundary condition for WC-Al₂O₃ (Reprint with permission from Ref. [20] Copyright© 2014 Imperial College Press) and (b) heat flux profile of AlN-5%TiN composites for thermal conductivity analysis. (Reprint with permission from Ref. [21] Copyright© 2014 Wiley Publishing Ltd).

The heat flux equation is solved by the conjugate gradient method with the help of a linear solver. The equivalent effective λ of composites is expressed by Fourier’s law:

$$\lambda =A\frac{Q_{\mathit{\text{avgdx}}}}{\mathit{dT}}$$E13

where Q_avg is average heat flux, A is the cross-sectional area, dT is temperature difference and x is width of the microstructure.

The λ of both composites WC-Al₂O₃ and AlN-5%TiN is predicated by OOF2 and found that results agree with various theoretical methods/models (effective-medium theory arithmetic and harmonic mean model, geometric mean model Lewis and Nielsen method, and Maxwell method) [20, 21]. A similar investigation is performed with various Al₂O₃ and TiN content and found that as the Al₂O₃ and TiN increase, λ of composites decreases and increases, respectively [20, 21].

Thus, the increase or decrease of composites λ is dependent on the filler’s λ. The Y, α and λ analysis shows that OOF2 can effectively predict composites’ thermal and mechanical properties. In order to look in-depth at the thermal and elastic stress/strain behavior under uniform and gradient temperature conditions, various contour plots are discussed in Section 4.

4. Microstructural stress-strain analysis

4.1 Stress-strain analysis under uniform temperature

A homogenous temperature is applied to the microstructure and thermal stress distribution is analyzed in the composites. Figure 4(a) shows typical thermal stress (σ_xx) contour plots when a homogenous temperature of 500°C is applied on NBT-6BT-4%ZnO composite. It is found that the stresses are highest at the sharp edges/interfaces of matrix and filler. At the interface, the stress distribution region is varied and nonuniformly distributed in the filler region. This is because the stress variation depends on the size of the filler particles/grains [50]. If the applied temperature is sufficiently large, local stress can reach above the yield strength, resulting in thermal cracking of the composite [17]. It is important to note that a critical radius is needed for the thermal or residual stress-based crack initiation in the composites. However, most of the case length of maximum stress is small compared to the critical length of crack initiation. Thus, the residual or thermal stress does not form a crack in the composites. Therefore, composites will not fail under such circumstances and can be utilized at higher temperatures. Moreover, when the composite is cool down to room temperature, the residual stresses play a vital role. The magnitude of the residual stress varies according to the applied temperature; however, the exact value of residual stress is hard to predict by OOF2. Further, the filler has more stress distribution compared to the matrix because filler/ZnO have lower α (4.3 × 10⁻⁶°C⁻¹) compared to matrix/NBT-6BT (7 × 10⁻⁶°C⁻¹) [2]. However, if the matrix has lower α compared to filler, then vice versa stress distribution can be observed. This type of stress behavior is reported for AlN-TiN composites where filler (TiN) and matrix (AlN) has an α of 9.1 × 10⁻⁶°C⁻¹ and 4.6 × 10⁻⁶°C⁻¹, respectively [21].

Figure 4.

Thermal stress (a) contour plots when a homogenous temperature of 500°C is applied on microstructure/image and (b) variation across the line scan as shown in (a), (c) thermal stress variation under different uniform temperature conditions, (d) thermal strain contours when a homogenous temperature of 500°C is applied. All plots are shown for the NBT-6BT-4%ZnO composite. (Reprint with permission from Ref. [2]. Copyright© 2020 Elsevier BV).

In order to look in-depth at stress distribution around the filler/matrix interface, a line scan is used (zoom part in the Figure 4(a)). The corresponding stress variation is depicted in Figure 4(b). At the interface, stresses increase to ∼10 and 3 times of average stress in matrix and filler regions. Moreover, the filler region has ∼3 times higher stress distribution compared to the matrix region. This is because of the difference between Y and α of filler and matrix of composites. It is found that filler consists of tensile stress of ∼60 MPa, whereas matrix has a compressive stress of ∼20 MPa. The compressive and tensile stress distributions are observed in different regions of composites. When the temperature is applied, one material expands more than others because of a mismatch of α, which compresses the other materials in some regions. However, in NBT-6BT-ZnO composites, the compressive stress value is lower than tensile stress.

In the same way, the homogenous temperature is varied from 100 to 1100°C and contours are plotted at each temperature. It is found that stress distribution is similar to Figure 4(a), whereas the magnitude of stress is varied according to applied temperature. Figure 4(c) shows the maximum stresses (tensile and compressive) when a homogenous temperature of 100 to 1100°C is applied. As the applied temperature increases, both tensile and compressive stress increases. This is because of that at higher temperature large mismatched between α of filler and matrix is obtained. A similar trend is also observed in other ZnO content compositions. Moreover, it was found that as the filler content increases magnitude of the stress (tensile and compressive) also increases. It is because the higher volume content of filler increases the larger surface area for the mismatch of α. A similar stress–strain simulation is performed on various ceramics composites such as Al₂O₃-SiC [51], Ni-Al₂O₃ [10], WC-Al₂O₃ [20], AlN-TiN [21], etc. They found that with an increase in temperature and filler volume, stress distribution also increases. The detailed discussion on the interface stress variation due to mismatch of the α and Y in composites is also discussed by Chawla et al. [13, 14] and Wang et al. [45].

Figure 4(d) shows the thermal strain contour when a homogenous temperature of 500°C is applied on the NBT-6BT-4%ZnO composite. The matrix has a high strain compared to filler due to the mismatch of α. Further, strain looks homogenous in both filler and matrix, resulting in constant stress in NBT-6BT and ZnO. However, Figure 4(a) shows that the stresses are nonuniformly distributed in the composites. Moreover, Figure 4(d) shows positive strain only, whereas negative stress is also observed in Figure 4(a). When a temperature is applied on the composites, each grain/interface regions of matrix/filler consist of different thermal strain distribution, which results in stress variation/distribution at the interface. At some point, the interface experiences large compressive stress that can be negligible or act on a small length compared to tensile stress. Further, it was found that as the applied temperature increases, the strain also increases. It is a result of an increase in α with temperature for both filler and matrix. The thermal strain variations also alter with an increase in ZnO/filler content in compositions. This is because the higher volume content of filler increases the larger surface area for the mismatch of α. In the same way, Al-MgO [19], WC-Al₂O₃ [20], Ni-Al₂O₃ [4, 10], AlN-TiN [21], etc., are also studied to look at the effect of uniform temperature and filler content. Apart from these, numerous other composites are also studied for better strain control via tailoring α of the filler/matrix [5, 10, 18, 19, 23, 41, 51]. Moreover, the effective α in composites analysis is potentially used in TBC and thermal resistance-based applications [7, 8, 9].

4.2 Stress-strain behaviors under gradient temperature

The OOF2 analysis is performed under a gradient temperature condition of 0–800°C for composites (AlN-5%TiN). For this purpose left and right edge is applied with a temperature of 0°C and 800°C. However, the bottom and top edges are considered adiabatic, as depicted in Figure 5(a). The resultant thermal stress contour under gradient temperature is shown in the Figure 5(a). It is observed from Figure 5(a) as the temperature increases from left to right, corresponding stresses also increase. The stresses are found to be maximum at the interface of the filler/matrix and nonuniformly distributed. Further, the stress distribution is localized at the filler region and increases composite failure/crack risk. As the filler content increases, the risk of failure is considerably higher due to clustering. The inclusion of the filler can decrease or increase the stress distribution in composites which again depends on the difference of Y and α of matrix/filler. Moreover, the possibility of crack generation depends on the length maximum stress distribution at the interface or critical length of the crack. The reason for the same is already discussed in Section 4.1 and the same is also valid here. Interestingly, it is found that the stress distribution largely varies according to the orientation of the filler grain. The inset of Figure 5(a) shows stress variation around two grains (i) and (ii). It can be seen that in the region (i) (Figure 5(a)), filler/grain is aligned with the applied temperature consists of higher stress compared to grain (ii) which is almost perpendicular to applied temperature. However, both the grain (i) and (ii) regions have similar temperature zones. Alike behavior is also observed in entire composites. Similar behavior is also observed in numerous compositions of Al-MgO [19], WC-Al₂O₃ [20], NBT-6BT-ZnO [2], Ni-Al₂O₃ [4, 10], etc. Further, with an increase in filler compositions, the magnitude of the stress increases, whereas the grain orientation behavior remains the same. It can be said that the various strength of composites can be fabricated via preferred grain orientation. Moreover, a large amount of thermal/residual stress can be increased/decreased in composites according to the direction of the grain (transverse and longitudinal). Hence, if composite is synthesized by oriented grain, i.e., filler as platelets/nanorods instead of power can show better performance. Recently, BaTiO₃:100xZnO composite ceramics are made with a micrometer and nanometer scaled power, supporting these observations [52].

Figure 5.

(a) Thermal stress contour when a gradient temperature of 0–800°C is applied, (b) maximum upper (tensile) and lower (compressive) thermal stress and (c) thermal strain contours when a gradient temperature of 0–800°C is applied for AlN-5%TiN composites. (Reprint with permission from Ref. [21]. Copyright © 2014 Wiley Publishing Ltd).

Similar to Section 4.1, tensile and compressive stress is observed in the gradient temperature conditions, whereas the magnitude is lower than that of uniform temperature conditions. The tensile and compressive stress increases with applied temperature, as shown in the Figure 5(b). The left side edge temperature is fixed as 0°C, whereas the right side temperature varies from 100 to 1100°C. The stresses linearly increase with the applied temperatures. Several researchers have reported alike trends in various composites. Thermal strain contour under a gradient temperature condition of 0–800°C for composites (AlN-5%TiN) is presented in Figure 5(c). It is observed that as the temperature increases from left to right, the corresponding strain also increases. The filler has a higher strain than the matrix because of filler has higher α compared to the matrix. Further, a large number of pores are observed in the AlN-TiN composites. When a temperature is applied, filler/matrix can easily expand at pores; hence, the variation in stress distribution is not observed. These types of strain contour and analysis are also discussed in the literature [2, 5, 13, 20]. It can be said the OOF2 can be an effective tool to see the effect of grain orientation or stress concentration at the interface.

4.3 Elastic stress-strain analysis

The elastic stress-strain is studied in the composite by deforming microstructure, as presented in Figure 6. Thus, the left edge is kept as fixed U|_left = 0 and the right edge U|_right = 0.3 μm is deformed and resultant stress-strain contours are depicted in Figure 6 for Al-10% MgO composites. Figure 6(a) shows that the matrix (Al) has lower stress as compared to filler (MgO). This is due to the fact that the matrix has lower Y than filler as 70 GPa and 294 GPa, respectively [19]. Further, maximum stresses are obtained at the interface of the filler and matrix in an entire composite. It gives crucial information about the composite that interfaces play a vital role in effective load transfer between filler and matrix. The sharp edges of the grain have enormous stress (stress intensification occurs) compared to smooth edges. These edges of filler act as stress concentration sites and are more pronounced when filler/grain is oriented in the direction of the applied load. Hence, if the uniform size of grain/filler is added to the matrix, then the strength of the composites can be improved. Similar behavior is also observed when the filler content is increased from 5–15%. However, the magnitude of stress increases with the inclusion of filler due to more frequent particle-to-particle contact. It is also observed that grain orientation also plays a vital role in stress distribution. The effect of filler orientation on elastic stress distribution is the same as discussed in Section 4.2.

Figure 6.

(a) Elastic stress (b) elastic strain contours for composites under μm displacement. Al-10%MgO composites. (Reprint with permission from Ref. [19]. Copyright © World Scientific Publishing Company).

The stress distribution can help to obtain the composite failure analysis. If the length of stress distribution (higher than yield stress) is more than the critical radius, crack generation will occur. Another software integration or analytical calculation is needed to find the exact failure stress analysis because OOF2 cannot provide direct failure analysis. The elastic behavior under various deformed conditions in different composites is previously studied [2, 5, 13, 20, 21]. They also found a similar type of stress–strain contour. However, the magnitude of the stress varies according to composite compositions. Effect of Y in matrix and filler interface, load transfer and orientation effect in composites are discussed by Chawla et al. [13] and Wang et al. [33], respectively [45]. They found that closely spaced filler encourages more deformation within the matrix. Moreover, the neighbor of weak filler zones evolves and, in the end, leads to mechanical failure. These low and high stresses can be calculated and used to forecast the composite’s failure displacement and temperature. Additionally, Figure 6(b) depicts the contour of strain variation under the elastic displacement condition. The matrix has a higher strain compare to filler and is concentrated at interfaces. A higher strain is observed at closely placed filler because of a large degree of constraint region. Other researchers in different ceramic composites previously observe a similar effect. This region can lead to failure of the composites; in this direction, matrix, filler and filler-matrix interface failure criteria are provided by Ha et al. [53]. It can be said that the geometric inhomogeneity in composites produces the nonuniform strain/stresses due to mechanical or thermal loadings. Thus, failure can initiate at any point within the composite, filler, matrix, or filler-matrix interface and consists of different failure processes according to the initiation of critical points. Recently, a similar boundary condition (shown in Figure 6) is also applied for NBT-6BT/ZnO composite to obtain the elastic stress–strain response [2, 38]. It was found that the stresses are concentrated at the sharp edges or interface of the composites and increase with the ZnO concentration. However, strain is uniform throughout the composites for particular compositions and varied with ZnO content [2, 38]. Further, the predicated stress–strain are compared with experimental data, which shows good agreement and validates findings by OOF2 [2].

In summary, ceramics composites have a wide variety of applications in aerospace, defense, energy, medical and transport sectors due to higher strength, toughness or fracture, etc. [54]. They are also extensively used for the TBC [7, 8, 9] and nuclear fuel cell [6] applications due to high α [10, 11, 12], thermal shock resistance [3, 6] and better λ [1, 3, 5, 6, 11]. However, predicting precise ceramics composites’ thermal and mechanical properties by numerical or simulation techniques has limited success so far. In this direction OOFEM based analysis only take care of the effect of size, shape, orientation, distribution and arrangement of individual matrix/filler phases during the analysis. Hence, considering these individual phases properties, it is used to predict ceramics composite λ, α, Y, σ_r, wear, cracking analysis and thermal/mechanical stress–strain contour on microstructure images. It can be said that OOF2 is an effective tool to obtain real microstructural analysis compared to other existing FE analysis techniques. Moreover, it provides the input meshing/parameter for various commercial software to achieve real-time thermal/mechanical/electric analysis. Thus, ceramic composites’ effective thermal/mechanical properties (λ, α, Y) should be analyzed with OOF2 to obtain the real interaction between filler/matrix or distribution of σ_r, wear and crack in individual phases.

4.4 Integration of OOF2 with other software

The OOF2 software is also integrated with other software where composite failure analysis is performed. In this direction, OOF2 is utilized for mesh generation because it uses the actual microstructure of composites. Then the mesh is imported into the other software such as ABAQUS for further analysis. In this direction, OOF2 is used for Al-Al₂O₃ composites, matrix, pores and particles selection in microstructure; then mesh is generated and converted to two dimensional-FE to analyze in ABAQUS [55]. Similarly, W/Cu composite is considered for preprocessing by OOF2; however, the commercial FEM code ABAQUS is used as a solver [11]. Moreover, SiC/SiC composite microstructure generates FE by OOF2 from simplified images of the composite cross-sections. The final FE mesh is exported directly into the ABAQUS to analyze effective material properties and stress distributions [56]. OOF2 and ANSYS are used to obtain the λ based on 2D and 3D calculations, respectively [57]. They found that a power-law function can accurately define the relationship between 2D and 3D results [57]. Further, they concluded that various microstructure images with pores are also valid for λ and 2D image-based model is easy to implement [57]. The YSZ based TBC image is imported in OOF2 and a mesh is constructed; then, a nonlinear elastic–plastic simulation model of micro-indentation is used for further simulation performed by ABAQUS software [58]. Finally, a three-dimensional (3D) OOFEM is also developed by NIST USA to simulate materials’ overall mechanical, dielectric, or thermal properties using actual or simulated micrographs [59]. In summary, it can be said that OOF2 can be used in varieties of microstructure mesh generation, which is further integrated into other software for detailed analysis.

5. Conclusions

OOF2 based FE analysis is discussed for various ceramics composites. The thermal (thermal expansion coefficient, thermal conductivity) and mechanical (Young’s modulus, residual stress) properties of the composites are predicated by OOF2 simulation and results are comparable with the analytical methods/models. A uniform temperature condition is applied to the composites and stress–strain distribution is analyzed. The maximum stress is found at the interface of the filler/matrix due to their mismatch of thermal expansion coefficient. The compressive and tensile stress distribution are observed in different composites with variation in their magnitude and localized at the interface. Similarly, a temperature gradient is also applied on the composites and found that the stress distribution depends on the orientation of the filler/grain. The thermal stress behavior is altered by fabricating composites via nano/micro-grain or making all grains parallel/perpendicular to the applied temperature field. The residual stresses are also varied according to the size, orientation and shape of the filler. Elastic stress also consists of a similar effect with filler size and orientation. A number of ceramics composites show analogous behavior at uniform/gradient temperature and displacement conditions. OOF2 mesh generation can be integrated with other software for failure analysis.

Acknowledgments

The authors want to thanks Mr. Manish Kumar Meena for providing one figure and write-up.

Conflict of interest

The authors declare no conflict of interest.