Multi-Scale Modeling of the Thermo-Mechanical Behavior of Cast Iron

This work presents a multi-scale modelling framework for thermo-mechanical behaviour of Compacted Graphite Iron cast iron. A general thermo-elasto-visco-plastic model is developed to describe the matrix (pearlite) behavior under thermo-mechanical cyclic loading, for which the parameters are identified from tests on pearlitic steel. The pearlite model takes into account the temperature dependent rate-dependency and kinematic hardening. The importance of properly accounting for the graphite anisotropy is emphasised, for which a numerical procedure for estimating the local anisotropy directions from the graphite particle geometry and experimental observations is proposed. A high quality conforming finite element mesh is generated on a representative volume element using discrete voxelized microstructural data in combination with signed distance functions from the interfaces. For fully constraint thermal cyclic loading conditions with different holding times, the capabilities of the developed multi-scale model are demonstrated at both scales: the macroscale, where the simulation results are in very good agreement with the experimental data, and the microscale, providing the evolution of local fields.


Introduction
Many components subjected to thermo-mechanical loads are made from cast iron. Examples include engine components such as the cylinder head, block and exhaust manifold, where the lifetime due to severe cyclic thermo-mechanical loads is a matter of concern for design and maintenance.
Cast iron is categorized as a ferrous alloy, which has a carbon content higher than 2.14% in weight. The high carbon content in the alloy, as well as the presence of silicon, lead to graphite dissolution and formation of graphite particles. Different graphite inclusion topologies provide different mechanical and thermal properties to cast iron. Based on the graphite shape, cast iron can be classified in three types [1][2][3]. Flake, or lamellar, graphite iron (FGI), also known as grey iron, has high thermal conductivity and low strength, whereas ductile, or spheroidal, graphite iron (SGI), on the contrary, has low thermal conductivity and higher strength. The need for materials with a more optimal property combination, has led to the development of another type, compacted graphite iron (CGI), which has intermediate characteristics between SGI and FGI regarding strength and thermal conductivity. This favorable combination of thermal and mechanical properties has made CGI an ideal candidate for application in engine and other thermally loaded components.
For thermo-mechanical analysis of cast iron, recently, several phenomenological models have been developed, e.g., [4][5][6][7]. Despite being computationally rather efficient, the disadvantage of these models is that they can not include the essential influence of microstructural features on the cast iron behaviour, entailing many parameters to be identified, requiring a lot of experimental data.
Multi-scale modeling of cast iron can provides, novel insights into the underlying micro-scale mechanisms, which are inaccessible in macro-scale phenomenological modeling. Accordingly, attention towards modeling cast irons at the micro-scale has increased recently. Most of these models were implemented in 2D [8,9] and are unable to give a comprehensive insight into the behaviour of the material with its intrinsic 3D microstructure. In other studies, an idealized 3D unit cell was constructed for the analysis of specific cases, e.g., SGI [10] or FGI [11]. However, for complex cases, in which the material morphology does not obey a generic pattern, such as in CGI, such an idealized unit cell is not representative. In [12], X-Ray tomography was used to examine the size, shape and spatial connectivity of graphite particles in high-strength CGI. A coral-tree-like morphology was revealed beneath the outer surface of the examined sample indicating the necessity for the development of 3D models, instead of 2D. In [13], 3D finite element models for SGI, FGI, CGI were constructed by interpolation between 2D X-ray tomography images. For the matrix behavior, a rate-independent Johnson-Cook model was used and graphite was modelled as an isotropic elastic material. The model was not examined under thermo-mechanical loading conditions; an experimental comparison for model validation was also not presented.
Given the wide-spread application of cast iron in severe thermo-mechanical loading conditions, an adequate multiscale model should contain the following ingredients, which will be addressed in this paper, representing the major novel contributions of this work: • A 3D representative volume element (RVE) which captures the real 3D microstructure; • Proper constitutive models for both the matrix and anisotropic inclusions; • Applicable to thermo-mechanical cyclic loading conditions.
With respect to the RVE generation, two prevailing approaches can be identified in the literature. One approach is based on synthetic statistics-based algorithms and another on direct digitalization of experimental data. The shortcoming of the first method is that it is usually applicable to simple geometries only, e.g., materials with spheroidal or ellipsoidal shape inclusions [14]. Examples of statistically generated RVEs for more complex, 3D interconnected, microstructures can be found in [15,16]. However, these synthetic statistical descriptions have so far only been tested for the prediction of linear macroscopic effective properties, which can be obtained as closed form expressions depending on statistical parameters describing geometerical characteristics. Exploitation of these methods in non-linear or localized (e.g., damage) regimes is yet to be investigated. The second approach, i.e., RVE construction based on experimental scans, gives a more precise geometrical description of microstructures. However, it needs dedicated tools to make the experimental data suitable for numerical analysis. To this aim, for complex geometries, a level-set based method was introduced in [17], which has no restrictions with respect to the inclusion shapes, allowing for an implicit mathematical description of the interfaces between microstructural constituents. This mathematical description of the geometry can be used as an input to mesh generation software. To generate a high-quality finite element mesh suitable for efficient nonlinear analysis, in [18,19] ad-hoc and force-based methods were used for constructing tetrahedral meshes. In the present work, this approach will be applied to a CGI microstructure, to the authors' best knowledge for the first time.
Another essential ingredient in setting up the multi-scale analysis of cast iron is the need for appropriate time and temperature-dependent models for the microstructure constituents, i.e. the matrix and the graphite. The matrix in cast iron can be ferritic or pearlitic, or a mixture of the two phases. In this work, focus will be put on a pearlitic matrix, often present in CGI cast iron for thermo-mechanical applications. In the literature, several fine scale models for pearlite have been presented, explicitly taking into account the lamellar ferrite-cementite substructure of pearlite colonies, e.g., [20][21][22], however, restricted to isothermal, rate-independent conditions. Another class of pearlite models, disregarding the fine scale structure and treating pearlite as a homogeneous phase have also been proposed for rate-independent [23] or rate-dependent conditions, representative for machining and tooling applications [20,24].
Only a few studies have so far been devoted to modeling the cyclic, time-and temperature-dependent behavior of pearlitic steels [25,26]. In [25], a model is presented to describe temperature-dependent creep and stress relaxation of pearlitic steels at finite strains. It, however, cannot properly describe the cyclic behavior since it includes isotropic hardening only, which does not account for the presence of back-stress effects emerging from the ferrite-cementite deformation incompatibilities in pearlite [20]. The model presented in [26] takes into account temperature effects and kinematic hardening, however, time-dependent effects are not addressed. Motivated by the literature, a new model will be presented here to describe the cyclic, thermal and time-dependent (i.e. creep and stress relaxation) behavior of pearlite. Although the model is developed and identified specifically for pearlite, it is general and can be applied to other metals.
Graphite has a hexagonal crystalline layered structure. The bonding between these layers (basal planes) is weak due to the van der Waals forces, while within the basal planes a strong covalent bonding is present. This layered structure gives graphite a pronounced anisotropy in its thermal, mechanical and other physical properties. Despite this fact, graphite has often been modelled as an isotropic phase in cast iron [27][28][29]. Although this might seem acceptable for nodular graphite, a large ambiguity exists with respect to its isotropic material constants, revealing a large spread between the upper and lower bounds [30], even including higher order corrections [31]. To our best knowledge, there are only a few papers where the anisotropic behavior of graphite in cast iron is addressed. The graphite anisotropy description for nodular graphite, where the arrangement of graphite crystallites in a nodule is well understood, has been investigated in [32]. The structure of the irregular shaped graphite flakes in FGI, and especially in CGI, is much less understood. Pina et al. [11,33] have modelled the anisotropy of flake graphite particles by assuming that the lamellar graphites grow parallel to basal direction. However, only a single anisotropy direction was assigned to the whole particle; this method thus cannot be applied to complex shaped particles with ambiguous growth directions. In this work, a systematic way of describing the local particle orientation based on geometrical features, along with the assignment of its local anisotropy directions, will be presented.
In this paper, a multi-scale model for CGI cast iron under thermo-mechanical loading conditions is presented. A general thermo-elastic-visco-plastic model for pearlite is developed in "Thermo-Elastic-Visco-Plastic Model for Pearlite" section. The pearlite model is aimed to capture the differences in ratedependency at low and high temperatures, typical of pearlite. In order to include the cyclic effects, kinematic hardening is introduced through a multi-component Armstrong-Fredrik type back-stress. The model parameters are identified from experimental data for a pearlitic steel. In "Graphite Modeling" section, the graphite anisotropy description of complex shaped particles is presented, which is based on the geometrical features of a graphite particle in combination with experimental evidence from the literature on the graphite crystalline a signed structure. Next, in "RVE Finite Element Mesh Generation" section, the CGI microscale RVE finite element model is generated from actual microstructural data. To this aim, the geometry is described by a signed distance function, which is used as the input for the mesh generation tool to produce a high-quality finite element mesh. Finally, in "Results and Discussion" section, after assembling all the above ingredients, the performance of the developed multi-scale model is assessed by comparison of the model predictions with experimental data for cyclic thermal loading conditions with varying holding times. The paper ends with concluding remarks in "Conclusions" section.

Thermo-Elastic-Visco-Plastic Model for Pearlite
In this section, several modelling choices will be made to describe specific features of pearlite in terms of its visco-plastic rate and temperature dependent behaviour. As explained in the introduction, pearlite is a two-phase composite material consisting of fine ferrite-cementite lamella colonies. In this work, however, pearlite will be considered as a homogeneous isotropic material, since the model developed here is intended to be used at a coarser scale. However, several features of pearlite originating from its composite nature will be phenomenologically included into the model, e.g., through the back-stress and kinematic hardening. Moreover, pearlite exhibits only a limited temperature sensitivity of the flow stress between 20 • C and 350 • C [25], see also section 2.3. In terms of its creep behaviour, pearlite is known to be highly creep-resistant at room temperatures, while it can exhibit significant creep at elevated temperatures [25]. This is due to the fact that cementite provides effective barriers for dislocation movement at room temperature, whereas it becomes highly deformable at elevated temperatures (above ca. 400 • C ) [34].

Visco-Plastic Constitutive Equation
The pearlite model is based on the multiplicative decomposition of the total deformation gradient tensor into its elastic and inelastic parts, with a subsequent split of the inelastic part into the thermal and visco-plastic parts. The general outline of the modelling framework is summarized in Appendix A, while the modelling choices specific for pearlite are presented in the following.
For temperature and rate-dependent models, the evolution equation for the visco-plastic multiplier ̇v p , which provides the magnitude of visco-plastic deformation rate, see eq. (A.11), can be described by a Zener-Hollomon relation given in a general form by where Z(Ξ eq , D) is the Zener parameter, which usually takes into account the strain rate sensitivity, defined in terms of an equivalent stress measure Ξ eq characterizing the current stress state, and the drag stress D, describing the resistance of the material to plastically deform under the applied stress; (T) is the thermal function, typically of an Arrhenius-type. The idea of having a thermal function instead of thermal dependency of the material parameters has been initially proposed in [35,36] and later adopted by many others [25,[37][38][39], since this approach allows separation of the temperature and rate-dependent effects in the model. The most typical choice for the Zener function is a power-law form, of either Perzyna type, i.e., introducing the yield stress, or fully viscous, i.e., without an explicit yield stress. However, using a single power law allows the description of only one active mechanism. When different mechanisms are active, for example a combination of dislocation climb and glide, a hyperbolic sine function has been suggested in the literature [37,40], which behaves as a power law for lower stresses, while for high stresses exponential response is observed. This allows capturing the transition of the creep mechanism from diffusion-controlled dislocation climb to obstacle controlled dislocation glide. For the same purpose, in [38] a product of the power and exponential functions for different stress ranges was used.
(1) vp = (T)Z(Ξ eq , D) 1 3 As mentioned above, pearlite exhibits distinctly different creep mechanisms at room and elevated temperatures, with a low strain rate sensitivity at room temperatures and a pronounced creep at elevated temperatures. To capture this effect, a visco-plastic model is introduced based on an additive combination of a Perzyna type model with a yield stress and a viscous power-law model, each with a distinct thermal function: where Ξ y is the yield stress at a specific temperature, < ⋅ > is the Macaulay bracket operator; and ′ are the fluidity parameters and n and n ′ are the rate sensitivity exponents. Through a specific choice of the thermal functions and ′ and the rate sensitivity parameters, the first and second terms in (2) can represent both the rate-independent and rate-dependent behaviour for different temperature ranges.
The expressions for the thermal functions and ′ are based on those proposed in [25,41]. The thermal function ′ is the classical Arrhenius-type function Since pearlite exhibits only a limited temperature sensitivity between 20 • C and 350 • C (as can also be seen from the experimental data in section 2.3), the thermal function is selected to be constant over this range (and below this temperature range, for simplicity, since modelling at low temperatures is not intended here) [25] ( In (3) and (4), T m and T c are the melting temperature and a selected cut-off temperature, respectively, Q is the activation energy for creep and R is the universal gas constant. Figure 1 shows these two thermal functions for T c = 320 • C and T m = 1371 • C . Note, that the vertical axis in Fig. 1 uses a logarithmic scale. The specific choice of T c = 320 • C is motivated by the observation in [34], where it has been stated that below 300 • C the deformability of cementite is confined to a single slip system, while above 400 • C new slip systems in cementite become active. Therefore, the transition temperature is expected to be in this range. In [25], 350 • C was used as the critical temperature for the ductility transition.
Inserting the thermal functions (3) and (4) into the evolution equation for the visco-plastic multiplier (2), allows plotting its first and second term in Fig. 2, which schematically illustrates the difference between the resulting rate-independent and rate-dependent behavior and their temperature sensitivity at a constant stress. The rate-independent term, Fig. 2a decays rapidly to a constant value and is essentially temperature insensitive, while the rate-dependent part, Fig. 2b, shows a significantly slower, temperature-dependent, decay. This is indeed the pearlite behaviour intended to be modelled here.

Flow Rule and Kinematic Hardening
As pointed out in [20], the presence of deformation incompatibilities between ferrite and cementite lamella's induces To take into account these effects and to be able to describe a cyclic behavior, a kinematic hardening model with back-stress is introduced. The flow direction ̄ and the equivalent stress Ξ eq are then defined through [42]: and where s dev is the deviator of the symmetric part of Mandel stress (defined in (A.7)) and is the back-stress. For more details on the large deformation formulations with kinematic hardening, see [42,43]. The evolution of the back-stress is described by the superimposed multi-component version of the Armstrong-Fredric model [43][44][45][46] where Here, H i is the temperature independent back-stress modulus, b i the dynamic recovery constant for the room temperature and R th,i the thermal recovery term describing the thermal dependency of the back-stress. At lower temperatures, the flow behavior of pearlite is essentially temperature independent. The thermal recovery thus only becomes relevant at elevated temperatures above a certain critical temperature T c . The thermal recovery function is here defined as: Here, a parabolic dependence of the thermal function on the temperature is used. Other formulations can also be found in the literature [37,47]. In the present model, two terms will be used for the multi-component model, thus n = 2 in (7). The first backstress term attempts to capture the initial stage of steep hardening in which the dislocations pile up near ferritecementite interfaces. The second back-stress term describes the reduced hardening rate, which can be attributed to the dynamic recovery leading to the reduction of dislocation densities and emergence of shear bands in pearlite colonies that entails the stress relaxation [48]. For both back-stress components, for simplicity, the same thermal function is used, i.e., R th,1 = R th,2 = R th .
Isotropic hardening will be neglected here, since it is known to be insignificant in pearlitic steels [26]. Hence, no evolution equation is introduced for the drag stress D.

Parameter identification
The elastic, thermal expansion and thermal function parameters of pearlite are listed in Table 1.
For the identification of the visco-plastic model material parameters, the experimental data from [25] has been used. In [25], uniaxial tensile tests at different temperatures and creep tests at 420 • C under different loads have been performed on a pearlitic steel with a chemical composition and pearlite spacing that are similar to those in pearlite in CGI cast iron. The experimental data from [25] is reproduced in Fig 3. The initial yield stress Ξ y at different temperatures has been taken directly from [25]. Next, the parameters for the visco-plastic multiplier evolution A schematic comparison of the behavior of the rate-independent (a) and rate-dependent (b) parts of the visco-plastic multiplier evolution law (2) equation (2) and kinematic hardening (8) and (9) are identified iteratively. First, n ′ and ′ are taken zero, while is fixed at unity, and n, D and the hardening parameters H i and b i are identified from a uniaxial tensile test at room temperature (20 • C ). Then, the thermal recovery parameters R (1) , R (2) and R (3) are identified from uniaxial test data at different temperatures. Next, parameters n ′ and ′ , mostly responsible for the creep response, are determined from the creep data for 420 • C . Finally, parameter values are iteratively improved to provide the simultaneous best match with all available experimental data. The final outcome is shown in Fig 3. The values of all pearlite model material parameters are listed in Table 1.

Graphite Modeling
The graphite morphology, together with its specific crystalline structure, has an important influence on the mechanical and physical properties of cast iron [52][53][54][55]. Graphite  [50]: ] Thermal functions and ′ , eqs. (3) and (4): Kinematic hardening (8): (9): has a layered structure, as schematically shown in Fig. 4(a). Each layer, known as a basal plane, called a-direction, has a hexagonal structure with strong covalent bonds between the carbon atoms. The layers are stacked upon each other in the so-called c-direction and connected through relatively weak Van Der Waals forces. This results in a profound anisotropic behavior of graphite along these two directions [56,57].
Considering the graphite structure, it can be modelled as a transversely isotropic material, whereby the basal plane is the plane of isotropy and the c-axis the anisotropy direction. The elastic stiffness tensor of the graphite can be given in matrix form with respect to the local a-and c-directions as [58] Measurement of the mechanical properties of graphite is an extremely challenging task and requires the use of advanced techniques, e.g., X-ray scattering [59] and ultrasonic resonance [60], which still do not always lead to the desired degree of accuracy. Moreover, experiments performed on different graphite types, i.e., hex-g graphite with a high degree of ordering and turbo-g graphite with random oriented stacking, typically result in large differences of the measured mechanical and thermal properties. It has been shown in [57] that the out-of-plane (c-direction) shear modulus has the greatest variation among different graphite types, being an order of magnitude higher in hex-g graphite compared to turbo-g graphite. To our best knowledge, measurements of the elastic properties of graphite in cast irons are not available in the literature. It can, however, be expected that the cast iron graphite will contain multiple imperfections and therefore have properties that are closer to turbo-g graphite. The elastic constants used in the simulations are summarized in Table 2. Similarly to the elastic constants, the coefficient of thermal expansion (CTE) of graphite also varies greatly between the a-and c-directions. The values used here have been taken from [61] and are given in Table 3. The material properties of graphite are assumed to be temperature independent in the temperature range considered here, i.e., 20 • C − 420 • C.
Graphite is modelled as a thermo-elastic material. The multiplicative decomposition of the deformation gradient tensor into its elastic and thermal parts is used The thermal deformation gradient tensor is given by where a and c are the linear coefficients of thermal expansion in the a-and c-directions, respectively. The stress-strain relation is expressed in terms of the second Piola-Kirchhoff (11) = e . th (12)    stress tensor and the elastic Green-Lagrange strain tensor e = 1∕2( eT ⋅ e − ): where the transversely anisotropic elasticity tensor is given by (10) and 0 < ≤ 1 is a reduction factor which will be discussed in the following.
As pointed out in [32], even with the proper account for the graphite anisotropy, the calculated effective Young's modulus of a SGI cast iron is over predicted compared to experimental values. This indicates that the graphite in cast iron is probably weaker than the stand-alone graphite samples. Indeed, graphite is a quasi-brittle material [62], and micro cracking inside the graphite particles has often been observed, even in an as-received material, especially along the basal planes. As explained in [63], the high in-plane stiffness of graphite can hinder the c-direction contraction during the cooling down process (during solidification or operation). This induces high tensile residual stresses in the c-direction which can easily separate the weakly bonded basal planes. Moreover, since the bonding between the graphite layers is weak and these layers are relatively thin (high slenderness ratio), they are prone to buckling under locally compressive loads, resulting in a reduction of the effective stiffness and possibly entailing more micro-cracking, as observed for composites [64]. Finally, the cohesion between the graphite particles and the matrix is often not perfect and (partial) decohesion is often observed. Since the explicit modelling of the above mentioned effects is beyond the scope of this work, their collective impact on the behaviour of the cast iron has been included into the model through the introduction of the reduction factor , effectively decreasing the stiffness of the graphite particles. A reduction factor = 0.3 has been chosen for the simulations presented in this work. The reduction factor is selected to best fit the experimental data as it will be discussed in "Results and Discussion" section.
The results of the cast iron model with the transversely anisotropic graphite will be compared to a model prediction using isotropic graphite. The isotropic graphite properties have been computed by applying either Voight (iso-strain, upper bound, UB) or Reuss (iso-stress, lower bound, LB) assumptions to the transversely anisotropic graphite averaged over all possible orientations described by the Euler angles 0 ≤ ≤ 2 and 0 ≤ ≤ ∕2 (the averaging over only two Euler angles, instead of three, is sufficient here due to the presence of one plane of isotropy): where the transversely isotropic elastic stiffness tensor is given by (10) and ℂ = −1 is the compliance. For the values of the elastic constants in Table 2, the isotropic Young's moduli computed according to (14) are E UB = 611. 9[GPa] and E LB = 1.63[GPa] (without considering the reduction factor ). A Poisson's ratio of 0.2, as suggested in [63], will be used in the simulations with isotropic graphite; it has been observed that Poisson's ratio only has a minor influence on the results.
The final step in setting up the graphite model requires the identification of the local anisotropy directions for each (integration) material point of a graphite particle. Since experimental acquisition of these directions for a 3D model would be prohibitively heavy and virtually infeasible, a novel approach is proposed here to approximate the directions of graphite anisotropy for a microstructural model. The approach is based on microstructural studies presented in the literature [65][66][67], where it has been shown that the nucleation and growth (crystallization) of graphite particles starts from the outer surface and results in either an 'onion'-like structure of graphite in spheroidal (SGI) cast iron, or a stack of graphite plates along the short direction of the lamellas in the lamellar (FGI) cast iron, as schematically shown in Fig. 4(b). The structure of the graphite in compacted graphite cast iron (CGI) is more complex, but in general can be assumed to be intermediate between that of SGI and FGI.
Based in the above observations, the proposed method is based on the assumption that the graphite layers at the surface of the graphite inclusions are positioned along the inclusion surface and they are stacked parallel to each other through the depth. Therefore, for the points on the surface of the graphite inclusion, the local c-direction is assumed to coincide with the local normal to the surface. For each graphite (integration) point inside the particle, an inverse distance weighting function is introduced, which averages the anisotropy c-direction based on the distance of the given point from all discretized facets at the inclusion surface. The c-direction in a point ⃗ x is obtained as where ⃗ n i is the normal to the facet i and n is the number of surface facets; d i is the distance from point ⃗ x to the i-th facet on the surface and the power p defines the sensitivity to the distance, here selected as p = 2 . Examples of the resulting local anisotropy directions for a few graphite particles computed using this algorithm are shown in Fig. 5.

RVE Finite Element Mesh Generation
Experimental microstructural images are usually provided in pixelized (2D) or voxelized (3D) form. Consequently, the most straightforward approach to create a microstructural geometry and discretization based on experimental data is by directly translating each pixel or voxel into a finite element (or other numerical discretization unit). This approach however typically leads to a uniformly fine discretization, which is often unnecessary. An even more severe drawback of pixel/voxel based discretizations is the jagged, 'stair-case', description of interfaces, which becomes particularly problematic in non-linear regimes and which is inadequate for proper interface modelling. Therefore, algorithms are required that, based on voxelized experimental data, allow generation of good quality finite element microstructural discretizations, with a conforming description of material interfaces and a spatially adaptive finite element size. A work flow for such an algorithm will be presented here and applied to generate an RVE of compacted graphite cast iron. The work flow is illustrated in Fig. 6.
In this work, the data set from [68] has been used as the starting point. The data set contains a voxelized representation of the microstructure of a compacted graphite cast iron obtained by serial sectioning EBSD. Here, only the information on the phase (graphite or pearlite) of a voxel was used, the (pearlite) crystal orientation was not exploited in the present analysis. The analyzed microstructural volume is 500 × 240 × 250 m 3 . The voxelized representation of the graphite particles within the considered microstructural region is shown in Fig. 6(a).
As the first step, the voxelized data is smoothed using a Gaussian filter with a volume fraction constraint, Fig 6(b). Next, the level set distance function description of the microstructure is generated by solving the Eikonal equation using the fast marching method [69]. The zero level set value corresponds to the interface, negative values lie inside the graphite inclusions, while the positive values are located in the matrix, Fig. 6(c).
The signed distance function is next used to compute the element size function that controls the mesh quality, i.e., to locally refine the element size based on geometrical microstructural features such as interface proximity, curvature of the interface and the distance between the inclusions [70,71]. The node distribution through the volume is generated using the octree partitioning method [71][72][73]. In this method, the domain is partitioned in small cubes. The cubes are subdivided into smaller cubes if the length of the cube edge is larger than the value of the size function at that location. The process of refinement is continued recursively until the desired sizes for all cube edges are achieved in accordance with the size function, Fig. 6(d). The obtained cube vertices will be used as the initial positions of the nodes for the subsequent mesh generation.
The process of triangulation for the finite element mesh generation starts from the inclusion surfaces and outer boundaries [71], Fig. 6(e). Next, the rest of the volume is discretized using a constrained Delaunay triangulation algorithm (with the position of the nodes on the interfaces and outer boundaries fixed) using TetGen software [74], Fig. 6(g).
The mesh generated in the previous step may in general still have poor quality, e.g., bad element aspect ratio's. Therefore, as the final step, the mesh is optimized using the 'truss analogy' proposed in [75]. In this method, the (tetrahedral) element edges are viewed as trusses that exert an axial  Fig. 6(f). If necessary, this process may be performed incrementally approaching the target element size. The final finite element mesh of the considered microstructural compacted graphite cast iron is shown in Fig. 6(i).

Results and Discussion
To demonstrate the performance of the developed microstructural model, the behaviour of compacted graphite cast iron is next simulated under thermo-mechanical cyclic loading conditions that are representative for heavy duty engine .
A thermal cyclic load under uniaxial total constraint conditions is considered, i.e., the sample is fully fixed in the axial direction, whereby the total axial strain is zero. Accordingly, the microstructural RVE is fully fixed in one direction, while the transverse sides are left traction free. The specimen is subjected to a cyclic thermal load whereby the temperature is assumed to be uniform over the whole microstructural RVE varying between 50 • C and 420 • C . Several holding times are considered, where the lowest and highest temperatures are held for 30s (short holding time), 1800s (intermediate holding time) and 18000s (long holding time); the long holding time is applied only at the highest temperature, in accordance with the experimental conditions that will be discussed later. Note, that the considered total constraint thermal cycling results in an out-of-phase thermo-mechanical loading condition with the lowest stress occurring at the maximum temperature and visa versa. Fig. 7 shows a sketch of the total constraint conditions and temperature variations during the heating up and cooling down process.  The results of the simulations for the three holding times are shown in Fig. 8. The simulation results are compared to experimental data, kindly provided by the authors of [76,77], where the experimental setup and procedure are described; the data used in this work corresponds to the measurements on unnotched tensile specimens. The simulation results adequately match the experimental data. These tests reveal several characteristic features. Due to the total constraint condition, at higher temperature, a significant amount of stress relaxation takes place, which saturates for longer holding times. Upon cooling, the accumulated inelastic strains result in tensile stresses. At low temperature, the rate dependency is limited, as expected from the behaviour of the pearlite matrix. All these effects are well captured by the developed model.
An advantage of the current multi-scale model based on a microstructural RVE, compared to single scale (phenomenological) models, is that in addition to the overall effective response it also provides information on the distribution of local fields at the scale of the microstructure. These can be of interest for e.g., the evaluation of the onset of damage. In this way, the microstructural model can reveal details in the local microstructural fields, which are very difficult, or even impossible, to quantify experimentally. The plastic strain fields in the microstructural RVE after the 1st and 20th thermal loading cycle for the simulation with the short holding time (30s) are shown in Fig. 9. The areas of high plastic strain initially concentrate in the pearlite matrix surrounding smaller, high curvature, graphite particles. With continued cyclic loading the plastic strain regions also appear near larger graphite particles and spread further away from the particles in the matrix bulk. Figure 10 compares the accumulated plastic strain after a total of 1 hour loading time for the tests with two different holding times, which corresponds to the 20th cycle of the short holding time conditions (Fig. 10(a)) and one fifth of the full cycle of the long holding time, (Fig. 10(b)). It is observed that for the short holding time, plastic strain areas form bridges between the particles since the material has already undergone multiple compression-tension cycles, while for the longer holding time the plastic field is more distributed since the material is still in the first part of the compression cycle.
Next, the impact of taking into account the anisotropy of graphite on the model predictions is investigated. As discussed in "Graphite Modeling" section, graphite is a highly anisotropic material by nature; however, in the literature, for modelling purposes, it is often assumed to behave isotropically, with a wide range of elastic moduli [78]. The prediction of the model developed here with Fig. 8 Microstructural compacted graphite cast iron model predictions versus experimental measurements for cyclic thermo-mechanical loading with three different holding times the anisotropic description of the graphite, as introduced in "Graphite Modeling" section, is compared to the case where the graphite is modelled as isotropic with the Young's modulus given by either the upper or lower bound values obtained in accordance with equation (14). Figure 11(a) shows the overall stress versus time response for the RVE model with different graphite descriptions. Clearly, the stress response varies significantly between the different graphite models; the upper and lower bound nature of the isotropic estimates is also well reflected in the results. Figure 11(b),(c),(d) demonstrate a striking difference in the local microstructural plastic strain fields in the matrix obtained using either the anisotropic or two isotropic graphite models. For the case of the lower bound Young's modulus, where the graphite modulus is significantly smaller than the elastic modulus of the matrix, a different plastic strain pattern is observed (Fig. 11(b)), compared to the anisotropic and upper bound isotropic models (Fig. 11(c),(d)). The latter two models, despite both including stiff graphite, also exhibit noticeable differences in the local plastic strain patterns, once more emphasising the important role of the graphite anisotropy. From these results, it can be concluded that the material behaviour of graphite has a significant impact on the model predictions, both in terms of the global stress response, as well as the local fields, and thus, for example, on the prediction of damage.
Finally, the local equivalent plastic strain field computed using the finite element mesh with a conforming description of the interfaces and the mesh refinement based on the signed distance function is compared to the field predicted by a voxel based finite element model, with approximately the same total number of degrees of freedom in the model. These fields are shown in Fig. 12. Clearly, significant differences in the equivalent plastic strain fields can be observed. The conforming mesh with local refinements provides an adequate description of the local plastic strain concentration areas near certain graphite particle features, while other interface regions do not exhibit significant plastic yielding. In the voxel based mesh, on the other hand, localized plastic strain concentrations are observed along almost the entire interface, which is mostly induced by the 'stair-case' interface description, rather than by the physical particle features. Moreover, on average, the extend of matrix plastic yielding is larger in the voxel based model, compared to the conforming mesh. It can, therefore, be concluded, that a good quality mesh is essential for obtaining an accurate description of local microstructural fields, especially in non-linear regimes.

Conclusions
In this paper a multi-scale approach for modelling the thermo-mechanical behavior of compacted graphite cast iron has been presented. The model is based on a detailed description of the microstructure of cast iron, consisting of complex 3D shaped graphite particles embedded in a pearlitic matrix.
To describe the material behaviour of each of the microstructural phases, suitable constitutive models have been developed and characterized. The thermo-elasto-visco-plastic model for pearlite takes in to account the transition between  the high creep resistance of pearlite at low (around room) temperature and significant creep at temperatures above 350 • C . Moreover, the model includes kinematic hardening, which is essential for a proper description of cyclic behaviour, which is the main interest in this work.
Concerning the graphite, the importance of properly taking into account the intrinsic anisotropy of graphite has been demonstrated. To this end, a numerical methodology to extract the estimated graphite anisotropy directions has been proposed, exploiting experimental observations of the graphite crystallization directions. The brittle nature of graphite, which often leads to local cracking and delamination, has been implicitly accounted for in the model through a stiffness reduction factor.
A systematic procedure has been elaborated for the generation of a high quality conforming finite element discretization on an RVE constructed on the basis of voxelized microstructural data. The microstructure description based on the signed distance function allows to control the mesh quality depending on the distance from the interface and the local curvature.
The model has been applied to thermal cyclic loading of CGI cast iron under total constraint conditions, representative for truck engine parts. An adequate correspondence between the model predictions and the experimental data has been demonstrated for several temperature holding cases, where the amount of stress relaxation and saturation is well described by the model. In addition to the overall stress response, the developed model provides rich information on the development of microstructural fields and related mechanisms. For example, the accumulation of the local plastic strain as a function of the number of cycles and the holding time has been illustrated. Moreover, the impact of the local graphite anisotropy on the plastic strain field in the matrix has been shown.
In the future, the developed model will be enhanced by including graphite-matrix decohesion and damage development in the matrix, which will allow the prediction of the (fatigue) failure of cast iron components under thermomechanical loading conditions. Finally, it should be pointed out that the modelling approach presented in this work is rather general and can be applied to multi-scale modelling of other materials with complex microstructures. The mesh generation procedure can be directly applied to a broad range of various materials, while the proposed thermo-elastovisco-plastic model, upon proper parameter identification, can also be used for other materials that exhibit a temperature dependence of the strain rate sensitivity.

Appendix A. General modelling framework
The model is based on the multiplicative decomposition of the deformation gradient tensor into its elastic e and inelastic ine parts, with a subsequent split of the inelastic part into the thermal th and visco-plastic parts vp This decomposition introduces two intermediate configurations between the reference configuration B 0 and the current configuration B. The first stress-free intermediate configuration B accounts for the pure thermal deformation. The second intermediate configuration B is the mechanical configuration which is free of any elastic deformation and is therefore also stress-free. Figure 13 shows the schematic representation of these intermediate configurations, the corresponding mappings and indicates some of the tensors defined on the different configurations.
The decomposition of the deformation gradient tensor (A.1) adopted here follows the work of Lion [79]. This decomposition implies that the mechanical deformation takes place isothermally after the pure thermal deformation. However, as argued in [79], this decomposition is theoretical and it occurs at a single time instance, i.e., it does not refer to any actual physical temporal sequence. On the other hand, from the implementation point of view, having e at most left instead of the thermal term th , enables pull-back and push-forward operations between elastic and inelastic configurations more conveniently through the elastic deformation gradient tensor without interference of thermal terms, as will be shown in the following.
The spatial velocity gradient tensor is defined on the current configuration and is given by The velocity gradient can be additively decomposed into its elastic, visco-plastic and thermal parts. The elastic part e only depends on the elastic deformation gradient, but the inelastic part vp + th involves the elastic as well as viscoplastic and thermal deformation gradients where Both ̄ th and ̃ th are the velocity gradients of temperature induced displacement but taken in different configurations, i.e., ̄ th is the push-forward of ̃ th from the thermal configuration B to the mechanical intermediate configuration B .
The strain measure used in the present model is a logarithmic elastic strain tensor ̄ e defined as where ̄ e is the elastic right Cauchy-Green tensor in the mechanical stress free configuration, which is related to the Cauchy-Green tensor = T . in the reference configuration through The current formulation is based on the Mandel stress tensor where is the second Piola-Kirchhoff stress tensor. It can be shown [43], that in the case of elastic isotropy, the antisymmetric part of the Mandel stress tensor is zero and the work rate s ∶̇̄ e gives the rate of the elastic stored energy. In this specific case, the symmetric Mandel stress tensor s can be related to the logarithmic elastic strain through the following elastic constitutive relation: where is an isotropic fourth-order tensor of elastic moduli, which is in general a function of the temperature. In the case of small elastic stretch ratios and isotropy it can be described by the Young's modulus E and Poisson's ratio . The thermal deformation gradient tensor for an isotropic material is given by with the second-order unit tensor. The scalar function (T) describes the linear thermal expansion induced by a temperature change (T − T 0 ) relative to a reference temperature T 0 and can be defined through: where is the coefficient of linear thermal expansion, which, in general, may itself depend on the temperature.
The visco-plastic rate of deformation ̄ vp = sym ̄ vp , which is conventionally assumed to take place in the direction of the normal ̄ to the yield surface with a magnitude given by the visco-plastic multiplier ̇v p , is defined through: It should be remarked that the use of the logarithmic strain furnishes an elasto-plastic correction that is similar to the one used in small strain frameworks, which eases the numerical implementation, i.e., where e * is the trial logarithmic strain defined by keeping the stress-free configuration from the previous increment frozen and assuming a purely elastic increment; The superscript t + Δt refers to the value of the tensor at the indicated discrete time instance.
Acknowledgements This work has been carried out under project number F23. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.