1 Introduction

Compacted graphite iron was first used in 1948 (Dawson 2008) and is considered an important metal matrix composite in the industry thanks to its beneficial physical and mechanical properties, thermal conductivity, and competitive price. CGI has wide applications in pipes and machinery, as well as in the automotive industry for disc brakes (Behera 2012), cylinder heads and blocks (Dawson 2008).

As a typical in-situ metal matrix composite, cast iron is composed of graphite, pearlite and ferrite (Chen and Chen 2011). The microstructure of CGI includes graphite inclusions and a metallic matrix, two phases with considerably different thermomechanical properties. Based on their circularity, graphite particles are classified as flake, vermicular and spherical (Fig. 1).

Fig. 1
figure 1

Microstructure of compacted graphite iron

Graphite-matrix debonding is experimentally observed as the main damage mechanism of compacted graphite iron (Di Cocco et al. 2010; Qiu et al. 2016a). Graphite is brittle and soft, generally separating from the metallic matrix at low stress levels (Di Cocco et al. 2014; Norman and Calmunger 2019) and resulting in plastic deformation of the matrix (Di Cocco et al. 2010). This interfacial-debonding (decohesion) phenomenon initiates microcracks in the matrix; the process of coalescence of such microcracks may cause main cracks or macrocracks (Di Cocco et al. 2010). Furthermore, under cyclic mechanical loads (Di Cocco et al. 2013), interfacial debonding can cause the initiation of fatigue cracks (Kohout 2001; Endo and Yanase 2014). Also, microcracks between the graphite particles and the matrix can result in yielding in cast irons (Hornbogen 1985; Wang et al. 2007; Qiu et al. 2016a).

In diesel engines that are manufactured in many cases of cast irons, the inner wall of the combustion chamber can reach 400–600 °C (Tang et al. 2022). Due to such high-temperature service conditions of automotive engines, CGI experiences intense thermal loadings. The mismatch in coefficients of thermal expansion between the graphite particles and the matrix can lead to interfacial debonding between these two phases, resulting in microcracks even without mechanical loads. The propagating microcracks can form a network and result in a main crack under high temperatures (Qiu et al. 2016a). Also, CGI softens when exposed to severe thermal loadings in engines enhancing debonding at lower stress. At room temperature, debonding occurs above 495 MPa (Wu et al. 2019). On the contrary, it was observed at stress levels below 50 MPa at 723 K (450 °C) (Qiu et al. 2016b). At the macroscale, the tensile strength of CGI decreases when the temperature surpasses 300 °C (Selin 2010). Therefore, the mechanical and thermal properties of graphite and matrix are significantly affected by thermal loads, which requires further investigation.

Morphology, size and distribution of the graphite phase are important factors in determining the damage mechanisms of cast irons (Di Cocco et al. 2014; Zhang et al. 2018; Salomonsson and Jarfors 2018; Zhan et al. 2022). The higher content and size of vermicular graphite particles cause easier initiation of cracks (Qiu et al. 2016a). Also, it was found that longer graphite inclusions reduced the fatigue life (Kihlberg et al. 2021). Both stiffness and toughness of composites changes with increasing volume fraction and aspect ratio of filler (Safaei et al. 2015). But in contrast to polymer-based composites, in cast irons the reinforcement is a weaker and softer constituent. Under increasing environmental temperature, the Young’s modulus of nodular cast iron remains constant whereas its ultimate tensile strength and yield stress decrease (Šamec et al. 2011).

Thanks to their various applications in industry, cast irons were researched for many decades with the development of science and technology. Advanced equipment for observation such as microscopy provided opportunities for researchers to examine and explore the surfaces of cast irons at microscopic scale in 1863 (Stefanescu 2019). The first macroscopic constitutive model for grey cast iron was developed in 1976, providing an accurate load to predict its failure (Frishmuth and McLaughlin 1976). In modern research, the thermomechanical behaviour of cast iron is numerically analysed mainly with phenomenological and micromechanical modelling schemes. Based on phenomenological descriptions, the yield surface and hardening parameters were modified to accommodate microstructural features (Frishmuth and McLaughlin 1976). The micromechanical methodology, on the other hand, focuses on predicting the effect of microstructure of a material directly from experimental observations (Andriollo et al. 2019; Yang et al. 2021). In addition, the constituents were usually assumed to be anisotropic (Andriollo et al. 2016) or elastoplastic (Andriollo et al. 2015).

So far, interfacial debonding was studied mainly under mechanical loads and there is not enough information in the literature about the thermal damage mechanism of CGI at the microscale under pure thermal loading (Palkanoglou et al. 2022). Due to differences in the coefficients of thermal expansion between its main constituents, CGI is vulnerable to high temperatures. However, because of its complex microstructure, it is difficult to study thermal debonding thoroughly only with experimental or microstructural studies. Therefore, this work assesses the effects of graphite morphology on high-temperature behaviours of CGI under thermal loading with three-dimensional numerical models. The inputs of into these models were obtained from the statistical characterisation of SEM images and mechanical testing.

2 Methodology

2.1 Microstructural characterisation of CGI

Statistical analysis of microstructure was performed for CGI specimens (EN-GJV-450). A set of 20 images was acquired using scanning electron microscopy (SEM) and analysed with ImageJ software. The image processes in ImageJ include scale setting, threshold adjustment, and outline shape adjustment of graphite (Fig. 2). The obtained results present the microstructure characterisation of CGI in terms of the perimeter, area, volume fraction, circularity, lengths of major and minor axes of graphite inclusions.

Fig. 2
figure 2

Procedure of statistical analysis of CGI micrographs

Based on the obtained 2D micrographs, the following assumptions about the geometry of both graphite and matrix in 3D numerical simulations were adopted:

  • nodular (spherical) graphite in 2D was simulated as a sphere in 3D;

  • vermicular graphite particles in 2D corresponded to regular ellipsoids in 3D;

  • the major axis of the vermicular graphite was considered equal to the diameter of a nodular graphite;

  • the rectangular matrix domain in 2D was assumed as a cube with equal length, width, and height in 3D.

A set of unit cells in 3D was created in Abaqus (Version 2021) using the methodology of representative volume element (RVE), where a graphite inclusion was embedded into a cubic metallic matrix domain (Fig. 3). An RVE with equivalent morphology could represent a constituent at the macroscale using a micro-unit cell (Zhan et al. 2022).

Fig. 3
figure 3

Geometry of FEA models: a model D1; b model D3 (see Table 3)

2.2 Numerical models

2.2.1 Geometry

Based on the assumptions for the microstructure characterisation results, the dimension of the spherical graphite was calculated with the following relations:

The 2D parameters of nodular graphite particle was obtained as follows:

$$A=\pi {\left(\frac{D}{2}\right)}^{2},$$
(1)
$$S={L}_{2\text{D}}^{2},$$
(2)
$${V}_{{\rm f}}=\frac{A}{S}.$$
(3)

The 3D dimensions of the numerical model follow the expressions:

$${V}_{\text{g}}=\frac{4}{3}\pi {\left(\frac{D}{2}\right)}^{3},$$
(4)
$$V={L}_{3\text{D}}^{3},$$
(5)
$${V}_{\text{f}}=\frac{{V}_{\text{g}}}{V}.$$
(6)

Based on Eqs. (1)–(3) and (4)–(6), the dimension of metallic matrix domain in the 3D numerical model was acquired with the following equations:

$${L}_{2\text{D}}=\sqrt{\frac{\pi {\left(\frac{D}{2}\right)}^{2}}{{V}_{{\rm f}}}},$$
(7)
$${L}_{3\text{D}}=\sqrt[3]{\frac{\frac{4}{3}\pi {\left(\frac{D}{2}\right)}^{3}}{{V}_{{\rm f}}}}.$$
(8)

The aspect ratio of vermicular graphite particle \({A}_{\text{R}}\) was acquired as follows:

$${A}_{\text{R}}=\frac{{L}_{\text{major}}}{{L}_{\text{minor}}},$$
(9)
$${L}_{\text{major}}=D.$$
(10)

The numerical models developed comprised a metallic matrix cube containing a nodular or a vermicular graphite inclusion. The dimensions of developed FEA models were established according to the range of spherical graphite diameters, as measured from the CGI micrographs. For nodular graphite, the diameter varied between 0.29–42.53 μm; hence, 15 μm was selected. For vermicular graphite particles, a major axis of 15 μm length was used; the maximum aspect ratio in the statistical analysis was 5, thus, the corresponding minor-axis length was 3 μm. Based on the 2D analysis, the dimension of cubic matrix in 3D was in the range of 25.74–32.23 μm; therefore, 30 μm was chosen in this research for comparability. The volume fraction of nodular graphite was 6.5%, in the range of statistical results (5.2–11.37%) (Palkanoglou et al. 2020). The dimensions of spherical and vermicular particles in the matrix are given in Table 1 for all the studied models. To assess the effects of the orientation of graphite inclusions on the mechanical behaviour of CGI, different orientations of graphite (at 0° and 45° with respect to the X axis) in the central XY plane of the cubic matrix were selected (Fig. 3). The length of diameter for the nodular graphite inclusion and the major axis for the vermicular graphite inclusion was 15 μm, and the distance between boundary of particle and the RVE face was 7.5 μm. Following the path AB of measurement for presentation of results, the region from 0 to 7.5 μm is inside the graphite domain, while the matrix region stretches from 7.5 to 15 μm in all diagrams (Fig. 3).

Table 1 Model dimensions

2.2.2 Constitutive behaviour

As a soft and brittle material, graphite was considered to display a limited plastic behaviour (Seldin 1966; Greenstreet et al. 1973; Andriollo et al. 2015). The graphite and matrix domain of the model had the properties of CGI used in the previous 2D study in Table 2 (Palkanoglou et al. 2020).

Table 2 Constitutive parameters for matrix and graphite

The mechanical behaviours of both graphite and the metallic matrix in this research were described with the J2-flow theory (Palkanoglou et al. 2020). For the graphite phase, beyond the elastic region, a damage model was employed. The initiation of damage was considered when a plastic-strain-based criterion was fulfilled. When the plastic displacement approaches the displacement at failure, the damage initiation was considered to make sure the debonding of graphite happens at the determined temperature matching the SEM experiments. After the onset of damage, the stiffness of the material degraded gradually to 0. The integral ductile damage criterion is given by the following equation (Hooputra et al. 2004)

$${\omega }_{\text{D}}=\int \frac{\text{d}{\overline{\varepsilon }}_{D}^{\text{pl}}}{{\overline{\varepsilon }}_{\text{D}}^{\text{pl}}}=1,$$
(11)

where \({\omega }_{\text{D}}\) is the state variable increasing monotonically. At each increment, the state variable increment:

$${\Delta \omega }_{\text{D}}=\frac{\Delta {\overline{\varepsilon }}_{\text{D}}^{\text{pl}}}{{\overline{\varepsilon }}_{\text{D}}^{\text{pl}}}\ge 0,$$
(12)

where \({\overline{\varepsilon }}_{\text{D}}^{\text{pl}}\) is the plastic strain at the onset of damage.

Damage was considered to evolve linearly with deformation after its initiation. Element deletion occurred when the value of damage exceeded a set threshold at all integration points of a given element. This damage criterion was only used in graphite inclusion and applied to every element in the graphite. The four neighbouring elements surrounding the node with maximum damage were selected. In simulations, the damage parameter was calculated as the average for the four neighbouring elements with maximum damage.

2.2.3 Boundary and loading conditions

Eight-node 3D stress elements (C3D8) were selected to mesh the models. After conducting a mesh convergence study, 1 μm was selected as the size of elements (Fig. 4). 26,400 elements were used leading to 28,957 nodes in model A.

Fig. 4
figure 4

Convergence study for developed model

Periodic boundary conditions (PBCs) and fully-fixed boundary conditions (FFBC) were considered in this study to represent the extreme cases of the in-service constraints. PBCs control the pairs of nodes on the corresponding surfaces with the same level of displacement in opposite directions. These conditions are widely applied in finite-elements analysis to simulate the microscopic and mesoscopic behaviours of the constituents with RVE (Garoz et al. 2019) and allow simulations of deformation of the RVE model using a small domain that represents the corresponding infinitely large system permitting the distortion of boundary surfaces (Omairey et al. 2019). Under PBCs, a pair of any two points \(x\) and \(x+d\) on the two corresponding boundary surfaces with distance \(d\) should meet the condition

$${\mathbf{u}}\left(x+d\right)={\mathbf{u}}\left(x\right)+\overline{\varepsilon }d,$$
(13)
$${\mathbf{t}}\left(x+d\right)=-{\mathbf{t}}\left(x\right),$$
(14)

where \({\mathbf{u}}\) is the displacement and \({\mathbf{t}}\) is the traction at \(x\), respectively; \(\overline{\varepsilon }\) is the average infinitesimal strain over the volume (Drago and Pindera 2007). In some simulations, fully-fixed boundary conditions (FFBCs) were applied at all the nodes along the boundary surfaces of the cubic metallic matrix, constraining all six degrees of freedom of these nodes. For each morphology of graphite unit cell, PBCs or FFBCs were applied for comparability. Fully fixed and periodic boundary conditions can give the upper and lower bounds of the realistic case in practise. The realistic case should be somewhere in between. Thermal loading was used for all unit cells with a linear increase from 25 to 500 °C. Specially, Model A2 was also followed by cooling down to 25 °C and then under further multiple thermal loading cycles (single cycle of 25–500–25 °C).

In summary, the details of FEA models are shown in Table 3:

Table 3 Summary of analysed RVE models

Thus, models A1 to A2 described nodular graphite particles, with all other models dealing with vermicular ones.

Complex 3D shapes of vermicular and flake graphite particles in CGI—as opposed to inclusions in nodular cast irons—prevent the establishment of general features of its thermomechanical behaviour at microscale. In contrast to direct introduction of such shapes (obtained, e.g., with computed tomography) into numerical simulations, this study focuses on the simplified–but controlled–morphology of the inclusions. Together with the concept of RVE, this allows the assessment of the geometrical features on evolution of thermal stresses and strains at microscale and distribution of the resultant damage.

2.2.4 Validation of numerical models

A tensile test was implemented on a CGI (EN-GJV-450) specimen at room temperature and the measured mechanical response was compared with results obtained with the developed FEA models for validation in Fig. 5. The mechanical data for the matrix domain in simulations were based on the literature (Niu et al. 2021). The tensile loading and PBCs were applied to the 3D numerical model at room temperature to emulate the experiments. The validation of the proposed models was implemented by the comparison between numerical simulations and tensile-test experiments since the interaction in microstructure was difficult to measure directly. Monitoring of three-dimensional particle debonding with SEM is not possible. Some in-situ SEM studies of microspecimens loaded with miniature load cells allowed 2d observations of separation between the nodular graphite particles and the matrix in ferritic–pearlitic ductile cast only on the specimen’s surface (Iacoviello et al. 2008).

Fig. 5
figure 5

Experimental and numerical stress–strain curves of CGI (model A2) at room temperature

3 Results and discussions

3.1 Effect of particle shape orientation

Results for the effect of graphite orientation were obtained with the respective models; the temperature of damage initiation in graphite in models with FFBCs is shown in Fig. 6. The graphite with a larger aspect ratio exhibited damage at lower temperatures since it is much sharper than the spherical particle (an aspect ratio of 1). These nodular particles are more difficult to damage. Damage appeared earlier in the thermal cycle as the inclusion aspect ratio increased, in a monotonic but non-linear way. The trend was similar for both graphite orientations, 0° and 45°. It can be concluded that damage in graphite initiates at lower temperatures in less spherical particles while graphite orientation does not affect this process significantly when FFBCs are used, to represent behaviour of constrained parts of the components exposed to the thermal load.

Fig. 6
figure 6

Temperature of damage initiation in graphite for different particle orientations FFBCs

3.2 Effect of graphite morphology

The evolution of the maximum damage in graphite particles with different aspect ratios is shown in Fig. 7 for FFBCs. As the aspect ratio of the inclusion increased, and the damage due to the application of thermal load appeared earlier. The cases with spherical graphite demonstrated higher strength, with the damage occurring at higher temperatures. The evolution of damage was affected by the inclusion shape with the sharper-shaped particles (i.e., higher aspect ratios) exhibiting a more rapid damage evolution. The trends of damage evolution differed significantly even for small variations in the inclusion’s aspect ratio.

Fig. 7
figure 7

Evolution of maximum value of graphite damage in CGI for inclusions with different aspect ratios for FFBCs

Under FFBCs, the unit cell was fully constrained; thus, at heating, the graphite and matrix were compressed by each other due to the mismatch in the coefficients of thermal expansion. The graphite with a smaller aspect ratio has higher circularity and larger volume so that it can withstand a higher compressive load. With an increase in the aspect ratio, the stress–strain behaviour of graphite particles became nonlinear due to the initiation of graphite-matrix decohesion. In contrast, the decreased aspect ratio reduces of the level of compressive stresses at FFBCs and made their mechanical response almost linear. Kinematically constrained graphite particles continue resisting compression even after local failures, under nearly hydrostatic loading; the current model does not account for the local cracking in compressed graphite.

Spatial distribution of damage in graphite inclusions with different morphologies is shown in Fig. 8 after heating to 500 °C. It was found that damage was mainly localised near the interface boundary between the particle and the matrix. For a higher aspect ratio, the volume of graphite inclusion in the direction of the Z axis in model E1 was smaller than that in model A1 (for the spherical particle) leading to a higher volume of the metallic matrix (Fig. 8). As a result, the higher compressive load from the matrix under FFBCs caused a more significant damage in the graphite inclusion in model E1 in Z direction.

Fig. 8
figure 8

Damage distributions in graphite inclusions with different morphologies after heating to 500 °C: a model A1; b model E1

3.3 Effect of boundary conditions

The effect of boundary conditions was studied by comparing the distributions of normal stress along the major axis of the graphite inclusion, from the graphite centre to the edge of the analysis domain (see path AB in Fig. 3a). In FFBCs, the mismatch caused compressive normal stresses in both constituents. The increased aspect ratio resulted in significantly higher stress levels. Transition to the PBCs that allow the thermal expansion of the unit cell fully changed the character of the normal stresses, generating tensile stresses both in the graphite particles and the matrix. Also, the charter of their distribution depending on the aspect ratio changed: stresses in the nodular particle were lower. With the increasing aspect ratio, the stresses grew, while the situation in the matrix was reversed: for the low aspect ratio, the normal stresses were smaller. For the nodular particle, for the PBCs, they were the highest in the matrix.

The calculated levels of the normal stresses in the graphite particles were rather high, potentially in excess of its strength and close to the maximum values observed in the literature (AZOM 2022). Hence, the future models should incorporate the account for the potential fracture of graphite.

The distribution of von Mises stress along the path AB in the models with FFBCs is presented in Fig. 9a. Under thermal loading, the thermal expansion of both phases is restricted by the fixed boundaries resulting in compressive stress. This interaction is also affected by the mismatch of thermomechanical parameters of the two constituents and the particle morphology. When the inclusion shape was closer to spherical, the stress state was effectively hydrostatic, i.e., with vanishing effective stresses. As a result, the von Mises stress in spherical graphite was rather low in all cases. In contrast, its magnitude in the matrix increased at the periphery of the unit cell with the aspect ratio.

Fig. 9
figure 9

Distribution of von Mises stress in CGI for different graphite aspect ratios at 500 °C: a FFBCs; b PBCs

The distribution of von Mises stresses in CGI changes for PBCs are in Fig. 9b; being able to expand, the unit cell did not develop excessive thermal stresses. As a result, the developed stresses in the matrix at the maximum temperature were far below the yield point even in the case of vermicular graphite inclusion. The trend of all curves for different aspect ratios was similar, characterised by considerable gradients in parts of the matrix close to the interface.

In FFBCs, the thermal stresses were still present after the cooling-down stage. The highest level of von Mises stress was 460 MPa at 500 °C and 553 MPa at 25 °C surrounding graphite particle (Fig. 10); this results in the plastic strain of the matrix.

Fig. 10
figure 10

Distribution of von Mises stress in model at different temperature under FFBCs: a 500 °C; b 25 °C

In case of PBCs, the mismatch of thermal expansion between the two constituents led to stress up to 160 MPa after the increase in the temperature (Fig. 11a). Residual stresses were observed in the model after each cycle; however, the level under PBC was negligibly small (up to 6.6 MPa) because of the relatively unconstrained thermal expansion after cooling down to 25 °C (see the difference in the scales at Fig. 11). The stresses in the models with PBC were lower than the yield point due to non-restricted boundary conditions. Thus, graphite damage was not observed in these PBC cases.

Fig. 11
figure 11

Distribution of von Mises stress in model at different temperature under PBCs: a 500 °C; b 25 °C

3.4 Outlook

The current model did not consider interfacial damage processes directly, although some analysis for formulations has been done by authors previously. This is related to the high intensity of the three-dimensional simulation using cohesive zone elements. The next step will be the introduction of these types of elements into simulations and to consider the interface between graphite and matrix.

The present paper approached 3D simulations and used the same damage criterion of graphite without interface. The simplified model without interface would be the basic stage of the series of model. The interface will be considered in the future models. However, direct measurement of interfacial damage properties is not possible and needs development of some additional experimental methods.

From the simplified models, damage was assumed to appear on the outer layer of graphite because of its brittle properties and inability to accommodate plastic deformation. Metallic matrix was generally considered isotropic and ductile. The ductile damage criterion has been used in ferrite matrix phase of cast iron because of the high ductility and moderate yield strength of matrix (Andriollo and Hattel 2016). However, the limitations of ductile damage criterion include its ability to properly describe the plastic volumetric strain (Murakami 2012).

Pure thermal load results in high levels of hydrostatic stresses, the J2 based plastic model is not fully adequate. Thus, the authors consider looking into other plastic potentials (such as Drucker-Prager) for their future simulations.

4 Conclusions

In this work, the effects of morphology and orientation of graphite inclusions as well as the boundary conditions of the studied domain on damage and thermomechanical behaviour of CGI were investigated by developing and employing a set of three-dimensional numerical models. The main conclusions are as follows:

  • The temperature of damage initiation in graphite were similar for both orientations (0° and 45°). Hence, orientation of graphite particles did not affect the damage onset significantly.

  • In contrast, the morphology of graphite influenced the onset of damage and plasticisation. When the aspect ratio of graphite inclusion increased, the initiation of damage in graphite started earlier and develops with increasing temperature.

  • Two types of used boundary conditions—PBCs and FFBCs—provide the lower and upper boundaries for graphite damage mechanisms, respectively. The results of models with FFBCs are more conservative than PBCs in FEA models (models A to E).

This work demonstrated a consistent trend for the graphite damage mechanism for various three-dimensional shapes with the change in the aspect ratios. The next step would be to consider the effect of deviation of real-life shapes of graphite particles from ellipsoids used in this study on thermomechanical behaviour of modern cast irons.