Title: A Digital Twin of the Powder Metallurgical Manufacturing Chain of High Strength Sintered Gears

This paper presents a digital twin of the powder metallurgical (PM) production chain of high-performance sintered gears based on an integrated computational materials engineering (ICME) platform. Discrete and finite element methods (DEM and FEM) were combined to describe the macroscopic material response to the thermomechanical loads and process conditions during the entire production process. The microstructural evolution during the sintering process was predicted on then meso-scale using a Monte-Carlo Model. The effective elastic properties were determined by a homogenization method based on modelling a representative volume element (RVE). The results were subsequently used for the FE modelling of the heat treatment process. Through the development of multi-scale models, it was possible on the one hand to obtain characteristics of the microstructural features and on the other hand the predicted hardness and residual stress distributions allow the calculation of the tooth root load bearing capacity of the heat-treated sintered gears.


Introduction
The rise of digital technologies in the industrial production, the fourth industrial revolution -Industry 4.0 -, requires interdisciplinary software solutions to create a platform for a crossdomain collaboration among the product development, production and usage. The motivation to utilize digital solutions in the production confirms the proof of the concept of the integrated computational materials engineering (ICME) within the field of material science. ICME describes the relationship among the manufacturing history, microstructural evolutions and the component properties on different length and time scales in a holistic system. Thereby, the main focus of the ICME concept lies on developing multiscale material models that describe the material response to inhomogeneous thermal, mechanical and further relevant process conditions (Allison et al., 2006). All decisions in product development, production and usage are based on the selection and understanding of materials. Hence, the use of ICME is essential for the optimizing manufacturing processes or the definition of alternative technologies in the modern production.
Powder metallurgy offers proven advantages as an alternative production route for conventional metal-forming technologies due to its high flexibility regarding the processes, materials, shapes and products (Danninger, 2018). One of the major applications of powder metallurgy is the production of sintered precision parts, used in automotive transmission. The increasing trends towards light weight applications, alternative (hybrid-) drivetrain systems and downsizing of combustion engines motivate the use of sintered gears in highly loaded applications (Leupold et al., 2017). Studies have shown a reduction in the use of materials in the entire production chain of sintered gears by approx. 50% and approx. 10% lower energy consumption (Kruzhanov & Arnhold, 2012). Furthermore, higher flexibility in shape optimization and better performance in noise vibration harshness are other potential advantages of the PM route (Leupold et al., 2017). However, sintered gears generally show lower strengths due to the remaining porosity from the process, which severely hinders their applications in automotive transmissions (Gräser et al., 2014). To improve the load bearing capacity, several techniques were suggested to densify highly loaded sintered gears at the functional surface (Kotthoff, 2003). Surface densified and heat-treated sintered gears can achieve a tooth flank load bearing capacity that is comparable to that of conventional gears at the flank. However, the tooth root bearing capacity still needs to be improved (Frech et al., 2018;Gräser et al., 2014). Computer simulations were used in this study to understand the correlation between the process parameters and the tooth root load bearing capacity of the sintered gears. The developed ICME approach comprises process and material models across scales that cover physical interactions in the entire process chain and describes the microstructure evolutions depending on the manufacturing history. Based on the simulation of the microstructure and residual stresses, the tooth root load bearing capacity of the sintered gear can be calculated.

The Process Chain
In powder metallurgy (PM), metallic powders are processed into semi-finished products and components such as high-strength sintered gears. The main process steps are the production and preparation e.g. by mixing a metal powder, the shaping mostly by mechanical pressing, sintering and, if necessary, a mechanical and/or thermal post-treatment as well as the final hard finishing, see Figure 1. While the shape of the later component is largely predetermined in the forming process, consolidation to a solid material takes place during sintering, which is understood as heat treatment below the melting temperature. Powder properties, shaping process and sintering parameters determine, in addition to the precision of the component, the porosity and grain size of the material. In particular, the mechanical properties can be specifically adjusted by subsequent mechanical densification and heat treatment. Material defects, which mainly affect fatigue strength and fracture toughness and thus the load-bearing capacity and service life of the components, can be created, modified or even eliminated in each of the manufacturing steps. A typical material for sintered gears is the water-atomized iron powder, pre-alloyed with 0.85 wt.% molybdenum and known as Astaloy 85 Mo. A lubricant is added to the powder to improve the compactibility and reduce the friction in the press die, and graphite is added to increase the hardenability. The relative density of sintered gears is typically about 90-95% without further densification (Beiss, 2013). During sintering, under the influence of temperature and time, the lubricant is expelled from the green compact and the strength is increased through diffusion processes. To achieve precision production of sintered gears, sintering process parameters are adjusted to minimize the sinter shrinkage (Altena & Danninger, 2005). The cold rolling process is an industrially widespread process for densifying the surface zone to increase the load-bearing capacity of sintered gears (Frech et al., 2018;Kauffmann, 2013;Kotthoff, 2003;Kruzhanov & Arnhold, 2012). After cold rolling, the heat treatment is carried out by case hardening. Figure 2 shows the microstructure of a tooth characterized by local pores and a density gradient as well as a transition from the martensitic edge to the bainitic core (part b). During carburizing of PM parts, the existing porosity can lead to undesired case hardening depths or carbide precipitations, which could impair the load-bearing capacity of the gear. In addition to porosity, hardness and residual stresses after the heat treatment are decisive factors influencing the gear performance. According to DIN 3979 (DIN, 1979), the main types of damage to gears are tooth breakage, pitting, scuffing and wear. In this work, only the tooth root load bearing capacity is considered. The fatigue fracture of the tooth root usually initiates on the tensile stresses side in the area of the contact point of the 30° tangent to the tooth root rounding, see Figure 2 (part c). The load bearing capacity is determined by the material selection (Doane, 1990;Sauter & Schmidt, I., Schulz, M., 1992), the tooth root geometry and manufacturing history (Feltkamp, 1967;Zuber, 2008) as well as the component size (Strasser, 1984).
To prove the validity of the ICME-approach, geometry, process parameters and gear testing results were adopted from (Hajeck, 2017). A spur gear with a modulus of 2 mm, a pressure angle of 20° and 27 teeth is produced by pressing Astaloy 85 Mo powder mixed with graphite, sintering at 1120°C for 20 min, surface densification by cold rolling, case hardening by gas carburizing and quenching, tempering at 200°C for 2 h and finally grinding. The initial carbon content and density were 0.25 wt.% and 7.1 g/cm 3 after sintering, respectively. More details about the geometry and process parameters are given in (Hajeck, 2017).

The ICME approach
The modelling approach comprises different methods that are applied individually to the corresponding process step to provide input data for the subsequent model. In each step, the "technological" side of the process is modelled by defining simple boundary conditions, which could be referred to as a "digital shadow" of the process. The main focus of the work is on the development of a material's "digital twin": a comprehensive model of "material processes" to quantitatively describe the material response during manufacturing. The process of die filling and pressing is modelled by means of the discrete element method (DEM) in the software LIGGGHTS-PUBLIC (DCS Computing GmbH). The result is the density of the green part on the macroscale. The macroscale sintering simulation is performed in the commercial finite element (FE) software ABAQUS and describes the density and shape of the sintered gear, considering the possible sintering shrinkage. The kinetic Monte Carlo (KMC) method is used to describe the diffusion process during sintering and the resulting morphology of the pores on a mesoscale, which provides the input information for the generation of a representative volume element (RVE) to obtain the effective elastic properties for the continuum mechanical modelling of the heat treatment process in ABAQUS. The heat treatment simulation describes the microstructure evolution, hardness and residual stresses. The density profile after cold rolling is determined by image analysis based on experimental results in (Gräser et al., 2014) and (Hajeck, 2017). Density, hardness and residual stresses are input data for the analytical calculation of the tooth root load bearing capacity.

Modelling of the die filling and powder compaction
Since the die filling and powder compaction process is based on a granular system and each particle in this system moves independently, it is difficult to predict the behaviour of the granular system using continuum mechanical models. In this context, the discrete approach developed for numerical modelling of granular materials at particle scale has become a powerful and reliable tool, which is generally referred to as discrete element method (DEM) (Cundall & Starck, 1979;Pöschel & Schwager, 2005). DEM is a particle method based on Newton's laws of motion. Particles can move with six degrees of freedom (three for translational and three for rotational movement). The particles are defined as soft particles that can only undergo elastic deformation. The simulation of the contact between particles considers a cohesive force (Eq. (1), which is calculated by contact model and cohesion model. Additionally, the external pressure and gravity are also considered (Luding, 2008).
The Hertz-Mindlin model as a widely used contact model describes the normal interaction by a nonlinear relationship representing the elastic contact behaviour between particles. The contact of two particles, and , is considered with radii and . The elastic deformation is adapted by the Hertz-Mindlin model for the normal and tangential contact interactions. In addition, damping in the normal direction and friction in the tangential direction are introduced into the model. A schematic description is shown in Figure 3. When the distance between two particles is smaller than their contact distance + , the model is used to calculate the Hertz-Mindlin contact force .
The first term is the normal force . between the two particles, and the second term . is the tangential force. The normal force has two terms, a spring force , , and a damping force , , . The tangential force also has two terms: a shear force , , and a damping force , , . and are the elastic constant and the viscoelastic constant, respectively. The lower limit of the tangential force . is given by the Coulomb friction where is the Coulomb coefficient of friction. Figure 3: Scheme of the contact model.
The simplified model of Johnson-Kendall-Roberts (SJKR) as the cohesion model adds an additional normal force contribution to the Hertz-Mindlin model. If two particles are in contact, the SJKR model adds an additional normal force ℎ , that tends to maintain contact. In addition to the contact forces, gravitational acceleration is also considered and the gravitational force is calculated.
Newton's equations of motion are used to determine the new positions of the particles. the equations of the interaction forces and moments in a time step for particle read as follows ̈= (4) where is the mass, ̈ is the acceleration, is the Inertia tensor, is the angular velocity, ̇ is the angular acceleration, and force and moment are the sums of all forces and moments that act on the particle.

3.2
Multiscale modelling of the sintering process

Continuum Mechanical Approach
FEM is used to simulate the sintering process in the ICME approach. Among different macroscopic models, the modified Abouaf model (Abouaf et al., 1988;van NGUYEN et al., 2016) , including the diffusion creep mechanism, is implemented to predict the densification process. The strain rate for sintering porous materials can be written as following: where ( ) = (− / ) and are the Dorn's constant and the creep exponent of the fully dense material, respectively. is the activation energy of the particular creep mechanism, is the gas constant, and is the absolute temperature. ′ is the deviatoric stress tensor and is the Kronecker delta. is the equivalent stress of creep deformation for a porous body, which considers the volumetric changes under hydrostatic and sintering stresses: where 2 and 1 are the second invariant of the deviatoric Cauchy stress tensor and the first invariant of Cauchy stress tensor. , is the sintering stress, which has been developed by many groups and reviewed in (Shinagawa, 1999). The sintering stress is derived from the change in the surface energy or from the chemical potential at the sintering neck. In this study, the sintering stress is calculated according Eq. (8) (Shinagawa, 1999): where 0 is the initial radius of the grain. 0 and are the relative density of the green body and the actual relative density, respectively. The rheological parameters ( ) and ( ) in Eq. (9) and (10) are determined from results of sintering experiments presented in .
The constitutive equations, Eqs. (6)-(10), are implemented by means of the user-defined material (UMAT) subroutine. Temperature dependent thermophysical properties such as thermal conductivity, thermal expansion, specific heat and Young's modulus are taken from literature (Desai, 1986;Rabald, 1958;Rayne & Chandrasekhar, 1961). In addition, the Young's modulus was considered to be dependent on temperature and relative density. The density dependent values are calculated according the equations from (Beiss, 2013).

Kinetic Monte Carlo Model
Kinetic Monte Carlo (KMC) models are often used to predict the microstructural evolution during sintering (Dudek et al., 1998;Olevsky et al., 2006;Raether & Seifert, 2018) . Unlike physically motivated approaches, the KMC method aims for the phenomenological description of the driving forces of sintering and grain growth. During the sintering process, grain coarsening as well as an increase in circularity of the pores can be observed, both leading to a reduction of interface area and consequently to a decrease of the internal energy. The KMC method predicts the microstructure by random modification of an RVE model. Experimentally obtained micrographs of green bodies are then incorporated in a KMC model implemented in Matlab by first assigning a random integer number to each grain, which lies in the range of the number of the grains within the respective image. In a subsequent conversion step, each pixel is replaced by the number of the respective grain, which is referred to as a 'state'. In this way pores are described by = 0. The model is initialized by computing the interfacial energy between each lattice site and each neighboring site , as follows (Zhang et al., 2019): where ( , ) is the Kronecker Delta and , the interfacial energy between and , which is assumed to be 1.93 J/m² at the surface and 0.54 J/m² at a grain boundary (Prokofjev, 2017). The microstructural evolution is simulated by selecting a random site within the domain of [1, ], where corresponds to the number of sites that contribute to the total energy of the system. When a surface site is selected, its state is randomly exchanged with a pore site in a predefined neighbourhood of and then the resulting change of energy Δ is determined. Grain growth is modelled in a similar way. An exchange of state is accepted depending on its transition probability , which is proposed as follows: where is a temperature dependent variable accounting for the activation energy. A Monte Carlo step (MCS) consists of such simulation trials. A common approach to derive a relationship between an MCS and a real time scale is to compare the experimental results with the numerical results, such as the grain size. However, a physical relationship can be derived from the mean square displacement 2 performed within a MCS (Takahashi et al., 2006) with the lattice size and the surface diffusion coefficient which is assumed as 10 −10 m²/s as an average value of literature data (Blundell et al., 2005;German, 1996;Seebauer, 1995).

Heat treatment
The continuum mechanical modelling of the steel heat treatment process represents a multiphysics problem, involving coupling interactions between different physical fieldsthermal, mechanical, metallurgical and further relevant fields (Simsir et al., 2014). The model applied in the present work is based on the early suggestions of Inoue (Inoue et al., 1985) for a coupled modelling approach that comprises the calculation of heat transfer, phase transformations, transformation strains and the elastoplastic material response. Similar approaches, each with a different focus regarding the observed mechanisms, already exist in the literature (Diemar, 2007;Eser et al., 2016;Nusskern, 2013;Şimşir & Gür, 2008).To model the case hardening process, carbon diffusion during the gas carburizing is considered as the first simulation step. In conventional carburizing of sintered components with open porosity, the case hardening depth is increased by gas penetration through the pore network under atmospheric pressure. Nusskern proposed an empirical approach in (Nusskern, 2013), which defines the activation energy of diffusion as a function of the density. Accordingly, a density dependent diffusion coefficient is defined here to model the gas carburizing of surface densified gears.
During quenching and tempering, phase transformations are triggered by local temperature changes. To calculate the temperature change during heat treatment, the Fourier's heat conduction equation is solved by the ABAQUS solver. The required thermophysical material data are adopted from (Nusskern, 2013), where the effective heat transfer coefficient is defined as a function of temperature and porosity. The heat generation due to phase transformations is described in the user subroutine "HETVAL" according to (Eser et al., 2016). In (Nusskern, 2013), it is reported that both carbon content and density have an influence on the kinetics of the transformations during quenching. The latter dependency was not confirmed in the dilatometric investigations of this work. For a continuum mechanical analysis of the heat treatment, most quenching models in the literature describe the microstructure quantitatively by calculating the overall phase transformation kinetics and the volume fractions of phases. Modified formulations and extensions of the Koistinen-Marburger equation (Koistinen & Marburger, 1959) and the Johnson-Mehl-Avrami (JMA) equation (Avrami, 1939(Avrami, , 1940(Avrami, , 1941 are suggested for the martensitic and the diffusion controlled transformations, respectively. The transformation kinetic models and parameters are adopted from (Nusskern, 2013) to model the transformation of austenite (A) into ferrite (F) and martensite (M). In the case of the diffusioncontrolled transformation of austenite into bainite (B) an implicit differential form of the JMA equation is implemented as in (14), which is numerically integrated for the calculation of volume fractions. The transformation kinetic depends on the carbon content and the cooling rate ̇. The rate of the austenite transformation into bainite ̇ changes with the amount of bainite . The kinetic parameters ̅ and are the maximum volume fraction of bainite and the time coefficient, which are determined based on continuous quenching dilatometric experiments.
The overall effect of microstructure evolution during tempering, including the precipitation of transition carbides and the decrease of the martensite tetragonality, is modelled as a diffusioncontrolled transformation based on tempering experiments. The phase transformation model is implemented in the user subroutine "USDFLD" and the volume fractions of the corresponding phases are defined as field variables that are used to determine phase dependent material properties.
The main constitutive law in the mechanical analysis describes the evolution of the strains and assumes that the total strain rate ̇ is the sum of independent strain rate components: ̇=̇+̇+̇ℎ +̇ (15) where ̇ is the elastic, ̇ is the plastic, ̇ℎ is the thermal and ̇ is the transformation induced strain rate due to the volume change and transformation plasticity. Elastic and plastic strains can be calculated by the ABAQUS solver itself. The Young's modulus depends mostly on density and temperature. The yield stress shows a strong dependency on density, temperature, microstructural phases and carbon content. The corresponding values can be found in (Nusskern, 2013). The overall yield stress of the multiphase microstructure is obtained using a nonlinear empirical mixture rule suggested by Leblond (Leblond et al., 1989). Thermal and transformation induced strains are modelled base on dilatometric investigations describing the continuous cooling transformation behaviour of samples with different densities and carbon contents. In the user subroutine "UEXPAN", the thermal and transformation strain components are defined. The thermal strain is determined based on the phase dependent thermal expansion of the phase mixture, which is obtained by a linear rule of mixture. Phase specific expansion coefficients and are adopted from (Nusskern, 2013). The transformation strain component has two parts, and is given by Eq. (16). The first part is the local isotropic strain, induced due to density and therefore volume changes as the result of microstructural evolutions during quenching and tempering. The terms and are the phase specific transformation strains of the product and the parent phases at room temperature with respect to the reference phase austenite. This formulation includes the transformation of austenite during quenching and the transformation of martensite into the tempered martensite (M′) during tempering. Obviously, the specific transformation strain of austenite itself as the parent phase during quenching is equal to zero. The second part is the transformation plasticity, that has a decisive impact on the development of residual stresses. The model to describe this effect is adopted from (Fischer et al., 1996). Transformation plasticity is caused by the deviatoric stress components during the phase transformation and therefore is defined as an anisotropic strain, due to its directional dependency. The transformation plasticity factor was determined for the bainitic and martensitic transformations during quenching and for the tempering of martensite. The corresponding values were found to be 5•10 -5 MPa -1 , 8•10 -5 MPa -1 and 2.5•10 -6 MPa -1 , respectively. In Eq. (16), and ∆ are the volume fraction and the incremental change in the fraction of the product phase k. Figure 4 summarizes the model parameters in Eqs. (14) and (16), determined by dilatometric experiments.

Effective properties by RVE-modelling
Porosity is a source of heterogeneity in the microstructure, which ultimately affects the overall properties of the final product. In (Nusskern, 2013), effective thermophysical properties are defined as functions of density, mostly based on experimental investigations. Porosity has an impact on elastic properties, yield stress and thermal conductivity. In this work, an RVE modelling was performed to obtain the elastic properties depending on the porosity. Elastic properties have the highest impact on the residual stresses and the stress distribution under loading, which determines the highly loaded volume. The first step is to investigate the actual microstructure and to extract the geometrical features such as grain size and distribution in a form of statistical distribution plots. Next, the obtained data set is inserted into the software Neper (Quey et al., 2011) , which generates the synthetic RVE based on Voronoi tessellation technique. The periodicity of the geometry and its mesh are crucial features for simulating the behaviour of the surrounding material that consequently imposes a realistic reaction force on the boundaries of the RVE. In order to apply the periodic boundary conditions to the RVE, all the periodic nodes of the opposite faces must be identified to couple their displacement to the corresponding corner nodes. The coupling of the node displacements is done through a set of linear equations. For calculating each component of the orthotropic elastic stiffness tensor, six separate displacement loading conditions are applied to the RVE (three tensile tests in each principal axis and three simple shear tests). Numerical homogenization of the effective elastic properties is based on volume averaging method so that the global stress ̅ and strain ̅ values are used for the calculation of the equivalent homogeneous Youngs' modulus ̅ and Poisson's ratio ̅ .
The parameter in Eqs. (18) and (19) is the volume of the RVE. and are the calculated stress and strain values of the integration points.

3.5
Load bearing capacity The calculation of the tooth root bearing capacity is based on the universal approach developed in (Keusemann et al., 2015) as a calculation method for designing cyclic loaded sintered parts with complex geometries. This approach assumes that sintered parts fail under cyclic loading by crack initiation at larger pores in highly stressed volume and applies the Weibull's weakest link theory for volume defects. To develop a model regardless of the part geometry, fatigue tests were conducted on samples of different geometries and densities. Different loading modes (torsion, bending and axial loading) and stress ratios were investigated. To obtain the highly stressed volume and calculate the fatigue stress for a complex geometry, the first principal stress is evaluated in an elastic FE analysis. The criterion for fatigue failure is based on the maximum normal stress valid for a proportional stress condition. Hajeck introduced the effect of the carbon content and residual stresses into the model and extended the method to predict the tooth root loading capacity of sintered gears (Hajeck, 2017). In this work, the model is implemented with a minor modification by replacing the carbon content with the hardness, so that the possible effect of the cooling rate is included. With the modified model from (Hajeck, 2017), the maximum allowable stress amplitude , =0.05 is defined in Eq. (22) for a 50% survival probability under cyclic loading with a load ratio of 0.05. The parameters V ref , ρ 0 and HV 0 are the reference volume, density and hardness: 1 mm 3 , 7.85 g/cm 3 and 835 HV. The highly stressed volume V 90 is defined as the material volume subjected to 90% or more of the maximum principal stress. Density ρ ̅ and hardness HV ̅̅̅̅ are average values within the volumeV 90 . The model coefficients are summarized in (Hajeck, 2017).

Simulation of powder compaction and sintering
The DEM modelling technique was used to simulate the die filling and powder compaction process prior to sintering. Due to the irregularity of the shape of the water-atomized Astaloy 85 Mo powder, multi-sized particles in different size categories were created in the DEM model. The size distribution of the multi-sized particles was adjusted to represent the experimental values obtained from real powder particles in different size categories (Table 3). To calibrate and validate the DEM model, die filling and powder compaction experiments were carried out by filling 9.5 g of powder into a cylindrical die (Figure 6(a)) with a diameter of D=9.5 mm. The filled powder was then pressed uniaxially with a pressure of 400 MPa by two punches from the top and bottom. To analyse the density distribution in the sample, the powder particles must be "frozen" in their positions after filling and compaction. Therefore, the cylindrical sample was cold embedded, ground, and polished. By applying a self-developed "Image Analysis" method (van Nguyen et al., 2015), the micrographs of the samples were captured, stored, and sorted . The resolution and contrast of all the images were adjusted to distinguish between porosity and powder particles, i.e. black and white contrast respectively in Figure 6(c). The obtained images were binarized, and the relative density distribution was calculated by adding the number of white pixels and dividing by the total number of pixels in small grid windows.   Figure 6(a)) determined by DEM and "Image Analysis". Results show high-density regions in the centre of the sample after compaction, while relatively low densities were measured the positions near the die wall. The difference between the local density in the centre of the sample and that in the peripheral region is approximately 10%. The average densities of the powder compact, measured experimentally and simulated by DEM, are 89.1% and 88.9%, respectively, which is also in good agreement. The sinter model was also validated on the cylindrical sample by comparing the simulated sinter shrinkage with the experimentally measured shrinkage from the literature (Selecka & Simkulet). As shown in Figure 7, the simulated length change of the sample was plotted and compared with the experimentally measured length change from (Selecka & Simkulet). Under the given sintering condition, almost no densification was observed in the cylindrical sample, while the density distribution was evened out by grain boundary diffusion. The validated DEM and FE models were used to predict the density distribution of the sintered gear from (Hajeck, 2017). For the simulation, a 3D model was created, and the initial boundary conditions were defined. The geometry was reduced to a single tooth segment by defining symmetry planes as boundary conditions to reduce the element number and simulation time.
In the simulation of the die filling process, approximately 58000 single particles were used to image multi-sized particles. Figure 8 (a)-(c) shows the simulated die filling by the free fall of particles from a powder reservoir into the die. After the die filling was completed, the powder was compacted to a height of 14 mm, as shown in Figure 8 (c) and (d). According to the Voronoi method (Aurenhammer, 1991;Deng et al., 2017), the initial local density was calculated and imported into the FE sinter model for each predefined element. The process parameters of pressing and sintering were kept as in the simulation of the cylindrical sample. Figure 9(a) shows the simulated final geometry and density distribution of a tooth in the PM gear. The local densities changed less than 2% in the entire geometry due to the low sinter shrinkage. A vertical section through a single tooth shows, as illustrated in Figure 9 (b). This corresponds to the density distributions of the horizontal sections at different heights (0, 3.5, 7, 10.5 and 14 mm) of the tooth as shown in (Figure 9 (c)-(g)). The low-density area in the centre of the powder compact, i.e. at the neutral zone, agrees well with the area predicted by the DEM simulation of the powder compaction.

Prediction of microstructure by kinetic Monte Carlo simulation
To simulate the evolution of pore morphology and grain growth, the initial state of the microstructure was defined using a longitudinal micrograph of a compacted sample, which is shown in Figure 10(a). It can be noted that the powder particles are mostly aligned transversely to the direction of compaction, and each consists of several grains. Another sample was compacted under the same condition as the previous, then reduced in a hydrogen atmosphere and sintered under vacuum for 20 min at a temperature of 1120°C. Further micrographs were taken to evaluate the microstructural evolution during sintering.
Based on the micrographs of the compacted sample prior to sintering, an RVE-model was created to predict the microstructure development during sintering and the final microstructure by KMC method. The change of energy caused by a possible state exchange was calculated using Eq. (11), whereas it was accepted according to the function described in Eq. (12). The activation energy was assumed to be 2.5, which leads to a transition probability of less than one percent during the entire sintering process in the case of energy constancy. This takes into account the comparatively low sintering temperature. According to Eq. (13), the number of MCS that correspond to the sintering time is dependent on the mean square displacement, which was defined as the average displacement within a square region. Here, its width was selected as 13 pixels, which leads to 2 = 28.16 ̅ and consequently to an approximate number of MCS of 2900.
A micrograph of an unetched sintered sample and the predicted microstructure of the sintered body are shown in Figure 10. To distinguish between the grains simulated with the KMC model, each grain was randomly assigned a colour, while the pores are shown in black. A qualitative comparison between the micrograph of the green body Figure 10 (a) with the simulated and the sintered microstructure (Figure 10 (b) and (c)) shows the formation of the sintering necks and an increase in circularity. Considering the evolution of the grain size, a steady decrease of the growth rate during the sintering process is predicted, which is mainly due to the interaction with the pores. Pores that are sufficiently large, cause a pinning effect that limits the grain growth, whereas small pores are mostly isolated and located within the grains and have no significant impact on microstructural evolution. These observed effects are consistent with pore-grain-interactions reported in numerous studies (German, 1996). The final grain size was estimated to be three times of the initial grain size, whereas only a slight increase of 30 percent could be observed in the experiments. This discrepancy is mainly due to the neglection of the phase transformation. Therefore, quantitative information on grain growth is not assessible so far and the accuracy of this prediction could be further improved by considering the subsequent cooling stage. In addition to the grain size, the shape factor Eq. (23) is an important characteristic that is decisive for the microstructural evolution.
is the projected area and the perimeter of the pore (Beiss & Dalgic, 2001). The average shape factor ̅ is calculated by weighting the shape of each pore with its area. The numerically predicted value of 0.35 agrees well with the experimentally derived (pore) shape factors of 0.3 to 0.4, which were calculated on the basis of several micrographs. Using the shape factor as an example, bending fatigue strength can be described as a function of , as stated in (Beiss & Dalgic, 2001).

Determination of the effective Young's modulus
Based on the analysis of several micrographs with the software ImageJ the geometrical features of the pores were determined, and the initial overall size of the RVE was defined. Using the statistical data obtained, such as total porosity, pore size and distribution, the RVE can then be generated as shown in Figure 11. The synthetic RVE were used for simulation studies under various loading conditions, e.g. tensile or shear load. The local notch effect of the pores leads to a heterogeneous strain distribution. Based on local stress and strain data, the Young's modulus and Poisson's ratio was then calculated for each specific load direction by homogenization. The stress distribution shown in Figure 12(a), corresponds to a sample with a porosity of 12.8% under uniaxial tensile test along the X-axis up to 5% total strain. The matrix was considered as a homogenous continuum with a Young's modulus of 210 GPa and a Poisson's ratio of 0.3. Figure 12(b) clearly shows an almost linear correlation between the Young's modulus and the density. The numerical results are in good agreement with the experimental data from (Nusskern, 2013). It is worth noting that the morphology of the pores can have a considerable impact on the nonlinear force-displacement behaviour of the RVE and homogenized mechanical properties, as they cause stress concentration around and at the tips of the sharp pores. These effects will be investigated more thoroughly in later studies.

Hardness and residual stresses after heat treatment
For the simulation of the heat treatment process, the geometry of the sintering model was further reduced to the half width of the gear. Thermal boundary conditions were defined to model different stages of the heat treatment process. The carburizing temperature, time and the carbon potential were set to T=940°C, t=35 min, and Cp=0.85 wt.%, to model the gears investigated in (Hajeck, 2017) as the reference variant, Gref. To show the effect of the case hardening depth (CHD), a second variant, Gh, was modelled for a two-step carburizing process with Cp1=1 wt.% for t1=50 min and with Cp2=0.85 wt.% for t2=35 min. To describe the quenching process, a rapid cooling in an oil bath at 80°C and a subsequent air cooling to room temperature was modelled for both variants. The temperature dependent heat transfer coefficient was adopted from literature (Steinbacher & Hoffmann, 2009). In the last simulation step, a tempering stage at 200°C for 2 h was modelled. Based on the measurements in (Hajeck, 2017), the density profile in the surface area after cold rolling was defined.
As the functional surface is densified, carbon is transferred to the surface and diffusion starts from the surface to the bulk. However, side surfaces are not densified and thus, the carburizing gas penetrates into the pore network and leads to a higher carburizing and hardening depth. This effect is incorporated in the model by defining a density dependent diffusion coefficient. Figure 13 (a) shows the simulated volume fraction of martensite. The martensitic region qualitatively corresponds to the carburized area (not shown in the Figure 13). According to the simulation, a retained austenite fraction of approximately 6% is given on the surface after quenching.
In Figure 13 (b), the simulated hardness profiles along the path A-B in the tooth root are shown. In sintered parts, the macro-hardness is influenced by the porosity and cannot be used to evaluate the microstructural state. Therefore, the calculated hardness profile of the reference variant Gref is compared to the measured micro-hardness values. Micro-hardness measurements usually show a higher scatter, due to local microstructural inhomogeneities. In a phase mixture, measured values can vary between the hardness levels of the "hard" and the "soft" phases. The local hardness is calculated in the simulation according to Eq. 19 by a linear rule of mixture. In fact, the martensite and bainite hardness values, HVM and HVB, determine the upper and lower limits of the measured hardness. The comparison of the simulation and experiments shows a very good agreement regarding the CHD, i.e. the depth at which a hardness of 550 HV is given. However, the simulation model underestimates the hardness in higher depths. For a better estimation of the hardness profile, further investigations of the transformation kinetics should be performed for lower carbon contents, 0.25 wt.%<cc<0.5 wt.%. In the second variant, Gh, a CHD of approximately 0.5 mm is achieved. Figure 10 (c) and (d) compare the tangential residual stresses in the tooth root (A-B) and the tooth flank (C-D) depths. Simulation results show a very good quantitative agreement with the measured stresses for the reverence variant Gref. As expected for case hardened parts, compressive residual stresses arise in the carburized surface due to the transformation strains by the formation of martensite. The depth of the compressive stresses and the level of tensile stress increase with the CHD. Surface stresses are not affected by the CHD in tooth root, but in the flank surface; a decrease of the compressive residual stresses is caused by increasing the CHD. The load bearing capacity of the tooth flank is assessed by considering both surface and volume failure. Higher compressive residual stresses on the surface increase the resistance to pitting on the surface. On the other hand, shifting the maximum compressive stresses to the depth of the maximum Hertzian pressure due to the rolling contact can reduce the risk of volume failure. According to the results presented in Figure 13 (d), it can be concluded that the CHD can have a decisive impact on the flank load bearing capacity. The tooth root load bearing capacity could be less sensitive to the CHD, since the highly stressed volume is located in the vicinity of the tooth root surface, where compressive stresses are not noticeably changed. Figure 13: Results of the heat treatment simulation and validation, (a) martensite volume fraction, (b) hardness profile in the tooth root, tangential residual stress profiles in the tooth root (c) and in the tooth flank (d).

Calculation of the load bearing capacity
For the assessment of the fatigue life of the tooth root according to the approach of Hajeck (Hajeck, 2017) a finite element analysis is required to obtain the highly stressed volume V90. For this purpose, a 3D-model was created by reducing the geometry to a 3-teeth segment with the half width of the gear. Symmetry planes and force direction were defined based on the plane bending pulsator tests carried out in (Hajeck, 2017), Figure 14 (a). The tested gears were loaded with a minimum clamping force of 1 kN and the pulsating force was applied over four teeth. With linear elastic FE analysis, a highly stressed volume of V90 = 0.5924 mm 3 was obtained for this geometry. Since the Young's modulus depends on the density, the density profile of the cold rolled gear tooth was introduced into the model as obtained in (Hajeck, 2017). The mean values of density and hardness were determined according to the results of the heat treatment simulation. The maximum allowable stress amplitude for a load ratio of R=0.05 was calculated by Eq. (22). Considering the minimum clamping force of 1 kN and the maximum force amplitude of 3.36 kN for a 50% probability of failure, a load ratio of 0.13 was given. The mean stress = 480 MPa was obtained by the FE analysis. As shown in Figure 14 (c) and (d), the case-hardened gears exhibit compressive residual stresses in the surface zone. Tangential residual stresses at the surface as calculated in the heat treatment simulation, were used to correct the mean stress and the load ratio for both variants Gref and Gh. The calculated stress amplitude , =0.05 was finally converted to the allowable stress amplitude for the obtained load ratios using the Haigh diagram in (Hajeck, 2017). The calculated stress amplitudes and corrected load ratios R* are summarized in Table 4. Due to small changes in the surface hardness and residual stresses, an increase of only 2.2% of the stress amplitude was predicted by increasing the CHD from 0.28 mm to 0.5 mm. According to the simulation results it can be concluded that the tooth root load bearing capacity is almost insensitive to the CHD in relevant ranges. However, it should be noted that the calculation method proposed in (Hajeck, 2017) assumes a proportional stress state, which requires an initial residual stress-free condition. Furthermore, residual stresses are considered as scalar in tangential direction, neglecting a possible effect of the axial stresses. Therefore, a reasonable assessment of the load bearing capacity with the ICME approach requires the integration of the analytic calculation method in the FE simulation. The local equivalent stress must be obtained by a multiaxial failure hypothesis.

Potentials of the digital twin
The presented ICME approach can be used as a platform for the digitalizing of process design for the realizing of an agile production concept. By utilizing mechanism-based models, the material response for different process conditions in the manufacturing chain can be described. The model considers the diffusion process during sintering, defines temperature and density dependent material properties and predicts phase transformations and stress evolution during the heat treatment process. The simulation tool is capable of quantitatively describing the correlation between manufacturing process, microstructure and mechanical properties of the component. Figure 15: Schematic representation of the development cycle including the ICME approach as a digital design tool Figure 15 shows schematically the vision for the integration of the ICME approach in the product development cycle. As illustrated in Figure 15, the greatest benefit of the ICME approach is the significant reduction of experimental investigations, which were previously required for process optimization mainly through "trial and error" to improve the component performance. Component testing can be replaced by laboratory material testing on small samples to determine material properties and model parameters. Furthermore, the simulation can provide results that that are not or only with difficulty obtained by experiments, such as a comprehensive measurement of the spatial residual stress state in the entire component.

Summary and outlook
The ICME approach presented in this paper covers most of the relevant process steps in the powder metallurgical manufacturing of sintered gears, including die filling, powder compaction, sintering, and heat treatment. The simulation tool integrates the digital shadow of the process chain and the digital twin of the material by comprehensive material modelling on meso-and macroscale to reflect characteristics of the real manufacturing process. The digital twin allows the prediction of key properties and parameters to determine the gear performance on the basis of the load bearing capacity of the tooth root. To combine the advantages of macroscopic and mesoscopic approaches, the DEM, KMC and homogenization methods were used at the mesoscopic scale to predict the pore morphology and the elastic properties of the porous material. On the component scale, the continuum modelling of sintering and heat treatment have achieved the aim of accurately predicting the density distribution, hardness and residual stresses by integrating results of mesoscopic simulations.
The simulation methods used and results obtained are summarized as follows:  The simulation with coupled DEM and FE modelling was used to calculate the powder density distribution after the die filling and powder compaction process and to predict the densification during the sintering process. Multi-sized particle models were implemented in the DEM modelling to consider the influences of powder morphology and size distribution on the initial density distribution after die filling and compaction. The simulation results of the density distribution in the cylindrical sample correspond well with the experimental measurements, which demonstrates the effectiveness of this approach in the sintering simulation.  The KMC model provides insights into the microstructural evolution during sintering.
Although the modelled pore morphology and the shape factor agree well with experimental data, further investigations are necessary to analyze the influence of sintering temperature in more detail and to derive a functional description of the dependence of the activation energy on the sintering temperature. In addition, an accurate measurement of the surface diffusion coefficient is necessary to correctly estimate the number of MCS. Furthermore, extending the model to a continuous description of the entire heat treatment process, which includes phase transformation during cooling, would provide the possibility to fully describe the microstructural evolution.  To determine the influence of porosity in sintered parts on the mechanical properties of the material, the real microstructure was analysed regarding shape and distribution of the pores. Based on the collected statistical microstructural properties, a synthetic RVE was created. The impact of the porosity on effective elastic properties was modelled with an FE-based homogenization technique. The predicted Young's modulus agrees well with experimental findings from literature. The influence of porosity on other effective properties such as heat conductivity and nonlinear mechanical properties is being investigated in ongoing studies.  The FE modelling of the heat treatment process covers the thermomechanical-metallurgical interactions. Carbon diffusion is modelled considering the effect of porosity during gas carburizing. Depending on the carbon content, local phase transformations and the associated volume change due to quenching are described to calculate the development of hardness and internal stresses. The tempering stage is modelled to predict the residual stresses after relaxation of hardening stresses. The predicted hardness and tangential residual stresses are in good agreement with experimental results from literature. The simulation model predicts compressive residual stresses in the surface zone and tensile residual stresses in the bulk. Predicted stress patterns indicate that the depth at which compressive stresses exists and the level of maximum tensile stresses increase with CHD. Surface stresses are affected by the CHD only in the tooth flank.  Hardness and residual stresses are input parameters for the analytical assessment of the tooth root load bearing capacity. The highly stressed volume is obtained by an elastic FE analysis, according to the plane bending pulsator tests from literature. The predicted residual stress state is considered as a scalar value in tangential direction and used to correct the applied mean stress. Due to small changes in the surface residual stresses, no significant influence on the tooth root load bearing capacity is expected by increasing CHD. The model is based on the assumption, that the subsurface porous material does not fail before the densified surface. The analytic calculation method cannot be used to predict the tooth flank load bearing capacity. The reason is that the basic assumption of proportional stresses under cyclic loading, is not valid for rolling contact conditions and the presence of considerable residual stresses. The fu-ture work would be to fully integrate the fatigue assessment into the continuum mechanical approach as a subsequent step of the heat treatment simulation. The numerical assessment of the fatigue life should be performed on the basis of valid failure hypotheses for non-proportional stress states for both the tooth root and tooth flank.

Availability of data and materials
Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.