An Eulerian thermomechanical elastic–viscoplastic model with isotropic and directional hardening applied to computational welding mechanics

An Eulerian thermomechanical elastic–viscoplastic model with isotropic and directional hardening is used to analyse the residual mechanical state resulting from the arc welding of a multi-pass weld. Details of the weld test plate, weld filler material, and numerical implementation of the model are provided, including integration algorithms and consistent tangent modulus. For the computational welding mechanics analyses, the austenitic ASME stainless steel grade 316L was considered so that no phase transformations of solid states needed to be considered. The maximum residual stresses were found to be about 500–600 MPa, which is of the order of the yield stress of the base material. Variations in the heat input and the resulting weld cooling time had a significant influence both on the residual stress state and on the resulting geometry of the weld. The predicted stress levels were compared to the experimental results. Overall, the proposed Eulerian framework seems to be a promising tool for analysing melting/solidification processes and residual mechanical states.


Introduction
Within the context of an Eulerian formulation of constitutive equations, consistent thermomechanical equations are proposed for modelling an elastic-viscoplastic material including isotropic and directional hardening. In particular, the constitutive equations model the phase transformation between solid and fluid response during melting and solidification. The implemented framework is applied to the challenging transient thermomechanical problem of simulating a 2D multi-pass weld in a narrow grove of a plate with emphasis on prediction of the residual mechanical state in the way of the weld. In welding mechanics simulations, the term Eulerian sometimes refers to a fixed mesh in the numerical simulation with material flowing through the mesh much like a fluid (cf. [38,50]). In contrast, here the term Eulerian refers to the constitutive formulation which depends only on state variables that characterize the current state of the material, without the need for a reference configuration.
The standard Lagrangian formulation of plasticity for finite strains utilizes a multiplicative decomposition of the deformation gradient (e.g. [6,24,26,33]), and use is made of reference and intermediate configurations.
That framework introduces an arbitrariness of reference and intermediate configurations as well as of total and inelastic deformation measures which does not exist in the Eulerian formulation of the problem. The Eulerian formulation adopted here is based on the works of Eckart [12] and Leonov [27], who proposed an evolution equation for the elastic deformation directly. Earlier versions of the present Eulerian framework have also been explored, e.g. Rubin [44,46], Rubin and Attia [47], and Kroon and Rubin [25].
In addition to rate dependency, the Eulerian model proposed here includes isotropic and directional hardening as well as the effect of thermal recovery. The model for isotropic hardening is adopted from Chan [9], who proposed a model of hardening for processes where there is sufficient time for thermal recovery (annealing effects) to influence the resistance to inelastic deformation rates.
Directional hardening is an alternative to kinematic hardening, which is usually included to account for such phenomena as the Bauschinger effect [4], and the literature contains a vast number of models for kinematic hardening (e.g. [2,8,57]). The model for directional hardening adopted here is one of the novelties of the present work, where directional hardening is modelled using a single scalar-valued state variable. Hence, for instance, the Bauschinger effect can be modelled using this state variable, which in a simplified way plays the role of the back-stress tensor in many kinematic hardening models.
Another contribution of the present paper is the application of a novel 2D weld heat source model. The 2D heat source is derived from Goldak's 3D heat source [15] and is applied to the simulation of a multi-pass weld joint, which is a complex and challenging transient thermomechanical problem. Since about 1920, arc welding is a very common manufacturing method and is nowadays often the most cost-effective 'off-theshelf' solution available for joining metallic materials. Wire arc additive manufacturing (WAAM) is another manufacturing process with similar physical phenomena involved. Arc welding and additive manufacturing constitute complex, transient, and thermomechanically coupled processes, which involve such phenomena as melting, solidification, and phase transformations between different solid states. During these processes, the temperature can change by thousands of degrees during a few seconds, causing dramatic changes in both thermal and mechanical properties of the materials involved.
In general, arc welding-based manufacturing processes leave the structure with significant geometrical distortions or residual stresses in the base and weld metal materials, which can be of great importance, for instance, with regard to the integrity of the structure. Accurate prediction of the residual stress state created in such processes requires material models that are able to predict the hardening behaviour of the material in a wide temperature range. In addition, material data for the same temperature range must be available, which is a great challenge by itself.
Computational welding mechanics (CWM) has a history that goes back to the 1970s [1,13,19,54]. In CWM, a few simplifications are usually made. For instance, the weld melt pool is, strictly speaking, in a flowing liquid form but is approximated to be a viscoplastic solid with low yield strength. Overall, the arc physics and the physics of the weld melt pool are simplified, such that only thermomechanical phenomena of the pool are considered, where the weld heat input (arc energy) of the arc welding process is modelled using Rosenthal's and Rykalin's solution [40,49] of the general heat conduction equation where the energy loss of each single arc welding process is taken into account by the use of an empirical arc efficiency factor [53]. Furthermore, rate of deformation effects (i.e. viscoplasticity) are often ignored, and a rate-independent plasticity formulation is usually adopted.
The arc welding process is, strictly speaking, a transient coupled thermomechanical problem. However, from a numerical point of view, it is sometimes convenient to ignore this coupling and instead solve the problem by a sequentially coupled thermomechanical analysis (e.g. [11,14,55,58]), where the thermal and mechanical problems are solved in sequence rather than in parallel. Lindgren [29] provides further information and an overview of CWM.
The material model proposed here is applied in the simulation of a welding process. More specifically, a gas tungsten arc welding (GTAW) process, applied to the manufacturing of a multi-pass weld in a narrow grove in an austenitic stainless steel plate ASME grade 316L, is considered. The austenitic steel is considered to avoid the need for modelling phase transformations between solid states. In addition, this is a material for which the experimental results are available for calibration (e.g. [10,30,31,57]). The outcome of the numerical analyses is the mechanical state resulting from the welding process.
The paper is organized as follows: In Sect. 2, the Eulerian theoretical framework and the constitutive model to be used are introduced. The basic features of the proposed model are illustrated in Sect. 3, and a calibration of the model and a comparison with the experimental results are performed in Sect. 4. The numerical example of the welding simulation is provided in Sect. 5. Finally, Sect. 6 contains a discussion and some concluding remarks.

Kinematics
The position of a material point in the present configuration at time t is denoted by x, and the velocity v is given by where( •) denotes the material time derivative. Also, the velocity gradient L and the rate of deformation tensor D are defined by

Evolution equations for elastic deformations
Using the Eulerian formulation of viscoplasticity developed in [47], the inelastic deformation rate is isochoric, so the elastic dilatation J e is determined by the evolution equatioṅ where Q : R = tr(QR T ), I is the second order unit tensor, and div(•) is the divergence operator relative to x. Also, the elastic distortional deformation is a symmetric, positive definite, unimodular tensorB e determined by the evolution equationḂ where Γ is a non-negative function that determines the rate of inelasticity and the direction of inelasticity is determined by the tensorĀ p , defined byĀ This tensor satisfies the conditionĀ which ensures thatB e remains unimodular (i.e. det(B e ) = 1). In contrast to the standard multiplicative Lagrangian formulation of plasticity (e.g. [6,24,26]), the Eulerian formulation of elastic-viscoplastic response is insensitive to arbitrariness of a reference configuration, an intermediate stress-free configuration, a total strain measure, and an inelastic deformation measure. In particular, it was shown in [45] that when this arbitrariness is removed from the multiplicative formulation, it must reduce to the formulation for elastically anisotropic response proposed in [44].

Governing equations
Within the context of the thermomechanical theory of Green [17], the conservation of mass and balances of linear momentum and entropy are proposed in the formṡ 2) ρη = ρξ + ρs − divp, (7.3) where ρ is the current mass density, b is the specific (per unit mass) body force, T is the symmetric Cauchy stress, η is the specific entropy, ξ is the specific internal rate of entropy production, s is the specific external rate of entropy supply, and p is the flux of entropy per unit current area. The balance of angular momentum requires T to be a symmetric tensor. Also, the local form of the balance of internal energy is where ε is the specific internal energy, r = θ s is the external rate of supply of heat energy, q = θ p is the heat flux vector, and θ is the absolute temperature. Next, the specific Helmholtz free energy, ψ, is introduced as and the rate of material dissipation, ρθξ , defined in [43], is given by where g = ∂θ/∂x is the temperature gradient. Then, the balances of entropy and energy in (7.3) and (8), respectively, can be used to write the rate of material dissipation in the form which is required to be non-negative for all processes.

Constitutive equations for an elastically isotropic material
For an elastically isotropic material, the Helmholtz free energy is taken in the form where α 1 is a scalar pure measure of elastic distortional deformation. For both rate-independent and rate-dependent material response, it is assumed that T and η are specified by where where ρ z is the zero-stress mass density at reference temperature θ z . Furthermore, it is noted that the tensorsB e andB e introduced here correspond to B e and B e , respectively, in [47].

Specific thermoelastic constitutive equations
The Helmholtz free energy in [41] is modified to take the form where f 1 (θ ) controls the purely thermal dependence of ψ, f 2 (θ ) controls thermal expansion, and K (θ ) and μ(θ) are the temperature-dependent zero-stress bulk and shear moduli, respectively. It then follows that Also, the function f 2 (θ ) is restricted by the condition which ensures that a zero-stress state at reference temperature is characterized by In addition, the von Mises stress σ e is defined by 1,2) where γ e is the effective elastic distortional deformation. For the numerical algorithm, the Kirchhoff stress, τ , is also introduced as 2.6 Specific constitutive equations for viscoplasticity An overstress model (e.g. [32,37]) is used to characterize the rate-dependent inelastic response with the inelastic distortional deformation rate controlled by the function Γ which is specified by the overstress form where a 1 is a temperature-dependent function that controls rate dependence (with higher rate dependence occurring for smaller values of a 1 ), and the Macaulay brackets g are defined by The yield function g is specified by where R is a function of the temperature that influences the temperature dependence of the yield strength and Z is a hardening variable that accounts for both isotropic and directional hardening defined by where κ and β are state variables accounting for isotropic and directional hardening, respectively. The variable U accounts for the mode of loading and is a variant of the Lode parameter: The stress states uniaxial tension, pure shear, and uniaxial compression correspond to U = 1, U = 0, and U = −1, respectively. Following the work in [20,21], let κ be a scalar measure of (isotropic) hardening determined by the evolution equationκ The first term in (26) accounts for work hardening. This term is proportional to the rate of plastic deformation, Γ , m 1 controls the rate of hardening, and κ s determines the saturated value of hardening. The second term accounts for thermal recovery at high temperatures, i.e. material softening, which is taken to be independent of the inelastic deformation rate. The rate of thermal recovery is governed by m 2 , and the target hardening (or rather softening) value for the thermal recovery process is κ a . Both the hardening and thermal recovery rates are taken to be functions of temperature.
Here, a new scalar form for directional hardening is proposed to determine the state variable β by the evolution equationβ The first term accounts for work hardening, and the second term accounts for thermal recovery. As before, the rate constants, in this case m 3 and m 4 , are taken to be temperature-dependent, and β s governs the magnitude or saturated value of β. Hence, the hardening variables are bounded according to In order to characterize directional hardening, a variant of the Lode angle, ϑ, is introduced as Since U varies between −1 and 1, it is sufficient to consider the range ϑ ∈ [−π/6, π/6]. Following Rubin [42], three special load cases are identified: -triaxial extension (TXE), U = 1, ϑ = π/6, -pure torsion (TOR), U = 0, ϑ = 0, -triaxial compression (TXC), U = −1, ϑ = −π/6.
The loading state is often characterized by the stress tensor or the stress deviator. In the present case, since use is made of a strain-based formulation of inelasticity, it is more convenient to consider the deviatoric deformation tensorB e , whose eigenvalues are denoted byB e1 ,B e2 , andB e3 , whereB e1 ≥B e2 ≥B e3 is understood. The three load states introduced above are then characterized by -TXC:B e1 =B e2 = 2γ e 3 ,B e3 = − 4γ e 3 .

Specific constitutive equations for heat conduction and convection
For Fourier heat conduction, the constitutive equation for the entropy flux is specified by where k(θ ) > 0 is the coefficient of heat conduction. The boundary condition associated with heat convection from a (hot) solid structure to the surrounding air is specified by the convection law where n is the unit outward normal vector to the boundary surface, q s is the normal component of the heat flux vector on the surface of the structure, α h is the coefficient of heat convection, θ s is the temperature at the surface of the structure, and θ ∞ is the temperature of the surrounding environment.

Numerical implementation of the constitutive model
The constitutive model above was implemented as a user subroutine in Abaqus, and details of the numerical implementation are provided in Appendix A.

Prerequisites
Several of the basic features of the present model have been described in previous works (e.g. [20,21]). For instance, the parameter R is intended to model the decrease in the yield strain with increasing temperature. Furthermore, a 1 accounts for the rate dependence of the material with lower values of a 1 causing a higher rate sensitivity of the material. The scalar-based description of directional hardening, however, is a new feature that will be explored in the present Section. In the following parametric study, the elastic properties of the material are given by Young's modulus E = 200 GPa and Poisson's ratio ν = 0.3, and the material is taken to be at room temperature, so that f 2 = 0. For the inelastic deformation rate and hardening, a 1 = 10 5 /s, R = 1 κ 0 = κ a = 0.0015, and κ s = 0.012 are assumed. The results are shown in terms of the Cauchy stress components T i j vs. the engineering strain components e i j = sym(∂u i /∂ x j ) relative to fixed rectangular Cartesian base vectors e i .

Evolution of the yield surface for different load cases
The yield function may be illustrated by its contour in the synoptic (or octahedral) plane, whose unit normal is given by the vector (p 1 + p 2 + p 3 )/ √ 3, where p i are orthonormal eigenvectors ofB e . Again, following Rubin [42], two orthogonal unit vectors in the synoptic plane are defined as The yield surface can be described by the vector b s , defined as where  The state variable κ governs the isotropic expansion (or contraction) of the yield surface in the synoptic plane, whereas evolution of β causes directional hardening, i.e. different hardening for different loading modes. Figure 1 shows contours of the normalized yield surface in strain space and in the τ 3 -τ 1 -plane for different values of β. From this Figure, it can be seen that loading in uniaxial tension (TXE) causes an increase in β and thereby increased hardening for loading modes with ϑ > 0 and softening for ϑ < 0. Loading in uniaxial compression (TXC) causes the opposite tendency. This normalized yield surface is pivoting around the yield point at ϑ = 0 (TOR), which remains unaffected by directional hardening.

Influence of directional hardening during cyclic uniaxial tension
In this Subsection, consider loading in uniaxial tension in the e 1 -direction of the material. A loading rate oḟ e 11 = 0.01/s is applied. Figure 2a-c shows graphs in terms of T 11 vs. e 11 . In Fig. 2a, the influence of the hardening rate m 1 is illustrated. In this case, directional hardening and thermal recovery have been deactivated. The material response is relatively intuitive with the hardening rate increasing with increasing value of m 1 . Figure 2b shows the influence of the rate of directional hardening, m 3 . The higher the value of m 3 , the faster the hardening curve reaches the saturated level of directional hardening, i.e. β → β s . For the highest value, m 3 = 0.5/s, directional hardening saturates already at a strain of about 0.01-0.02. It should also be noted how the model clearly has the potential to model a Bauschinger effect, since yielding during reversed loading occurs at a smaller effective stress in compression compared to the previous yielding point in tension.
In Fig. 2c, variations in β s are considered. It can be clearly seen how an increasing value of β s adds more potential for hardening to the material. The value of m 3 chosen in these simulations again implies that the hardening saturates after a strain of about 0.01 − 0.02.
Finally, Fig. 2d demonstrates the influence of thermal recovery, that is the influence of the rate constants m 2 and m 4 . In these simulations, m 2 = m 4 is applied. Since thermal recovery is a process that depends explicitly on time, T 11 is plotted vs. time, t. In this case, the material is loaded with a strain rate ofė 11 = 0.01/s up to a strain of e 11 = 0.05. After this, the strain is held constant, implying that thermal recovery is the only potentially active process. The simulation continues up to the time t = 100s, during which time thermal recovery causes the hardening to decrease and thereby the stress to decrease. As can be seen from Fig. 2d, thermal recovery causes a relaxation of the stress in the material. For the highest recovery rates, the hardening variables κ and β have sufficient time to relax to their annealed values κ a and 0, respectively, which is associated with the stress stabilizing at about 220 MPa.

Preliminaries
The purpose of the present Section is to calibrate the thermomechanical model, so that it can be used for prediction of the residual mechanical state in the simulation of a multi-pass weld in the next Section. Temperaturedependent thermoelastic and viscoplastic properties are determined on the basis of the experimental results from literature.
In the remainder of this paper, the symbol T is used for the temperature in degrees centigrade ( • C), while θ is used for the absolute temperature.

Model calibration: temperature-dependent thermal and elastic properties
The density of the material in a zero-stress state at room (reference) temperature, T z = 20 • C, is ρ z = 7, 966 kg/m 3 . The temperature dependence of Young's modulus, E, is illustrated in Fig. 3a and is adopted from Lindström [31]. The bulk modulus was taken to be constant and independent of temperature. Assuming that Poisson's ratio at room temperature is ν(T z ) = 0.3, the bulk modulus was calculated as The temperature-dependent shear modulus was then given by Thermal expansion was quantified by the function f 2 (T ), illustrated in Fig. 3b. The function f 2 (T ) vanishes at T = T z .

Model calibration: inelastic yielding during reversed loading
Information about the present material's response during reversed loading is only available for room temperature (as far as the current authors have found). Hence, the Bauschinger effect of the austenitic steel at hand is considered for room temperature. The inelasticity in the model is compared to data from Choteau et al. [10],   Table 1 Model parameters for reversed loading at room temperature data from Choteau et al. [10]. The model was calibrated manually, and the model parameters thus determined are listed in Table 1.

Model calibration: temperature-dependent inelastic flow rate properties
The temperature-dependent inelastic flow rate properties were determined from uniaxial stress compression tests performed at different strain rates and different temperatures [30]. The flow response was quantified in terms of true (Cauchy) stress component T 11 and the associated true (logarithmic) strain component ε 11 . The tests were performed at constant engineering strain rates, and the plasticity model was calibrated using data forė 11  For Model B, a true distinction between isotropic and directional hardening was only possible at room temperature, according to the model calibration in the previous Subsection. To generalize the results for higher temperatures, it is assumed that the ratio m 1 /m 3 , i.e. the ratio between the hardening rates of isotropic and directional hardening, is constant and independent of temperature. This means that the ratio m 3 = 0.2/0.0014m 1 ≈ 143m 1 , determined at room temperature, is taken to hold for higher temperatures as well.
It is also assumed that the thermal recovery rates for isotropic and directional hardening are the same functions of temperature, i.e. m 2 = m 4 . Finally, it is assumed that the maximum value of the state variable β that governs directional hardening, i.e. β s , is constant and independent of temperature. The last assumption is supported by the fact that the saturated value of isotropic hardening, κ s , was found to be fairly independent of temperature, see below.
In Figs. 5 and 6, the experimental data from Lindgren et al. [30] are shown together with the calibrated model responses at the different temperatures. For the temperatures 20-700 • C, it was assumed that there was not sufficient time thermal recovery to be active so m 2 and m 4 were set equal to zero.
In the calibrations, an accurate fit in the strain range ε 11 ∈ [0, 0.3] was prioritized, since the strains in the present applications were essentially confined to this regime. Lindgren et al. [30] also emphasized that in the tests performed at the highest strain rates the relative temperature increase was significant in some of the tests at the lowest temperatures. For this reason, in cases where the rate dependence in the experimental data is not  consistent with the behaviour predicted by the model, priority is given to the data for the lowest strain rates, since these data can be expected to be least affected by temperature variations. In principle, the initial value κ 0 of κ needs not be equal to the fully annealed value κ a . However, during the fitting procedure it was concluded that taking κ 0 = κ a enabled a good fit to the data, so the same value was used for these two variables. Also, it was concluded that κ a and κ s could be taken to be independent of temperature. For both models A and B, κ s = 0.012 was adopted, whereas κ a was set to 0.002 and 0.0015 for models A and B, respectively. As indicated above, β s = 0.3 was specified for all temperatures.
The remaining model parameters of the inelasticity rate were fitted at each test temperature, and the estimated parameters are listed in Tables 2 and 3.
The temperature changes continuously both in time and space, and for a specific temperature, the associated model parameters were attained by linear interpolation between the discrete values as listed in Tables 2 and  3. The melting temperature was taken to be 1, 500 • C, and between 1,300 • C and 1,500 • C no data for the flow properties were available. In the range 1,300 • C to 1, 500 • C, the properties corresponding to 1,300 • C were therefore applied, except for the parameters R, m 2 , and m 4 . The melting process was accounted for by letting R decrease to a low value, i.e. R(1, 500 • C) = 0.01. Also, the thermal recovery process was taken to be accelerated at temperatures approaching the melting temperature, and m 2 (1, 500 • C) = m 4 (1, 500 • C) = 100/s was therefore applied.
The following comments are based on the values recorded in Tables 2 and 3. The parameter a 1 governs the influence of the strain rate on the flow stress, with lower values increasing the strain-rate sensitivity. The strain  rate has a stronger influence on the flow stress at higher temperatures, which corresponds to lower values of a 1 at higher temperatures. The parameter m 1 governs the rate of work hardening. Although the variation is not completely monotonic, the hardening rate does tend to decrease at high temperatures. Thermal recovery, quantified by m 2 and m 4 , kicks in at 800 • C and seems to increase exponentially with temperature. The parameter R starts at unity and decreases monotonically with increasing temperature. This parameter controls the temperature dependence of the yield strain, and R → 0 at 1,500 • C models melting. For numerical reasons, R = 0.01 was adopted as the minimum value.

Problem formulation
A numerical example is provided in the form of a transient, thermomechanical analysis of a multi-pass weld. Figure 7a shows a full 3D image of the problem of arc welding. The geometry that was actually modelled and analysed is the 2D profile, shown in Fig. 7b, with a close-up view of the weld joint region in Fig. 7c. The dimensions are included in each of these Figures. The base material occupies the domain 0 . The weld joint is produced in two steps/runs. In the first run, the root weld, domain 1 , is produced, and in the second run, the second weld, domain 2 , is produced, as indicated in Fig. 7c.
The welding process was simulated by introducing a prescribed heat source in the domains 1 and 2 during runs one and two, respectively. The two runs were simulated, and the output from the analysis yielded 2D distributions of the field variables.

Prerequisites for the numerical analysis
The finite element code Abaqus [18] was used for the numerical analysis. The elastic-viscoplastic model was implemented as a user material (UMAT), using the implicit integration of the evolution equations outlined above. Instead of solving the fully coupled energy equation (8), use was made of the approximate formulation for heat conduction in Abaqus. Specifically, the heat conduction equation was solved with c p (θ ) being the specific heat capacity at constant pressure. The problem was analysed as a 2D problem, which seems to give sufficiently accurate results in terms of the resulting residual stress state [29]. The computations were performed in parallel on 8 cores. In contrast to many other CWM simulations, the coupled thermomechanical problem was solved in parallel, i.e. no staggered or sequential approach was used. For the out-of-plane dimension, generalized plane strain was assumed. (The 2D structure is then essentially constrained in the out-of-plane direction between two rigid blocks which are themselves free to slide in the out-of-plane direction.) A constitutive thickness, H = 1 mm, for the out-of-plane dimension was specified.
Both the base and weld materials were taken to be the austenitic stainless steel 316L, and the material model was calibrated using the experimental results from the literature [30,31], as shown in the previous Section.
The mechanical state at a material point was characterized by the state variable vector [B e,11Be,22Be,33Be,12Be,13Be,23 J e κ β], whereB e,ij are the rectangular Cartesian components ofB e , and the initial value of this vector at room temperature was specified by [1 1 1 0 0 0 1 κ 0 0].

Discretized geometry
The geometry was discretized, and the standard mesh used in the simulations is shown in Fig. 8. Symmetry was utilized, such that only half of the 2D model was analysed. The smallest elements in the refined zone close to the boundary between the weld and the base material had a size of l rz = 0.2 mm, and the maximum size of the elements in the areas without any significant gradients was l cz = 5 mm. Quadrilateral, hybrid elements were used in the finite element analysis, i.e. quadratic and linear shape functions were used to interpolate the displacement and temperature fields, respectively, and a linear variation in the pressure was enabled.
The mechanical boundary conditions applied are indicated in Fig. 8. Apart from the displacement constraints, indicated in Fig. 8, the boundaries were traction-free. Also, influences of inertia and gravity were ignored.
Along the red dashed line, the boundary heat flux was governed by the convection law in (31). Heat convection along the weld area was not modelled explicitly. This loss of heat from the weld pool was instead accounted for by modifying the power input from the welding heat source, see below.

Implementation of the welding process
The welding process was implemented as a series of computational steps, as illustrated in Fig. 9. Heat was supplied to the weld regions in two ways: first, the weld pool was given an elevated, prescribed temperature, T wp = 2000 • C, which corresponds to a liquid state of the material; second, heat was supplied by an internal heat source, r w , during a subsequent, transient analysis.
Before the beginning of the welding process, at t < 0, the whole geometry was in a zero-stress state at the temperature T rm = 20 • C. At t = 0, the process for the root weld was initialized by instantly elevating the temperature of the domains 1 and 2 to T wp . At the same time, the material in 1 and 2 was initialized by assigning the values J e = 1/(1 − f 2 (2000 • C)), κ = κ 0 , and β = 0, which approximates a thermally expanded liquid at T = 2, 000 • C. In the subsequent transient analysis (0 < t < t 1 ), the root weld was simulated and a heat source, r w , was applied. The time t 1 was chosen large enough, so that at t = t 1 the whole model (including 1 ) had virtually cooled down to T rm . During 0 < t < t 1 , the domain 2 should, strictly speaking, be a void. This was approximated in the following way: during the simulation of the root weld (0 < t < t 1 ), the domain 2 was ascribed the properties of a liquid by enforcing R = 0.01 in 2 (independent of the temperature). In this way, 2 was prohibited from enacting any significant mechanical influence on the root weld or the base material. As indicated above, heat  convection from the weld surface to the surroundings was not modelled explicitly but was instead accounted for by applying an effective (i.e. reduced) amplitude for the heat source.
The processing of the root weld was followed by a short time period t 12 , during which the conditions for the second weld were set. During this period, the temperature of 2 was raised to T wp and the material in 2 was re-initialized by (again) assigning the values J e = 1/(1 − f 2 (2000 • C)), κ = κ 0 , and β = 0.
During the last part of the analysis, the second weld was simulated. A heat source, r w , was again applied, but this time to 2 , as indicated in Fig. 9.
In the present analysis, the values t 1 = t 2 = 10 5 s, and t 12 = 0.001s were specified.

Heat source
The welding process was simulated by introducing a prescribed heat source, r = r w (X 1 , X 2 , t), in the domains 1 and 2 , see Fig. 10 where an orthogonal coordinate system with the coordinates X 1 , X 2 , and X 3 has been introduced. The heat source was confined to a region with the radius l w , as indicated in Fig. 10.  The heat source applied in the present analysis is a 2D version of the 3D heat source proposed by Goldak et al. [15]. Specifically, where the function r G2D (t) is defined by In this expression, P w denotes the (effective) electric power input, and the function r G3 (X 30 , t) is defined as where X 30 is a constant, v is the velocity of the heat source (in the X 3 -direction), t denotes time, and c f and c a define the profile of the heat source (in the X 3 -direction). In the full 3D analysis, the heat source would move in the (negative) X 3 -direction with the velocity v. Since we consider a 2D analysis in the X 1 -X 2 -plane, the shape of the heat source in the X 3 -direction determines the amplitude as a function of time for the 2D heat source. In Fig. 11, an example of a heat source amplitude is illustrated for the case c f = 1, c a = 3 mm, and X 30 = −c f √ ln 0.001/(−3). The value of X 30 is chosen such that for t = 0 the amplitude function has a value that is a fraction 0.001 of the maximum amplitude value.
Additional explanation of these expressions for the heat source is provided in Appendix B. As indicated above, at the beginning of the analysis (t = 0), the weld pool is assigned the initial temperature T wp , whereas the base material has the initial temperature T rm . The effective power amplitude of the heat source may therefore be expressed as  where P w0 is the nominal electric power supply (without any preheating), e in is the energy (per unit mass) required to raise the temperature of the weld pool material from T rm to T wp , and A * denotes the area in the model that is given an elevated temperature at the beginning of a welding step. The areas of the domains 1 and 2 are specified by A 1 = 13.4 and A 2 = 17.765 mm 2 , respectively. This means that for the root weld in the two-step welding process A * = A 1 + A 2 , whereas for the second step A * = A 2 . The nominal power input, P w0 , is a reduced value in the sense that it also models the energy loss due to heat convection from the weld surface. Additional discussion of (41) can be found in Appendix A. Hence, the standard power values for the two-step analysis were P w = 615W and P w = 711W for the root weld and for the second weld, respectively. These input values are consistent with the experimental conditions. Table 4 lists the standard parameters used for the welding process in the simulations.

Numerical results: prerequisites
The results are shown in terms of temperature and stress distributions. The geometry of the analysed problem is shown in Fig. 12, where the weld region is indicated by the black line. In most cases below, the results are shown for the region in the vicinity of the weld, i.e. the blue dotted box in Fig. 12.
The predictions of stress entities along three specific paths are also considered, i.e. lines 1-3, as shown in Fig. 12. Lines 1 and 3 are horizontal lines along the surface and at a distance of 3 mm below the surface, respectively, whereas line 2 is the vertical, interior symmetry line. Stress distributions along these lines may be compared with the experimental results presented in [56]. The coordinates x 1 and x 2 , associated with the deformed configuration, are used to indicate the position along the different lines.
A convergence study was made for the spatial discretization, and the mesh adopted in the analyses below enables the converged results, i.e. further mesh refinement does not significantly affect the outcome of the analysis.

Numerical results: weld analysis with only isotropic hardening
In Fig. 13, the residual stress state from the two-step welding procedure is shown. The stress state is represented by the effective von Mises stress, and the T 11 -and T 33 -components. The weld domain is indicated by the dashed black line.
The von Mises stress reaches a peak of about 700 MPa in the base material close to the welding area. This is higher than the initial yield stress, and it is clear that the material close to the weld has undergone plastic deformation during the solidification process.
The stress in the horizontal direction, T 11 , does not reach as high a magnitude as the effective stress. The magnitude is actually the highest for the compressive stresses that appear in the welding zone and above all at some distance below the welding zone. These stresses are of the order of −300 MPa. The distribution of T 33 shows that the out-of-plane stress component dominates the residual stress state and is the main contributor to the residual effective stress. The T 33 stress reaches its peak at the bottom of the weld zone and just below it, with a peak stress of about 600 MPa. Also, at a distance of about 20-30 mm from the weld zone there is a large area of material under compression at about −300 MPa. Figure 14 shows some sample images from the simulation for the evolution of the temperature field and the associated effective stress. For clarity, only the region closest to the weld is shown.
At t = 0, the whole weld region ( 1 and 2 ) is given the initial temperature T = T wp , and the structure is stress-free. After a few seconds, the weld region is essentially stress-free but it has been heated by the heat source. At this stage, the region surrounding the weld experiences a temperature gradient which introduces distortions that cause an effective stress. At the end of the first step, the whole structure has cooled down to room temperature, but there is a remaining residual stress state. As mentioned above, during this first welding step, the domain 2 is (artificially) assigned liquid-like properties independent of its temperature, and for this reason, 2 remains stress-free at the end of the first step.
At the start of the second welding step (t = t 1 + t 12 ), the temperature of the domain 2 has been raised to T wp . The structure then undergoes a heat cycle similar to the first step with heating and cooling. At the end of the second welding step, there is a significant residual stress state as seen at the bottom of Fig. 14. The peak stress in the residual field exceeds 600 MPa. It is clear from Fig. 14 that the effective residual stress after the second step is higher than that after the root weld.

Numerical results: influence of heat input
The influence of changes in the heat power input, P w0 , is now investigated. The standard value, P w0 = 837 W, was both increased and decreased by 25%, yielding the alternative power inputs 628 and 1046 W. Simulations of the two-step welding process were performed using these alternative power inputs, and the resulting stress states are shown in Fig. 15.  Figure 15 shows the resulting distributions of the effective stress, σ e , for the three different values of P w0 . Only the region around the weld is shown. In and in the vicinity of the weld, the residual effective stress seems to decrease somewhat with increasing power input. The reason for this tendency seems to be that when a plate of finite size is infused with heat energy, the temperature distribution in the plate becomes more homogeneous the more heat energy is infused into the plate. Hence, more heat input means that the intricate expansion and contraction processes in the material and the cooling process take place under the influence of weaker temperature gradients, which suppresses the formation of residual stresses. The maximum effective stress in the vicinity of the weld is about 750, 700, and 630 MPa for P w0 = 628, 837, and 1,046 W, respectively.
It is also noted that the resulting deformed geometry of the structure depends on the power input. The highest power input produces a remaining bulge, whereas the lowest power input results in a small valley at the weld. For the standard power input (P w0 = 837 W), a small bulge also remains at the end of the welding process.

Numerical results: influence of hardening model
The solutions from the present analysis are now compared with the experimental results from the literature [56]. In addition, some results from analyses using a Lagrangian constitutive model [31] have been included for comparison. Lindström [31] performed FE simulations of the same problem as has been studied here. In  that study, the welding problem was analysed using a thermomechanical, staggered coupled approach, and the mechanical material model used was a nonlinear, thermoelastic-plastic model that used a mixture of isotropic and kinematic hardening. Lindström [31] did not publish predictions for all stress contours that are discussed below, but for the cases where they are available, they are included in the presentations below.
In Figs. 16, 17, 18, 19, and 20, numerical predictions of T 11 and T 33 along lines 1-3 are shown together with the experimental results. When performing the experimental measurements, different laboratories were involved, and different experimental techniques were used for the measurements, see Wohlfahrt et al. [56]. Therefore, different symbols have been used to denote the experimental results from each of the laboratories.
It is evident that the scatter in the experimental data is significant and that the differences in measurements from different laboratories are sometimes dramatic. As discussed in [56], these differences may partly be attributed to surface treatments performed in the laboratories before the residual stress measurements were made. Figure 16 shows the distribution of T 11 along line 1 in Fig. 12. The peak stress in the experiments varies between 300 and 400 MPa. The difference between the models is not dramatic, but the models underestimate the peak stress by almost 200 MPa. Also, the main stress peak in the experiments appears closer to the symmetry plane of the weld compared to the peak in the numerical predictions.  Figure 17 shows the distribution of T 33 along line 1 in Fig. 12. The peak stress in the experiments ranges between 250 and 500 MPa, and the models predict a peak stress of about 500 MPa (model A) and 400 MPa (model B). The peak stress in the solution from Lindström [31] is between that predicted by models A and B. Also, the location of the stress peak is predicted fairly well by the simulations, i.e. between 5 to 10 mm from the symmetry line. In general, the models are able to predict the stress level in the weld region relatively well. In the region closest to the weld, the model with both isotropic and directional hardening predicts a stress level that is about 100 MPa lower than the model with only isotropic hardening.
In Fig. 18, stress contours for the vertical symmetry line (line 2) are shown. Starting with the results for T 11 , it is noted that all of the model predictions agree fairly well with the experimental data. For this case, the results are available from Lindström [31]. In particular, all three numerical simulations predict a peak stress of about 150 MPa, which agrees well with at least one of the experimental sets. The location of the peak is also relatively well predicted by the numerical models. In general, the solution from Lindström [31] seems to follow more closely to model B, which is natural, since both of these models account for directional hardening.
With regard to the distribution of T 33 , the model predictions clearly overestimate the peak stress compared to the experimental measurements. Model A overestimates the peak stress by several hundred MPa, and the peak stress in model B is about 150 MPa lower. Hence, in the predictions of T 33 there is a significant influence of directional hardening.
The distribution of T 11 is shown in Fig. 19 along line 3 in Fig. 12. The results from the models have been extended to the negative side assuming symmetry. The experimental results are ambiguous. However, it may be concluded that there seem to exist stress peaks at a distance of about 5 mm from the central line of the structure, and the magnitude of these peaks seems to be of the order of 100 MPa. The magnitude of this peak is well predicted by model A, whereas model B predicts a somewhat lower peak stress. The location of the peak is predicted by the models to be about 15 mm from the central line, which is about 10 mm from the experimental results.
Finally, the distribution of T 33 is considered in Fig. 20 along line 3 in Fig. 12. Both of the models from the present study as well as Lindström's solution overestimate the stress level in the vicinity of the weld region. Both models and experimental estimates suggest the existence of a peak stress at about 5-8 mm from the central line, the magnitude being somewhat lower in the experiments (250-300 MPa) compared with the numerical predictions (480-580 MPa). There is again a significant difference between the predictions from models A and B, where the peak stress in model B is about 150 MPa lower than in model A. The predictions of model B and from Lindström both have a peak stress of about 450 MPa. Again, there is a significant influence of directional hardening.

Discussion and concluding remarks
The process of welding is relatively complex one from a thermomechanical point of view because it includes such phenomena as melting, solidification, and phase transformations. The large changes in temperature during the process cause dramatic changes in the elastic and viscoplastic properties of the materials involved. The combination of these factors determines the resulting residual stress state, which is often of great interest for predicting the risk of fatigue in welded structures.
The present work uses an Eulerian formulation of elastic-viscoplastic response which is thermomechanically consistent. In contrast to standard Lagrangian formulations for large deformations, the Eulerian formulation proposes an evolution equation directly for a second-order unimodular tensor that measures elastic distortional deformations. In particular, this Eulerian formulation is not sensitive to arbitrariness of choices of the reference configuration, an intermediate configuration, a total deformation measure, or an inelastic deformation measure. Since there is no need for an inelastic deformation measure, there is also no need for a plastic potential. However, it is possible to introduce a scalar measure of equivalent plastic strain from the initial state for engineering analysis. Also, the constitutive equation for the rate of inelastic deformation Γ allows for a smooth transition from viscoelastic fluid-like response in the melted region to elastic-viscoplastic response of the solidified material. It is expected that this feature of the model is particularly appropriate for modelling residual stresses in the weld material, which has experienced melting and solidification.
Overall, the Eulerian elastic-viscoplastic formulation adopted here seems to be able to handle the complex processes in welding well. However, the numerical simulations are not fully able to predict the experimental measurements in terms of residual stresses. Quantitatively, there were discrepancies in peak stress levels and the locations of these peak stresses. There may be several reasons for these discrepancies. First, it must be noted that there are notable differences in measurements between different laboratories, which suggests that there are significant uncertainties in the measurements themselves. Second, it might simply be that the 2D representation of the problem employed here is too simplistic and that a full 3D representation is needed. For example, it is evident that the out-of-plane stress component T 33 is the dominant stress component in the residual stress state, and the generalized strain boundary condition may not represent the 3D situation well enough (cf. [5,34]). Third, the present study is the first attempt to apply this Eulerian viscoplasticity framework to this type of problem. The constitutive model was calibrated using a relatively detailed set of the experimental results. The simulation of the welding process and the resulting residual mechanical state was truly a pure prediction without any attempts to fine-tune the framework in order to optimize the outcome with regard to the experimental measurements of residual stresses. Fourth, the modelling of the heat source and the thermal process may also introduce errors in the numerical predictions. Unfortunately, separate measurements for the thermal process are not available to validate this part of the simulation.
Previous attempts at predicting the residual stress state in the present type of material and for the same type of welding process (e.g. [35,52]) also show the difficulty of accurate predictions. The study by Smith and Smith [52] points to the need for fine-tuning of the modelling framework for accurate numerical predictions of the residual stress state. In particular, the out-of-plane stress component is difficult to accurately reproduce in simulations.
It is clear that introducing directional hardening is of significant importance, and the difference between the model predictions of stress levels for isotropic hardening only and combined isotropic and directional hardening was about 150-200 MPa in the most extreme cases. The same trend was observed by Lindgren [28] for kinematic hardening. The influence of the hardening model seems to come into play primarily during the cooling phase of the welding cycle, as the material in the vicinity of the weld is subjected to reversed plastic yielding. Proper treatment and modelling of reversed plastic loading thus seem to be very important. In general, it seems that introducing combined hardening brings the model predictions closer to the experimental results. Additional new features can be added to this Eulerian framework to improve the numerical predictions further.
The importance of the choice of the hardening model (isotropic, kinematic) has been investigated in other studies, for instance by Bammann and Ortega [3] and Muransky et al. [35], who investigated the effect of assuming isotropic and/or kinematic hardening. These studies also confirm that the choice of the hardening model has a very strong influence on the predicted residual stresses close to the weld.
The predictions of peak stresses and the associated measurements suggest that the residual stresses are of the order of the yield stress of these materials at room temperature. The residual stress state is therefore expected to strongly affect such phenomena as fatigue crack growth. Furthermore, variation in heat power input was considered, and it was demonstrated that a relative change of 25% in heat power input may significantly affect both the residual stress state and the resulting weld geometry. In this regard, it is noted that varying the heat power input is equivalent to varying the welding speed, see (41), so the cooled residual stress state and shape of the weld might be optimized by controlling the welding speed.
Within the context of CWM, there is no consensus as to how to apply the heat to the weld area. In the present study, heat was supplied by a combination of a prescribed temperature and a prescribed heat flux, using a 2D version of the heat source proposed by Goldak et al. [15]. It should be mentioned that other models of heating during the welding process have been proposed, see for example [7,16,23,39].
In many material models used in CWM, residual stress release at elevated temperature (annealing) is accounted for by a reduction or removal of accumulated plastic strains. The elastic-viscoplastic model used in this study does not use a plastic deformation measure. Instead, it accounts for thermal recovery and annealing by allowing the state variable κ (which models hardening) to decrease continuously at increased temperatures.
In summary, an Eulerian elastic-viscoplastic formulation has been applied to the analysis of a welding process. The residual stress state was quantified in terms of the normal Cauchy stress components and the effective stress, and the maximum residual stresses were found to be of the order of the yield stress of the material at room temperature. Variations in the heat input had a significant influence both on the residual stress state and on the resulting geometry of the weld. It was also clear that the choice of the hardening model has a significant influence on the residual stress state, especially on the peak levels of stress.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Numerical implementation of the thermomechanical model
Consider an arbitrary time step which starts at t = t 1 , ends at t = t 2 = t 1 + Δt with time increment Δt. At t = t 1 , the set of values {J e (t 1 ),B e (t 1 ), κ(t 1 ), β(t 1 )} is known. The numerical integrator proposes values of these quantities at the end of the time step, i.e. {J e (t 2 ),B e (t 2 ), κ(t 2 ), β(t 2 )}. In addition to the necessary state variables, the auxiliary variable ε e p is updated from ε e p (t 1 ) to ε e p (t 2 ). Following the work in [36,48,51], it is possible to develop an implicit robust and strongly objective numerical algorithm for integrating the evolution equations by introducing the relative deformation gradient, F r , which is introduced as Also, the relative dilatation is defined by J r = detF r , andF r = J −1/3 r F r is the unimodular part of F r . In the following, the notation '(t 2 )' is omitted in most cases for clarity. Hence, the values of state variables at the end of the time step, i.e. t = t 2 , are intended, unless otherwise specified. It can be shown that the deviatoric part of the evolution equation (4) becomes (cf. [46]) Then, the exact solution of (3) is given by Furthermore, it is convenient to introduce the elastic trial solutionsB * e andB * e by the expressions The evolution in (43) is then approximated bẏ which guarantees an exact solution for elastic response (Γ = 0). Using a backward Euler approximation of the derivative, the solution of (47) is written in the implicit form with λ = 1 1 + ΔΓ , ΔΓ = ΔtΓ, (49.1,2) where Γ = Γ (t 2 ) [20,21]. Then with the help of (19.2) and (48), it follows that The evolution equation for κ in (26) is discretized, and an implicit estimate of κ at the end of the time step yields Also, discretization of (27) and assuming implicit integration yields where Then, insertion of into Eq. (49.2) gives a nonlinear function that implicitly defines ΔΓ , which can be solved by an iterative solution procedure. In the solution procedure, trial values of the state variables are first computed, i.e. values assuming elastic deformation (ΔΓ = 0). In addition toB * e , the other trial values are evaluated according to κ * = κ(t 1 ) + Δtm 2 κ a 1 + Δtm 2 , Hence, ΔΓ and the updated values of the hardening parameters κ and β are determined by the numerical scheme below: -Evaluate trials:B * e , γ * e (B * e ), κ * , β * , Z * -If γ * e ≤ R Z * then -ΔΓ = 0,B e (t 2 ) =B * e , κ(t 2 ) = κ * , β(t 2 ) = β * -Else -While |δ(ΔΓ )| > tolerance -Evaluate κ, β, Z for ΔΓ I where dκ Once ΔΓ has been determined,B e (t 2 ) is determined using (48). Finally,B e (t 2 ) is expressed in the form where α 1 is determined by the condition detB e (t 2 ) = 1, which is a cubic equation with a closed form solution [47]. The Kirchhoff stress is given by (20) and is a function of F r and θ . The variation of τ is given by δτ = (∂τ /∂F r ) : δF r + ∂τ ∂θ δθ = (∂τ /∂F r )F T r : δF r F −1 see Jabareen (cf. [22]). The consistent tangent moduli, C and c, can then be identified as For the mechanical part, the stiffness C = 1 J e K (1 − f 2 (θ ))I ⊗ ∂ J e ∂F r + μ ∂B e /∂F r F T r (64) needs to be determined. Recalling that the current value of J e is given by J e = J e (t 1 )J r , it follows that The second derivative needed is given by Differentiation of (45.2) yields where ∂ F rBe (t 1 )F T r /∂F r = (I F r + F r ⊕ I)B e (t 1 ), and the notation (Q R) i jkl = Q ik R jl and (Q ⊕ R) i jkl = Q il R jk has been introduced for convenience. Differentiation of (49.2) and some restructuring then yields where N 2 is defined by the entities where I is the fourth-order unit tensor. Furthermore, a fourth-order tensor D is defined as Finally, the total mechanical stiffness is given by Above, v is the speed of the heat source that is moving in the (negative) X 3 direction, t denotes time, and a, b, c f , and c a determine the shape of the heat source. The total power input into the X 1 -X 2 -plane (for a given value of X 3 and t) may then be computed as In the present 2D analysis, the heat source is located within a half-circle with radius l w , as indicated in Fig. 10. Within this half-circle, the intensity is taken to be constant and to be given by where X 30 is a position where initially (i.e. at t = 0) the intensity is very low. With regard to the initial temperature of the weld pool, it is noted that raising the temperature of the weld pools in the domains 1 and 2 from room temperature, T rm , to the weld pool temperature, T wp , requires an energy input of where A * is the area of the preheated domain that was defined in Sect. 5, H is the thickness of the 2D model, and The energy supplied from the heat source to the weld pool during the subsequent, transient simulation is Since the total energy input in the simulation must correspond to the energy input in the experiments, the amplitude of the heat source is decreased to account for the initial energy of the weld pool. Hence, energy considerations require that where E tr0 = P w0 H/v is the reference heat supplied by the heat source during a transient simulation with no preheating of the weld pool. From (87), it may be deduced that the power amplitude, P w , must take the form