A thermo-viscoplasticity model for metals over wide temperature ranges- application to case hardening steel

In this contribution, a model for the thermomechanically coupled behaviour of case hardening steel is introduced with application to 16MnCr5 (1.7131). The model is based on a decomposition of the free energy into a thermo-elastic and a plastic part. Associated viscoplasticity, in terms of a temperature-depenent Perzyna-type power law, in combination with an isotropic von Mises yield function takes respect for strain-rate dependency of the yield stress. The model covers additional temperature-related effects, like temperature-dependent elastic moduli, coefficient of thermal expansion, heat capacity, heat conductivity, yield stress and cold work hardening. The formulation fulfils the second law of thermodynamics in the form of the Clausius–Duhem inequality by exploiting the Coleman–Noll procedure. The introduced model parameters are fitted against experimental data. An implementation into a fully coupled finite element model is provided and representative numerical examples are presented showing aspects of the localisation and regularisation behaviour of the proposed model.


Introduction
During the manufacturing of parts from case hardening steel high temperature differences occur. The first step, forging, is usually followed by the first chipping step, together defining the geometry of the parts. These intermediate parts are then exposed to a carbon rich environment at elevated temperature which promotes diffusion of carbon into the material and therefore leads to an increase of the carbon content close to the surface. Thus, these case hardened parts gain surface hardness while still having a ductile interior.
Experimental data provided by [4,11,13,25,29,32]  Institute of Mechanics, TU Dortmund University, Leonhard-Euler-Straße 5, 44227 Dortmund, Germany a coupled temperature and rate dependency of the plastic deformation behaviour. As elaborated in, e.g. [39], other thermomechanical properties, such as the thermal expansion, the heat capacity and the heat conductivity, are also highly temperature dependent. To obtain reasonable results in numerical simulations, these effects have to be covered. Commonly, this is done by, e.g. interpolating between yield curves experimentally obtained at constant temperatures and deformation rates, as applied by, e.g. [32], or by using direct approximations of the yield curve, as discussed in, e.g. [15,18]. However, the temperature evolution either remains unconsidered or is approximated by semi-analytical methods.
Another method is to consider the temperature and rate dependency a priori and to formulate the material model within the framework of rational thermodynamics established by [9]. This framework is used for thermoplasticity in, e.g. [3,31,34]; see also references cited in these works. The related numerical treatment is proposed by [34] along with an implementation into a coupled finite element framework. In [5] thermo-viscoplasticity is discussed and non-linear models for the treatment of rate-dependent yield stress are proposed. The works [6][7][8] demonstrate the advantages of considering temperature dependent material properties along with nonlinear kinematic hardening. The special focus on the amount of plastic work transformed to heat and calibration to experiments is emphasised in [12,30].
The localisation and regularisation behaviour of a Perzynatype viscoplastic model with strain softening is investigated in [21] and the influences of the applied loading rate and the choice of the model parameters are studied.
Thermoplasticity and the related localisation effects are numerically investigated in, e.g. [27,28,42] and gradientenhancement by an additional differential equation introducing an internal length scale is proposed as a solution strategy to overcome mesh dependency issues.
The key focus of this paper is the elaboration of a thermo-viscoplasticy model which is applicable to various ductile thermo-viscoplastic materials, especially metals. The material model includes temperature-dependent material parameters in a thermodynamically consistent manner and covers various effects related to thermo-viscoplasticity, such as non-linear hardening, thermal softening as well as heat generation due to structural cooling and plastic deformation. The choice of the ansatzes for the temperature-dependencies makes the model applicable over wide temperature ranges. The thermo-viscoplasticity model is embedded into a fully coupled finite-element framework. Moreover, aspects of the related localisation and regularisation behaviour are discussed based on numerical examples. This paper is organised as follows: In Sect. 2 the theory of thermo-viscoplasticity, involving temperaturedependent material properties, is introduced and is specified for 16MnCr5 (1.7131) case hardening steel. The implementation into a fully coupled finite element framework is outlined in Sect. 3 followed by the specification of the approximation functions for the temperature-dependent model parameters for steel grade 16MnCr5 and their identification in Sect. 4. Results of representative numerical examples, demonstrating the performance of the proposed model, are presented and discussed in Sect. 5. Finally, this paper is closed by a summary in Sect. 6.

Thermo-viscoplasticity
In this section, the essential balance relations for the closed thermo-mechanical system are presented in a general framework using thermodynamics with internal variables. Thereafter, the model used for the implementation as this work proceeds shall be specified.

Essential balance relations
The thermomechanical system considered shall fulfil the balance of linear momentum here given for the quasi-static case, as well as the balance of internal energy ρu = σ :ε + ρr − divq (2) and the Clausius-Duhem inequality describing the first and the second law of thermodynamics, respectively, cf. [31]. Therein, T denotes the absolute temperature, σ the stress tensor, the small strain tensor, u the displacement vector, and q the heat flux vector. Since these local forms of the balance laws are given per unit volume, all quantities per unit mass, i.e. the body force b, the heat source r , the entropy s, and the internal energy u, are transformed by the mass density ρ, whereby conservation of mass is assumed. The operator• indicates the time derivative d dt of the quantity •. In addition, for any Cauchy-Boltzmann continuum the balance of angular momentum is fulfilled a priori by the symmetry of the stress tensor, i.e. σ = σ t .

Dissipation inequality and heat equation
It is common to rewrite the energy balance and entropy inequality into more convenient formats. Therefore, the Helmholtz (free) energy is obtained from the Legendre transformation based on ψ = u − T s and inserted into (2), which results in the balance of Helmholtz energy Inserting this relation into (3) yields the dissipation inequality which defines the dissipation rate per unit volume, D. For the Helmholtz energy of the thermo-mechanical system, the ansatz is assumed, depending on observable quantities, i.e. the temperature and the total strain tensor as well as internal variables, see [9], here denoted as (second order tensor) z. Inserting (7) into (6) specifies the dissipation inequality Therein, the entropy and the stress tensor are identified as entities, energetically conjugated to the strain tensor and the temperature and, similarly, following the framework of rational thermodynamics, also the dissipative stress, ζ , is defined according to With these relations, the terms in brackets in (8) vanish and the reduced dissipation inequality takes the compact form Following the Coleman-Noll procedure [10], restricts this condition by separating the two terms into independent conditions, i.e.
which are referred to as the mechanical and the thermal dissipation inequality. By inserting the constitutive ansatz (7) and the constitutive relations (9) into the energy balance (5), it is transformed into the heat balance where denotes the heat capacity at constant strains and where represents the mechanical heat production rate. While a heat conduction model is used to fulfil the thermal part of the dissipation inequality, the mechanical dissipation inequality must be satisfied and, moreover, is assumed to fulfil the postulate of maximum plastic dissipation. Therefore, following the framework of generalised standard materials [14], the existence of a yield function f (T , ζ ) is assumed limiting the elastic region as the admissible space, i.e.
The optimisation problem is formulated, here using the penalty method with the penalty factor τ > 0, commonly interpreted as relaxation time, and the penalty function From the optimality condition the associated evolution equation follows for the internal variableṡ The definition of the viscoplastic flow rateλ vp can be rewritten into the residual format that provides the rate-dependent yield surface and, accordingly, f vp can be denoted as the viscoplastic yield function. In addition, this formulation resolves the definition gap for τ → 0 and thus also allows the rate-independent case, see [33].

Model specification
In the case of elasto-plasticity in an infinitesimal strain setting, the total strain tensor decomposes additively into the reversible thermo-elastic strain tensor, ε r , and the irreversible plastic strain tensor, ε p , as Thus, the internal variables contain the plastic strain tensor as well as, e.g. hardening variables. The reversible strains split into the thermal strains, ε ϑ , and the elastic strains, ε e , according to For these effective thermal strains, related to the absolute temperature of an initial stress-free state T 0 , the isotropic relation is assumed with the mean coefficient of linear thermal expansion, α m , as defined in Remark 1, and T 0 • C = 273.15 K providing the Celsius temperature scale.
This model will restrict to isotropic cold-work hardening with the related internal variable k, so that z = {ε p , k}. Inserting the linearised kinematics into (11) and (14) leads for the mechanical dissipation inequality and the mechanical heat production rate to and respectively, with κ = −ρ ∂ψ ∂k denoting the isotropic coldwork hardening related stress. The mechanical heat production rate can be split into an reversible and irreversible (plastic) contribution according to where the reversible part considers the Gough-Joule effect and the plastic part provides the heat production due to plastic yielding and hardening effects. It is assumed that the Helmholtz energy can be decomposed according to into a purely thermal part, depending only on the temperature, a thermo-elastic part, depending on the temperature and elastic strain tensor, and into a purely plastic part, depending only on the hardening variable. For the respective contributions the following specifications are introduced, i.e.
Therein, K and G denote the elastic moduli for bulk and shear response, respectively. Parameters c 0 , and w c i (i = 1, . . . , n − 1) are model parameters for the temperaturedependency of the heat capacity with respect to a reference temperature T r . For the thermoelastic contribution to the free energy, a volumetric deviatoric split is assumed, with ε dev e = ε e − 1 3 trε e I denoting the deviatoric elastic strain tensor. The chosen plastic part of the free energy is able to describe exponentially saturating hardening with the isothermal saturation yield stress increment Y ∞0 and the saturation parameter k 0 , as well as linear hardening with the isothermal hardening modulus H 0 . The mix parameter r mix ∈ [0, 1] specifies the fraction of linear hardening. The stress tensor and the cold-work hardening stress (9) read and The isomorphic heat capacity as present in the energy balance (12) results in Therein, γ = 3 α m + dα m dT [T − T 0 • C ] denotes the coefficient of volumetric thermal expansion. In experiments usually the heat capacity at constant stresses is determined, i.e.
as discussed in Remark 2, where ϕ denotes the Gibbs (free) energy of the thermo-mechanical system obtained from the Helmholtz (free) energy by the Legendre transformation ϕ = ψ + 1 ρ σ : ε. Based on the ansatz specified in (28), (29), (30) and (31) one obtains Therein, p = − 1 3 tr(σ ) and σ dev = σ + p I denote the hydrostatic pressure and the deviatoric stress tensor, respectively. Combination of (35) and (36) yields It is remarked that for σ = 0 the isobaric heat capacity reduces to a polynomial in T . Similar approaches of second order have been proposed by, e.g. [36,41].
Since the plastic deformation of metals is assumed deviatoric, i.e. trε p ≡ 0, and, moreover, considering isotropic plastic response, the use of the von Mises yield function as yield criterion and plastic potential is appropriate. Therein, y 0 denotes the initial yield stress and ϕ κ accounts for the thermal softening of the contribution related to cold-work hardening, κ. For the penalty function, the power-type ansatz is chosen, with the temperature-dependent parameter m, and d 0 , the stress-like scaling parameter, cf. [5]. Here, d 0 = √ 2/3 y 0 is assumed. The associated evolution equations of the internal variables (19) can be derived from these quantities, to be specifiċ wherė denote the Perzyna-type viscoplastic flow rate and the deviatoric yield direction, respectively. The formulation of the viscoplastic yield function (20) in terms of stresses, allows the identification of viscous overstress which describes a power-law relation between the overstress and the viscoplastic flow rate. The thermal part of the dissipation inequality is fulfilled by Fourier's law of heat conduction with the positive semi-definite heat conductivity tensor K . As this work proceeds, the isotropic relation K = Λ(T ) I shall be assumed.

Remark 1
In the literature, definitions of two different coefficients of linear thermal expansion can be found: The real coefficient of linear thermal expansion and the mean coefficient of linear thermal expansion where L(T ) denotes the length of a specimen at temperature T and T 0 • C is the reference temperature, commonly chosen at T 0 • C = 273.15 K, see, e.g. [19,37,41]. From both relations, the thermal strain between an initial temperature T 0 and the current temperature T , i.e.
can be obtained. From (48) one may conclude that the relation between the real and mean coefficient of linear thermal expansion holds. Due to its simple definition, α m (T , commonly provided in the literature, as in, e.g. [20,35,37].

Remark 2
The isobaric heat capacity is defined by where Q specifies the heat in a body of mass m. For the experimental determination of c σ (T ), the heat in a specimen of mass m is increased by the portion Q and the resulting temperature change from T 1 to T 2 is measured, see, e.g. [ see, e.g. [19,36]. If the temperature difference T 2 − T 1 becomes sufficiently small, the mean heat capacity approaches the real heat capacity, i.e. c σ m ( Thus, values of the mean isobaric heat capacity c σ m (T i , T 0 • C ) at different temperatures T i with respect to a reference temperature, commonly T 0 • C = 273.15 K, are often found in the literature and interpolating functions for c σ m (T , T 0 • C ) are provided as in, e.g. [20,35,36]. In analogy to the relation between the coefficients of thermal expansion, the relation holds between the real and mean isobaric heat capacity.

Remark 3
The particular model proposed in the present paper introduces temperature-dependencies of cold-work hardening within the evolution equations of the related internal variables rather than within the related driving stress-type quantities. This ansatz implies that κ shows no direct dependencies on the temperature T so that Q p = D m ≥ 0 (27) can be satisfied straightforwardly. The similar ansatz for the plastic contribution to the free energy and the yield function including the temperature dependency of the cold-work hardening directly in the free energy by the temperature dependent hardening moduli y ∞ and h, implies that κ(T , k). Although, the contributions of the cold work-hardening ϕ κ (T )κ(k) in (38) and κ(T , k) in (54) are identical for both models, i.e. κ = ϕ κ κ holds and thus, the resulting dissipation rates are equivalent, the rates of plastic heat production (27) are not. Therefore, such model introduces the additional constraint Q p ≥ 0 ∀ k ≥ 0 to the ansatz for the temperature dependencies and the choice of the related model parameters.

Numerical treatment
In this section, the numerical treatment of the balance equations in a finite element framework and of the presented material model in a predictor-corrector scheme are outlined. Table 1 Algorithmic box of the presented material model. The subscript K denotes the use of Kelvin notation for symmetric tensors

Finite element framework
The balance equations (1) and (12) of the thermo-mechanical initial boundary value problem describe its mechanical and its thermal sub-problem. A related graphical illustration is depicted in Fig. 1. The mechanical sub-problem, Fig. 1a, is built up by the surface of the body ∂B being subjected to displacement (Dirichlet) boundary conditions u =ū on ∂B u and traction (Neumann) boundary conditions σ · n =t on ∂B t with ∂B = ∂B u ∪ ∂B t ∧ ∂B u ∩ ∂B t = ∅ and outward surface normal unit vector n. In the interior domain B volume forces ρ b are active.
The thermal sub-problem, Fig. 1b, consists of heating of the body by volumetric heat sources, ρ r , and mechanical heat production Q m in the interior domain, as well as by surface heat fluxes (Neumann), −q · n =q n on ∂B q , and by prescribed temperatures Following the standard procedure, as outlined in, e.g. [1], and discretising the body by n el non-intersecting finite elements B e , with n sur el surface elements ∂B e , the weak forms of the balance equations (1) and (12), are obtained for all admissible virtual fields, δu and δT . Within the finite elements B e , the approximations δu(x) ≈ n en A=1 δu e A N A (x) and δT (x) ≈ n en C=1 δT eC N C (x) are applied, where n en denotes the number of element nodes of element e and where δu e A and δT eC are the values of the respective global field at node A, respectively C, of the particular element e. The element shape functions are denoted as N A (x) and an isoparametric approach is adopted. The residual format of the balance equations in global form reduces to Therein, the operator A The constitutive equations for c ε , σ , q, b, r , Q m , t and q n may, in general, depend on the state variable fields, i.e. Fig. 4 Approximation of the mean isobaric heat capacity c σ m , experimental data from [35], resulting real isobaric heat capacity c σ , and isomorphic heat capacity c ε . For c σ and c σ m a purely volumetric stress state with the atmospheric standard pressure p = p at = 1013.25 hPa is assumed and, moreover, ε = 0 is assumed for c ε Fig. 5 Approximation of the heat conductivity Λ(T ) and experimental data from [35] displacement u and temperature T , as well as on internal variables z, and on related gradients and rates.
According to the Bubnov-Galerkin formulation, the displacement u(x) ≈ n en A=1 u e A N A (x) and temperature field T (x) ≈ n en C=1 T eC N C (x) + T 0 as well as their rates are approximated by using the same ansatzes as for their virtual equivalents. Therein, T 0 indicates the temperature of the initial stress-free state, here assumed to be constant over the entire body. Global vectors for the global degrees of freedom (DOFs) of the n np nodes of the system can be defined as which similarly applies to the rates,ḋ. These are approximated within a time integration scheme referred to discretised time intervals with the increments t = t n+1 −t n as function of the current and previous values of the degrees of freedom and initial conditions such thatḋ(d, d n ,ḋ n ). Therein and in the progression of this work, the index n + 1 is omitted for better readability. In the material subroutine, discussed in Sect. 3.2, the internal variables z and their ratesż are computed implicitly from the constitutive assumptions, which depend only on the current values of the state variables d as well as on the internal variables of the previous time step z n .
Since in the numerical examples in Sect. 5 only surface loads independent of the nodal state variables and their rates are discussed, the further derivations will be restricted to 'dead' loads. To solve this discrete spatio-temporal non-linear residual system of equations by means of the Newton-Raphson scheme, the global tangent matrix J is required. It is defined by Therein, the tangent matrix of element e is defined by For the implementation of the proposed material model, it is convenient to use the AceGen-Toolbox [22,23] for Wolfram Mathematica [43]. It generates exact derivatives and : applied on the residuum vector of the evolution equation of the internal variables r z as defined in Table 1 in Sect. 3.2 at the converged state, i.e. r z ≈ 0.

Algorithmic treatment of the material model
Based on the finite element discretisation, the values of the temperature field T , its gradient ∇T , the displacement gradient ∇u, its rate ∇u, and the time increment t = t − t n are provided. In addition, the internal variables, i.e. the plastic strain tensor ε p n and the isotropic hardening variable k n of the previous time step are known. The implementation of the material model presented in Sect. 2.3 follows a predictor-corrector scheme as outlined in [33]. Therefore, in the predictor step all variables independent of the internal variables as well as trial values for those variables that depend on the internal state variables are computed based on which the yield condition can be evaluated. If the yield condition is fulfilled, the trial values are used as an update for the variables. If not, the internal variables are updated by using a Newton-Raphson scheme until the viscoplastic yield condition is fulfilled. The algorithmic box in Table 1 states the necessary steps for the implementation of the presented material model.

Parameter identification
In this section, the approximation functions for the temperature dependent parameters of the presented material model are given and the identification procedure of the model parameters based on experimental data is shown along with a table of the results for 16MnCr5 (1.7131) case hardening steel.
For the temperature dependency of the elastic bulk and shear moduli, K and G, as well as for the initial yield stress y 0 , the relaxation time τ and the viscoplastic exponent m, the representation is assumed which includes a continuous function ϕ X . Therein, X and X 0 denote the respective temperature dependent quantity and its value at reference temperature T = T r , respectively. For the continuous functions ϕ K , ϕ G , ϕ y0 , and ϕ κ , the arctangent ansatz is assumed, where φ X defines the initial 'angle' of the arctangent curve at reference temperature T r = 0. Moreover, the function ϕ τ is modelled by the exponential ansatz Both functions (70) and (71) include scaling parameters ω X and temperature dependencies are introduced by χ X . The function ϕ m is assumed to be linear according to the recursive formula for an n-th order polynomial with n ≥ 1 and ϕ p0 X = 1. In general, the coefficients w Xi with i = 1, . . . , n, are included, whereas subsequent elaborations are restricted to n = 1. The temperature dependencies of the mean coefficient of thermal expansion α m and of the heat conductivity Λ are described as piecewise continuous functions Moreover, n-th order polynomial relations, as defined in (72), with degree n = 1 are assumed for the functions ϕ α1 , ϕ α2 , n = 2 for ϕ c1 , ϕ c2 , and n = 4 for ϕ Λ1 . The contribution ϕ Λ2 is modelled by the exponential relation given in (71). For the identification of the material parameters, all parameters p X describing the temperature dependent quan- 12 Comparison of the temperature distribution T (x) at initial temperature T 0 = 20 • C and averaged shear ratesγ xy tity X (T , p X ) are identified by minimising least squares functionals of the form where X exp (T i ) denotes experimental (literature) data of the quantity X at discrete temperatures T i . For the minimisation, the MATLAB [26] implementation fminsearch of the Nelder-Mead simplex algorithm [24] is used. Since the simplex algorithm is an unconstrained optimisation method, special care has to be taken for those parameters which are defined only on an interval, such as φ • ∈ [− π 2 , π 2 ] and {Y 00 , Y ∞0 , k 0 , τ 0 , m} ∈ (0, ∞]. Dur-ing the identification procedure, these have been replaced by parameters p • according to φ • = arctan( p φ • ) and

Thermo-elastic properties, G, K, and˛m
This subsection presents the identification of the thermoelasticity related parameters of the model. To identify the parameters for the elastic moduli, the experimental data points for the Young's modulus E and Poisson's ratio ν are converted into values for the bulk modulus K and shear modulus G, respectively, according to The obtained values of the parameters of the arctangent model (70), Table 2. To show the quality of the approximation, the resulting curves for K (T ), G(T ), and E = 9K G 3K +G are presented in Fig. 2 along with the experimental data from [13,38]. Therein, the solidus temperature of pure iron T Fe s = 1538 • C [20] provides an estimate about the feasible temperature range of the model. For all moduli the model shows a good agreement with the experimental data but does not catch the relatively steep decline around room temperature of the bulk modulus, and thus as well of the Young's modulus.
The mean coefficient of thermal expansion α m is modelled by piecewise continuous linear functions with a jump at the austenitisation temperature T A . The obtained parameters p α = {α 01 , w α 1 , α 02 , w α 2 , } are given in Table 2. The good agreement of the resulting curve α m (T ) with the experimental data from [35] is shown in Fig. 3.

Thermal properties c and
The model parameters for the heat capacity cannot be identified directly from experimental data, since in experiments the mean isobaric heat capacity (51) is determined, from which the real isobaric heat capacity (50) can be obtained by (52). For the specimen subjected to standard conditions ( p = p at = 1013 hPa) in an ideal setting with σ = −p I, the heat capacity at constant hydrostatic stress (37) with n = 1 results in and thus depends on K (T ) which has already been identified. Since the material considered, i.e. 16MnCr5 case hardening steel, undergoes a phase transformation, values for the model parameters c 0 , w c1 and w c2 are identified independently below austenitisation temperature T ≤ T A as c 0 1 , w c1 1 and w c2 1 with respect ot the reference temperature T r , and above i.e. T > T A as c 0 2 , w c1 2 and w c2 2 with respect to T A  w c1 1 , w c2 1 , c 0 2 , w c1 2 , w c2 2 } are given in Table 2. Inserting the condition ε = 0 into (34) leads for the heat capacity at constant (zero) strain to  Fig. 4. Therein, the real isobaric capacity c σ , the resulting mean isobaric heat capacity c σ m , and the isomorphic heat capacity c ε are plotted against the temperature. The plots of the real and mean isobaric heat capacity show good agreement between the model and the experimental data from [35]. Moreover, the differences between the values of the different heat capacities become visible.

(b) (a)
The heat conductivity is approximated by piecewise continuous functions. To be specific, a fourth order polynomial below the austenitisation temperature T A and an exponential ansatz above T A are chosen. The obtained parameters p Λ = {Λ 01 , w Λ 11 , w Λ 21 , w Λ 31 , w Λ 41 , Λ 02 , ω Λ 2 , χ Λ 2 } are given in Table 2. To show the quality of the approximation, the resulting curve Λ(T ) is depicted in Fig. 5 together with the experimental data from [35]. -viscoplastic properties y 0 , ' Ä , r mix , H 0 ,  1Y ∞0 , k 0 , , and m The identification of the parameters for the visco-plasticity related part of the model follows a staggered scheme: In a first step, exponentially saturating hardening, i.e. r mix = 0, is assumed. From yield curves obtained at different temperatures T i and different strain rates, found in [32], values of the initial yield stress y 0 (T i ), the saturation yield stress

Thermo
In a second step, the parameters p y 0 = {Y 00 , ω y 0 , χ y 0 , φ y 0 } as well as updated values of y ∞ (T i ), k 0 (T i ), τ (T i ), m(T i ) are identified from the data obtained in step 1 as well as from experimental data from [25].
In a fourth step, linear hardening, i.e. r mix = 1, is assumed and the hardening modulus H 0 is identified. Finally, the mix parameter r mix is identified in a fifth step. It appears that purely exponential hardening fits best and thus r mix converges towards zero. All parameters of the presented model are provided in Table 2.
In Fig. 6, the resulting curve for the initial yield stress y 0 (T ) and the experimental data are plotted. Both are in good agreement. The depicted saturation yield stress y ∞ = y 0 (T ) + ϕ κ (T ) Y ∞0 also shows distinct thermal softening.
In Fig. 7, the viscous over-stress σ is depicted for selected effective plastic strain ratesε p = √ 2/3 ε p in a logarithmic scale over the temperature range of interest. Due to the power-type ansatz for the viscous over-stress, a tenfold increase of the effective plastic strain rate leads to a constant spacing between the curves. It can be observed that for plastic strain rates > 1/s, the viscous over-stress may dominate the total yield stress with increasing temperature, whereas for low plastic strain rates it does not.
In Fig. 8, the stress-strain curves for a uniaxial stress state predicted by the material model at constant strain rates and temperatures are compared with experimental data from [32]. At room temperature T = 25 • C the initial yield stress as well as the saturation behaviour are captured well by the model for all strain rates. At T = 800 • C the initial yield stress overestimated by the model but the correct stress range is met. This over-estimation of the initial yield stress can be explained by comparably large derivations in the mechanical properties when comparing one batch to another and by a lack of experimental data performed on specimen from one batch. In addition, the increasing share of the viscous overstress with increasing temperature is visible. The comparison of the curves obtained at T = 25 • C and T = 800 • C with the curves obtained at T = 600 • C shows the non-linearity of the thermal softening: Although the temperature difference between 25 • C and 600 • C is approximately three times larger as between T = 600 • C and T = 800 • C, the curves obtained at T = 600 • C differ less to the ones obtained at T = 25 • C than to those obtained at T = 800 • C.

Numerical examples
In this section, aspects of the localisation and regularisation behaviour of the representative initial boundary value prob-

Shearing of plate with perturbation
To investigate the effect of the model parameters on the localisation behaviour of the model, a perturbed shear problem as depicted in Fig. 9 is analysed. To be specific, simulations of a block with edge lengths h = 1 mm subjected to shear loading are performed. The shear loading is applied in vertical directionū y (t) on the left and right edges in opposite directions, the displacement in the remaining directions is inhibited on all nodes,ū x = 0 andū z = 0. Since localisation is expected when hardening does not compensate thermal softening effects in this analysis, the applied magnitudes of shear loading are in general not within the usual range of a smallstrain model. The boundaries are assumed to be adiabatic. For the investigation of the mesh-dependency of the problem, the block is discretised by n el = 20, . . . , 2000 eight-noded hexahedral Lagrangian elements along the horizontal xdirection and by 1 element along the other directions. The influence of the loading-rate, the viscous over-stress and the internal heat flux on the results is of special interest. Therefore, the problem is simulated at different loading rates, u y = 0.05 mm/s,u y = 0.5 mm/s andu y = 5 mm/s, for different ratios of the heat conductivity Λ 0 /Λ 16MnCr5 0 = 0  In Fig. 10, the resulting force-displacement curves are depicted for both problems, i.e. with, Fig. 10a, and without internal heat flux, Fig. 10b. In both cases the initial temperature is T 0 = 20 • C. Whereas the problem with internal heat flux does not show mesh dependency in the simulated regime, the problem without heat conduction does: At certain loading points in the considered range of applied displacement, the force-displacement curves observed at average shear rates oḟ γ xy = 1 /s andγ xy = 0.1 /s drop, where their negative slope increases with increasing element-number n el , i.e. decreasing mesh-size h el = h/n el . Also, with decreasing loading rate, the differences between the force-displacement curves of the different discretisations initiates at smaller magnitudes of the boundary displacements, e.g. atū y ≈ 0.2 mm forγ xy = 1 /s and atū y ≈ 0.1 mm forγ xy = 0.1 /s.
The comparison of the averaged shear strain distributions γ xy (x) at the integration points of the different discretisations depicted in Fig. 11, shows only minor localisation with decreasing loading rate for the simulations considering heat flux within the specimen, and all discretisations converge to the same solution. The relative differences of the averaged shear strainsγ xy at the centre of the specimen between the coarse (n el = 20) and the finer discretisations are only ≈ 6.5 % for the slowest loading rate (γ xy = 0.1 /s) and ≈ 2 % for the fastest loading rate (γ xy = 10 /s). For the simulation neglecting heat flux, the comparison shows an obvious localisation with decreasing mesh size for the slowest loading rateγ xy = 0.1 /s. Its maximum values atγ xy = 0.1 /s are 4.8 times the average for the coarse discretisation with n el = 20, 17 times for n el = 200, and 48 times for n el = 2000. For the fast loading rate,γ xy = 10 /s, the shear strain distributions of the two fine discretisations do not show obvious differences in magnitude and converge to the same solution.
In Fig. 12 the resulting temperature distribution is shown. For both loading rates, the results of the simulations of considering heat flux converge to the same results. Due to the influence of the heat flux, a decreasing loading rate leads to a more uniform temperature distribution. Comparing the two plots of the simulations neglecting the influence of heat flux, a localisation in the temperature increase T (x) caused by the localisation in strain and thus by the localisation of the mechanical heat production can be observed for the problem at slow loading rateγ xy = 0.1 /s. There, the maximum temperature change rises with the increasing number of elements from 122 K for n el = 20 up to 1330 K for n el = 2000. For the problem at fast loading rate, which in Fig. 10b does not show any divergence of the force-displacement curves obtained for the different discretisations in the investigated displacement range, the temperature increase converges also for the fine discretisations, n el = 200 and n el = 2000 to 213 K. Figure 13 compares the strain rate dependency of the width of the shear zone b at different initial temperatures for simulations run with Λ 0 /Λ 16MnCr5 0 = 1, Fig. 13a, and Fig. 13b. The width of the shear zone b is identified by the relation x ∈ [−b, b] | γ xy (x) ≥γ xy . It can be observed that for both cases, i.e. Λ 0 /Λ 16MnCr5 0 = 1 and Λ 0 /Λ 16MnCr5 0 = 0, with increasing strain rate, b seems to converge against values which are only little below the dimension of the imposed cosine-shaped perturbation b * = 0.25 mm. For low loading rates, however, the curves of the shear zone decline towards b ≈ 0.14 mm for Λ 0 /Λ 16MnCr5 0 = 1 and b ≈ 0.08 mm for Λ 0 /Λ 16MnCr5 0 = 0. This decline appears to show stronger saturation for Λ 0 /Λ 16MnCr5 0 = 1 than for Λ 0 /Λ 16MnCr5 0 = 0 which may indicate a regularising effect of the heat conduction. The tendency of an increasing width of the shear zone b with increasing initial temperature T 0 can be observed. Figure 14 depicts the resulting force-displacement curves for different values of the heat conductivity for the case τ 0 /τ 16MnCr5 0 = 0 atγ xy = 1 /s for two different discretisations n el = 200 and n el = 2000. Up to u y ≈ 0.06 mm all curves are congruent. Beginning at that point, the curves for low or no heat-conductivity show stronger softening behaviour than the curves at Λ 0 /Λ 16MnCr5 0 = 1. This is caused by heat being more efficiently transported away from the centre of the specimen in the latter case (Λ 0 /Λ 16MnCr5 0 = 1), see Fig. 16). The curves for the different discretisations for Λ 0 /Λ 16MnCr5 0 = 1 and Λ 0 /Λ 16MnCr5 0 = 10 −3 remain congruent in the sequel or show only minor deviations, whereas for lower or no heat conductivity, Λ 0 /Λ 16MnCr5 0 ≤ 10 −4 , the force-displacement curves for both discretisations diverge from each other starting at u y ≈ 0.15 mm.
In Fig. 15, the resulting distributions of the normalised shear angle at u y ≈ 0.16 mm are compared for both discretisations. Both figures show an increase of the maximum of the normalised shear angle γ xy /γ xy at x = 0 with decreasing heat conductivity. The comparison of both discretisations shows an increase of the maximum values from γ xy /γ xy ≈ 6.9 (Λ 0 /Λ 16MnCr5 0 = 10 −5 ) and γ xy /γ xy ≈ 7.4 (Λ 0 /Λ 16MnCr5 0 = 0) for n el = 200, Fig. 15a, to γ xy /γ xy ≈ 54 (Λ 0 /Λ 16MnCr5 0 = 10 −5 ) and γ xy /γ xy ≈ 75 (Λ 0 /Λ 16MnCr5 0 = 0) for n el = 2000, Fig. 15b, which is an increase of the order of the reciprocal mesh size. For higher heat conductivity only minor deviations appear between the two discretisations. Figure 16 shows the resulting temperature distribution T (x). In contrast to the curves obtained for Λ 0 /Λ 16MnCr5 0 = 1, which appears to be almost a homogeneous distribution, the curves at reduced heat conductivity (Λ 0 /Λ 16MnCr5 0 ≤ 10 −3 ) show a distinct peak at the centre of the specimen (x = 0), which increases with decreasing heat conductivity. Similarly to the shear angle, the maximum values for low or no heat conductivity increase significantly with decreasing mesh size from T ≈ 119 K (Λ 0 /Λ 16MnCr5 0 = 10 −5 ) and T ≈ 126 K (Λ 0 /Λ 16MnCr5 0 = 0) for n el = 200 , Fig. 16a, to T ≈ 300 K (Λ 0 /Λ 16MnCr5 0 = 10 −5 ) and T ≈ 690 K (Λ 0 /Λ 16MnCr5 0 = 0) for n el = 2000, which is of the fac-tors ≈ 2.5 and ≈ 5.5, respectively, larger for the finer mesh (n el = 2000), Fig. 16b, than for the coarser mesh (n el = 200), Fig. 16a. Figure 17 compares the influence of the viscoplastic parameters τ 0 (Fig. 17) and m 0 (Fig. 17b) on the stability of the solution. The simulations have been performed at an averaged shear rateγ xy = 1 /s starting at the initial temperature T 0 = 20 • C and adiabatic conditions Λ 0 /Λ 16MnCr5 0 = 0 have been assumed. It can be denoted that a larger characteristic time τ 0 has a stabilising effect, making the solution meshindependent for a larger strain-range, whereas a smaller viscous exponent m 0 has the opposite effect, although both increase the current yield stress. For both cases, a variation of τ 0 and a variation of m 0 , the force displacement-curves obtained for both discretisations, n el = 200 and n el = 2000, diverge as soon as certain critical points have been reached and softening occurs.
The resulting distributions of the normalised shear angle γ xy (x)/γ xy are shown in Fig. 18. The comparison of the distributions depicted in Fig. 18a gated examples, a continuous differentiable perturbation has been imposed, and the distribution of the shear angle also seems to be continuous differentiable. This is in good agreement with analytical derivations found in, e.g. [2,17].

Plate with a hole under tension
To investigate aspects of the regularisation behaviour of the model, a non-homogeneous problem of a plate with a circular centric hole, as depicted in Fig. 20, is considered.
Starting at an initial temperature of T 0 = 20 • C, the plate is subjected to vertical displacementū y in opposite directions along its upper and lower surfaces. The response of the system is investigated for three different constant loading ratesu y = 2 mm/s,u y = 20 mm/s andu y = 200 mm/s, which are equivalent toε yy = 0.1 /s,ε yy = 1 /s, anḋ ε yy = 10 /s, respectively. The remaining displacement components at both loaded surfaces are equal to zero, i.e.ū x = 0 andū z = 0. All boundaries are assumed to be adiabatic. Based on the multiple symmetry of the problem, only the part of the plate located in the positive octant is considered. It is discretised by 10272 eight-noded hexahedralB-elements, cf. [16, ch. 4.5.2], where four elements are considered in thickness direction z. The resulting force-displacement curves, compared in Fig. 21, increase to maximum values which are between 5.1 and 5.2 kN for the three different loading rates considered. The maximum of the curve obtained atε yy = 1 /s at u y ≈ 2 mm, F y ≈ 5.12 kN, is slightly below the maximum of the curve obtained atε yy = 0.1 /s at u y ≈ 1.8 mm, F y ≈ 5.13 kN. The curve obtained forε yy = 10 /s shows its maximum of F y ≈ 5.2 kN at u y ≈ 2.3 mm. For displacements larger than at the force maxima, all three curves show a decreasing force F y , where the magnitude of the negative slopes of the force-displacement curves increases with increasing strain rate. Figure 22 shows the resulting distribution of the isotropic hardening variable k atε yy = 0.05 which, for moderate temperature changes and monotonically increasing loading, gives a good approximation of the effective plastic strain. For all loading rates compared, the distribution of k shows a localisation at the hole under 45 • with respect to the overall loading direction. Since the maximum value of k is highest forε yy = 1 /s, the localisation seems to reach a maximum for a loading rate betweenε = 0.1 /s andε = 10 /s. Figure 23 compares the resulting temperature distribution for the considered loading rates atε yy = 0.05. Like the hardening variable, the temperature increase T also concentrates at the hole, whereas, with increasing loading rate, a more distinct localisation in direction of 45 • is observable.

Summary
In this paper, a material model for thermo-viscoplasticity has been established in the setting of thermodynamics with internal variables and implemented into a fully coupled finite-element framework. The model incorporates thermal softening of the thermal and thermo-elastic properties, as well as of the thermo-plastic properties. Mechanical heat production due to plastic deformation as well as the Gough-Joule effect are covered. The model parameters have been identified based on experimental data for 16MnCr5 (1.1731) case hardening steel.
The regularisation effects arising from heat conduction and strain-rate hardening (viscosity) on the localisation phenomena related to thermal softening and saturating strain hardening included in the model have been investigated for different strain rates and values of the heat conductivity, viscoplastic characteristic time and viscoplastic exponent. The examples show a regularising effect of the heat conduction as well as the viscoplastic characteristic time and the strain rate, whereas a decreasing viscoplastic exponent has a nega-tive effect on the convergence of the finite element solution. It can be concluded that for the given parameter set and boundary conditions considered, the regularisation due to the heat conduction is more pronounced than due to viscoplasticity.
Since the model is formulated in a small strain setting, it does not cover geometric effects like, e.g. necking. Thus, for the investigation of such effects, a formulation in a finite strain framework is deemed necessary.