Progress of Physics-based Mean-field Modeling and Simulation of Steel

The progress of mean-field modeling and simulation in steel is presented. In the modeling, the focus is put on the development and application of a physical modeling base, including Calphad, diffusion assessment, nucleation and growth of precipitates, and dislocation dynamics. This leads to an improved prediction of the materials response after different thermo-mechanical treatments in terms of microstructure evolution and mechanical properties. The presented case studies represent the success of the integrated computational materials engineering approach.


Introduction
Materials research has evolved from the early days' phenomenological characterizations of microstructures, driven by extensive trial and error experimental work, towards nano-to mesoscale analysis of phase stabilities and the in-depth understanding of physico-chemical interrelations among phases and of the evolution of different microstructures. This has been accompanied by the formulation of physics-based materials models with integration of, e.g. thermodynamic, crystal structural, phase stability, diffusion data, and interfacial energies.

From Calphad to Simulation of Materials
Properties-an Integrated Modeling Base The Core of Calphad Is Constituted by 1. physically appropriate phase descriptions, including crystallographic structure, stoichiometry, and stability range of a phase, coupled with the thermodynamic state variables that define the temperature-and compositiondependent molar Gibbs energy of a phase, 2. the optimization of model parameters taking into account selective, weighted/assessed input data from all available experimental data that are directly or indirectly connected to the model, i.e. thermodynamic properties and phase diagram data, 3. the extension from low-order systems to high-order systems, with (optional) non-ideal mixing added to the "frozen" combinatorial base of all low-order descrip-tions. In case of quaternary and higher-order phases occurring, their individual thermodynamic descriptions are further required.
One key towards Applied Calphad for multi-component and multi-phase alloys is an iterative validation process of 1) and 2), which results in the refinement of 3) and the associated improvement of the Calphad-based modeling application, such as the prediction of thermokinetic phase evolution, precipitation strengthening, or recrystallisation. Inspired by this powerful route, the Calphad approach is extended beyond the conventional, pure thermodynamics use in that sense that the "highest-order" system is represented by the integrated, predictive kinetic modeling base. Consequently, kinetic experimental results, such as precipitate sizes and compositions after different heat-treatments, are considered in this Calphad-like iterative bottomup and top-down validation with simultaneous assessment of model parameters. Fig. 1 shows a (partial) scheme of model integrations into a computational materials simulation package.
A predictive simulation of materials behaviour requires a physically appropriate integrated modeling base, consisting of thermodynamics, kinetics consisting of diffusion mobilities, nucleation and growth modeling, phase transformations modeling, and modeling of dislocation dynamics. We believe that different modeling activities, as shown by the representative examples as follows, constitute puzzle pieces towards the realisation of a fully integrated, predictive computational materials design. We show current modeling approaches of bainitic transformation and microstructures, dislocation evolution, and recrystallisation in microalloyed steel including carbo-nitrides (MX) precipitation kinetics, bake hardening simulation, and multiscale modeling of micro-mechanics.

Modeling
Modeling Aim Numerous models on the kinetics of bainite transformation have been proposed based on the displacive theory, which we accept here [1][2][3][4][5][6][7]. Most of the current displacive models for bainite transformation kinetics apply the following concepts initially developed by Bhadeshia [8,9]. A fundamental question for the bainite transformation kinetics concerns the completeness of reaction. When carbide precipitation is suppressed, as well known for high-Si steels, the bainitic transformation can be incomplete and carbon enrichment is enhanced in the surrounding austenite, strongly affecting the transformation kinetics. Thermodynamics of the bainitic transformation have not been rigorously used in previous modeling. Namely, the change of critical temperature T0', where the Gibbs energy of austenite Gγ M equals that of ferrite Gα M plus a certain amount of strain energy, has not been taken into account. T0' rapidly decreases with the carbon enrichment in austenite and determines whether the growth of bainite can proceed or whether it will come to a stop. This further means that T0' evolves indirectly with the volume fraction of bainite, f. Moreover, T0' and the isothermal transformation temperature T are not independent from each other. T must be lower than T0', which is actually the growth criterion of bainite. The most complete and physical description of the bainitic transformation to date [6] includes a rigorous description for the autocatalytic nucleation and the partitioning of carbon between austenite and bainite. However, the carbon concentration of bainite, Xb, in that work has been a constant which is not sufficiently precise. Actually, the formulation of the evolution equation of Xb involves answering the most challenging and difficult Modeling Technique In the work of Wei et al. [10], an evolution equation of Xb with two adjustable parameters, α and β, has been chosen, which captures the main mechanisms of carbon redistribution, trapping of carbon in bainite and carbon enrichment in film-like austenite. Carbon atoms can be trapped inside of ferritic subunits or within film-like austenite between ferritic subunits. In the modeling, this is described by the parameter α, which is related to the probability of carbon partitioning from ferrite into austenite. The scheme in Fig. 2 shows two sheaves, one of which already ejected all possible carbon into the austenite and one that still holds carbon trapped in bainitic ferrite and in austenite films. The equilibration time necessary for a "carbonsupersaturated" sheaf depends on transformation temperature and morphology of the bainite sheaf. Based on the experimental observation for the evolution of the carbon content of bainitic ferrite with time, an empirical equation to capture the evolving feature of Xb (Eq. 1) can be formulated, With fitted α values below 0.5 and β in the range from 0 to 1 (0≤ β < 1), we observed good agreement of the predicted bainite evolution with experiments, as shown in the following. The parameter β acts as the theoretical potential of the final carbon enrichment of austenite caused by the transformation of bainite. β is dependent on several materials depending factors, such as, the steel composition or the configuration of bainite subunits.

Results and Discussion
To date, predictions of bainite transformation kinetics have remained problematic in particular for high-Si steels at high C alloying around 0.4 wt.%. This problem has The consistency between calculated and experimental kinetics for steel grades with varying C-contents at high Sialloying [10] indicate the soundness of our assumption that the evolution of Xb with time can be described by the semiempirical α-β equation.
The correct temperature dependence of the simulated bainitic transformation is strongly governed by α, which needs to increase with temperature. This trend is in line with increasing entropy at higher temperatures, thus facilitated carbon partitioning. Further, the strong sensitivity of the parameter α to temperature indicates that α is heavily related to the carbon diffusion process in bainite.

Conclusions and Outlook
The proposed semi-physical α-β model is capable of capturing the dynamic evolution of carbon originating from bainite sheaves and its enrichment in austenite, considering the evolution of Xb. Our semi-physical α-β modeling approach contributes to the solution of the important and long-term remaining question of how the carbon atoms redistribute during bainite transformation. Current research is conducted to understand the exact physical background of α and β, aiming on the realisation of an explicit formulism for direct determination of α and β.

Modeling
Modeling Aim A prediction of the bainitic ferritic subunit size, the thickness of the residual austenite films, their cor-responding C-enrichment and the accompanying stabilization of the residual austenite is obtained.

Modeling Technique
The austenite film thickness is controlled by the diffusion of C around a bainitic ferrite platelet into the surrounding residual austenite [11]. In our modeling approach, we combine the T0-temperature concept with the numerical simulation of C-diffusion profiles, utilizing the 1D cell diffusion module of the thermokinetic software package MatCalc. Fig. 4 schematically shows the evolution of the C profile with time for one subunit in contact with the austenite matrix. This methodology gives the opportunity to predict the C-distribution under consideration of consecutively forming subunits, which is necessary to estimate the C-content of austenite films. The simulations also take into account APT measurements from ref [13] the effect of C trapping at the dislocations formed inside the ferritic platelets due to plastic deformation. Different geometrical configurations with a number of bainitic subunits can be studied, and, finally, each cell can perform a fully coupled precipitation kinetics calculation to evaluate the simultaneous long-range diffusion of C and precipitation kinetics of carbides inside the ferritic platelet and the surrounding austenite. Fig. 5 shows the simulation results of the described approach and outlines the importance of C trapping at low bainite transformation temperatures. More details concerning this approach and further results can be found in [12]. The locations of the critical carbon content (horizontal dashed line), di, correspond to the simulation times, ti. While the predicted film thickness increases from t0 to t3, it decreases afterwards and eventually reaches zero again after longer simulation times.

Modeling
Modeling Aim The strain dependent dislocation density evolution represents a key factor for recrystallization during thermo-mechanical treatments as deformation induced dislocations mainly control the driving pressure for recrystallization. In this contribution, the dislocation density evolution and the recrystallization behaviour of microalloyed steel are experimentally examined via single-hit and double-hit compression tests, respectively. The obtained results are subsequently reproduced using the simulation software MatCalc. A new set of physics-based, experimentally determined strengthening parameters is introduced for the dislocation density evolution. Additionally, the applicability of the parameters for recrystallization simulation of microalloyed steel grades is assured.

Modeling Base-Dislocation Density Evolution
Stressstrain (σ-φ) curves of face-centered cubic metals can be divided into five different strain hardening stages. This work is focused on step III, which is strongly dependent on the dislocation density development. At this stage, the strain hardening rate θ = dσ / dφ decreases and the rising yield stress σ can be separated into an initial yield stress σ0 as well as a deformation dependent plastic stress σP [14]. The initial strain hardening rate θ0 and the saturation stress σ∞ can be determined via the corresponding Kocksplot, which depicts θ as a function of σP [15]. For the strain dependent dislocation density evolution modeling, an extended Kocks-Mecking model is employed, which comprises interactions at elevated temperatures [16]: where ρ describes the dislocation density, ϕ the true strain, M the Taylor factor, b the Burgers vector, A the free travelling path of dislocations until their movement is impeded by immobile dislocations, B the dislocation annihilation coefficient for dynamic recovery, dcrit the critical annihilation distance between two dislocations, C the dislocation annihilation parameter for static recovery, Dd the diffusion coefficient, G the shear modulus,φ the strain rate, k the Boltzmann constant, T the temperature, and ρeq the initial dislocation density. The experimental determination of the physical parameters σ0, θ0, and σ∞ is essential for the evaluation of the strengthening parameters A, B, and C, the applied model is based on [16]: Modeling Technique Single-hit compression tests are performed at varying test temperatures (913-1473 K) and strain rates (0.01-1 s -1 ) for five different microalloyed steel grades. Using the attained flow curves, the needed input variables σ0, θ0, and σ∞ for the calculation of the coefficients A, B, and C are identified. The initial yield stress is individually defined at Rp0.2, while the initial strain hardening rate and the saturation stress are determined via the Kocks-Plots for every recorded stress-strain curve [15].

Results and Discussion
Flow Curves The determined individual strengthening parameters at varying temperatures and strain rates are applied for the dislocation density evolution model implemented in MatCalc. The experimental strain dependent plastic stress σP is consistent with the simulative dislocation yield stress contribution using the strengthening parameters as determined for the specific material, temperatures and strain rates (Fig. 6). Flow curve modeling with an averaged ABC set resulted in slight deviations, but a sufficient accuracy is achieved for a general estimation of the needed input parameters.
Recrystallization To ensure the applicability of the new averaged ABC parameter set to recrystallization modeling of alloys with different compositions, recrystallization calculations are performed and compared to the experimental results from double-hit compression tests. Thermo-mechanical tests of a Nb-and V-microalloyed steel grade are carried out at three different test temperatures (1273-1373 K). The specimens are deformed twice to 0.2 strain at a time with 0.5 s -1 strain rate and varying holding times between the compressions. The recrystallized fraction is subsequently evaluated using the 5% true strain method [17,18].
The applied MatCalc recrystallization simulation model is based on the research by Buken et al. [19,20], combining simultaneous precipitation pinning by Calphad-described fcc-structured (Nb,V)(C,N) carbo-nitrides (MX) and solute drag with the impact of dislocation density evolution. Fig. 7 demonstrates that the experimental recrystallized fractions of the examined steel grade can be successfully reproduced via simulation using the averaged parameters of A, B, and C for the dislocation density evolution. At temperatures be- Fig. 7: Effect of deformation temperatures on recrystallized fractions of Nb/V-microalloyed steel grade derived from experimental (symbols) as well as simulated (lines) isothermal compression tests with 0.5 s -1 strain rate low 1373 K, both methods agree that the strain induced formation of MX leads to retardation of static recrystallization and a plateau is formed.

Conclusions
The evaluated physics-based parameters for the dislocation density evolution in combination with precipitation kinetics allow for an accurate modeling of flow curves and recrystallization kinetics of microalloyed steel.

Modeling of Bake Hardening Kinetics and
Carbon Redistribution in Dual Phase Steel DPW 600

Aim of Simulations
The bake hardening response of the dual phase (DP) steel DPW600 with a martensite volume fraction of 0.2 is studied by means of computational simulation based on experimental austenitization and quenching, followed by isothermal treatment at 100 and 220°C for times up to 10000 min. The increase of yield strength during bake hardening follows the well-known two-stage characteristics. The kinetics of bake hardening is analyzed on basis of computational models for the formation of Cottrell atmospheres and precipitation, which represent the cause of strength increase, as well as the long-range diffusion of C, which is partitioning between the martensite and ferrite phases during aging. The model describes the shape of the experimentally observed bake hardening curves as a function of temperature, time, and degree of pre-strain [21].

Modeling Base
The dislocation density for pre-strained samples is evaluated with an extended Kocks-Mecking model [15], which considers the deformation-induced dislocation generation as well as dynamic and static recovery. The redistribution of C between martensite and ferrite, as well as across the phase interface is accounted for in terms of the multicomponent diffusion kinetics module of MatCalc and is shown in Fig. 8a. We use the Cottrell atmosphere formation model [22], Classical Nucleation Theory [23], and the SFFK (Svoboda-Fratzl-Fischer-Kozeschnik) growth model [24] to simulate the precipitation processes that accompany the bake hardening step. In the simulation, both, Cottrell atmosphere formation as well as carbide precipitation, are solved simultaneously on basis of the same set of state parameters. Accordingly, they are fully thermodynamically interconnected. To predict the strength increase during bake hardening, a linear relationship to describe the strength increase in low-alloyed steel based on an Ashby-Orowan relationship is used [25].

Results and Discussion
The observed two-step increase of the bake hardening response as function of time is reproduced by the present model, as seen in Fig. 8b. Our simulations suggest that the decrease of the bake hardening response after long-term annealing is due to the depletion of C in ferrite that leads to a decomposition of existing Cottrell atmospheres and carbides and an accompanying reduction in yield strength.

Conclusions
The present work represents a fully coupled simulation framework for the evolution of all relevant state parameters, such as the size and distribution of carbide precipitates, the degree of C segregation to dislocations, and the amount of dissolved elements in the DP matrix. In addition, due to long-range diffusion of C from ferrite into martensite, the characteristic decrease of the bake hardening response after long-term annealing is explained.

Modeling
Modeling Aim The effect of a material's constituting phases on its macroscopic elastic properties can be described using continuum micromechanics modeling. The complexity of the steel microstructure makes the question of transition between scales critical and the use of a multiscale homogenization approach necessary, extending Hellmich's group work [26,27] on biological or construction materials to polycrystalline metals. This approach is tested on hypereutectoid steel constituted of pearlite colonies and grain boundary proeutectoid cementite. Pearlite can be either lamellar (LP) or globular (GP), and grain boundary cementite may either form a continuous (CC) or fragmented (FC) network.

Modeling Technique
We define a statistically representative volume element (RVE) [28] following the rule of separation of scales. The homogenized stiffness tensor C hom [29] is basically defined by the phase fractions and stiffness tensors of the constituting phases. Moreover, the implemented Hill tensor accounts for the shape of the inclusions formed by the phase p within a matrix described by the stiffness tensor C 0 . More details of the modeling can be found in our recent publication [30]. The homogenization of an individual pearlite colony is first performed from the properties of ferrite and cementite, followed by the homogenization of the macroscopic material from the properties of pearlite and cementite. The classical self-consistent scheme (CSCS) [31]-suitable for aggregates and polycrystals-and Mori-Tanaka (MT) [32,33]-adapted to matrix-inclusion composites-are respectively used for the first homogenization step in the LP and GP case. The generalized self-consistent scheme (GSCS) [34]-suitable for aggregates of isotropic inclusions surrounded by a homogeneous film-and CSCS are respectively used in the second homogenization step in the CC and FC case.

Results and Outlook
Our model delivers values of the stiffness tensor and of the different elastic moduli similar to experimental data for pearlitic steels ( [35] ; Fig. 9). The influence of microstructure parameters such as grain boundary cementite thickness or pearlite colony size is small. The difference between LP/GP and CC/FC is also small due to similar elastic properties of Fig. 9: Calculated Young's modulus as a function of grain boundary cementite thickness from the models with LP and CC or FC at grain boundaries, with material input data for ferrite [36] and cementite [37] ferrite and cementite. Moreover, this result agrees with the fact that the pearlite and proeutectoid cementite morphologies influence mainly the material's plasticity and toughness. This directly leads to our current research focus, aiming on the implementation of simulated elastic behaviour within an elasto-plastic modeling routine.

Summary and Outlook
The presented integrated modeling approach in the meanfield is universal and-in combination with careful model validations based on key experiments-predictive allowing for an in-depth understanding of phase evolutions and microstructural trends and their consequences for macroscopic materials properties with relatively manageable computational time and costs.
Various examples of coupling of materials at different scales are found in technological applications. Simulative answers to questions concerning efficiency, sustainability and degradation of multi-material technological systems need to consider, e.g. effects of thermal cycling on phase stabilities and microstructures and of high-temperature oxidation of alloy parts. For the time being, we are still in the fledgling stages on the way to predictive computational materials engineering. A major research task refers to the assessment of the role of local chemical and microstructural heterogeneities as well as defects and interfaces and their dynamics for materials properties and materials behaviour during thermo-mechanical treatment.
The providential combination of ever increased resolution of characterization techniques and computational power gives us a great possibility for predictive, physicsbased materials simulations of complex, multi-component, multi-phase materials.
Technology and Development and the Christian Doppler Research Association is gratefully acknowledged.
Funding. Open access funding provided by TU Wien (TUW).
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://creativecommons.org/licenses/by/4.0/.