Modelling of temperature and strain rate dependent behaviour of pearlitic steel in block braked railway wheels

Block braked railway wheels are subjected to thermal and rolling contact loading. The thermal loading results in high temperatures and thermal stresses which cause slow time dependent processes such as creep, relaxation and static recovery of the wheel material. At the same time, the rolling contact loading implies a very fast mechanical load application. This paper is focused on material modeling of pearlitic steel for a wide range of loading rates at elevated temperatures. The starting point is a viscoplasticity model including nonlinear isotropic and kinematic hardening. The Delobelle overstress function is employed to capture strain rate dependent response of the material. The model also includes static recovery of the hardening to capture slower viscous (diffusion dominated) behaviour of the material. Experiments for the pearlitic wheel steel ER7 in terms of cyclic strain-controlled uniaxial tests with hold-time, uniaxial ratchetting tests including rapid cycles and biaxial cyclic tests with tension/compression and torsion are used to calibrate the material model. These experiments were performed under isothermal conditions at different temperatures. In the ratchetting tests, higher loading rates are obtained and these have been used to calibrate the high strain rate response of the viscoplasticity model. The paper is concluded with a numerical example of a block braked wheel where the importance of accounting for the viscoplasticity in modelling is highlighted.


Introduction
The material in railway wheels is, during regular railway operations, subjected to very high stresses and sometimes also elevated temperatures. The stresses and temperatures are caused by the combination of rolling contact between the wheel and the rail together with frictional heating due to block braking. This may result in wheel tread damage which limits the capacity of wheels from reaching desired mileage performance. The main damage mechanisms are rolling contact fatigue (RCF) and thermal cracking [1,2].
The behaviour of near-pearlitic railway wheel steels at elevated temperatures have been studied in, e.g. [3][4][5]. It is important to recognize that due to mechanical loading and elevated temperatures, the cementite lamellas tend to break up. This microstructural degradation is quantified by the degree of spheroidisation by Nikas et al. [5]. Microstructural degradation rates and other diffusion controlled processes increase with temperature elevation causing a viscous behaviour of the material. This viscous behaviour is observed in [5] to be significant in cyclic tests with hold time in terms of stress relaxation for temperatures above 200 C.
In [6] numerical simulations to predict damage mechanisms in railway wheel treads have been performed for the combined mechanical loading from the wheel-rail contact and the thermal loading from block braking. A methodology to simulate brake rig experiments has been developed in [6,7]. The methodology includes a thermal axisymmetric finite element (FE) analysis to find the temperature field caused by the braking and a mechanical 3D FE analysis of the wheel-rail contact. The material modelling is based on uniaxial cyclic strain-controlled experiments [5] at temperatures from 20 to 625 C. The choice of material model in [6,7] is a (rate independent) plasticity model of Chaboche type [8].
However, the significant viscous behaviour of the wheel steel has not been accounted for in previous work. The main goal of this paper is to improve modelling of the wheel material subjected to mechanical and thermal loadings in terms of cyclic plastic and viscous behaviour. To do so, the experimental data of the steel ER7 from [5] that includes uniaxial strain-controlled cyclic experiments with hold time are used to calibrate the material model. A viscoplastic material model including combined nonlinear isotropic and kinematic hardening (see, e.g. [8]) has been chosen for this purpose.
A main challenge in the present application is the large difference in loading rates of the very fast load application introduced by the rolling contact and the slow thermal loading caused by braking. In this regard, experimental data of the steel ER7 from uniaxial ratchetting experiments with rather fast loading at elevated temperatures are used to characterize the viscous behaviour of ER7 at high strain rates. The experiments were partly described in [9].
Different approaches have been suggested in literature to capture wide loading rate ranges and associated physical mechanisms in viscoplasticity models. In Liang and Khan [10], different constitutive viscous models and their performance to predict the response of BCC (body-centered cubic) and FCC (face-centered cubic) metals are examined. The formulations of the models correspond to overstress functions of exponential type. The models are capable of capturing the response for a wide range of loading rates. Furthermore, in Becker-Hackenberg [11] and Wang-Steinmann [12] models, the inelastic strain is additively decomposed into a viscoplastic and a creep strain. The viscoplastic strain in the Becker-Hackenberg model is governed by a sinh overstress function, whereas in the Wang-Steinmann model, it is governed by an exponential overstress function. Different overstress functions and their behaviour are examined by Chaboche in [13] where it is concluded that the exponential and sinh overstress functions (proposed by Delobelle) yield similar response. The appearance of a saturation effect of the viscosity for high strain rates regimes in some materials justifies the use of hyperbolic sine function, see, e.g. discussion in [13]. To limit the stress values at very high strain rates, a possibility is to use a limit/dynamic yield surface [11,14] and thereby for such strain rates obtain a reduced rate dependence.
In the current paper, we adopt a viscoplasticity model with Delobelle sinh overstress function. In addition, to further improve the viscous modelling, static recovery of the hardening variables is included in the model to capture slow (diffusion dominated) processes in the material, cf. [15].
Also, the Armstrong-Frederick model [16] has been adopted as the nonlinear kinematic hardening rule. This model has a tendency to overpredict the ratchetting in many metallic materials under both uniaxial and multiaxial loading conditions (see [17][18][19][20]). As simulation of the ratchetting greatly depends on the modeling of kinematic hardening, modifications of the dynamic recovery term in the Armstrong-Frederick model were proposed by Burlet and Cailletaud [21], Chaboche and Nouailhas [17], Chaboche [22], Ohno and Wang [18].
Also, in order to predict the behaviour of the material when it is subjected to multiaxial loading, similar to wheel-rail contact loading, some work has been done, in the literature, on rate-independent and rate-dependent plasticity modelling of railway steels (see [23][24][25]). In this study, to account for ratchetting under multiaxial loading, the description of the hardening is extended by including the Burlet-Cailletaud model [21]. In this regard, experimental data of the steel ER7 from biaxial (tension/compression and torsion) experiments at elevated temperatures from [26] are used to identify the multiaxial hardening parameters of the model. The paper is organized as follows: in Sect. 2, the viscoplastic material model is introduced together with the numerical procedure to incrementally solve the equations. In Sect. 3, the experiments and the calibration of the material model are presented. The type of experiments are cyclic uniaxial strain-controlled tests, uniaxial ratchetting tests and non-proportional cyclic biaxial strain-controlled tests at different temperatures. FE models to simulate brake rig tests are introduced and analyzed in Sect. 4. In particular, the influence of the viscoplasticity modelling is examined and the obtained results are briefly discussed. In this section, the adopted viscoplasticity model is described. The model is formulated in the small strain framework and the strain is assumed to be additively decomposed into an elastic strain e , a viscoplastic strain p and a thermal strain th : The thermal strain is assumed to be isotropic and is given by th ¼ a ðT À T 0 Þ I, where a is the thermal expansion coefficient, T is the temperature, T 0 is the initial stress-free temperature, and I is the 2nd-order identity tensor. The elastic response is assumed to follow the isotropic Hooke's law such that the stress r can be decomposed into a volumetric r vol ¼ I : r and a deviatoric r dev ¼ r À r vol =3 I part: where K b is the bulk modulus and G is the shear modulus. This can also be written using the 4th-order stiffness tensor E e : where I dev denotes the 4th-order deviatoric projection tensor. The evolution of viscoplastic strain is assumed to be of associative type and given by a yield function f together with an overstress function g: since the plastic multiplier _ k describes the amount of plastic strain rate and m ¼ signðrÞ. Moreover, t Ã is a material parameter. The yield function f is assumed to be of von Mises' type, i.e.
where X is the kinematic hardening stress, R is the isotropic hardening stress and k the initial yield stress. This choice of yield function f results in m, introduced in (4), as follows: In this paper, the overstress function g is described by the proposal in Delobelle [27] using a sinh function, i.e.
gðf Þ ¼ sinh f D n ; ð7Þ where D denotes the drag stress and the Macauley brackets Using Macauley brackets \ [ in the sinh function will result in an elastic stress regime where no viscoplastic strain will develop. The evolution of the isotropic hardening R is assumed to be of nonlinear type [8]: where Q is the saturation hardening and b affects the initial hardening slope. Clearly, this equation can be solved analytically (with initial condition Rð0Þ ¼ 0) resulting in the expression R ¼ Q ð1 À expðÀb kÞÞ.
The kinematic hardening X is assumed to be additively decomposed into several hardening stresses X i , i.e. X ¼ P n X i¼1 X i . This allows the material model to capture more sophisticated cyclic hardening characteristics. The evolution equations for the kinematic hardening stresses X i are based on the Armstrong-Frederick model [16] with a linear and a nonlinear (dynamic recovery) part. However, we further extend the hardening by including the Burlet-Cailletaud model [21] for multiaxial ratchetting and the static recovery effect [13,15]. The resulting evolution equation is formulated as where C i controls the initial hardening, c i controls the saturation hardening, b i is the multiaxial loading parameter, and parameters s i , m i and M i control the static recovery of the hardening. Note that for proportional loading (and the special case uniaxial stress loading), then the parameter b i does not affect _ X i . The static recovery term represents slow diffusive relaxation processes in the material.
An important restriction of material models is that they should fulfill the 2nd law of thermodynamics in order to be physically sound, see, e.g. Lemaitre-Chaboche's book [28].

Time integrated formulation
The numerical implementation of the model has been conducted by adopting a time integration scheme of backward Euler (BE) type from time t n to time t nþ1 . The assumed input to the constitutive model in each time step is the strain and thermal strain th together with old state variables ( n k, n X i ). For simplicity, we denote quantities at time t n by n and quantities at time t nþ1 without any superindex. By adopting BE on the time integrated form of Hooke's law in Eq. (2), we obtain the unbalance R r : where D denotes change between time t n and t nþ1 . Furthermore, BE integration of the definition of _ k in Eq. (4) gives where it has been used that k ¼ n k þ Dk and its relation to the isotropic hardening R ¼ Q ð1 À expðÀb kÞÞ. Note that 0 in Eq. (10) represents a null vector, i.e. (½0; :::; 0 T 1Â9 ), while 0 in Eq. (11) is a scalar value. Finally, also the evolution law for kinematic hardening stress is BE integrated. Here, we make a simplification, and the static recovery term that represents slow physical processes is integrated using a forward Euler integration scheme. The purpose is to reduce the nonlinearity of the resulting equation system and thereby increase the robustness and efficiency of the numerical implementation. The resulting unbalance equations become The total equation system with unknowns (for a given D and DT) can be summarized as where 0 represents a null vector (½0; :::; 0 T 1Âð9þ1þ9n X Þ ). The adopted algorithm to solve these equations in each time step is summarized as follows: The stress r and hardening stresses X i can be computed as The algorithmic tangent stiffness becomes (equal to standard Hooke's law stiffness) If f ðr; X; 0Þ\0 then STOP else GOTO 2.
(2) Solve the full equation system R tot ðX tot ; D; DTÞ ¼ 0 with a Newton iteration scheme. After that the solution is obtained, the algorithmic tangent stiffness can be derived using (see, e.g. [14]): Hence, the algorithmic tangent stiffness can be computed as part of We used the automatic differentiation software Acegen [29] to produce the code for the expressions of the derivatives oR tot =oX tot and oR tot =o. The algorithm has been implemented as a user material subroutine (UMAT) to the commercial FE software Abaqus [30].

Uniaxial fatigue experiments
To characterize the cyclic behaviour of the ER7 wheel material, low cycle fatigue tests (using the standard ASTM E606-12 [31] as a guideline) have been performed in terms of uniaxial strain-controlled experiments with hold-time and uniaxial ratchetting experiments at elevated temperatures. The specimens were extracted from unused wheels from around 15 to 20 mm below the wheel tread surface. The strain-controlled experiments are conducted with fully reversed (R ¼ À1) strain control with a triangular time history in [5]. In the tests, the strain rate is 5 Â 10 À3 s À1 except for the two initial cycles which are conducted at the strain rate 5 Â 10 À4 s À1 . The tests that we will use were performed at isothermal conditions from room temperature 20 C up to 600 C with the strain amplitude 0:6% at 20 C and the strain amplitude 1% for the remaining temperatures. All the tests, except the test at 20 C, included a 30 min hold time at constant compressive strain after a number of cycles to show the slow viscous (diffusion dominated) behaviour of the material in terms of stress relaxation. More details and experimental results can be found in [5,32].
While the cyclic strain-controlled experiments with hold-time are used to identify the cyclic and the low loading rate behaviour of the ER7 material, in order to characterize the ratchetting behaviour of the ER7 material at medium to high loading rates, isothermal ratchetting tests were performed in [9]. The specimens were subjected to cyclic loading with r m AE r a ¼ 100 AE 500 MPa where r m is the mean stress and r a is the amplitude stress. A sinusoidal waveform was used for the stress control. After 5 initial slow cycles at 0.5 Hz, each block of hundred cycles contained 15 rapid cycles at frequency 5 Hz and then 85 slower cycles again at frequency 0.5 Hz (which result in _ % 3 Â 10 À3 À 9 Â 10 À2 s À1 ). The tests were performed in isothermal conditions at elevated temperatures including 200, 300, 400 and 500 C, see Figs. 1 and 2 .

Multiaxial fatigue experiments
In order to enhance the understanding of the material behaviour under more complex loading conditions, similar to multiaxial loading conditions occurring during wheelrail interaction, biaxial strain-controlled experiments with alternating tension/compression and torsion loading are conducted in [26] on thin walled test specimens. A non-proportional loading scheme with a 90 phase lag was applied during the experiments. Our tests were performed in isothermal conditions at 20, 200, 300 and 400 C with a constant total equivalent von Mises strain amplitude loading of D t =2 ¼ 1% (see Figs. 3 and 4 ). A sinusoidal waveform with the frequency 0.25 Hz was used in all the tests except for the room temperature test which was run with the frequency 0.1 Hz.

Identification of uniaxial model parameters
It is observed from Figs. 1 and 2 and in [5] that the material exhibits several phenomena such as cyclic hardening, the Bauschinger effect, ratchetting, dynamic strain aging (associated with dislocation motion hindered by interstitials) and other diffusion processes including material degradation (mainly spheroidisation of cementite lamellae). In order to identify the parameters of the material model associated with the observed phenomena, an optimization procedure is performed. The optimization (minimization) is conducted using the matmodfit program ( [33]) by employing the Nelder-Mead algorithm [34] to minimize an objective function E obj . The objective function for a straincontrolled experiment is defined by the difference between experimental stress r exp and simulation stress r sim from time t 0 to t N r . The explicit expression of the objective function is in such a case formulated as follows: where the experimental scaling factor l r is chosen as ðmax i ðr exp ðt i ÞÞ À min i ðr exp ðt i ÞÞÞ À1 , the number of experimental data points is N r þ 1 and the time increments Dt i are chosen as t iþ1 À t i . For a stress-controlled experiment, the objective function E obj; is defined analogously, and when both strain-and stress-controlled experiments are considered in the parameter identification, then the objec- To increase the probability of finding the global minimum of the objective function and to devise a systematic approach, the model calibration strategy is divided into two main steps. In the first step, parameter values are identified for selected experimental data with the purpose to find good initial guesses for the second step of the strategy. In the second step, the parameter identification process is continued using all the available experimental data. However, as can be seen in Table 1 Modelling of temperature and strain rate dependent behaviour of pearlitic steel... 367 at 325 C, while the ratchetting and biaxial tests are conducted at 300 C. Despite this we have, for simplicity in the parameter identification, assumed that the behaviour of the material at 325 and 300 C for the strain-controlled tests are sufficiently similar. Hence, the strain-controlled experimental data at 325 C are used for parameter identification of the model at 300 C as is written in Table 1.
To find out initial guesses of Young's modulus (E), the yield stress (k) and the hardening parameters (C i ; c i ; b; Q) for each temperature, the cyclic part of the strain-controlled test data with strain rate _ ¼ 5 Â 10 À3 s À1 (excluding the two initial cycles and the hold-time) together with (when available) the 5 initial slow cycles with frequency 0.5 Hz (which results in _ % 3 Â 10 À3 À 9 Â 10 À3 s À1 ) of the ratchetting test data were used. In this study, two kinematic hardening stresses are employed, where the evolution of the second kinematic hardening is at this stage chosen to be linear which implies that c 2 % 0. This kinematic hardening is used to achieve a better agreement with the strain-controlled test data in terms of maximum and minimum stress, while the elastic-plastic transition of stress-strain curves is well captured by the first kinematic hardening. Also, it was found at this stage when studying a limited number of cycles, the cyclic hardening-softening expressed by isotropic hardening is relatively low (see also [7]). Consequently, the isotropic hardening is assumed to be linear for all temperatures, i.e. b ! 0 (in the code 10 À6 ) and Q ¼ 1.
The linear isotropic hardening is controlled by the product b Q (cf. Eq. (8)). The results in terms of parameter values are shown in Table 2 and examples of model predictions (after calibration) together with experimental results at 200 C are shown in Fig. 5.
Next in the procedure to find good initial guess of the parameter values, we identify the parameters t Ã and exponent n in Delobelle overstress function. To do so, the strain-controlled test data are extended by including the two slow initial cycles ( _ ¼ 5 Â 10 À4 s À1 ) as well as the hold time, and the ratchetting test data (when available) are extended by including the 15 rapid cycles (5 Hz). By these extensions, the strain rate range in the test data covers 5 Â 10 À4 À 10 À1 s À1 . It was found difficult to identify both t Ã and D in the Delobelle overstress function since they influence the response in a similar way, and hence, the value D ¼ 1 MPa was chosen for simplicity, whereas t Ã was identified. Furthermore, it was found that also the yield The mark U denotes that the experimental data are available for a specific temperature and -denotes that the experimental data are unavailable stress k had to be re-calibrated at this stage to achieve a good fit with experimental results. The obtained parameter values are given in Table 3. Note that since the ratchetting test data are not available at 600 C, the values of t Ã and n cannot be identified against medium to high strain rate test results. Instead, the t Ã and n values at 600 C were kept the same as at 500 C and only the parameter value of the yield stress k was identified. With these parameter values, the predictions of the (slow) relaxation in the strain-controlled experiment become poor. Therefore, the static recovery is also activated in the model. From the model formulation, it can be concluded that not both s i and M i can be identified.
We have chosen M i ¼ 1 MPa and identified s i . During calibration, it was found that changing the parameter value for m i did not influence the quality of the predictions. Instead, it was found sufficient to vary the s i parameters, and hence, m i ¼ 2 was chosen for simplicity. The obtained values of s i are given in Table 4, and numerical results of the strain-controlled experiments for 500 and 600 C are shown in Fig. 6. It can be observed that the relaxation is now well captured for these temperatures.
In the second main step of the identification strategy, all the available experimental data are used to identify the model parameters at the different temperatures. More specific, the difference from the previous step is that all the cycles in the ratchetting experiments are considered. To be able to capture the ratchetting behaviour for the increased number of loading cycles satisfactorily, the hardening modelling was extended to allow for nonlinear isotropic hardening as well as nonlinear hardening of both kinematic hardening variables. The results in terms of parameter values are shown in Tables 5 and 6 and model predictions after the calibration together with experimental results are shown in Figs. 7 and 8 . Although a more complex modelling of the hardening has been used, it should be noted that the calibration against all cycles in the ratchetting experiments has the drawback that the accuracy of the simulations of the initial ratchetting cycles decreases. So, if the focus of the simulations is the initial loading cycles then it can be well motivated to adopt the model parameters in Tables 2, 3 and 4 . However, in the brake rig experiments in Sect. 4, the focus is on the material behavior during many loading cycles, and therefore, the parameter values in Tables 5 and 6 are used. For temperatures in between the temperatures with available   experimental data, linear interpolation of the parameter values is adopted.

Simulations of strain rate dependence
It is worth mentioning again that the material model parameters are identified from experiments with strain rate in the range _ ¼ 5 Â 10 À4 À 10 À1 s À1 . However, the simulations of brake rig experiment using the plastic material model in [7] suggest a range of strain rate _ ¼ 10 À7 À 100 s À1 resulting from the applied thermal and mechanical loadings. Figure 9 shows the predictions of the viscoplasticity model in terms of monotonic stress strain curves at the temperatures 200, 400 and 500 C, employing the parameter values from Tables 5 and 6 . The predictions are assessed for the strain rates _ ¼ 5 Â 10 À4 s À1 (the lowest cyclic strain rate that the viscoplasticity model is calibrated against) and _ ¼ 100 s À1 (the highest expected strain rate in the simulation of brake rig experiments). We can observe how much the flow stress is increased depending on the strain rate, i.e. Dr u ð _ Þ, and that the rate dependence increases with increasing temperature. It can be noted that the plasticity model assumes the same material behaviour as in the slow LCF strain-control tests. Because of extrapolation to higher strain rates than the model was calibrated for, too high stresses are predicted for the higher temperatures 400 and 500 C at the strain rate _ ¼ 100 s À1 . A comparison can be made with results from high strain rate tests at elevated temperatures for different steels in [35]. To improve the reliability of the predictions for higher strain rates, the model has to be calibrated against representative experimental data or estimates.

Identification of multiaxial model parameters
As the final step of the parameter identification, the experimental data obtained from the biaxial tests are employed to identify the multiaxial hardening parameters b i in the Burlet-Cailletaud model, cf. Eq. (9). The parameter identification is performed by minimizing a straightforward extension of the objective function in Eq. (14), i.e. E obj ¼ E obj;r þ E obj;s . It was noticed during   Modelling of temperature and strain rate dependent behaviour of pearlitic steel... 371  Table 7 and examples of model predictions together with the experimental results at 200 and 400 C are shown in Fig. 10. We adopt linear interpolation also for the temperature dependence of b 1 and assume that it is temperature independent above 400 C.  Finally, since uniaxial tests do not give any information about the transversal behaviour of the material, Poisson's ratio m in the model cannot be identified from the available experimental data. In this paper, we assume that Poisson's ratio m and its temperature dependence are similar for ER7 as for mild steels. The values for mild steel from [36] are given in Table 8.
It is noteworthy to mention that the material parameter values are significantly temperature dependent, as can be seen in Tables 5, 6, and 7 . To get as accurate simulations of the brake rig tests (in Sect. 4) as possible, we have chosen to use all experimental data for calibration of the model. If more experimental data would have been available to us, then it would be suitable to use some of those for model validation.

Results and discussion
In this section, the influence of viscoplastic material modelling on simulations of RCF brake rig experiments is investigated. The brake experiments have been conducted to study repeated stop braking using a full-scale brake dynamometer featuring a rail-wheel, cf. [37]. Each stop braking test starts by enforcing the prescribed wheel-rail contact force after which the rig is accelerated to the chosen test speed. The stop braking tests have been performed with a sinter brake block with initial velocities of 160 and 130 km/h as well as with an organic composite brake block with initial velocity of 160 km/h. The brake block is then pressed towards the wheel tread with a prescribed brake force until the wheel stops. Finally, this is followed by cooling combined with a low rail-wheel contact force. Videos acquired by thermocamera reveal that the rise in temperature on the tread in some cases is not uniformly distributed over the block-wheel contact surface and it is found that hot bands occur on the tread which can be explained by the phenomenon known as frictionally induced thermoelastic instability. The experimental set-up is described in detail in [37] and the corresponding simulation methodology in [7]. The main components in the simulation methodology are as follows: (1) Axisymmetric thermal analyses of wheel taking into account heat partitioning between brake block, wheel and rail. In this paper, we limit the study to the case of one hot-band (50 mm) between the wheel and brake block with initial velocity 160 km/h for a sintered brake block (case S160).
(2) Detailed 3D contact FE-simulations (without friction) between wheel and rail to obtain contact patches and pressure distributions for different temperatures and strain rates. (3) Given temperatures (from step (1)) and contact loading (from step (2)), 3D FE simulations of braking cycles with a representative temperature history and a chosen number of wheel-rail contact load passages.
In steps (2) and (3) of this simulation methodology, the temperature and strain rate dependence in the viscoplasticity model are important features. In the results below, we particularly focus on the difference in terms of equivalent von Mises strain between the plasticity model presented in [7] and the viscoplasticity model with Delobelle overstress function presented in Sects. 2 and 3. Figure 11 shows the evolution of ratchetting for the plasticity model, in terms of equivalent von Mises plastic strain p eff , during one braking cycle simulation at a point 1.4 mm below the center of wheel-railwheel contact. It is noteworthy that based on the investigations in [7] the brake case S160 with 1 band 50 mm width results in highest temperatures and ratchetting strains among all 9 combinations of brake cases and banding patterns. The corresponding results for the viscoplasticity model are shown in Fig. 12.
From the results in Figs. 11a and 12a, it can be seen that the ratchetting predicted by the viscoplasticity model during the first 10 rolling passages at 20 C is lower than the ratchetting predicted by the plasticity model. The overprediction of the ratchetting by the plasticity model is expected since it does not account for strain rate hardening. The decrease in predicted ratchetting is even more clear during the braking cycles (see Figs. 11b and 12b) which shows the necessity to account for viscous effects in modelling.
To further highlight the influence of the viscoplastic material model in the simulations, an assessment of fatigue life based on ratchetting failure (RF) can be performed. In order to account for RF, the empirical model proposed by Kapoor [38] is adopted which accounts for fatigue crack initiation due to ratchetting. In this method, the ratchetting life N r is obtained by dividing a critical ratchetting plastic strain c that the material can sustain by the total ratchetting strain for the entire braking cycle D r (for details see [7]): The value of total ratchetting strain employing the plastic material model gives D r ¼ 0:13, while this value decreases to D r ¼ 0:01 employing the viscoplastic material model. Note that the fatigue life prediction in the Kapoor model is controlled by critical strain c which in [39,40] was assessed experimentally by twin disc tests. Since such data are not available for ER7 material, in [7] we estimated this value as c ¼ 1:6 by comparing the fatigue life predictions from numerical analyses with the results from brake rig experiments. If we repeat that procedure for the viscoplastic material model the value of critical strain decreases to c ¼ 0:1 which is in the same order of magnitude with the experimental value obtained for mild steel in [38,41] and the experimental results for the uniaxial ratchetting presented in Figs. 1 and 2 . However, from these experimental results, it can be noted that the critical strain c is in general temperature dependent which is not accounted for in our estimations. Finally, we evaluate the sensitivity of the multiaxial Burlet-Cailletaud parameters b i . In this regard, the same simulations with the viscoplasticity model are now repeated but with the identified values of b i in Table 7 against the multiaxial fatigue tests (instead of b i ¼ 1). The results from the simulation are given in Fig. 13. One can observe, when comparing with Fig. 12, that the results differ at the room temperature cycles but that difference is decreased during the braking.

Concluding remarks
To simulate the material behaviour of railway wheel steels ER7 subjected to combined tread braking and rolling contact loading, a viscoplastic material model is proposed and calibrated against experimental data. A main challenge is the large difference in loading rate between the very fast rolling contact loading and the slow thermal loading caused by braking. In this regard, the viscoplasticity model was calibrated against LCF uniaxial strain-controlled tests with hold-time (to capture slow viscous diffusion dominated behaviour) at elevated temperatures and also against uniaxial ratchetting tests including rapid cycles at elevated temperatures. It is observed that for strain rates in the range ( _ ¼ 5 Â 10 À4 À 10 À1 s À1 ), the rate dependent behaviour can be successfully modelled by the Delobelle overstress function. To also capture the slower relaxation behaviour of the material, it was found necessary to include a static recovery term of the kinematic hardening in the model. Biaxial strain-controlled cyclic tests of ER7 at elevated temperatures have been conducted with the purpose to find its behaviour at loading situations more similar to rolling contact loading. In this context, the kinematic hardening is modified by using the Burlet-Cailletaud model, and the (a) (b) Fig. 11 Plasticity model. Case S160, 1 band 50 mm. Evolution of the equivalent von Mises plastic strain p eff at 1.4 mm below the center of wheel-railwheel contact a during the first 10 rolling passages at room temperature and b during the braking cycle. Solid circles indicate the beginning instance of every rolling passage (a) (b) Fig. 12 Viscoplasticity model when biaxial parameter is set to b i ¼ 1. Case S160, 1 band 50 mm. Evolution of the equivalent von Mises plastic strain p eff at 1.4 mm below the center of wheel-railwheel contact a during the first 10 rolling passages at room temperature and b during the braking cycle. Solid circles indicate the beginning instance of every rolling passage biaxial test data were used to identify the corresponding parameters. To summarize, for the temperatures and strain rate range which the model was calibrated against it can mimic the viscous, cyclic inelastic and multiaxial behaviour of ER7 steel very well.
The capabilities of the material model are challenged by simulating a full-scale test rig experiment reported in [7]. It is observed from the simulation results that the predicted ratchetting strain using the viscoplastic material model during a braking cycle is significantly lower than the corresponding results from using a plasticity model in [7]. This highlights the importance of accounting for viscous effects in the material model. However, for this demanding case including very high strain rates (in the order of 100 s À1 ), the viscoplasticity model does need additional calibration with supplementary sets of experimental data. Further, it is observed that including the effect of multiaxial hardening in the material model does not show a visible effect in the predicted ratchetting during a braking cycle.
The simulations with the viscoplasticity model could in future work be further exploited and validated by, e.g. repeating the simulations of the test-rig tests in [7] for different brake block materials and initial velocities. The utilized fatigue life criterion based on the critical ratchetting strain needs to be further investigated. In particular, its influence of temperature, loading rate and multiaxial conditions need to be addressed in order to ensure the quality of the life predictions with the criterion.  Table 7. Case S160, 1 band 50 mm. Evolution of the equivalent von Mises plastic strain p eff at 1.4 mm below the center of wheel-railwheel contact a during the first 10 rolling passages at room temperature and b during the braking cycle. Solid circles indicate the beginning instance of every rolling passage