Computationally efficient simulation methodology for railway repair welding: Cyclic plasticity, phase transformations and multi-phase homogenization

Abstract The in-situ railway repair welding process consists of multiple weld passes, which makes it significantly different from other rail welding processes. In this study, finite element simulations of repair welding are performed to predict the resulting microstructure and residual stresses. To accurately simulate the material behavior, the modeling includes phase transformation kinetics, cyclic hardening plasticity, transformation induced plasticity, and multi-phase homogenization. More specifically, four different homogenization methods are investigated: isostrain, isostress, self-consistent and linear mixture rule. The performance of the material modeling is demonstrated by simulating multiple weld passes using a classical three-bar welding experiment. Based on the results, the self-consistent method and linear mixture rule are used in a 3D full-scale railhead repair weld simulation, in which the former generates a more realistic mechanical response. The immense computational cost associated with 3D full-scale, full-detail multi-pass welding simulations is addressed by exploring different model reduction schemes. From this study, a 2D generalized plane strain model, extended with out-of-plane axial and bending stiffness, is found to replicate the full-scale model at a mere fraction of the computational cost. Finally, the longitudinal residual stress distribution obtained from the reduced model is shown to correlate well with experimental measurements.


Introduction
In the last century, progress in railway technology has led to a significant decrease in train wheel and axle failures, but the occurrence of rail failures has, despite devastating consequences, not decreased to the same extent, see Smith [1].It is argued that, depending on the types of rail and traffic, between 25% and 70% of rail failures occur at welds, see e.g., Cannon et al. [2], Kondo et al. [3] and Meissner [4], confirming that welds act as a weak link in the rail system.Welding in railway rails is typically performed for two reasons: joining the rails when laying the track and repairing damaged parts of the railhead.In Sweden, when laying the track, rail sections are usually constructed in a stationary plant by joining several rail pieces using flash butt welding.These sections are then transported to the field and joined by alumino-thermic (thermite) welding, see Skyttebol and Josefson [5].Once in the field, repairing surface damage on the railhead is generally performed using fusion welding to, layer-by-layer, fill the gap after the damaged section has been removed.Common for all three welding procedures is that the localized heating generates rapid cooling rates in the weld and surrounding material, which may cause phase transformations and result in varying hardness and toughness of the rail surface as well as considerable residual stresses, see e.g., Porcaro et al. [6] and Allie et al. [7].Important to note is that these residual stresses differ significantly from the preceding manufacturing residual stresses, see e.g., Mutton and Alvares [8] and Tawfik et al. [9], and may have detrimental effects on rail performance.Moreover, rail welding may also introduce additional complications: surface irregularities that magnify dynamic loads, see e.g., Wen et al. [10] and Li et al. [11]; subsurface porosity and inclusions defects which deteriorate the mechanical fatigue properties, see e.g., Feddersen et al. [12]; and a microstructure with varying grain size, see e.g., Mansouri and Monshi [13].With these aspects taken into account, Tawfik et al. [14] points out that continuously welded rail is discontinuous in terms of microstructure, mechanical properties and residual stresses.Altogether, there is consensus in research that welds are the Achilles heel of rails, and therefore rail welding has been a prominent research topic for the last decades.As demonstrated by the rail failure statistics, it is still highly relevant.
In improving the quality of railway welds, numerical simulations have proven to be a valuable tool for understanding the driving mechanisms of the residual stress field and how residual stresses contribute to crack nucleation and crack growth.For example, numerical analyses of the residual stresses from flash butt welding, see e.g., Skyttebol and Josefson [5] and Tawfik et al. [14,15], and thermite welding, see e.g., Josefson et al. [16] and Liu et al. [17], indicate high tensile residual stresses in the rail web for both the base material and the weld, balanced by compressive stresses in the head and foot of the rail.Perhaps more importantly, numerical studies have also shown that the tensile stresses in the web are not redistributed during typical operational loads, see e.g., Fang et al. [18], and may therefore contribute to the growth of subsurface horizontal fatigue cracks initiated by internal defects, see e.g Skyttebol et al. [19], Josefson [20] and Lee et al. [21].To mitigate residual stress driven defects, the accuracy of numerical tools has improved by including aspects such as phase transformations, see e.g., Cal et al. [22] and Ma et al. [23], as these strongly affect the material properties and cause cycling straining despite monotonic heating or cooling, see e.g., Rammerstorfer et al. [24].
However, when simulating the repair welding process, the requirements of the numerical tools are higher than for thermite or flash butt welding, since the evolution of the stress field is more complex.Instead of one intense heat pulse as in thermite or flash butt welding, the repair welding process involves multiple, weaker heat cycles from several weld passes, each affecting a smaller region, which gives higher cooling rates.Further adding to the complexity, previous weld layers and the surrounding region are reheated and possibly remelted by each additional weld pass, causing more pronounced cyclic straining and a greater risk of martensite formation and varying surface hardness, see e.g., Zheng et al. [25].Although efforts to numerically simulate rail repair welding mechanically are far less common in the literature than for thermite or flash butt welding, three noteworthy examples are Jun et al. [26], Kabo et al. [27], and Lee et al. [21].The former illustrates the risk of martensite formation, and the latter two indicate that operational rolling contact loads redistribute the residual stresses at the surface and magnify the stresses beneath the surface of the repaired rail.To complement these studies, we present a repair welding simulation methodology combining three main improvements: (i) a more advanced material modeling presented in previous works [28,29], which includes explicit modeling of phase transformation kinetics, multi-phase homogenization, cyclic hardening plasticity and transformation induced plasticity (TRIP), (ii) explicit modeling of a moving Goldak heat source [30], (iii) and modeling continuous addition of weld filler material.Compared to previous studies, explicit modeling of phase transformation kinetics, including solid-to-solid and solid-to-liquid-like transformations, and using cyclic hardening plasticity allow for in-depth focus on the cyclic nature of the residual stress evolution.Furthermore, investigating how different homogenization methods account for the multi-phase stages during the simulation adds further insight into the complex evolution of the residual stresses and highlights the importance of the choice of homogenization method.
The improvements of the numerical tools presented in this study add needed details to the weld simulation in terms of physical phenomena, material modeling and practical weld process aspects.However, these improvements also add computational cost, seemingly to the point where practical use (with today's performance of computers) can be questioned if additional measures are not taken.It is well established that detailed simulations of wire arc multi-pass welding require modeling of multiple overlapping, interacting physical phenomena, see e.g., Runesson et al. [31].Equally established is that many of these interactions are negligible compared to the interactions which couples thermal, structural and metallurgical processes, Lindgren [32] refers to these as classical computational welding mechanics (CWM).Nevertheless, simulations in typical CWM level of detail still require tremendous computational power due to the large scale of the full rail repair welding process.Most likely, this is why previous studies have used less detailed approaches compared to CWM simulations in other research fields.In an attempt to circumvent the issue of computational cost, this work presents a model reduction study which focuses on obtaining accurate residual stress results for the center cross-section of the weld repair.The simulation results are validated using the residual stress measurements by Jun et al. [26,33].Using reduced simulation models, the three main improvements, i.e., more advanced modeling of material behavior, heat source movement and addition filler material, enable the interacting phenomena in CWM to be accounted for at reasonable computational cost and with accurate results.Altogether, the numerical tools presented in this study can be used as viable aid in optimizing repair weld process parameters, thus adding further insight into how to achieve a successful and effective repair.

Preliminaries
The study presented in this paper regards railway repair welding, i.e., multiple heating cycles, of an initially fully pearlitic rail steel R260 [34] with an approximate chemical composition Fe-0.72C-0.3Si-1.0Mn.Thus, the material has a slightly hypo-eutectoid composition, but still exhibits a fully pearlitic microstructure after production.When simulating repair welding for this material, the following possible temperature driven phase transformations are considered: 1. Pearlite into austenite during heating 2. Austenite into martensite during high cooling rates 3. Austenite into pearlite and/or ferrite during cooling 4. Austenite into bainite for moderate cooling rates 5. Tempering of martensite and bainite during heating 6. Martensite into austenite during heating 7. Melting and solidification Note that this paper considers hypoeutectoid steels only, whereby cementite phase transformations are not considered.Throughout the modeling in this the paper, pearlite, bainite and tempered martensite are referred to as phases even though they are in fact microstructures.

Phase transformation kinetics
The adopted modeling of the kinetics of the phase transformations 1 to 6 is described in detail in previous works, see [28,29], and will here only be summarized briefly.To adapt the modeling to the repair welding process, a liquid-like phase is introduced and also phase transformations to and from this liquid-like phase, see phase transformation 7.Moreover, using the same nomenclature as in previous works [28,29], the following indices are used to denote the phases: austenite (a), pearlite (p), ferrite (f), bainite (b), martensite (m), tempered martensite (tm) and liquid-like phase (l).During the heating cycles, the sum of all phase volume fractions, p x , must equal 1: To describe the kinetics of the first phase transformation, austenitization of pearlite, an IT-diagram is used.Where for a constant temperature, the transformation kinetics of the decreasing pearlite volume fraction is described by the Johnson-Mehl-Avrami-Kolmogorov (JMAK) equation see e.g., [35,36].To handle the varying temperatures, the Scheil's additive rule [37] is adopted see e.g., [38].
The kinetics of transformations 2 and 3 are based on a CCT-diagram.Transformation 2 is the formation of martensite during rapid cooling of austenite, i.e., quenching.It is assumed to be diffusionless and is modeled by using the purely temperature dependent Koistinen-Marburger [39] equation.To improve the numerical performance of the finite element solver, a smoother evolution of very low martensite volume fractions is adopted by replacing Koistinen-Marburger equation with a hyperbolic tangent function to compute the volume fraction during the initiation of martensite.
Under low to moderate cooling rates, austenite transforms into ferrite, pearlite, and/or bainite which is the fourth phase transformation considered in this paper.As these transformations are of diffusive nature, we again use the JMAK equation and Scheil's additive rule to describe the transformation kinetics.
The solid to liquid and liquid to solid phase transformations are modeled by linear development of the liquid-like phase p l ðTÞ as follows, cf.[40][41][42]: where the temperatures T l0 and T l1 denote start of melting and when a fully liquid phase is obtained, respectively.Experimental data of the mechanical behavior at near melting temperatures is scarce (compared to data available at lower temperatures), whereby a simulation cutoff temperature can be motivated and generally does not impair the fidelity of the welding simulation outcome, see e.g., Lindgren [43] and Dong et al. [44].Therefore, the governing temperature limits T l0 and T l1 in Eq. (2) are set to temperature levels below the physical melting temperature for the simulated liquid-like phase to serve the same purpose as a simulation cutoff temperature.By this modeling and for this purpose, rapid heating events result in solid to liquid phase transformations without preceding full austenitisation.The initial melting temperature is set to T l0 ¼ 1200 � C and the fully melted temperature to T l1 ¼ 1500 � C. The motivation for implementing a liquid-like phase explicitly in the material model is twofold; i) a stress free state is obtained for p l ¼ 1 (pure liquid-like phase), ii) the internal state variables for all phases are reset (i.e., set to zero) without affecting the stability of the numerical finite element solver.Whereby the resulting effect is that a virgin material state is obtained during the liquid-like to solid phase transformation.This is an important aspect of the simulation model to capture the correct stress-strain behavior during melting and cooling, see e.g., Dong et al. [44].However, the internal variables are not reset for a partially liquid like state, i.e., 0 < p l < 1:

Homogenization methods
The constitutive behavior of the multi-phase steel is assumed to be obtained from homogenization of the behavior of the individual phases.Four homogenization methods will be adopted: isostrain (Voigt assumption), isostress (Reuss assumption), self-consistent and the linear mixture rule.The implementation of these are presented in a previous work, see [29], and will only be summarized briefly in the following section.For all methods we assume that the temperature T is the same in all the phases.Further, the methods are implemented by using an incremental straindriven algorithm.The state from the previous time step n t is assumed to be given in terms of homogenized strain n � � and stress n � r as well as strain n � x , stress n r x and state variables for all the phases.Then, a time increment with a strain increment d� � is applied whereby the updated homogenized strain can be computed as � �¼ n � � þ d� � but the increments d� r, d� x and stress dr x are determined by the chosen homogenization method.
During a time increment of the mechanical problem it is assumed that the temperature and phase fractions remain constant.Hence, the relation between the homogenized strain increment and strain increment of the phases can be written as: and similarly between the stress increments d� r ¼ Using the isotrain (Voigt) assumption, the strain increments of all the phases d� x are assumed to be equal to the homogenized strain increment, i.e., d� x ¼ d� �: Whereas for the isostress (Reuss) assumption the stress increment dr x is assumed to be the same in all the phases and thereby equal to the homogenized stress, i.e., dr x ¼ d� r: However, since the model will be used together with FEM, our numerical implementation of the cyclic plasticity model is based on a strain controlled algorithm (given d� �).This means that the isostress method requires an additional Newton iteration scheme to find each phase's individual strain increment, d� ¼ d� 1 , :::, d� n x ½ � T , such that the isostress criteria is fulfilled.This procedure has several similarities to a 2D plane stress algorithm.Moreover, by adopting the self-consistent framework described in [45], the strain increment in each phase d� x is governed by a fourth order concentration tensor A x which is computed using the global algorithmic stiffness tensor � E, the algorithmic stiffness tensor of the current phase E x and the Eshelby tensor P: These are computed using the fixedpoint iteration technique described in [46,47] and thus requires a pre-processing step at the start of each time increment to determine the individual strain increment of each constitutive phase, d� x : However, the most used and most straight forward way of homogenizing multi-phase states in steel is to use the linear mixture rule, see e.g., [48] and [49].Using this method, the constitutive phases are not treated by individual material models, but as one nominal model where all thermal and mechanical properties and parameters are computed as the volume fraction governed average of phases.However, the underlying constitutive model for the linear mixture rule is the same as for the individual phases in the other homogenization methods.
In this study, the primary purpose of using intricate multi-phase homogenization can be exemplified by the case of a phase transformation with phases of different strengths.By homogenizing individual phase material models in the self-consistent method, the weaker phase can deform plastically while the other phase behaves elastically.This is not the case using the linear mixture rule, where the model has one nominal hardening, which can cause less physical plastic inheritance during the phase transformation.

Preliminaries
For the isostrain, isostress and self-consistent homogenization methods, the constitutive models of the individual phases are run in parallel.As mentioned in the previous section, the homogenization methods form the total constitutive model response for the mixture of the phases using the phase volume fractions of the current time step.Phases which have yet to materialize, i.e., p x ¼ 0, are assumed to not accumulate hardening or plastic strain, thus behaving linear elastic.As a new phase starts to nucleate during transformation from its parent phase, the stress-strain state enforced by the homogenization method is imposed and the behavior follows the elasto-plasticity model described below.However, as discussed in Subsection 2.1.2,during the liquid-like state (p l ¼ 1), all internal variables for the individual constitutive models are reset and the following solidification generates a virgin material state.Internal variables for phases reborn without melting having taken place, e.g., austenite transforming back into pearlite, the internal variables are, for simplicity, not reset.This explicit treatment of the internal variables for the individual phases does not apply to the linear mixture method as this method simply has one nominal constitutive model.

Kinematics and elasticity
For each phase x in the material, we assume that the total strain � x is additively decomposed as (see e.g., [50]): where � e x is the elastic strain, � th x the thermal expansion strain, � tv x the transformation strain, � p x the plastic strain and � tp x the strain due to transformation induced plasticity (TRIP).For each time increment, the strains � th x and � tv x and their increments d� th x and d� tv x are given as the temperatures at the start and end of the current increment is known.Note also that the TRIP strain, � tp x , is implemented only for the austenite, martensite and bainite phases (x ¼ a, x ¼ m and x ¼ b).
The elastic strain governs the stress and we assume linear isotropic elasticity by adopting Hooke's law: with the fourth order deviatoric identity tensor I dev ¼ I − 1=3 I � I, the fourth order identity tensor I and the second order identity tensor I 1 .Furthermore, the material parameters for elasticity are the shear modulus G x and bulk modulus K x, b : The stress can be decomposed into a deviatoric and a volumetric part r ¼ r dev þ 1=3 r vol I with: x, dev and r x, vol ¼ I : The thermal expansion strain is assumed to be linear isotropic: where DT is the temperature increase from a reference temperature T 0 and a x is the thermal expansion factor.To compute the transformation strain, � tv x , the initial density of the material is q 0 while for each phase the density is assumed to be q x .The initial density of the material is that of pearlite until the material melts and solidifies into a virgin material state, then the initial density is set to that of austenite.Using the densities and the conservation of mass, a phase transformation will lead to a change of volume which defines the transformation strain: The TRIP strain � tp x is caused by phase transformation under applied stress.The physics behind the mechanism is described in e.g., [51][52][53].In this paper, we allow for TRIP to occur during the phase transformation of austenite into martensite or bainite by adopting the formulation proposed in [49,54], with stiffness parameters K tp acquired by [50].To obtain the correct TRIP strain output after homogenization, the TRIP strain � tp x is implemented in each individual material model of the phases active in the transformation, see previous work [29].For the linear mixture rule, this is obtained by including the transformation stiffness parameters in the mechanical properties for austenite, martensite and bainite.In the martensite constitutive model, � tp m is adopted as follows: with For better numerical stability when the constitutive model is implemented and used in an FEsolver, we modify the saturation function f ðp m Þ used by [50,55] as follows (satisfying f ð0Þ ¼ 0 and f ð1Þ ¼ 1): Together with the smooth evolution of martensite volume fraction (as discussed in Subsection 2.1) during initial martensite nucleation (p m � 1) this saturation function allows for a gentle initiation of the TRIP strain evolution.In these expressions we have introduced the equivalent von Mises stress r m, e ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 3=2 r m, dev : r m, dev p ¼ ffi ffi ffi ffi ffi ffi ffiffi 3=2 p jr m, dev j: For the austenite during martensitic phase transformation, the TRIP strain model is chosen as follows: Note that the expression includes the martensite volume fraction p m and its evolution _ p m : This prevents TRIP strains to develop during austenitization.

Plasticity and hardening model
For each phase we adopt an individual Chaboche plasticity model, as proposed in e.g., [56], that includes the von Mises yield function, non-linear isotropic and kinematic hardening.The von Mises type yield surface is here defined as: where X x is the kinematic hardening stress (back-stress), R x is the isotropic hardening stress and r x, Y is the initial yield stress.This yield function is used to distinguish elastic and plastic response via the conditions: where _ k x is the plastic multiplier.It is when _ k x > 0 (and U x ¼ 0) that plastic strain and hardening variables evolve.The evolution equation for the plastic strain is assumed to be of associate type: The total kinematic hardening stress X x is obtained from adding n x kinematic hardening stresses X x, i each following [57] type of evolution law: Finally, the evolution of isotropic hardening is adopted as follows: The material parameter values are assumed to be temperature dependent and chosen according to [28], linearly extrapolated to the temperature range considered in the present study.

Liquid-like phase and virgin material solidification
For a material in the liquid-like state to be stress free and to impose a virgin material during solidification, a strain component � l x is added to the total strain as follows: The liquid strain is set to zero initially and then it is controlled by the phase fraction liquid p l according to: Hence, the strain component � l x produces a zero stress state during the simulated liquid-like state, p l ¼ 1, when used in Eq. ( 6).As this strain component is fixed during the solidification process, the solidification starts from a stress free virgin material state, similar to the proposal by Chiumenti et al. [58].Also, the reference temperature, T 0 , for the thermal strain, see � th in Eq. ( 5), is changed from room temperature to 1200 � C during this simulated solidification process.In this way contraction of the cooling material are accounted for and spurious stresses are avoided during the liquid-like state.
We model the liquid-like phase as a soft, elastic-ideal plastic solid, similar to Draxler et al. [59], using the constitutive model described in Section 2.3, hence the governing material parameters are Young's modulus, Poisson's ratio, and the yield limit.For melting metals, Poisson's ratio approaches � !0:5, cf.e.g., Greaves et al. [60], whereby we use � ¼ 0:48 as a compromise for numerical stability.For the same reason, to get a soft material response of the liquid-like phase, Young's modulus is set to 0.1% of its room temperature value, much like e.g., Banik et al. [61] and Derakhshan et al. [62].Similarly, the yield stress is set to 1% of that at room temperature, cf.e.g., Elcoate et al. [63].Moreover, to ensure that the liquid-quid like phase does not affect the surrounding material structurally, its thermal expansion coefficient is set to zero, i.e., a l ¼ 0: The material properties during the solid to liquid and liquid to solid phase transformations, i.e., the so called" mushy" state (0 < p l < 1), are obtained by the different homogenization methods described in Section 2.2.Note that it is only during the mushy state that the liquid-like phase has an elastic-ideal plastic material response, as the introduction of the strain component � l x produces a zero stress state for p l ¼ 1:

Three-bar experiment
The material models presented in Section 2 are demonstrated using a three-bar model, see Figure 1a, and thermal and mechanical properties from [28,29].This is a classic example to conceptually model the residual longitudinal stresses in a butt-welded plate, see Figure 1b [64,65].It is here used to demonstrate the ability of the material models to predict residual stress during the multiple heating cycles of multi-pass welding.The weld joining the plates and the HAZ are represented by the center bar, whereas the plates are represented by the two outer bars.The welding process is simulated by a prescribed temperature cycling in the center bar.The resulting expansion and contraction of the heated bar are influenced by the outer bars, similar to how the plates restricts the HAZ from expanding freely.
The simulation is performed using the commercial software Abaqus [66] together with user defined subroutines for the material models presented in Section 2. The simulation model consists of second order 20 node brick elements.All bars are constrained vertically (z-direction) at the bottom end and their vertical displacements at the top end are connected via a rigid body coupling, marked by green lines in Figure 1a.The heated center bar is free to expand in the xand y-direction.Based on the analytical reasoning presented by Dong [44], the areas of the nonheated bars are set to five times that of the heated bar.The prescribed temperature history of the heated bar is a series of linear pulses from room temperature to 1600 � C and back to room temperature in 10 s.Even though this cooling rate is more reminiscent of quenching than of welding operations, it is used for demonstration purposes to induce martensite formation.To demonstrate pearlite formation, the last cooling cycle is extended to 30 s. Two bottom graphs of Figure 2 present the temperature history prescribed to the center bar and the time history of the material phase volume fractions.These graphs show how the cyclic temperature load causes phase transformation from pearlite into austenite during the initial heating, how the material melts and solidifies, and then how the austenite is quenched to form martensite.The graphs also shows how the material undergoes identical phase transformations during the following heat cycle; martensite is tempered, tempered martensite is transformed to austenite, the material then melts and solidifies, and the austenite is quenched during the cooling process.The slower rate of the last cooling results in a mixture of pearlite and martensite.
The two top graphs of Figure 2 presents the axial (z-direction) strain and stress history of the heated bar.The distinct kinks of the curves are due to the limited selected temperatures used when linearly interpolating in the parameter data of the material models, see previous work [28].The graphs show the general behavior of the HAZ during welding; The initial temperature increase causes the heated bar to expand, generating positive axial strain.The cool bars prevent free axial expansion of the heated bar, which thereby are subjected to compressive stress.As the temperature increases, the heated bar softens, whereby the strained cool bars start to compress the heated bar slightly.At about 2 s into the simulation, before the first phase transformation takes place, the effect of the thermal expansion again becomes dominating and the axial strain increases.As the heated bar melts, the material is in a stress free state.In the following the solidification process, the cool bars prevent the heated bar from contracting and generates tensile residual stresses in the heated bar.Following the solidification transformation, the effect of the virgin material state is clearly seen as each new heat cycle generates an identical residual state.Typical longitudinal residual stress levels found in the welding mechanics literature, see e.g., Masubuchi [64], for this three-bar example are at the order of the yield limit.This is also what is seen in the results presented in Figure 2, although the resulting microstructure complicates the comparison.
Other than side by side comparison of the different homogenization methods, the simple three-bar simulation also illustrates the nature of the different methods, see Figure 2.During the quenching of austenite, the isostrain method imposes the parent phase strain of the austenite onto the stiffer martensite phase, generating high stresses.Whereas for the isostress method, the imposed stress of the austenite parent phase is within the elastic region for the martensite, which generates a drop in strain.Using the self-consistent method, the behavior is presumed to be more realistic as the individual phase strains are governed by the stiffness and concentration tensors of the transforming phases.Whereas for the linear mixture rule, the quenching of austenite, and pearlite transformation, reveal an inherent drawback as the plastic strain accumulated in the austenite phase is inherited by the martensite or pearlite, causing a sudden leap in its stress-strain response curves.In addition, the top two graphs of Figure 2 also illustrate how all homogenization methods generate identical results in the initial single phase state and in the virgin material states, but significantly different stress and strain histories during the phase transformations.Moreover, it also shows how the different plasticity and hardening behavior of the phases lead to that the homogenization methods generate different stress and strain histories also in the succeeding single phase material states.Thus, the homogenization methods result in different residual states.
To further demonstrate the variation in the structural response obtained using the different homogenization methods and to motivate the modeling of a liquid-like phase, the three-bar weld experiment is rerun with the prescribed temperature load lowered to 1000 � C. The heated bar axial stress obtained in this simulation is presented Figure 3. Here, the effect of not obtaining a virgin material state in the cycling heating causes the homogenization methods to accumulate different internal variable histories and the axial stresses starts to deviate.Concluding the three-bar experiment, we find the isostrain and isostress method unfit for the multi-pass welding simulations presented in this study.The former, mainly due to the unrealistic stress-strain response upon virgin material cooling, and the latter mainly due to the additional computational cost and numerical instability.Therefore, in the simulations presented in Section 4, we use the self-consistent homogenization method for its realistic behavior and the linear mixture rule as it is the most frequently used method.

Railhead repair welding simulation
We adopt a decoupling of the thermomechanical finite element analysis into a transient thermal analysis and a quasi-static structural analysis.The temperature field from the thermal analysis drives the metallurgical processes and the thermal expansions, which together with mechanical loading drive the response in the mechanical analysis.This setup does not account for plastic energy dissipation, latent heat during transformations, or stress induced phase transformations.Omitting these is commonplace in CWM as the effects are negligible in comparison the weld heat input and temperature induced phase transformations.The presented rail head repair welding simulations replicates the experimental work presented in [26], both in terms of process parameters and experiment setup.In accordance, the rail material is a pearlitic steel and the filler material a C-Mn-Si steel with a chemical composition given in [26].

Moving heat source
From the three-bar experiment, the next step toward a 3D welding simulation is to model the moving heat source of the wire arc welding process in the heat transfer simulation.Explicitly modeling the movement of the heat source is found to be significantly more accurate than modeling an instantaneous heating of the entire weld pass, see e.g., Kiyoshima et al. [67] or Pu et al. [68].When modeling the heat source movement, two common techniques are found in the literature; prescribed nodal temperatures, see e.g., Lindgren et al. [69], Seles et al. [70] and Lee et al. [71]; or prescribed heat flux, see e.g., Deng and Murakawa [72,73] and Taraphdar et al. [74].Modeling of the actual arch physics of the heat source weld arc is a complex research field and it is deemed out of scope for the macro-scale purpose of this study, though it is possible, see e.g., Hu and Tsai [75].However, to accurately predict of the temperature field, the shape of the heat distribution from the arc is important, several shapes discussed in the literature, see e.g., Kik [76].Common choices are the Gaussian distribution, see e.g., Eeagar and Tsai [77], a double ellipsoid distribution, see Goldak et al. [30,78], and a conical distribution, see Spina et al. [79].Farias et al. [80] compared the results from the three methods and found that all give valid predictions of the temperature field, with the two latter giving very similar results.In this study, we use a double ellipsoid heat distribution with the specific heat flux defined as: where qðx, y, zÞ is the specific heat flux and the parameters defining its size.The adopted numerical values for the size and heat flux parameters of the heat source are set such that full penetration of the simulated weld melt pool, for each weld pass is obtained, i.e., such that the liquid-like phase extends into base material or previous repair weld layers; depth a ¼ 6 mm, half of the width b ¼ 3 mm, front extent c i ¼ c f ¼ 3 mm, trailing extent c i ¼ c r ¼ 6 mm and total length c ¼ c f þ c r : The heat input fractions f i are set to f i ¼ f f ¼ 0:6 for the front section and f i ¼ f r ¼ 1:4 for the trailing section, in accordance with Goldak et al. [30,78].Furthermore, Q 0 is the nominal energy input which is set to Q 0 ¼ 4 kW, implemented together with a heat input cutoff temperature of 2500 � C.This is not to be confused with the mechanical material response cutoff obtained by the explicit modeling of the liquid-like phase.Moreover, the coordinates x, y, z in Eq. ( 21) are the distance from the center point of the heat source, which moves with the simulated welding speed of 5 mm/s.The moving heat source is implemented in Abaqus [66] by a user subroutine.

Addition of filler material during the weld process
In wire arc welding, filler material is added continuously during the welding process, and the supplied heat melts both the filler material and the base material, fuzing the two constituents.Therefore, to accurately simulate the welding process using FEM, the thermomechanical simulation model must be able to account for elements being continuously added during the simulation.
To address this modeling challenge, two common approaches are found in the literature; the quiet (or silent) element technique, see e.g., Lindgren et al. [69] and Kollar et al. [81], and the element birth (or inactive) technique, see e.g., Chen et al. [82] and Lee et al. [21].Studies comparing these two approaches show that both techniques generate similar temperature and stress fields for thermal and mechanical simulations, see Lindgren [83].However, the quiet element approach is more convenient to implement in the commercial software Abaqus used in this study and is, therefore, our method of choice.
To use the quiet element approach, we activate the filler elements individually as the moving heat source passes over the elements.This approach is more computationally expensive than lumping the filler elements into larger segments, so called macro beads, see Kik [76], and activating segment by segment, see e.g., Kiyoshima et al. [67] and Lundb€ ack and Lindgren [84].Further adding to unwanted computational cost, we simulate each weld pass of the multi-pass welding process explicitly, in contrast to lumping several layers of weld passes together, see e.g., Liu et al. [85].However, despite the added computational cost, this modeling allows for full utilization of the explicit modeling of the phase transformation kinetics presented in Section 2.1.2and thereby for a more detailed analysis of the residual stress evolution.
For the thermal simulation, elements in the quiet state are assigned a low thermal conductivity and an initial temperature of 1500 � : When activated, normal thermal properties are assigned, allowing the elements to participate in the heat transfer.To isolate quiet state filler elements form the heat transfer, the filler elements do not share nodes with the base material elements or the adjacent weld pass elements.Instead, the interfaces are assigned a user defined contact conduction condition, cf.gap conductance elements used by Michaleris [86].This contact condition is activated in the same instance as the quiet elements.In this way, the room temperature base material nodes are unaffected by the neighboring 1500 � C quiet nodes and fictitious temperature gradients avoided.
For the mechanical simulation, the filler elements are activated in the liquid-like state as the activation takes place by the passing of the weld torch.Similar approaches are found in the literature, e.g., a study by Brust and Rybicki [87] where they used a temperature limit of 1150 � C for filler elements' activation.To mechanically isolate quiet filler elements, elements in the quiet state are assigned material properties of the liquid-like phase.Using the FE mesh of the heat transfer simulation, the filler elements do not share nodes with the base material, therefore tie contacts are used to constantly connect the surfaces.Nodal positions of the filler elements are not adjusted upon element activation, in contrast to the work by Lindgren and Hedblom [88] and Lundb€ ack and Lindgren [84].

FE-model setup
The basis for both the thermal and mechanical simulations is a full scale rail FE-model consisting of a 1000 mm long rail section of European rail profile 60E1 [89].In the rail head, a damaged section of 80 mm length and 10 mm depth is cut out and is to be filled by repair welding, see Figure 4a, in accordance with the experimental work by Jun et.al. [26].The repair welding is done in four layers with 23 weld passes for the first two layers, 21 passes for the third layer and 19 passes for the final layer, see Figure 4b.Excess weld material of 1.5 mm thickness is included in the final layer.This excess weld material is removed after the material is cooled by using the quite element technique, i.e., it does not induce additional stresses.Moreover, the mesh of the weld layers and HAZ has characteristic element length of approximately 1.0 mm, which gives each weld pass cross-section a width and height of three elements.To constrain the model, the end surfaces of the rail are each connected to a rigid body coupling, for which the master node is fixed in all three translational degrees of freedom.This allows for a bending deformation mode as the top surface of the rail is heated and cooled.Note that the heat transfer and mechanical simulations use the same mesh and (linear) element type.
To validate the presented simulation methodology using experimental results from [26], the mechanical and thermal properties of the pearlitic rail material are obtained from [28,29].For the weld filler material, the TTT-diagram data is obtained from JmatPro [90], see Figure 5. Here, it is seen how the filler material is prone to produce a bainitic microstructure, even for slower cooling rates.However, a bainitic microstructure is not explicitly included in the modeling in the previous works [28,29] upon which the presented study is based.Therefore, we have included bainite and adapted its mechanical properties to data from JmatPro.Furthermore, the density of bainite is estimated to be that of martensite, cf.[91].Based on this estimation, the TRIP-related properties of bainite are approximated as that of martensite.The tempering of martensite and bainite is modeled as the diffusive transformations described in Section 2.1.2,cf.[92] and [93].The mechanical properties of tempered bainite is set to that of tempered martensite, obtained from [94].Whereby these tempered phases will be illustrated as one single tempered material phase.

Heat transfer simulation
Before the repair welding procedure starts, the rail head is preheated using a gas torch on the cutout and on a region 100 mm to either side of the cutout until the temperature of the entire surface is above 350 � C. To simulate this, a uniform heat is applied to the surface.It is linearly ramped up to 100 kW=m 2 for the first 100 s and is then held constant for 800 s.After this  preheating, the minimum temperature of the heated region is 350 � C and the repair welding process starts.The heat transfer simulation also accounts for heat radiation and convective heat transfer from the outer surfaces of the rail, see Eq. ( 22): where T 1 is the ambient temperature of 20 � C and T s the surface temperature.The heat radiation coefficient is set to b ¼ 0:7 [95] and natural convection heat transfer coefficient to a ¼ 25 W=ðmKÞ [96], using the Stefan-Boltzman constant of c SB ¼ 5:67 � 10 −8 W=ðm 2 KÞ: For the filler material to fuze with the base material at the start and end points of each weld pass, the start and end points for the moving heat source extends 2.5 mm into the base material, resulting in a total welding time of 17 s for each weld pass.Each new weld pass starts 1 s after the previous ends, to account for the repositioning of the weld torch.After a full layer is completed, the model is cooled for 60 s to allow the temperature to cool down to the initial preheat temperature level, following [26].Once all layers are completed, the model is cooled, without any forced cooling, for 2500s.The post welding corrective grinding is not simulated explicitly.As previously mentioned, the excess material is simply removed using the quiet element method.Figure 6 presents the temperature history of points below or at the surface of the center of the repair cutout section, highlighted in Figure 4a.The figure clearly shows the cyclic nature of the heat loads during the repair welding process.Worth noticing is that for all measured depths, the cycling temperature in one point passes the eutectoid temperature several times, demonstrating the cyclic nature also in terms of phase transformation kinetics.The maximum depth of the austenitization as indicated by the heat transfer simulation is approximately 5 mm.

Model reduction
The computational time for a mechanical simulation of the complete weld repair process using the full-scale 3D FE-model would by far exceed the limits for being of any practical use.Therefore, we here investigate different model reductions and compare results against the full 3D model for the pre-heating and the initial weld pass of the first weld layer.The temperature field imposed to the reduced models is from the full-scale 3D heat transfer analysis.Four reduced models are considered; three vertical cutouts with different thicknesses, taken from the full scale 3D FE-model at the center of the weld repair; and a 2D generalized plane strain model.The thicknesses of the cutouts that we consider are; a 30 mm, a 10 mm and a 5 mm section, illustrated in Figure 7. Three different boundary conditions (BC) for the cutout interfaces of the reduced models are also considered, see Figure 8; (a) zero longitudinal displacement u x and zero vertical traction t z on all the nodes, hereon referred to as fixed BC.(b) constant longitudinal displacement u x with zero total longitudinal force F x and zero vertical traction t z on all the nodes, referred to as sliding BC.(c) boundary displacements obtained from a full-scale simulation model using a coarser mesh and linear-elastic material model without phase-transformations, illustrated by the orange model in Figure 8 and referred to as submodel BC.
The fourth reduced model that is investigated is a 2D Generalized Plane Strain (GPS) model of the rail cross-section that uses Abaqus' CPEG elements [66].Using this generalized plane strain theory, it is assumed that the cross-section lies between two bounding rigid planes which can move in axial translation and rotation around the pitch and yaw axes of the cross-section, with respect to a defined pivotal point.The translation and rotation of the bounding planes allow for the 2D model to simulate both membrane and bending stresses in the cross-section.Moreover, the translation and rotation degrees freedom for the bounding planes are constrained by assigning axial and torsional stiffnesses which have been tuned to match the stiffnesses of the full scale 3D FE-model and validated using first order beam theory.
To evaluate which of the proposed model reductions that most accurately replicates the behavior of the full scale model, vertical, longitudinal and transverse stress histories obtained from the reduced models are compared to those obtained using the full scale model.Figure 9 presents the stresses histories from the different models, measured at the center of the repair cutout section (see highlighted point in Figure 4a) during the pre-heating and initial weld pass of the first weld layer.Note that the colors of the lines represent the different FE-models and the boundary conditions are indicated by the line type.As expected, the results show how the choice of boundary condition has greater impact than the thickness of the reduced model.Perhaps less expected, the GPS model is the reduced model that best replicates the full-scale model.Moreover, the results show how the sliding boundary condition setup fail to give consistent results.
An additional advantages of the 2D GPS model is its outstanding performance in terms of computational time, presented in Figure 10.Here, comparing the computational time for the different models again shows the inconsistency of the sliding boundary conditions and how, apart  from the sliding boundary conditions, the computational time decreases in proportion to the size of the FE-model.In conclusion, with both accuracy and computational cost considered, the GPS model, with tuned axial and bending stiffness, will be used when simulating the full repair welding procedure in Section 4.5.2.Moreover, the improved computational efficiency allows for the use of a finer mesh.Therefore, the full repair welding simulation uses an FE-mesh with an element length of 0.5 mm, which gives each weld pass cross-section a width and height of six elements, i.e., four times the total number of elements shown in Figure 4b.

Repair weld simulation
Figure 11 presents results of the simulated rail head repair welding process.To illustrate the process, three different time instances are presented column-wise; (i) the first weld layer, (ii) the third weld layer, and (iii) the final state after cooling and grinding.For each time instance, the following four results are presented row-wise; (i) the temperature field, (ii) a rudimentary illustration of the phase fractions, (iii) the longitudinal stress obtained using the linear mixture and (iv) using the self-consistent homogenization methods.The temperature field plots show how heat in each weld pass is applied to a small region and that there are strong temperature gradients to the surrounding material.This localized heat concentration is also shown in the phase fraction illustration, where a small melt pool is surrounded by austenitic material.Moreover, the illustrations show how the decelerated phase transformations of the filler material, see Figure 5, produce a bainitic-martensitic microstructure and how this is tempered by subsequent weld passes.In addition, note that the final microstructure becomes asymmetric.This is due to the heat from the pre-heating stage diffusing during the repair process, generating faster cooling rates during its later stages.The effect of the asymmetric final microstructure is also observed in the longitudinal stress plots, where the stress state differs not only from top to bottom but also from left to right in the repaired section of the rail.Furthermore, the plots show how the residual stress states generated by the two homogenization methods match in single-phase regions and differ significantly in multi-phase regions.This inconsistency is revealed by noticing how the stress fields from the methods match at the first time instance and diverge during the prolonged multi-phase states of the latter time instances.Note also how the results correlate to the three-bar results presented in Figures 2 and 3, i.e., how the linear mixture rule generates more compressive longitudinal residual stresses than the selfconsistent method.Further analysis of the in-process stress field shows how the decelerated phase transformations affect the HAZ during the welding process; a zero stress state in the melt pool, low stresses in the "mushy" zone, compressive stresses in the bainitic-martensitic microstructure preceding the melt pool and tensile stresses in the austenitic-bainitic microstructure succeeding the melt pool.Finally, concurring with the results of the three-bar experiment, by accounting for the intricate interaction of different phases in multi-phase stages, the stresses obtained using the self-consistent method are deemed more realistic than those obtained using the more simplistic linear mixture rule.
Figure 12 presents longitudinal residual stress plotted along a vertical axis through the center of the rail cross-section.The figure shows how simulation results obtained using the linear mixture rule and self-consistent homogenization methods give similar overall stress distributions, except in the repaired region where the linear mixture gives compressive stresses and the self-consistent method gives tensile stresses, as also seen in the third column of Figure 11. Figure 12 also presents experimental measurements found in the literature [26], which show that the longitudinal residual stress is tensile in the rail foot, compressive in the web, tensile in the rail head and compressive at the rail surface.Comparing the measurements to the simulation models shows how the general behavior of the residual stress field is replicated in the foot, web and head of the rail but not at the top surface and how the self-consistent method correlates more closely, overall, than the linear mixture rule.However, the discrepancy at the surface of the rail is of great significance since tensile residual stress, as suggested by our simulation models, can promote fatigue crack initiation and propagation.Large tensile longitudinal residual stresses can be expected on the rail head surface as only the upper part of the rail cross-section is heated and cooled to room temperature, which is also shown by Kabo et al. [27].They also show that the tensile stresses are redistributed during operational conditions.Further contemplating the surface stress discrepancy, the effect of post-repair corrective surface grinding in the experimental procedure [26] should be considered.Surface grinding is known to produce compressive stresses, see e.g., [97], which likely explains the compressive surface stresses in the experimental results shown in Figure 12.However, this procedure is not modeled explicitly in the presented simulations, as removing the excess weld material elements using the quiet element method does not induce any stresses.

Figure 12.
Longitudinal stress in the rail cross-section, comparing experimental results from [26] to simulations using the linear mixture rule (blue) and the self-consistent homogenization method (orange).

Conclusions and outlook
The study presented in this paper simulates rail repair welding by extending the material modeling approach formulated in previous works [28,29].The previous material modeling includes cyclic plasticity, phase transformation kinetics and homogenization of individual constitutive models, and is in this paper extended by explicit modeling of a liquid-like phase and a virgin material state when liquid-like material solidifies.The importance of these modeling extensions is demonstrated by simulating a classic three-bar welding experiment.This simple experiment is also used to study the effect of using different homogenization methods to account for the multi-phase stages of the welding process, where the self-consistent homogenization method is found to give the most realistic behavior.However, since the linear mixture method is the most frequently used in the literature, we compare these two methods in the succeeding rail weld repair simulations.
For the railhead repair welding simulation, we use a one-way coupled thermomechanical finite element analysis.By comparing our simulations to previous works in the literature, e.g., Jun et al. [26], Kabo et al. [27], and Lee et al. [21], we present a detailed simulation methodology using more advanced material modeling, explicit modeling of a moving Goldak heat source [30] and modeling of continuous addition of weld filler material.However, for the simulations to be of practical use, a model reduction study is performed to develop a model that is computationally efficient and also accurate.We use the full-scale 3D heat transfer simulation and different model reductions for the mechanical simulations.To ensure that simulation accuracy is not lost, the simulation results from the reduced models are compared to the 3D full-scale mechanical simulation.
The model reduction study shows that results obtained from the 2D general plane strain (GPS) model, extended with axial and bending stiffness tuned against the 3D model, correlate surprisingly well to results from the full-scale 3D model at a fraction of the computational cost, allowing the use of a finer mesh.In further successful validation of the 2D GPS model, the longitudinal residual stress results correlate well with experimental measurements by Jun et al. [26,33].In this validation, the self-consistent method correlates more closely than the linear mixture rule.Furthermore, the characteristics of the two homogenization methods are highlighted during the decelerated phase transformations of the weld filler material, causing multi-phase material states during the entire repair welding process.In these states, the difference between the homogenization methods increases significantly as the softer austenite and harder martensite (and bainite) are combined differently by the two methods.Using the self-consistent method, the individual constitutive models of the phases and the analytically based modeling of their mechanical interaction allow for a more realistic material response.Moreover, this in-process evaluation also shows the pronounced cyclic nature of the stress history, where material may be melted up to five times during the repair process.This remelting during the multiple weld passes further demonstrates the importance of accounting for both a virgin material state and the melted (liquid-like) phase.
Other than typical simplifications used in computational welding mechanics (CWM), such as not accounting for plastic energy dissipation or stress-induced phase transformations, the welding simulations in this study use a rudimentary discretization of the shape (square) of the weld beads in each weld pass.This is deemed to be sufficiently accurate for the overall simulation outcome, as the heat flux input and the amount of material deposited in each weld pass govern the material response during the repair process to a larger extent than the bead shape.Also, as in all weld simulations, uncertainties in the input heat flux limit the deterministic nature of the study to some extent.However, these uncertainties do not undermine the simulation methodology's capability to perform intricate parameter studies or its ability to provide insight into the complex evolution of the residual stresses during the repair welding process.
Overall, the simulation methodology presented in this study enables detailed repair welding process simulations to be performed iteratively in future process parameter studies.Even though the main focus of this study is demonstrating the methodology, the simulation results highlight the importance of regulating and monitoring the operational temperatures during the repair process.Process parameters, such as heat input, inter-pass and inter-layer cooling times, should be adjusted so that the temperature and temperature rate do not allow for excessive spheroidization of the pearlitic microstructure or cooling rates causing martensite formation.The risks involved are dire as both have detrimental effects on the mechanical performance of the rail.Therefore, the presented simulation methodology could aid in avoiding such consequences and facilitate a truly effective repair, defined by Workman and Kral [98] as a repair that performs similarly to the parent rail and the loss of service time kept to a minimum.To this end, further studies using the presented numerical tools should be carried out as support when determining optimal process parameters of the repair welding process, and if needed, update the welding procedure specifications (WPS) for rail surface repair.

Figure 1 .
Figure 1.(a) A three-bar FE model with a heated bar in red, cold bars in blue, and vertical constraint in green, cf.[44, 64, 65], (b) schematic illustration of two plates butt-welded together.

Figure 2 .
Figure 2. Simulation results for the heated bar in the three-bar model.Bottom to top graph descriptions; prescribed temperature load, volume fraction, axial stress, and axial strain.Temperature load magnitude is 1600 � C.

Figure 3 .
Figure 3. Simulation results for the heated bar in the three-bar model.Axial stress in heated bar when temperature load magnitude of 1000 � C.

Figure 4 .
Figure 4. (a) FE model of a 1000 mm long rail section with the European rail profile 60E1 [89] which is used in the welding simulations.The temperature and stress measure point is shown with a white cross.(b) Cross-section view of the FE model used in weld simulation, mesh size, weld passes and excess material are shown by different colors in highlighted circle.

Figure 5 .
Figure 5. TTT diagrams for the base and filler material used in the simulations, obtained from JmatPro [90].

Figure 6 .
Figure 6.Temperature history at the center point of the rail head cutout highlighted in Figure 4a, measured at different depths.

Figure 9 .
Figure 9.Comparison of the longitudinal stress component obtained using different FE-models, measured (a) at the surface and (b) 8 mm beneath the surface of the point highlighted in Figure 4a.

Figure 10 .
Figure 10.Comparison of computational time using different FE-models.The computational time is normalized using the total time of the full scale 3D model.

Figure 11 .
Figure 11.First row: Temperature field.Second row: Rudimentary phase fraction field illustration; red -martensite, white -austenite, blue -pearlite, green -tempered martensite and/or bainite, yellow -melt and black -bainite.Third row: Longitudinal stress field using the linear mixture rule.Forth row: Longitudinal stress field using the self-consistent homogenization method.