State Parameter-Based Yield Strength Model for Integration in Finite Element User-Material Routines

: A new state parameter-based user-subroutine for finite-element software packages, which can be used to simulate microstructure-dependent stress–strain relations, is presented. Well-established precipitation kinetics, strain hardening and strengthening models are brought into a condensed form to optimise computational efﬁciency, without losing their predictive capabilities. The framework includes main strengthening mechanisms, such as, precipitation strengthening, solid solution strengthening, the cross-core diffusion effect and work hardening. With the novel user-subroutine, the microstructure evolution of various thermo-mechanical treatments on the full integration point grid of the ﬁnite element (FE) mesh can be calculated. The validation of the simulations is carried out by mechanical testing as well as microstructure characterisation of an Al-6082 alloy, including transmission electron microscopy (TEM) investigations after various annealing times at 180 ◦ C.


Introduction
The knowledge of residual stresses and distortions, which can be introduced into a workpiece during the manufacturing process, is essential for the optimisation of both, the design of the components and the processing parameters. Since plasticity is dependent on the thermo-mechanical history of the material, the entire process chain, starting with solidification, must be taken into consideration for the simulation of temperature and strain rate-dependent stress evolution. An accurate material model is mandatory to relate flow stress to plastic strain. Many constitutive models were developed in the past decades, such as the Ludwik approach [1], the Voce type approach [2], the Johnson Cook model [3], Zerilli-Armstrong model [4], or the model of Khan and Huang [5]. Empirical material models or data table-based methods, which often describe experimental results but with limited physical meanings, are standard nowadays and they are included in most finite element analyses (FEA) tools. In contrast, advanced microstructure evolution models, which are, for instance, implemented in the MatCalc software package (http://matcalc.at, accessed on 12 July 2022), would lead to very long computing times and consume extensive computer memory, when included in FEA tools. FEA material models are presented in [6,7], which include strain hardening, recovery and recrystallisation, but do not consider precipitation kinetics or the evolution of excess vacancies. In the present paper, a mechanism-and state parameter-based model is introduced, which allows simulating the microstructure evolution in a metallic material during the manufacturing process and in-service conditions. In addition, the model allows for an extrapolation of these operating conditions, such as strain rate and temperature, while maintaining stable convergence behaviour. Although the computation time can be up to twice as long as with typical standard models, a comprehensive microstructure description is provided, which gives insight into the underlying metallurgical processes that occur during the thermo-mechanical treatment.
A new user-subroutine, suitable for incorporation in commercial FEA-software packages, is presented in the first part of this work, which is called 'simple MicroStructure Evolution' (sMSE) model in the following. The sMSE model is based on the mechanical threshold concept and an extended Kocks-Mecking approach, which allows to calculate the evolution of the average dislocation density during thermo-mechanical treatments. With the temperature and strain rate profiles given from the finite element software, the sMSE materials subroutine returns local yield stresses as well as the corresponding derivatives with respect to strain and strain rate back to the finite element software.
The second part of the paper illustrates the application of this new framework to the simulation of flow curves in a precipitation-hardenable Al-6082 alloy in dependence of the particular heat treatment condition. For this task, efficient models of precipitate nucleation and growth are included, which are an important component of the microstructure evolution beside the average dislocation density evolution. The work is, finally, experimentally validated by compression tests and electron microscopy investigations. The new framework conveniently balances low calculation time, which allows to simulate residual stresses of complex components, such as cylinder heads, and the evolution of physically based material properties on the full integration point grid of the FE-mesh.

The Model
The strength model consists of an athermal stress contribution σ ath (see Section 2.1) and a thermally activated stress contribution σ th (see Section 2.2). σ ath is caused by forest dislocations and is expressed as an average dislocation density as formulated in an extended Kocks-Mecking approach [8]. σ th includes strengthening mechanisms, such as, solid solution strengthening, precipitation strengthening, the cross-core diffusion effect and grain boundary strengthening in the mechanical threshold framework. The total stress contribution σ is given as a function of temperature T, strain rate . ε and the vector of state parameter-based coefficients χ, which can include the concentration of the alloying elements in the matrix, the radius and number density of precipitates or the dislocation density, for instance.
Adaptions of original models from literature are introduced for a suitable integration into the user-subroutine in FE simulations. Grain growth and recrystallisation effects are not included, because precipitation strengthening, and work hardening are the dominant mechanisms in the present experimental set up. Since a 6082 Al alloy is investigated, the model is introduced for an Al matrix with the two main alloying elements Mg and Si. Extension to higher-order systems is straightforward but not explicitly elaborated, here.

Athermal Stress Contribution
An extended one-parameter model of Kocks and Mecking is used to calculate the temperature and strain rate-dependent dislocation density evolution [8,9] M is the Taylor factor, b is the Burger's vector, k B is the Boltzmann constant, G is the shear modulus, A, B and C are temperature and strain rate-dependent coefficients. Kreyca et al. [8] related the coefficients A, B and C to the initial hardening rate θ 0 and the saturation stress σ ∞ . d crit is the critical annihilation distance between two dislocations of opposite sign of the Burgers' vector [10] and ρ eq is the equilibrium dislocation density.
dρ + dε represents the dislocation generation rate, which is inversely proportional to the current mean distance between dislocations and the total dislocation density, ρ. dρ − dε and dρ − s dε take into account the annihilation of dislocations by cross slip processes and vacancy-assisted climb, respectively. The first annihilation term describes dynamic recovery processes at low and medium temperatures, whereas the second term represents static recovery at elevated temperatures. The latter quantity is marked by the subscript s. For each time integration interval ∆t, the average dislocation density ρ is correlated with the athermal plastic stress contribution by the Taylor equation [11,12] with where α is the strengthening coefficient. Further details on the implementation of the athermal strength framework are described in [13].

Thermal Stress Contributions
Many state parameter-based yield strength models capture the influence of temperature and strain rate on plastic deformation of polycrystalline materials in various ways [14][15][16][17][18][19][20][21][22][23]. One example is given in [8], where a low temperature and a high temperature part can be distinguished, labelled by 'lt' and 'ht' in the following. When dislocation motion is characterised by glide processes, the stress can be expressed as where the subscript 'σ 0 ' refers to the initial yield stress, . ε 0 is a constant andτ is the mechanical threshold, defined by the sum of a basic stress, solid solution hardening (see Section 2.2.1), cross core diffusion hardening (Section 2.2.2), grain size hardening, sub-grain size hardening and precipitation hardening (Section 2.2.3) in the absence of thermal activation.
When dislocation climb becomes dominant at high temperatures, the following stress contribution can be used [8] c is the speed of sound, and the exponent n of the power law equation varies between 3 and 10 [24]. The true strain rate is modified in the present framework by the exponent n . ε to reduce the strain rate dependency with ε . The activation energies ∆F lt σ 0 and ∆F ht σ 0 depend on the effective solute concentrations in the matrix [25][26][27][28], which vary due to the cross-core diffusion effect and by the nucleation and growth of precipitates (see Section 2.3). The total thermal stress σ 0 is evaluated by the summation rule according to where n c is a coupling coefficient. The overall temperature-and strain rate-dependent stress response at an applied constant strain rate . ε is given by the summation of the athermal stress contribution σ p of Equation (3) and the thermal stress contribution σ 0 of Equation (6).

Solid Solution Hardening
To incorporate the solid solution strengthening contribution into the present framework, the following equation is utilized, based on the Labusch approach [29] The subscript 'SS' refers to solid solution and the sum is taken over all alloying elements, i. The contribution of the alloying elements with the concentrations c i are calculated separately and summed up, using the exponent n SS . k i is defined by f max is the maximum interaction force between the solutes and the dislocations, which is defined by [29] f max = w = 5b is the interaction distance and the dislocation line tension E L = 1 2 ·G·b 2 [30]. ν is the Poisson's ratio and ε m is the misfit strain between solute and matrix atoms.

Cross-Core Diffusion Hardening
Dynamic strain aging (DSA) phenomena occur by the interaction of diffusing solutes with the stress fields of dislocations and retard the ongoing dislocation motion. Since the possible negative strain rate sensitivity can lead to plastic material instabilities due to local material softening, it plays a key role in material processing [31]. A well-known example for DSA is the Portevin-LeChatelier effect, which manifests itself in serrated stress-strain behaviour. In the sMSE framework, the cross-core diffusion effect is included as developed by Curtin et al. [31]. The model is based on single atomic jumps of solutes from the compression side to the tension side in the core of a dislocation, leading to the additional strain rate-dependent strengthening [31] as α = 0.56, c 0 is the bulk solute concentration in each iteration step, ∆W is the average binding energy difference between the core compression and tension sites, Ω is a constant in this framework, Γ c is the reference core transition rate with where ν 0 is the attempt frequency and ∆H c is the average activation enthalpy for transitions from tension to compression and vice versa. Generally, higher temperatures and smaller strain rates lead to more diffusion and, consequently, to higher cross-core diffusion strengthening. The effective concentration of solutes on the tension side of a dislocation core is calculated as [31] Continuous changes of the effective concentration c eff influences the strengthening contribution ∆τ s , which is linearly superimposed to the mechanical threshold. The crosscore diffusion effect is directly included to the sMSE framework without any simplification.

Precipitation Hardening
Obstacles within a material, such as precipitates, can hinder dislocation motion and thus lead to an increased strength of the material. In general, small and coherent precipitates can be sheared by dislocations, whereas, above a certain critical size, a transition to bypassing and formation of dislocation loops around the precipitates occurs. In the following, the stress contribution based on the Orowan dislocation looping mechanism [32], as well as the shearing mechanism is outlined and the equations, which are implemented in the sMSE framework, are presented. The mechanism, which delivers the least contribution to the total stress, is assumed to be the operative one and it is labelled σ prec in the following. For simplification, no distinction between weak and strong precipitates is included, as described in detail, for instance, in [33].
The interaction of dislocations and non-shearable precipitates was first described by Orowan [32] and later modified by Ashby [34], Brown and Ham [35], and Ardell [36]. The implemented equations are taken from Ahmadi [33], based on the original Orowan model. The subscript 'O' is used subsequently and the Orowan strength is taken as where C O is the precipitation strengthening coefficient and L S is the mean distance between two equally sized spherical precipitate surfaces with [37] N is the number density and r is the radius of the precipitate; the evolution equations are given in Section 2.3. The equivalent radius R eq in Equation (13) describes the precipitatedislocation interference area with R eq = π 4 ·r (15) To avoid negative values for the Orowan contribution, the equivalent radius has a minimum value of 4·b.
The shearing mechanisms involve the coherency effect, the modulus effect, the stacking fault effect and the interfacial effect. For simplicity, only the coherency effect, which provides the largest strengthening contribution in many practical applications, is included in the present framework. The coherency effect is based on the interaction of the elastic strain field, which is caused by the difference between the molar volumes of matrix and the precipitate, with a moving dislocation [33]. The stress contribution for coherency strengthening is given as where θ is the angle between dislocation line and its Burgers vector. θ = 0 for pure screw dislocations and θ = π 2 for pure edge dislocations. In the present treatment, θ = π 4 is taken as a mean value. The linear misfit ε is approximated as 1 3 of the volumetric misfit. C Coh is the precipitation strengthening coefficient, which is a calibration parameter.

Precipitation Kinetics Model
Supersaturated states are unstable since a driving force exists to minimize the Gibbs free energy by the formation of a new phase, rearrangement of existing phases or redistribution of alloying constituents. For the evolution of the precipitate microstructure, the SFFK model [38][39][40] can be utilized, which is, for instance, implemented in the MatCalc software package (http://matcalc.at, accessed on 12 July 2022). However, these models need to be simplified to minimise both, the calculation time as well as the memory resources, within an FE framework. To determine the driving force for nucleation of precipitates, the knowledge of the equilibrium concentrations of Mg and Si in the Al matrix can be taken as a starting point. X 0 Mg and X 0 Si are the nominal Mg and Si mole fractions in the system and X p Mg and X p Si are the Mg and Si fractions inside the precipitates. The following description is given for one specific precipitation phase, e.g., clusters or the β" phase, but the extension to a multi-phase system is straightforward. The Mg and Si concentrations in equilibrium within the fcc Al matrix is given by mass conservation.
f eq is the equilibrium phase fraction of the precipitate, which is calculated by solving the following solubility product by numerical methods where ∆G is the energy of (precipitate) dissolution normalized with respect to RT C and D are input parameters within the sMSE framework. In contrast to classical precipitation calculations based on, e.g., the CALPHAD method, no thermodynamic databases are used in the sMSE framework, and the driving force is calibrated by the parameters C and D. The molar driving force for precipitation nucleation is approximated with The nucleation of new precipitates can be evaluated on the basis of the steady-state nucleation rate, which is defined as the number of newly formed precipitate nuclei per unit volume and unit time as [41,42] where N 0 is the number of available nucleation sites, Z is the Zeldovich factor, β * is the atomic attachment rate and G * is the critical nucleation energy. The Zeldovich factor is expressed as [41,43], ν α represents the molar volume of the precipitate. In this simplified framework, an average molar volume of 10 −6 m 3 /mol is assumed for all phases. γ is the specific interfacial energy. The critical nucleation energy is given by The atomic attachment rate reads with the critical radius The effective diffusion coefficient D eff is taken as D 0 is the pre-exponential factor, Q is the activation energy for diffusion, X Va is the current vacancy concentration and X Va,eq is the equilibrium vacancy concentration. The evolution of quenched-in vacancies is described on the basis of the FSAK-model [44], which considers the formation and annihilation of vacancies at grain boundaries, dislocation jogs or Frank loops. The change of mean radius in each time step due to precipitate growth is evaluated from the original SFFK treatment by solving the evolution equations for a single diffusing species and under the assumption that the precipitate is a stoichiometric compound. The rate of the radius due to particle growth, . r g , is then obtained with After the growth of precipitates has seized due to decreasing supersaturation, further growth of large particles commences at the expense of the smaller particles. This so-called coarsening process occurs with continuously increasing mean radius of the precipitates and a simultaneous decrease of their number density. The radius change due to coarsening in the classical LSW mean-field approximation reads where the subscript 'LSW' refers to the original work of Lifshitz and Slyozov [45] and Wagner [46]. The coarsening rate constant is obtained as [47], where the effective diffusion coefficient D eff is used in this framework and η LSW_fact is a fitting parameter.

State-Dependent Variables
For each time integration step dt, the evolution of the vector of state parameters, χ, is calculated within the subroutine. Table 1 assigns all state dependent variables to either thermal or athermal stress. The final stress σ(T, . ε,χ) is then calculated by the sum of σ th and σ at .

Material and Heat Treatment
The process of parameter identification and calibration of the present model is divided into two parts: (i) analysis of the nucleation and growth kinetics of precipitates during artificial aging, and (ii) characterisation of the work hardening behaviour evaluated with compression tests. The chemical composition of the commercial AA6082 aluminium alloy, which is used in the present experiments, is given in Table 2. After die-casting, the material is homogenised at 560 • C for 2 h, which represents the as-received state for the investigations. The specimens are then solution heat treated at 560 • C for 50 min in a circulating air furnace (Carbolite Type 3508), water-quenched and artificially aged in an oil bath between 0.5 h and 8 h for the first experimental setup (i). Between the solid solution heat treatment and the artificial aging (ii), the material remained at RT for 20 min.

Mechanical Testing
The compression tests are carried out in a high-speed quenching and deformation dilatometer DIL 805 A/D. Cylindrical specimens with lengths of 10 mm and diameters of 5 mm are produced from the as-received state. Prior to the deformation step, the specimens are solution heat treated at 530 • C for 5 min and helium-cooled with a cooling rate of 50 K/s to the deformation temperatures (25,100,200, 300, 400 and 500 • C), where they are held for 10 s to achieve sufficient thermal equilibration. Each deformation test is repeated at least twice with applied true strain rates of 0.1 s −1 and 1 s −1 .
Brinell hardness measurements (HBW 1/10) are carried out in an EMCO-Test M1C 010 unit, where at least eight measurements are performed for each aging time.

Electron Microscopy
Transmission electron microscopy (TEM) is applied for microstructure characterisation. The specimens are ground with SiC-papers until a thickness of approximately 0.1 mm is achieved. Subsequently, discs of 3 mm in diameter are punched out and electrochemically etched using a Struers Tenupol-5 in order to obtain electron-transparent regions. This process is carried out in a 2:1 mixture of methanol and nitric acid at temperatures between −20 • C and −10 • C applying a polishing voltage of 10 V. The analysis is performed on an FEI TECNAI F20 microscope equipped with a field emission gun and operated at 200 kV acceleration voltage. Bright field images in <001>Al zone axis reveal the lengths of the needle-shaped precipitates, whereas precipitate diameters are estimated using highresolution transmission electron microscopy (HRTEM). The thickness of the observed regions is measured by electron energy loss spectrometry (EELS) using the t/λ (log-ratio) methodology [48]. Upon data collection, average lengths, diameters, and number densities of the precipitates are estimated.

Simulation
To demonstrate the potential of the present state parameter-based concept, all simulations are conducted on the basis of one single set of input parameters. Table 3 summarises all input parameters for the precipitation kinetics calculation as described in Section 2.3. The parameters are calibrated on the precipitation statistics determined by TEM. In this setup, only one precipitate type is listed, since only β precipitates are detected in TEM during the deformation tests. Although a second precipitate type is available for simulation in the sMSE framework, only the observed β is accounted for in the following example according to experimental evidence.  Table 4 lists the input parameters for the strengthening models, which are described in Sections 2.1 and 2.2. The parameters are calibrated on basis of the flow curves obtained in the compression tests. Since precipitation strengthening significantly impacts the yield stress, calibration of the precipitation kinetic models is carried out prior to the calibration of the work hardening models.

Hardness Tests
The experimentally observed hardness curve for artificial aging at 180 • C is shown in Figure 2. The hardness rapidly increases until a peak is reached after approximately 8 h. Afterwards, the hardness decreases again due to overaging (coarsening) of the β precipitates. In addition, β , which is assumed to be the main hardening phase in the 6xxx series alloys [49], transforms into β . The latter process is not considered in the present analysis, though.

Precipitation Evolution
N v is the number of precipitates of the field of view (FOV) of the TEM image, l is the mean length of the precipitates, t is the thickness of the specimen, which is measured by EELS, and A FOV is the area of the field of view. Since the evolution of the mean radius r of a spherical particle is calculated within the sMSE framework, a conversion to the needle-shaped β is necessary. This is done by Equation (32), with the shape parameter h [57]. The length l of the precipitates increases in good agreement with the measured values as shown in Figure 4c.  It should be emphasised that all precipitates are assumed to be β in this simplified framework, because β is known to be the main strengthening phase during a T6 heat treatment. The measured number density is almost constant at low annealing times and increases slightly after two hours. The simulations show that the number density of β rapidly increases when the material is heated up to 180 • C and remains constant during the holding time.  Tables 3 and 4. The stress simulation comprises all strengthening contributions, which have been discussed in Sections 2.1 and 2.2. The simulations show that, the lower the deformation temperature is, the higher is the work hardening contribution. The stronger hardening behaviour at 200 • C and 300 • C at the lower deformation rate is due to strain-induced precipitation hardening. This effect is accounted for by an (artificially) adjusted specific interfacial energy of γ = 0.085 J/m 2 . To demonstrate the influence of precipitation strengthening on the work hardening behaviour, simulations without precipitation strengthening, which are marked with *, are included. To calculate internal stresses using the FEM mesh, an accurate simulation of the yield strength is required. Figure 6a shows the values for R p0.2 as a function of temperature for both strain rates, 0.1 s −1 and 1 s −1 . The stress values are higher for higher strain rates, except for a temperature of 100 • C and 200 • C. Although the quenching time and the holding time on the deformation temperature are very short, an explanation is given by the cross-core diffusion process on the one hand, and precipitation strengthening on the other hand. The calibration for the precipitation nucleation process and precipitation growth is discussed in Section 2.3. In Figure 6b, the experimental results of the yield stress are compared to the results of the simulation for both strain rates and six different temperatures. Although simplified models are used in the sMSE framework, experiment and simulation are in good agreement for different temperatures and strain rates with a single set of input parameters.

Conclusions
A state parameter-based framework for plastic deformation modelling is introduced, which is suitable for implementation in a user-subroutine in commercial FE-software packages. Well established nucleation and growth models for precipitation kinetics, as well as models for strengthening mechanisms, such as solid solution strengthening, precipitation strengthening and work hardening, are brought into a condensed form to reduce the calculation time for the stress simulation of possible complex structures. In this framework, two alloying components are considered, as well as two precipitate types, which are defined by constant stoichiometry, are included. Although severe simplifications are adopted, flow curves at various temperatures and strain rates are successfully reproduced and the trend of nucleation and growth of precipitates are reflected.