Residual stresses in gas tungsten arc welding: a novel phase-field thermo-elastoplasticity modeling and parameter treatment framework

The fusion welding process of metallic components, such as using gas tungsten arc welding (GTAW), is often accompanied by detrimental deformations and residual stresses, which affect the strength and functionality of these components. In this work, a phase-field model, usually used to track the states of phase-change materials, is embedded in a thermo-elastoplastic finite element model to simulate the GTAW process and estimate the residual stresses. This embedment allows to track the moving melting front of the metallic material induced by the welding heat source and, thus, splits the domain into soft and hard solid regions with a diffusive interface between them. Additionally, temperature- and phase-field-dependent material properties are considered. The J2 plasticity model with isotropic hardening is considered. The coupled system of equations is solved in the FE package FEniCS, whereas two- and three-dimensional initial-boundary-value problems are introduced and the results are compared with reference data from the literature.


Introduction
In fusion welding processes, two or more metallic parts are permanently joined by applying a suitable heat source, e.g. arc, laser, or electron beam sources. Such a way of joining is of particular importance, especially for metallic structures, such as pressure vessels and piping systems. Due to the highly-localized nature of the heat source, the highest temperatures occur in the zone near the welding torch and decrease with increasing distance from the weld centerline. As a result, welded parts undergo higher nonuniform thermal expansion followed by contraction due to the heating-cooling cycle, which in turn leads to the generation of plastic deformations within and around the weld zone. The expansion and contraction result in high-level tensile residual stresses B Yousef Heider heider@iam.rwth-aachen.de Baharin Ali baharin@iam.rwth-aachen.de Bernd Markert markert@iam.rwth-aachen.de 1 Institute of General Mechanics, RWTH Aachen University, Aachen, Germany and distortions when the weld structure cools to the ambient temperature. The presence of weld-induced tensile residual stresses has detrimental effects on load-carrying capacity, fracture toughness, stress corrosion resistance, and fatigue crack initiation and propagation when subjected to cyclic loading [2,11,14,54,66]. Therefore, understanding and predicting weld-induced residual stresses and deformations in the welded structures are of high significance and should be accounted for in the early design stages to avoid possible reductions in their performance and reliability.

Thermo-elastoplastic models of residual stresses computation
The presence of residual stresses in the welded structures is complex to be only calculated by performing experiments or applying analytical formulations. For instance, the experimental measurements of residual stresses and distortions involve some shortcomings related to the complexity, cost, and difficulty of obtaining complete and detailed values. For these reasons and due to the advancements in numerical techniques and computational capacities, researchers have focused on developing and implementing numerical models. The developed models are then validated and cal-ibrated with the experimental measurements using, e.g., the Finite Element Method (FEM). These FEM approaches have shown acceptable capabilities to perform computational analyses and to predict the residual stresses for different fusion welding processes in multi-dimensional space, e.g., [12,22,27,79,80,110,[116][117][118]. Other works have intensively studied the effect of finite and small strain theory on the numerical results of the welding deformations, e.g., [28,31,76,116].
In welding computational analysis, it is a common practice to use temperature-dependent material properties. However, at high temperatures, the complete temperature-dependent properties for many materials could be very difficult to measure and are not available in the literature [42,115]. The parts being welded undergo complex multi-physical processes that involve heat transfer, material flow, and other metallurgical changes. As a result of these coupled phenomena, the physical and chemical properties of the final welded parts are usually changed. Therefore, the need for considering the dependency of material properties on temperature has been intensively investigated in the literature, where assumptions and simplification approaches could also be found. In most cases, numerical modeling approaches related to welding are evaluated by their accuracy and performance to predict the transient temperature field and the welding-induced residual stresses and permanent deformations [7,8,13,15,42,55,61,67,84,115]. For instance, Little and Kamtekar [67] have found that varying the thermal properties, namely the heat capacity and the thermal conductivity, highly influence the computed temperature fields in the weld and its surrounding regions. Zhu and Chao [115] performed an analysis on welding of an aluminum plate using three different sets for the thermo-mechanical properties. The study showed that the temperature distribution, residual stresses, and distortions can be predicted with sufficient accuracy by taking constant room-temperature properties, except for the yield stress. Armentani et al. [8] showed how the residual stresses of butt-welded joints can be affected by thermal properties and concluded that varying the thermal conductivity has a significant influence on the calculated residual stresses. Barroso et al. [13] applied a simplified material model in their FEM simulations and found that the residual stresses in the longitudinal direction can be predicted with reasonable accuracy when the mechanical properties are taken as constant values. Joshi et al. [55] performed a FE analysis of two overlapping beads, where the chemical compositions of the base and deposited metal were used in the Weldware package to derive the equivalent material properties in the weld pool. Bhatti et al. [16] investigated the influence of thermo-mechanical material properties of various grades of steel on the predicted residual stresses and distortions. They found that considering constant properties, except for the yield stress, can provide good results for residual stresses. Additionally, they concluded that acceptable angular distortion can be achieved when the yield stress, the heat capacity, and the thermal expansion coefficient are functions of temperature. Similarly, Perić et al. [84] investigated the influence of the temperature-dependent properties on the FE model outputs, and suggested that using constants properties, except for the yield stress, is suitable for the prediction of the deformations. Moreover, temperaturedependent properties can be used for the assessment of the temperature distributions, the heat-affected zones, and the residual stresses. Kong and Kovacevi [61] reported that assigning a small value of the elastic modulus at the weld pool for temperatures exceeding the melting point can reduce the computational time and increase the solution convergence rate. Anca et al. [7] studied the melting-solidifying behavior on the residual stresses calculation during the liquid/solid and solid/liquid phase changes. For the mechanical analysis, the weld pool was assumed as a soft region with no further changes in the material properties above a specific cut-off or zero-strength temperature. He et al. [42] introduced an additional temperature-dependent internal state variable to capture the irreversible welding interfaces and accounted for the effects of the melting state on the thermomechanical properties during the lap joint welding process of two sheets. The aforementioned research works showed the importance of considering state-or temperature-dependent material parameters in the FE simulations to get reasonable numerical results, as for the residual stresses. They also showed that these considerations are associated with many challenges concerning the parameter values measurement or estimation. This leads to the conclusion that further investigations and reduction of model complexity are badly needed to achieve more accurate numerical computation of, e.g., residual stress calculations, at moderate efforts and assumptions.
While utilizing the phase-field method to capture the melting of metal during welding, the underlying work considers the phase-change effects in both the thermal and the mechanical analysis. In particular, it considers phase-fielddependent material properties, such that the model domain is divided into a hard (unmelted) and a soft (melted) region. These regions are separated by a diffusive transition interface, where the material properties can smoothly change. Hence, soft material properties are assigned to the fully melted region in the weld pool, whereas hard material properties are assigned to the material below the solidus temperature.

Phase-field modeling of thermally-induced phase transition
In the literature, several approaches exist to describe the process of phase transition of phase-change materials, such as during melting, and to determine the position of the corresponding interface between melted and unmelted phases.
These methods, such as the level-set [81,104], the volumeof-fluid [51], and the phase-field method (PFM) allow to describing the phase change within continuum mechanical framework as needed in the current research work. The developments and uses of the PFM have been discussed by many research groups [18,23,34,58,98,109,111,114]. In this, the PFM is originally applied in the mechanics of materials to describe the complex microstructural evolution of the interfaces in non-equilibrium states [18,58,59,98].
Here, we utilize the energy-based PFM to describe the melting of metals as a phase-change phenomenon between a hard state for temperatures below the melting point and a soft state for temperatures above the melting point, whereas the thermal energy during welding is the driving force for the phase change. To this, a phenomenological phasefield variable is used to indicate the state of the metal, which allows for a diffusive state change across the interface between the two states. Unlike the sharp-interface approaches for describing the phase-change materials, the diffusive-interface approaches allow a much easier and stable FE implementation, whereas the finite thickness of the interface has its origin in the lower-scale phase-change description, see, e.g., [18,86]. In this regard, the PFM is applied in many research works to describe material phase change on the continuum scale, see, e.g., [101,102,[112][113][114]. In this case, the macroscopic interface thickness, which is dictated by the FE mesh size, could be orders of magnitude larger than the physical interface thickness. However, studies related to the variation of the interface thickness parameter in the PFM on the macroscale show that below a certain limit, decreasing the value of this parameter will have an insignificant effect on the overall solution or the interface position [19,103]. To this, references and overview of the origin of this method in the context of continuum mechanics can be found in, e.g., [57] and [102,103]. Apart from PFM of phase-change materials, the PFM has recently been intensively applied within fracture mechanics for solids as well as for porous media as in, e.g., [3,5,[44][45][46]48,82,83], among others. It is also worth mentioning that several research works have recently employed the PFM to study the microstructure evolution, such as the investigation of the growth of the columnar grains accompanied by the solidification process in welding and additive manufacturing process, for example, the works done in [36,63,90,107,108]. Therefore, the underlying modeling framework together with the developed algorithms will be implemented in future works in the context of modeling direct metal laser melting (DMLM) in the additive manufacturing process. Additionally, future works will address the possible embedment of machine learning approaches, like in [35,49,50,60,99,100], to capture the microscopic thermo-mechanical processes in the continuum modeling.

Highlights and content overview
In summary, the main objective of this work is to gain new insight into the challenging coupled inelastic thermomechanical processes that occur during gas tungsten arc welding (GTAW) via applying an advanced numerical study and comparison with experimental and numerical data from the literature. The proposed modeling methodology allows an accurate estimation of the spatial and temporal distribution of the residual stresses and the permanent plastic deformations. This approach includes: (1) Presenting and implementing a continuum mechanical framework with unified kinematics of the melted and unmelted phases of the solid metal. (2) Applying a thermodynamically consistent phase-field approach to describe study the effect of the phase transition during melting/solidification processes on the size and the shape of the weld pool. (3) Proposing a new ansatz for the phase-fielddependent material parameters, which reduces the modeling complexity and gives a realistic description of the properties' change across the diffusive interface. (4) Investigating two different plasticity models and describing their implementation in the open-access/open-source FEniCS package to solve two-and three-dimensional initial-boundary-valueproblems (IBVPs).
To give an overview, Sect. 2 describes the essential kinematics and the formulation of the coupled phase-field thermo-elastoplasticity within a thermodynamic framework. In Sect. 3, the material properties are properly formulated to be functions of the temperature field as well as the phasefield variable, so that their values change smoothly across the diffusive interface. The finite element implementation is discussed in Sect. 4, which includes writing the weak formulation and discussing the discretization with mixed finite elements and the time-stepping via the backward finite difference method. Two-and three-dimensional numerical examples are presented in Sect. 5, where the proposed models are implemented to solve IBVPs. These implementations allow to verify and validate the numerical model through comparison with reference data from the literature. This is followed by the conclusions in Sect. 6, where also a summary and an outlook of future works are presented. As the focus of the current work is to develop a themomechanical model for the prediction of residual stresses, currently small strain theory is considered and geometrical nonlinearity will be involved in future works.

Mathematical formulations of the phase-field thermo-elastoplasticity model
In this work, a coupled phase-field thermo-elastoplasticity problem is considered for the modeling of the GTAW welding process and the occurring residual stresses, where the material modeling is addressed within a thermodynamically consistent framework. The three primary fields considered in the following treatment are functions of the space x and time t. These are the vector-valued displacement field u = u(x, t), the scalar-valued temperature field θ = θ(x, t), and the scalar-valued phase-field variable φ = φ(x, t). A brief description of the phase-field modeling of temperatureinduced melting of metals is presented in ("Appendix A"). However, in the following treatment, the phase-field evolution equation along with the other constitutive formulations are derived based on the restrictions imposed by the Clausius-Duhem inequality.

Kinematics of the thermo-mechanical problem
The formulation of the thermo-elastoplasticity in the current contribution is restricted to the small deformations assumption. Therefore, the total strain tensor ε can be expressed in terms of the displacement gradient and its transpose as For the thermo-elastoplasticity treatment, ε can additively be decomposed into an elastic strain tensor ε e and an inelastic (plastic) strain tensor ε p as

Constitutive formulations of the phase-field thermo-mechanical problem
In analogy to the work presented in [4] in the context of gradient thermo-plasticity, the focus in the constitutive formulations will be on the following set of independent state variables where α is the equivalent plastic strain, which is chosen to characterize the state of the material hardening and φ is the phenomenological phase-field variable (see, "Appendix A"). Therefore, the free energy function Ψ (V) for the coupled problem can be expressed as The above definitions of the free energy densities (per unit volume) can be found, e.g., in [4,6,78,95]. In this, Ψ e is the elastic stored energy, Ψ p is the plastic contribution due to hardening, Ψ e-th is the energy term due to thermoelastic coupling, Ψ φ is the phase-field contribution, and Ψ th is the purely thermal energy. H is the linear hardening modulus, which is taken (H ≈ 0.01E) in analogy to [92], c ε is the volumetric heat capacity, and K, μ are the bulk and shear moduli, respectively. Moreover, e = tr ε e andε e = (ε e − 1 3 e I) are, respectively, the trace and the deviatoric (isochoric) part of the elastic strain tensor ε e . Furthermore, α θ is the temperature-dependent coefficient of linear expansion, θ ref is a specified reference temperature, and I is the identity tensor. In connection with the temperature-related phase change of the metallic material, the phase-field related energy function Ψ φ is expressed in terms of a bulk free energy contribution f (φ, θ ) and the interfacial energy terms. This can further be expressed, in analogy to [18], as where W , g(φ), p(φ) andθ M are the energy barrier height at the interface, the double-well potential, the interpolation function, and the material melting temperature, respectively. Moreover, ρ, L and l are, respectively, the material density, the latent heat of fusion, and the gradient-energy coefficient that affects the width of the diffusive interface. In the current work, we consider a constant reference material density at room temperature. For the heat capacity, thermal conductivity, Young's modulus, yield stress, and thermal expansion coefficient will be either temperature-or phase-depended as discussed later in Sect. 3 and given in Table 2 and Figs. 1 and 2.

Plasticity-related response and yield functions
We proceed in the following discussion with the dissipative force fields that represent the dual fields to the internal state variables of the plastic response, i.e. σ D dual to ε p and q α dual to α . For the J2 plasticity with isotropic hardening, the von Mises yield criterion can be defined according to, e.g., [70,94,96], as where σ D , q, and y 0 (θ ) are the deviatoric part of the stress tensor, the internal stress-like dynamic variable, and the initial yield stress, respectively.

The free energy and the entropy evaluation
With the given definitions of the energy contributions in Eqs. (4) and (5), the total free energy per unit volume can now be expressed as The time derivative of the energy formulation, needed for the entropy inequality evaluation, can be expressed aṡ Using the divergence theorem in analogy to [32], the last term in the above derivative can be reformulated as follows Following this, we rely on the 2 nd law of thermodynamics (entropy inequality) in order to generate thermodynamically consistent constitutive formulations for the coupled problem under consideration. This inequality can be expressed in the Clausius-Duhem general representation as where η and q θ are the entropy and the heat flux vector. Substituting the formulation in (8) into the entropy inequality (10) and applying reformulation yields The evaluation of the entropy inequality (11) towards deriving thermodynamically consistent constitutive formulations is applied following the procedure presented in [24]. In this, the fulfillment of the inequality is guaranteed if each term is greater or equal to zero. Thus, we conclude the following equilibrium relations: Analogous to the continuity equation in fluid mechanics, relation (12d) describes the continuity of the vector term in the parentheses. Following this, integration of (12d) over the volume of the domain and applying the Gaussian integral theorem as presented in [32] yields where n is the outward unit surface normal. The resulting equation (13) can be realized in the numerical implementation through applying the following Neumann boundary condition: Additionally, the evaluation of (11) results in the following three non-equilibrium (dissipative) restrictions representing the thermal, phase change, and mechanical dissipation contributions, respectively. In the case of the thermal dissipation or the so-called Fourier's inequality in (15a), this is related to the existence of the temperature gradient with the restriction such that the heat must flow from warm material points to colder ones [26,62]. In this regard, the constitutive equation for the heat flux vector may be given by Fourier's law q θ := −κ grad θ .
To satisfy the dissipation inequalities in (15b) and (15c), proper evolution equations should be formulated. First, an evolution equation for the phase-field variable can be introduced by assuming the following natural proportionality [32,71] where M > 0 can be interpreted as the interface mobility parameter, which is discussed in "Appendix A" and expressed according to [18] in Eq. (49). It is worth mentioning here that the phase-field evolution equation can be derived based on other approaches, such that presented by Gurtin [41] and proceeded from the microforce equilibrium laws. This approach is however beyond the scope of the current work. In addition, due to the time history of the plastic deformation, the first and second terms in (15c) do not vanish [29]. The thermodynamic variable q α conjugate to α can be expressed as Having the latter definitions, an expression for non-negative reduced dissipation, or the so-called mechanical dissipation, given by the Clausius-Planck inequality, can be written as

Evolution and governing balance equations
The starting point in the derivation of the evolution equations forε p andα is the application of the principle of maximum plastic dissipation, which is widely used in connection with a plasticity variational formulation [96]. In this regard, the yield surface (6) of the associative plasticity model together with the plastic dissipation relation (18) are employed to transfer the maximum plastic dissipation into a minimization problem by introducing the Lagrangean functional [29] whereλ ≥ 0 is the Lagrange (or plastic) multiplier, which enforces the restriction χ ≤ 0 on the yield surface potential. Finally, the evolution equations can be obtained by taking the derivative with respect to σ and q α as Following this, the governing balance equations in their local form will be introduced. The balance of linear momentum can in general be defined as whereü is the acceleration vector and b is the body force vector, defined per unit volume. Making use of the balance of internal energy, the scalar-valued local energy balance equation can be expressed aṡ In this,Ė and r are the internal energy rate and the internal energy supply per unit volume. A relationship between the energetically conjugate variables θ and η is obtained via the application of Legendre transformation for the internal energy Recalling the reduced local dissipation relation (18), the definition of the time derivative Eq. (8) can be expressed as followṡ With the latter definition of the internal energy densityĖ, Eq. (23) can be rewritten as Following this, we derive an explicit relationship to the entropy rateη in analogy to the procedure presented in [4]. In nutshell, we start from the general definition of the entropy η := − ∂Ψ (V) ∂θ with V = {ε − ε p , α, θ, grad θ, φ, grad φ} and carry out the time derivative using the chain rule to obtaiṅ η. This leads to the following relation By subtitling Eq. (27) into (26), we finally obtain the governing energy balance equation in its local form with where c and H ep represent the volumetric heat capacity and the structural elastoplastic heating effect. The term, H pc , represents the phase-change contribution. For linear elastic isotropic materials with temperaturedependent elastic modulus 4 C, the stress given in Eq. (12a) can be expressed using the additive strain decomposition in Eq. (2) as in this, ε th can be expressed within the isotropic small strains framework as Following (30), the rate form can be expressed aṡ To account for plastic strains, a rate-independent associative plasticity model with J2 von Mises plasticity and isotropic strain hardening law is used. The derivation of the plasticity model is presented briefly in the following, whereas a detailed description can be found in, e.g., [1]. According to Eqs. (20) and (21) with the definition of the yield function in (6), the evolution for both plastic strain and equivalent plastic strain can be established as followṡ Incorporating the consistency conditionχ = 0 and using the rate of the yield stress function, one finds another relation for the rate of the equivalent plastic straiṅ Combining equations of the equivalent plastic strain (33) 2 with (34) and using the stress rate Eq. (32) augmented with the flow rule (33) 1 , the plastic multiplier evolution can finally be written aṡ Onceλ is calculated, the rates of stressσ and equivalent plastic strainα can be computed, as will be discussed in details in Sect. 4 within the numerical treatment.

Summary of the governing balance relations
Following the previous discussion, Table 1 summarizes the local system of coupled nonlinear partial differential equations (PDEs), which are needed to solve IBVPs of phasechange thermo-elastoplasticity.

Temperature and phase-field dependent material properties
In the numerical simulation of the welding process using the FEM, the choice of the material properties and their dependencies on the temperature or the material state (melted/unmelted) plays a vital role in the thermal and mechanical analysis [40,64,93]. For the transient temperature field computation, the material density (ρ), the thermal conductivity (κ ), and the heat capacity (c) are of particular importance. The Young's modulus (E), the yield stress (y 0 ), the hardening parameter (H ), the Poisson's ratio (ν), and thermal expansion coefficient (α θ ) play the major role in the calculation of elastoplastic deformations and, thus, the residual stresses. However, Poisson's ratio has only a slight effect on the final state of the residual stress as discussed in [21,105]. In the literature, it is a common practice in the thermomechanical analysis of the welding process in metals to consider the values of ρ and ν to be independent of the Heat capacity mJ/(mm 3 Yield stress MPa 8 · 10 −6 Soft/hard-state-related values are adapted from [31] temperature or the material state, whereas their values are usually taken equal to that measured at room temperature. This simplification will also be applied in the underlying work. However, the consideration of the dependency of the other material properties, i.e. E, c, y 0 , H , κ, and α θ , on the temperature field or the material state via the phasefield is crucial for the accurate computation of the residual stresses. In this case, the dependency on the temperature can experimentally be figured out through conducting special measurement techniques as presented in, e.g, [85,93]. Alternatively, to capture the variation of the parameter's values during the welding process, we propose in this work an alternative and easier-to-implement approach, in which the phase-field variable is utilized to generate state-dependent material properties. In particular, depending on the state of the material, i.e. soft or hard, a state-dependent value is assigned to each parameter (one value for the soft state and one value for the hard state). In this case, a smooth transition between these two values is achieved via the utilization of a phase-field-dependent interpolation function. Table 2 summarizes the proposed expressions of the material properties for the nickel-based alloy IN718 with their soft/hard-staterelated values.
To show the effect of the gradient-energy coefficient l on the resulted parameter curves given in Table 2, a transient one-directional melting analysis is performed. In this, the right-end of the rectangular FE domain is kept at a constant temperature of 1700 K above the melting point, while the left-end is kept at 293 K as illustrated in Fig. 1a.
In this way, the melting interface moves from the right to the left-hand side of the considered domain. Figure 1a shows the contour plots of the temperature distribution and the phase-field variable at time equals 50 s. The spatial variation of the phase-field-dependent properties with change of the gradient-energy coefficient is given in Fig. 1b. On the other hand, Fig. 2 shows the counterpart contours and curves when the temperature-dependent properties are used [31].

Finite element formulation
The following discussion focuses on the numerical solution of IBVPs related to residual stresses in gas tungsten arc welding using the finite element method (FEM). The primary variables in connection with the applied phase-field thermo-elastoplasticity modeling approach are {θ, φ, u} and the governing PDEs are summarized in Table 1. Let Ω ⊂ R 3 be the homogenized spatial domain under consideration with ∂Ω being its boundary, ∂Ω is split for each variable into two parts, i.e. Dirichlet or essential (Γ D ) and Neumann or natural (Γ N ) boundaries. These boundaries fulfill the conditions ∂Ω = Γ D ∪ Γ N and Γ D ∩ Γ N = ∅. More specifically, the Dirichlet and Neumann boundary conditions for each variable together with the initial values at t = t 0 can be specified as follows: In this, n is the outward unit normal vector to ∂Ω andq is a prescribed external heat flux.  Fig. 2 Temperature and phase-field contours with l = 10 (mJ/mm) 1/2 (a) and temperature-dependent properties adapted from [31] (b)

Governing weak formulations
In connection with the implementation in the open-access/ open-source FE package FEniCS, one needs to derive the weak or variational form of the governing PDEs, given in strong form in Table 1. The required trial spaces V θ , V φ , and V u for the primary variables and the related weighting functions can be defined as with H 1 (Ω) denoting the Sobolev space of order one. For the numerical examples of welding process in Sect. 5, the energy balance equation will be simplified by ignoring the terms D red loc and H ep from the energy balance equation in Table 1. This is justified as the amount of heat produced by the mechanical deformations is too small in comparison with the heat provided by the welding heat source, see, e.g., [25,89,92] for analogous discussion. In this regard, investigations performed in [77] for entropic thermoelasticity showed also that the structural elastoplastic heating term has a negligible contribution. Following this, the variational problem after multiplying the governing equations in Table 1 by the test functions and applying integration by parts over the entire domain can be posed as:

Temporal and spatial discretization
In the numerical solution of time-dependent IBVPs in FEn-iCS, a time-stepping scheme based on the finite difference approximation over a time interval Δ t = t n+1 −t n is applied. This results in time-discrete (stationary) equations, where the time stepping is followed by constructing the variational formulation towards the spatial FE discretization, see, e.g., [68] for more details. In this work, an "isothermal split"-like strategy is applied in accordance to the descriptions presented in, e.g., [9,72,73], which results in a robust solution of the volumetrically-coupled problem. In particular, at each time step t n+1 , the solution of the coupled problem is carried out in two sequential steps: (i) The thermal and the phase-field problem are solved monolithically under a constant mechanical state to get θ n+1 and φ n+1 . (ii) The mechanical subsystem of elastoplasticity is solved under an isothermal state, which allows computing the updated displacement u n+1 . For the spatial-discretization using the FEM, the continuous homogenized domain Ω is transformed into a discrete domain Ω h subdivided into N e finite elements. For an abstract representation, the primary unknowns of the problem are summarized in a vector v := v(x, t) with v = [ θ, φ, u ] . Thus, the trial and test functions read In this,v h are the approximated Dirichlet boundary conditions and N v denotes the total number of FE nodes. Additionally, N v(i) and M v(i) denote the global basis functions of the trial and test functions at node i, whereas v (i) are the nodal degrees of freedom and δv (i) represent the nodal values of the test functions. Moreover, V h andV h are the discrete, finite-dimensional trial and test spaces, respectively. In this work, we follow the Continuous Galerkin (CG) approach with first-order basis functions, which provides standard linear Lagrange elements, e.g., a triangle with nodes at the three vertices of each element [68]. For the FE treatment of the last term in Eq. (29) 3 , i.e. div(l 2 grad φ), we carried out numerical studies employing second-order basis functions to check the importance of this term. The results showed that this term has a negligible value in comparison to the contribution of other terms in H pc . Thus, it is neglected in the later numerical treatments. Furthermore, the application of the Bubnov-Galerkin procedure, in which v and δv are approximated using the same basis functions N v(i) ≡ M v(i) , implies that V h andV h coincide [74]. Moreover, for the implementation in FEniCS, a common space V h is used for both test and trial function asv h are not specified as part of this function space [68]. For further details on the finite element formulations, the interested reader is referred to, e.g., [10,38,43,47,52,53,68], among others.

Integration algorithms for plasticity
For the numerical solution of IBVPs of rate-independent plasticity with linear isotropic hardening, proper integration algorithms have to be specified to integrate the evolution equations of the plastic strain and the related hardening parameters. Adopting the FE package FEniCS in this contribution, a special focus is laid in the following on the implementation and comparison between two plasticity integration algorithms, which are the "return mapping algorithm" and an "incremental plasticity" approach.

The return mapping method (RM)
Proceeding with a known material state {θ, φ} n+1 at t n+1 and {ε, σ , ε p , α} n at t n , the RM applies a strain increment Δε and aims to figure out the updated current state at t n+1 as summarized in Algorithm (I).

Algorithm (I): Solution procedure using the return mapping method (RM)
Step 1: while t n+1 < t end do Step 2: Apply Δε and compute trial state: Step 3: If χ trial n+1 < 0 then −→ elastic loading step Compute the tangent modulus (pure elastic) Continue to step (1) Step 4: If χ trial n+1 ≥ 0 then −→ plastic loading step Solve nonlinear problem and find Δλ with j = 1 Update quantities: Compute the consistent algorithmic tangent Assemble the global Jacobian and residual Check global convergence: If Res < Tol then (·) n+1 −→ (·) n else j −→ j + 1 iterate within step (4) Increment time t n+1 + = Δt and continue to step (1) In this, the updated values have to satisfy the yielding condition and the global equilibrium within a predefined tolerance. For this purpose, an iterative procedure is used to drive the residual (Res) between the internal forces and the external force vector. Furthermore, in the FEniCS implementation of plasticity, the solution of the internal variables is carried out at the quadrature points, which are interpreted as degrees of freedom when using the so-called "quadrature element". This element type is developed within the FEniCS framework to avoid a suboptimal convergence rate for Newton's method [17,68]. The suboptimal convergence arises due to the linearization procedure followed by FEn-iCS Form Compiler (FFC) when nonlinear continuous weak forms are linearized using Newton's method. Thus, the implementation of RM in FEniCS requires the definition of quadrature elements with the application of Newton's method to compute the stress and its linearization at the quadrature points.

The incremental method (IM)
The incremental method (IM) is presented here as an alternative procedure to solve the plasticity problem. In the IM, the plasticity problem can be solved by avoiding the quadrature element definition and the nonlinear problem in step (4) given in Algorithm (I). The method simply uses the stress state σ n from the previous step t n to calculate the plastic multiplieṙ λ n+1 and update the internal state variables. Then updated σ n+1 and ε p n+1 at the current time (t n+1 ) can be computed as where (•) n denotes the value computed at the previous time step. The method provides an accurate numerical solution for small time increments [1]. The solution steps for the IM are presented in Algorithm (II).

2D modeling of welding-induced residual stresses
In this section, the numerical simulations are carried out for an IBVP problem of GTAW spot welding in 2D plane strain settings. The focus in the phase-field thermo-elastoplastic model is laid on the performances of the two different algorithms for plasticity, i.e. the return mapping method (RM) and the incremental method (IM) previously discussed in Sect. 4.3. The phase-field-dependent material properties, depicted in Fig. 1, are employed in the analysis. Figure 3 illustrates the geometry of the symmetric 2D problem together with the mechanical and heat flux boundary conditions. The left symmetry surface is thermally insulated, while the bottom and right surfaces are subjected to heat losses by convection according tō with θ amb being the surrounding ambient temperature 293 K and h c is the convective heat transfer coefficient of the surrounding air. The heat flux at the top surface represents the heat exchange between the top surface and the GTAW Melting temperatureθ M 1514 (K) Latent heat of fusion L 219 × 10 9 (mJ · tonne −1 ) torch, i.e.
where q max = 81280 mW being the welding power and r 0 = 1.4 mm is the Gaussian radius of the heat source. Additionally, a constant value of the material density ρ = 7737 × 10 −12 tonne·mm −3 is considered, and the value of convective heat transfer coefficient is chosen to be h c = 700 mW·mm −2 K −1 . For the 2D analysis case, the gradientenergy coefficient is set to 3 (mJ/mm) 1/2 and the rest of the parameters used for solving the phase-field variable Eq. (38) 2 are listed in Table 3. Thermal analysis results are presented in Fig. 4, which depicts a contour plot of the phase-field variable after 5 s heating time with the melted weld pool as a soft region φ = 1 and the unmelted hard region with φ = 0.
At each time step, the computed values of the temperature and phase-field variables are passed to the mechanical problem for the computation of the residual stresses. The distribution of the residual stress components in x 1 and x 2directions are shown in Fig. 5 at the end of the analysis, i.e. at t = 50 s. For quantitative comparisons of the results obtained by RM and IM, Fig. 6 shows curves of σ 11 and σ 22 along the sections indicated in Fig. 5. It can be seen that a good agreement is obtained between the solutions by RM and IM. This indicates that both methods are reliable to be used for performing the mechanical analysis. However, they differ in the convergence behavior and computational costs, as will be discussed in the following.
Within the model sensitivity study, the implementation is tested with different time steps and mesh refinements. The results of the sensitivity study are summarized in Tables 4, 5 and 6. The computations are performed using a workstation of an Intel Core™ i7-8700 CPU@3.20 GHz ×12 and 16 GB RAM. The results presented in these tables show that both RM and IM have good stability and convergence behavior. However, IM attains convergence for a wider range of time steps and mesh refinements over RM, whereas the total time needed for the IM computation for a given time-step size is longer. Due to its stability and convergence characteristics for a relatively large time-step size, the IM will be used in the subsequent section for the simulation of residual stresses within a 3D IBVP.

Problem setup
In the following, a 3D problem of the GTAW process will be studied. It represents bead-on-plate welding as illustrated in Fig. 7 (left). The boundary conditions, the geometry, and the material parameters are adopted from [31] for the case study on an autogenous welding process of a plate made from Nickel-based alloy IN718.
For the considered bead-on-plate welding IBVP and due to the symmetry of the surface heat source (illustrated in Fig.  8) along the welding central-line, only half of the domain is considered in the FE model. In this way, the computational time is reduced and, thus, the employed domain has the dimensions of 50×200×2mm 3 , as shown in Fig. 7 (left). The unstructured FE mesh, which consists of four-node tetrahedral elements, is created by the mesh generator Gmsh [37] and illustrated in Fig. 7 (right). The mesh density is higher towards the welding zone and gradually decreases far away from this zone. The side length of the smallest elements utilized in the FE model is (0.6mm), while the largest allowed element side has a size of 5mm. The FE mesh consists of 107,151 elements and 27,147 nodes. The same finite element mesh is used for the thermal and mechanical analyses choosing first-order Lagrange elements. Moreover, thermal and mechanical boundary conditions are defined on the considered domain. In the thermal analysis, the moving heat source model is applied as a surface heat flux at the top surface of the domain. For the plate with a small thickness of 2 mm, a 2D heat source with circular Gaussian distribution is a suitable choice for giving a good approximation of the heating process with few parameters [31,92]. As indicated in Fig. 8, the maximum power density mW/mm 2 will be experienced at the center of the top surface of the plate and continuously decays as the distance increases. The surface heat flux is given as a Neumann boundary condition at the top surface of the considered FE model [33,39,106], which can be written as In this, Q, η s , and r 0 , represent the input power provided by the moving welding torch, the arc welding efficiency, and the radius of the flux distribution, respectively. x 01 and x 02 refer to the initial position of the heat source and v and t are Initial position x 01 , x 02 , x 03 (0, 10, 2) (mm)  Table 7.
On the top of the plate, the heat source travels along the x 2direction for 180 mm, where the start and the end locations are 10 mm distance from the edges. The heat losses due to the interaction with the surroundings are considered to occur through heat convection and modeled by Newton's law. All external boundary surfaces, except for the symmetry plane, are assumed to interact with the surroundings by applying the Neumann boundary conditions given in Eq. (41). To prevent rigid body motion, artificial boundary constraints are applied during the mechanical analysis, as depicted in Fig. 7 (left). This includes a mechanical symmetry condition on the plane of symmetry (u 1 = ∂u 2 /∂ x 1 = ∂u 3 /∂ x 1 = 0 mm), constrained degree of freedoms (u 2 = u 3 = 0 mm) at point D, and (u 3 = 0 mm) at point E. Furthermore, at the beginning of computation, the primary variables are initialized with the following values (φ init = 0, θ init = 293 K and u 0 = 0 mm). The total simulation time is 300 s, which includes the heating time of 113 s. The time steps employed in the thermal and mechanical analysis are Δt = 0.08 s and Δt = 0.16 s, respectively.

Discussion of the thermal analysis results
The results of the 3D welding problem related to the temperature field and phase change are presented in Figs. 9, 10, 11 and 12. Figure 9 shows temperature profiles assembled  at three different locations corresponding to points A, B and C, which are, respectively, at a distance of 6 mm, 8 mm, and 10 mm from the weld centerline. As shown in Fig. 9a, when temperature-dependent properties are used, the predictions of the thermal analysis have a good agreement with experimental measurements conducted in [31] with some deviations related to the peak temperatures. The deviations of the predicted peak temperatures are expressed by the relative errors 1 of ERR θ,A ≈ 6%, ERR θ,B ≈ 7% and ERR θ,C ≈ 7% at points A, B and C, respectively. The corresponding temperature profiles at these points for the analysis based on phase-field-dependent properties are shown in Fig. 9b. In comparison with the experimental ones, higher peak temperature relative errors of ERR θ,A ≈ 19%, ERR θ,B ≈ 13%, and ERR θ,C ≈ 10% are observed. The increase in deviations is attributed to the lowered thermal conductivity at the hard state of the material, which in turn affects the flow of heat from the soft region (the weld pool) to the surrounding regions. Bhatti et al. [16] have also reported that more heat accumulation occurs at lower thermal conductivity values and thus it takes a longer time for the specimen to conduct heat to the surrounding material. As shown in Fig. 10, this effect can be mitigated if we use an averaged value for thermal conductivity of the hard phase, i.e., κ H = (κ H + κ M )/2. Reduced relative errors ERR θ,A ≈ 10%, ERR θ,B ≈ 8%, and ERR θ,C ≈ 8% can be achieved for the profiles at the specified points. Thus, the averaged value for thermal conductivity will be considered in the forthcoming thermal and mechanical analysis.
When it comes to the thermomechanical modeling of welding processes, several research works ignore the effect of the phase change term, H pc , occurring in the energy equation. Investigations of this term exhibits a sound effect on both size and shape of the weld pool, where a wider and shorter weld pool is observed if this term is neglected as shown in Fig. 11. Thus, by utilizing the phase-field approach, it is possible to account for the effect of energy release or absorption in the shape of latent heat during the phase change. Furthermore, this term has a prominent effect on the temperature evolution in the weld pool, which was also observed in, e.g., [56,97]. The impact on the temperature profile is shown in Fig. 12 for the point (0, 100, 2)mm within the weld pool.  Fig. 9 Temperature profiles of points A, B, and C compared to the experiment in [31] for temperature-dependent properties (a) and phase-field-dependent properties (b) 25 [31] for temperature-dependent properties using an averaged value of thermal conductivity

Discussion of the mechanical analysis results
The results of the 3D mechanical analysis are presented in Figs. 13 and 14. In this, the final states of the residual stress distributions for both transverse σ 11 and longitudinal σ 22 stress components are presented in Fig. 13 on a plane at x 3 = 1 mm. To closely investigate the stress state, Fig. 14 shows line plots for σ 11 and σ 22 at the mid-length and mid-depth of the FE model along the path P1 indicated Fig. 7 (left). The residual stress components are compared with the results in [31] as shown in Fig. 14a for the case of a small deformation assumption. From a qualitative point of view, the resulted stresses using the temperature-dependent properties are compared to the benchmark results by calculating the Weighted Mean Absolute Percentage Error (WMAPE). The definition of WMAPE can be expressed as [75] WMAPE : where A t and F t being the actual (reference) values and the predicted values by the proposed model. The transverse stress deviates by 16% while the longitudinal stress shows 17% deviation. Moreover, Fig. 14b Fig. 12 Effect of the phase-change term on the calculated peak temperatures. Neglecting this term leads to a relative error of ERR θ ≈ 6% ison of the results obtained from the model and benchmark temperature-dependent analyses with the ones from the phase-field-dependent properties analysis. It is shown that, in the region of the weld bead, the longitudinal stress σ 22 starts decreasing earlier than in the case of the temperaturedependent analysis. The WMAPE deviation error is 24% for the transverse stress and 8% in the case of the longitudinal stress component. From both Fig. 14a and b, one observes that the region close to the weld centerline experiences tensile stress of the yield stress magnitude, which gradually decreases and turns to compressive stress as the distance from the weld centerline increases. However, in this region, it is noticed that the predicted longitudinal stress is lower than the one predicted in [31]. This can be attributed to the treatment of the plastic strain evolution in our model. We involve the annealing effect in the proposed model. The annealing effect accounts for the hardening reduction due to the change in the dislocation structure caused by the solidliquid phase change [65]. This is modelled such that for any material point going through melting/re-melting effect, the artificial accumulated plastic strains are eliminated [25,30]. From a quantitative point of view, the residual stresses predictions by the numerical model are in good agreement with the small deformation analysis given in [31].

Conclusions and future works
The presented work proposes a modeling methodology for the phase-field thermo-elastoplasticity problem of gas tungsten arc welding in a thermodynamically consistent framework. The involvement of the phase-field method promotes the use of phase-field-dependent material properties for the melted (soft) and the unmelted (hard) material states. In this way, the thermomechanical properties are driven by the melting/solidification state with the phase-field variable. Two-and three-dimensional IBVP numerical examples are introduced to investigate the performance of the proposed model. In the two-dimensional example, two approaches (algorithms) are tested for the computation of the residual stresses using the von Mises isotropic hardening plasticity. The first approach relies on the conventional return-mapping method by solving a nonlinear problem using Newton's method to correct the stress state. The second approach applies simplified "incremental plasticity", in which the stress from the previous step is used to compute the plastic multiplier. The numerical models are implemented in the open-source finite element package FEniCS, where for the return-mapping algorithm the quadrature element type developed with the FEniCS framework is used. The capability of the proposed model in capturing realistic scenarios of the gas tungsten arc welding process was investigated and a three-dimensional bead-on-plate analysis was performed. The numerical results of temperature profiles and residual stresses have good agreement with results available in the literature. The model shows that using phasefield-dependent parameters can provide reliable predictions for both thermal and mechanical analyses. Furthermore, as a matter of implementation in FEniCS, employing the described incremental plasticity approach gives accurate numerical results when compared to those obtained from the return-mapping method.
To this end, the proposed numerical treatment and modeling framework can serve as a base for future research works in related fields to GTAW and the prediction of residual stresses. As future works, the approach can be extended to consider thermo-elastoplasticity within finite deformations, which allows for more accurate description of the coupled processes if large deformations occur. One of the important extensions would also be the consideration of the material in a melted state as a fluid, which can be treated within a unified kinematics framework together with the solid phase. Another future aspect in connection with the proposed phasefield thermo-elastoplasticity approach is the application to predict residual stresses induced by additive layer manufacturing processes of metal parts. Moreover, modeling the weld-induced cracks can also be achieved by incorporating the phase-field fracture approach.  Fig. 7 (left). Solid lines represent the results obtained from the 3D mechanical analysis and dashed lines for small-deformation results in [31], wherein case a temperature-dependent properties are used and case b phase-field-dependent properties are used   Spatial distance φ unmelted metal. In particular, putting the material parameters as functions of the phase-field variable allows for a smooth variation of their values between the two metal states and across the diffusive interface. For illustration, Fig.  15 illustrates the concept of a diffusive interface in phasechange materials in comparison with the sharp-interface concept. The first step in the PFM approach is the definition of the phenomenological phase-field variable that governs the smooth change across the interface. This is defined as : melted (soft) state , φ = 0 : unmelted (hard) state , 0 <φ < 1 : diffusive interface .
The treatment in the following assumes that the phasechange process occurs merely due to changing of the thermal energy and ignores the contribution of the mechanical energy in this process. Thus, in analogy to the thermodynamically consistent approach in [18] and following the Ginzburg-Landau theory, the global potential energy function F over the spatial domain Ω is expressed as the sum of the bulk free energy term f (φ, θ ) and the interfacial energy terms as with l being a gradient-energy coefficient that influences the width of the diffusive zone. At equilibrium, the variational derivative of (46) results in the following restriction: The phase-field evolution and the resulting migration of the diffusive interface is a time-dependent process. To describe this evolution, the application of the Allen-Cahn reactiondiffusion equation is a common practice. It can be expressed as In this, M represents the interface mobility parameter, which can be expressed according to [18] as M = μ φθ M 6 δ ρL (49) with μ φ ,θ M , and δ as the kinetic coefficient, the material melting temperature and the interface width, respectively. In absence of the mechanical or chemical free energy contributions, the free energy density of the system is defined by where the first term W g represents the energy barrier related to the interface between the soft and hard phases of the metal. Moreover, W , g(φ), and p(φ) are the energy barrier height at the interface, the double-well potential, and the interpolation function, respectively: with σ φ being the surface tension. For melting/solidification phase-change problems, the main thermal force driving for the evolution of the phase-field parameter, Q(θ ), can be expressed as Using Eqs. (50) and (51), the evolution of the phase-field variable, Eq. (48), can be rewritten aṡ More details, references, and applications can be found in [20,69,87,88,90,98] among others.