A New Nodal-Integration-Based Finite Element Method for the Numerical Simulation of Welding Processes

: This paper aims at introducing a new nodal-integration-based ﬁnite element method for the numerical calculation of residual stresses induced by welding processes. The main advantage of the proposed method is to be based on ﬁrst-order tetrahedral meshes, thus greatly facilitating the meshing of complex geometries using currently available meshing tools. In addition, the formulation of the problem avoids any locking phenomena arising from the plastic incompressibility associated with von Mises plasticity and currently encountered with standard 4-node tetrahedral elements. The numerical results generated by the nodal approach are compared to those obtained with more classical simulations using ﬁnite elements based on mixed displacement–pressure formulations: 8-node Q1P0 hexahedra (linear displacement, constant pressure) and 4-node P1P1 tetrahedra (linear displacement, linear pressure). The comparisons evidence the efﬁciency of the nodal approach for the simulation of complex thermal–elastic–plastic problems.

on the modeling of the various interactions that exist between several physical phenomena and features: heat transfer, metallurgical transformations, and mechanical stresses and strains [2,4] (Figure 1). In the present work, metallurgical transformations are disregarded; the application presented in fine will be about the AISI 316 stainless steel devoid of such transformations. With such a simplification, welding simulations are generally performed in two steps. A heat transfer simulation is first carried out to determine the temperature distributions during the process. Then, these distributions are imposed as a variable thermal loading in a mechanical simulation; complex thermoelastoplastic or elastoviscoplastic constitutive laws are used to predict residual stresses and distortions. Performing simulations in such a way, without any feedback of the mechanical calculation upon the thermal calculation implicitly means neglecting the energy dissipation due to plastic straining, as compared to the power of the welding source.
Solving elastoplastic or viscoplastic problems by the FEM leads to difficulties connected to the plastic incompressibility condition resulting, via the normality of the flow rule, from the insensitivity of the yield criterion (generally that of von Mises) upon the hydrostatic stress. The condition of quasi-incompressibility which results from there-compressibility arising only from elasticity-is enforced at all points where the constitutive equations are solved and plasticity (or viscoplasticity) occurs; that is, at all integration points except those rare ones which remain purely elastic. If the number of integration points is too large compared to the number of nodes (where the degrees of freedom are defined), volumetric locking occurs: The number of degrees of freedom becomes insufficient to satisfy both the equilibrium and quasi-incompressibility equations in the whole structure.
A measure of the potential ability of a given element to deal with such problems was proposed by Nagtegaal et al. [5], who defined the following constraint ratio of a given mesh: Constraint ratio = number of degrees of freedom number of constraints .
The "optimal" constraint ratio, considered to characterize the best possible behavior with respect to volumetric locking is given by Optimal constraint ratio = number of equilibrium equations number of incompressibility equations . (2) This optimum ratio is equal to 2 in 2D and 3 in 3D. If we consider, for instance, a 3D mesh made of first-order tetrahedral elements, the number of tetrahedrons is approximately 5 times larger than the number of nodes. Since there is one integration point per element and 3 degrees of freedom per node, the value of the constraint ratio is 3/5, which is way too small compared to the optimum value.
This means that such an element cannot be used in practice for elastic-plastic (or elastic-viscoplastic) simulations.
Several finite element formulations have been proposed to overcome this difficulty, based on either reduced or selectively reduced integration schemes, or mixed displacement/pressure formulations [6,7].
In practice, formulations based on reduced or selectively reduced integration techniques are easily applied to hexahedral elements only; one example is the so-called B-Bar element [6]. However, in the case of complex geometries, hexahedral elements demand heavy meshing operations-in contrast to tetrahedral elements which can be used for all kinds of geometries, even quite complex, using the robust automatic meshing tools now available.
Mixed formulations [8] can be applied to both hexahedral and tetrahedral elements, using Lagrange multipliers and considering extra degrees of freedom. However, the coexistence of degrees of freedom of different physical nature (displacement and pressure), and for some elements-such as the P1+/P1 element [9,10]-the presence of additional internal degrees of freedom, can lead to a significant increase of computation cost for nonlinear problems. Likewise, high-order elements combined with reduced integration techniques can be used to circumvent the locking problem, but they unavoidably lead to a significant increase of the size of the discretized problem, and consequently of the CPU time.
Choosing nodes as integration points emerges as a natural idea to solve the locking problem, since the constraint ratio then takes its optimal value. This idea has been developed within the context of finite elements during the past 20 years, mainly for 4 node tetrahedral meshes. Bonet et al. [11] first developed a tetrahedral element with average nodal pressure, for applications to dynamic problems with an explicit algorithm. The method was extended by Dohrmann et al. [12], who proposed a new linear tetrahedral element using an averaging procedure applied to all strain components, for applications to small strain elasticity. In their approach, the strain tensor was defined, at each node, as the average of the strain tensors in the various elements containing it (nodal averaging technique). Several authors developed similar approaches for dynamic problems in a large strain with an explicit algorithm [13], and quasistatic elastoplastic problems [14]. Other types of elements were also considered [15,16]. The nodal smoothed finite element method rests upon similar ideas [17,18].
The idea of defining the strain tensor at points, rather than over volumes, is also standard in meshless methods. In this context, numerical integration techniques have been developed to calculate the integrals involved in the weak formulation of the problem directly from a cloud of points. Among these, the stabilized conforming numerical integration (SCNI) technique proposed by Chen et al. [19][20][21] seems of particular efficiency. The approach basically consists of splitting the structure represented by a cloud of points into subvolumes containing one point each, and considering, at each point, a strain tensor defined as the average of the strain tensor over the subvolume containing it. However, application of the SCNI technique is not in fact limited to the context of meshless methods, and its combination with the FEM was recently discussed [22,23].
However, although finite element computations using nodal integration techniques have been considered extensively for purely mechanical problems, their application to thermomechanical problems does not seem to have attracted attention up to now. Such problems are bound to raise new difficulties, since they may involve very localized loadings induced by thermal strain fields.
In this paper, we propose a new finite element approach based on linear tetrahedral meshes, nodal integration, and the SCNI technique, for the numerical solution of 3D thermoelastic-plastic problems. The aim is essentially to assess the robustness of the method and its ability to provide accurate results in welding simulations, as compared to more classical finite element approaches.
The paper is organized as follows: • In Section 2, the nodal-integration-based finite element approach proposed is developed. This presentation is also taken as an occasion to clarify the relations between several classical nodal integration techniques.
• Then, Section 3 is devoted to the application of the method to some welding simulations. The example treated is provided by Task Group 4 of the European Network on Neutron Techniques Standardization for Structural Integrity (NeT). The results are compared to those obtained with more classical finite element analyses involving elements with mixed displacement/pressure formulations: Q1P0 hexahedrons [5,8] (linear displacements and constant pressure in each element) and P1P1 tetrahedrons (linear displacements and pressure in each element) [24].

General Principles
For the sake of simplicity, the approach is developed in the context of the geometrically linearized theory (infinitesimally small displacements and strains). This hypothesis is harmless here, since the focus is on welding simulations, which generally involve small displacements and strains.
Let Ω denote the domain of study and ∂Ω its boundary. This boundary is the disjoint union of a part ∂Ω u where the displacement u is prescribed, and a part ∂Ω t where surface forces t are imposed.
Let us first shortly review the FEM in its classical form. The method permits calculating an approximation of the displacement field u, using the principle of virtual work which reads, dynamic effects being disregarded: In Equation (3), u * is an arbitrary virtual displacement field satisfying u * = 0 on ∂Ω u , b denotes the body force, ε * = 1 2 grad u * + grad T u * the virtual linearized strain associated to u * , and σ the Cauchy stress tensor.
The FEM consists in looking for a subdomain-based nodal approximation. The domain of study Ω is split into disjoint subdomains or finite elements Ω e (meshing step). These elements are connected through nodes which carry the degrees of freedom of the problem-in mechanics, the components of the displacement field. The nodal values of the displacements are interpolated inside the elements using the so-called shape functions, to obtain an approximation of the displacement field over the whole domain of study.
The integrals appearing in Equation (3) are decomposed over the various finite elements and then evaluated through some numerical integration technique-generally that of Gauss; for instance, for the integral Ω ε * : σ dΩ, one writes: where ε * g and σ g denote the virtual strain and stress tensors at Gauss point g, and ω g the weight associated to this point. Equation (4) makes it clear that the constitutive equations of the material must be used at each Gauss point g to calculate the local stress tensor σ g . This calculation is carried out at each discretized instant from the local mechanical state (stresses plus internal variables) at the previous instant, plus the local strain increment tensor ∆ε = 1 2 grad ∆u + grad T ∆u where ∆u denotes the increment of displacement. Note that this procedure makes it compulsory to store all necessary data (strains, stresses, internal variables) at each Gauss point.
Let us now come to the nodal-integration-based FEM. The general principle of nodal integration is different in the sense that integrals are not decomposed over finite elements, but directly evaluated as sums over the various nodes. Thus, for the term Ω ε * : σ dΩ, for instance, one writes: where Ω p denotes the subdomain containing node p, ε * p and σ p the virtual strain and stress tensors at this node, and ω p the volume of the subdomain. The various subdomains Ω p must be disjoint and cover the domain of study Ω. Equation (5) implicitly relies on a suitable definition of the nodal strains, which is discussed next.

Definition of Nodal Domains and Nodal Strains
In the approach adopted, one uses a first-order triangular (in 2D) or tetrahedral (in 3D) mesh of the domain of study. The subdomain associated to a given node p is then defined as the union of some subsubdomains Ω e p , embedded each within a finite element e containing this node. The construction is illustrated in Figure 2. In 2D, we define, following [20], the subsubdomain Ω e p as the union of the two triangles with vertices at the node p, the midpoint of one of the two edges containing it, and the centroid of the triangular element. In 3D, we naturally extend the 2D definition by defining Ω e p as the union of the six tetrahedra with vertices at the node p, the midpoint of one of the three edges containing this node, the centroid of one of the two triangular faces containing this edge, and the centroid of the tetrahedral element. (The subsubdomain Ω e p may be characterized in a more abstract, but quicker way as follows. Let (α, β, γ) in 2D, or (α, β, γ, δ) in 3D, denote the triplet or quadruplet of barycentric coordinates defined within the triangular or tetrahedral element e, with α being the coordinate associated to node p. Then, Ω e p is the subset of the finite element e, which consists of those points having a barycentric coordinate α larger than the other coordinates β, γ or β, γ, δ.) In the context of the SCNI technique, Chen et al. [19,20] propose to define the strain tensor at node p as the average of this tensor over the subdomain Ω p containing this node: Application of Green's formula to the second integral here permits expressing it as a surface integral: where n ≡ (n i ) denotes the outside unit normal vector to ∂Ω p . Equation (7) is advantageous to evaluate the nodal strains in the context of meshless methods, because it permits circumventing the tricky calculation of the derivatives of the displacement field. However, in the context of the FEM, the situation is different, because computation of the strain components within each finite element, from the gradients of the shape functions, is a completely standard operation: Equation (6) is then preferable because it eliminates difficulties tied to the somewhat complex definition (especially in 3D) of the boundaries of the subdomains Ω p .
In addition, the calculation of nodal strains may be simplified using the following remark: within each finite element e, the volumes of the various subdomains Ω e p containing the different vertices p of the element are equal. To establish this property, one may first note that it is trivial for an equilateral triangle or a regular tetrahedron, as a consequence of the "geometrical equivalence" of the three or four vertices in this case. Next, for an arbitrary triangle or tetrahedron, there exists a unique affine transformation F mapping some "reference" equilateral triangle of regular tetrahedron onto that actually considered. Then, with obvious notations, vol(Ω e p ) = detF vol((Ω e p ) ref ). Since detF is constant within element e, equality of the volumes of the reference subdomains (Ω e p ) ref containing the various vertices p of the element implies equality of the volumes of all actual subdomains Ω e p . This remark may be put to use to calculate the nodal strains by noting that as a consequence, within each element e, for every p, vol(Ω e p ) = 1 N vol(Ω e ), where N denotes the number of vertices (3 in 2D, 4 in 3D) of the element, Ω e the domain it defines, and vol(Ω e ) its volume. It then follows from Equation (6) plus the homogeneity of the strain ε e ij within each element e (first-order triangle or tetrahedron) that where S p denotes the set of elements e containing node p. (Note that the factors 1 N in the numerator and the denominator finally cancel out!). This means that to calculate the nodal strains, one need not in fact worry about the subdomains Ω p containing the various nodes p; it suffices to average the strain tensor over the elements containing each node. In other words, the nodal strains calculated by the SCNI technique, using the nodal subdomains defined above, are identical to those resulting from the simple nodal averaging technique used by Dohrmann et al. [12] and many other authors.
One may also think of defining the nodal strains ε p ij so as to ensure identity, in a weak sense, of the strain field ε nod ij defined by these nodal strains via the shape functions of the elements), and the strain field ε FE ij calculated using the standard FEM. This would lead to the following conditions: where N p classically denotes the shape function associated to node p. Since N p vanishes in all elements not containing node p, this yields: where N e p denotes the value of the shape function N p in element e. Lumping the "left-hand side matrix" of components Ω e N e p N e q dΩ (this approximate procedure is classical and justified for such "mass matrices") and using the homogeneity of the strain ε e ij within each element e, one gets from there: Noting that Ω e N e p dΩ = 1 N vol(Ω e ), (the argument leading to this conclusion is the same as that used above to establish the property vol(Ω e p ) = 1 N vol(Ω e ) for every node p in element e) one finally gets: Thus, again, the nodal strains calculated in that way turn out to be identical to those resulting from the nodal averaging technique of Dohrmann et al. [12].

Definition of Nodal Thermal Strains
To get the thermal strains at the various nodes, the first, simplest option is to directly calculate them from the nodal temperatures resulting from the previous thermal analysis.
However, this option is not logically consistent with the SCNI approach, in which the nodal (total) strain at node p is calculated on average over the nodal subdomain Ω p containing this node (Equation (6)). Therefore a second, more logical option consists in similarly calculating the thermal strain at node p on average over Ω p : A third option, consistent with the nodal averaging technique of Dohrmann et al. [12], consists in calculating the thermal strain at a given node as the average of the thermal strain in the elements e containing this node. This option is similar to the second one, except that for any node p and any element e containing it, in the averaging procedure at this node, the subsubdomain Ω e p is replaced by the whole element e, which is 3 (in 2D) or 4 (in 3D) times larger; this should result in a slightly smoother thermal strain field.

Calculation of Stresses and Internal Variables
Once the nodal (total) strains and thermal strains are known, the latter may be subtracted from the former to get the "mechanical part" of the nodal strains. The stresses and internal variables may then be calculated at every node, in exactly the same way as they are calculated at every Gauss point in the classical FEM; the calculation (classically termed the correction of the elastic predictor) is based on the constitutive equations of the material: elasticity law, plasticity criterion, (visco)plastic flow rule, evolution equations of hardening parameters. Any elastic-plastic or elastic-viscoplastic behavior available in the computer code may be considered. In the application considered in Section 3 of this paper, however, an elastic-plastic model based on the von Mises criterion, the associated plastic flow rule, and an isotropic type of hardening are considered.

Benefits and Drawbacks of the Nodal Approach
The ensemble of methods proposed in the former sections will simply be termed collectively the nodal approach in the sequel. Its advantages are as follows: 1.
The restriction to linear triangular (in 2D) or tetrahedral (in 3D) meshes is compatible with the use of the standard automatic meshing tools available nowadays; 2.
Volumetric locking, for incompressible or nearly incompressible media, is eliminated by the very principle of the method, which places the integration points at the nodes carrying the degrees of freedom; 3.
Transfer of quantities from one mesh to another (if remeshing is required) also becomes easier. Indeed, with the classical FEM, such a transfer requires, just like the plotting of contours, a preliminary extrapolation of these quantities from the Gauss points to the nodes of the first mesh, which again is made unnecessary by the nodal approach; 4.
If the elements become too distorted near one or several node(s), local remeshing may be achieved by redefining the elements without moving the nodes, thus circumventing any transfer of quantities; 5.
The size of the results files is greatly reduced, thus permitting to consider larger meshes. This advantage arises from the storage of physical quantities at the nodes instead of the Gauss points, combined with the fact that at least for big structures, the total number of nodes is generally much lower than the total number of Gauss points. (In the case of first-order tetrahedral meshes, the ratio is of the order of 1/5.) However, the nodal approach is not free of drawbacks: 1.
One shortcoming arises from the very definition of the nodal strains (Equation (8)). The strain tensor at node p depends on the displacements at all nodes belonging to elements containing node p (first-neighbor nodes). This introduces couplings between first-neighbor nodes in the so-called "B-matrix" connecting the strain tensor at a given point to the displacements of nodes, and thus couplings between second-neighbor nodes in the global tangent matrix (Figure 3). In contrast, in the classical FEM, couplings in the tangent matrix arise only between first-neighbor nodes. This must induce an increase of the computation time within the nodal approach. In practice, however, this disadvantage is more than compensated by the solving of constitutive equations at the nodes, which are much less numerous than the Gauss points where these equations are solved in the classical FEM; Figure 3. Couplings between a node and its first-and second-neighbors.

2.
Another drawback arises from the possible presence, in some cases, of "hourglass" modes-that is, ensembles of nodal displacements generating a strain field with zero elastic energy, although not defining a rigid-body motion. To overcome this shortcoming, Puso et al. [25,26] proposed to modify the weak form of the equilibrium equations in the following way: where α is a small, positive "stabilization parameter", and σ elast denotes the stress tensor calculated assuming a purely elastic behavior. In the expression between big parentheses, the first integral is calculated using nodal integration, while Gaussian integration is used for the second.

Application to Some Welding Simulation
The finite element method based on the nodal integration technique discussed above has been implemented into the SYSWELD commercial finite element code. All calculations presented below were performed using this special version of SYSWELD [27].
An application to some welding simulation, involving very strong thermomechanical aspects, is presented. Some comparisons of the results with those obtained with more classical FE approaches are made.

Benchmark Model
Welding simulations have become more and more common in the past 30 years. For a company such as FRAMATOME, the French manufacturer of nuclear plants, such simulations have become a decision-supporting tool used to speed up and improve the development and qualification of welding and repair techniques.
The example considered here is a benchmark proposed by Task Group 4 of the European Network on Neutron Techniques Standardization for Structural Integrity (NeT). The general aim of the NeT is to develop experimental and numerical techniques, and standards for the reliable characterization of residual stresses in welded structures [28]. Task Group 4 deals with the estimation of residual stress fields in a 3-pass slot weld obtained with the tungsten-inert-gas (TIG) process, in an 18 mm thick plate made of AISI 316L austenitic stainless steel [29]. Dimensions are illustrated in Figure 4; the bead length and depth are both somewhat variable, the bead length varying from 74 to 82 mm and the bead depth from 2.1 to 2.7 mm. The welding process parameters are presented in Table 1.  A preliminary remark is in order concerning what is expected from numerical simulations of such a benchmark problem. Of course, the anticipated benefit lies in a validation of the numerical method and results. However, two types of validation should be carefully distinguished in such an instance:

1.
Validation of numerical results versus experiments. What is essentially validated in this way is the quality of the physical model(s) used in the simulation: for instance, for elastoplastic problems, the relevance of the modeling of strain hardening phenomena (isotropic, kinematic or mixed istropic/kinematic model); 2.
Validation of the numerical results through comparison with results obtained with other numerical methods (including, if possible, one acceptable as a reliable "reference"). What is tested in that way is the quality of the numerical solution obtained with the method investigated, the physical model(s) being given and fixed.
What is in question in this paper is the second type of validation. However, this does not mean, of course, that the quality of the theoretical model used to simulate the benchmark problem is unimportant; on this topic, see [29].
The simulations are carried out in two steps. A thermal analysis is first performed; it provides the temperature distributions during the whole process. Then, the calculated temperature distributions are applied as a thermal loading in a thermomechanical analysis, which provides the resulting stresses and strains.
The thermal analysis is carried out using the classical (Gaussian) FEM, using meshes made up of 4-node tetrahedra or 8-node hexahedra. For the thermomechanical calculation, three types of formulation are considered. The first one uses the nodal-integration-based finite element formulation depicted above, with linear tetrahedral meshes. The results are compared to those obtained with the classical FEM, using both mixed displacement/pressure Q1P0 hexahedra (linear displacements, constant pressure within each element) and P1P1 tetrahedra (linear displacements and pressure within each element). The comparison of the three types of simulations is made on the residual stresses obtained after the first pass.

Meshes
The meshes were created using the VISUAL-MESH TM software [30]. Thanks to the symmetry (see the D-D line in Figure 4), only half of the geometry needs to be discretized. The three different meshes considered are presented in Figure 5. The first mesh (A) is intended for element type Q1P0, and the other two (B and C) for the stabilized P1P1 element [24] and the nodal approach; this is summarized in Table 2. Mesh A was generated manually, but both meshes B and C were obtained with an automatic mesh generator. Note that mesh C is more refined than the other two in the weld pass (colored green in Figure 5) and in the heat-affected zone.

Thermal Simulations
Thermal simulations have been performed for the three different meshes, using the parameters provided in Table 1. The initial temperature is 20 • C. A two-offset-double-ellipsoid heat source moving at a constant speed v is introduced within the elements of the first weld pass, in the same way as in the work of Xu [29]. (An element activation/deactivation option is used to simulate the material input in the weld pass.) The volumic density of power q in the front and rear parts of this heat source is given by Equation (15) below, where these regions are denoted by the subscripts 1 and 2, respectively: where f 1 = 50.63, f 2 = 25.31, a = 1.799, b = 1.299, c 1 = 1.5, c 2 = 2. These parameters have been adjusted to find by the simulation, the boundary of the experimental melted zone [29]. The volumic density of power given by Equation (15) is applied at the Gauss points of the finite elements of the first weld pass. The total power really input in the model may therefore depend on the local mesh where the heat source applies. Therefore, the Q value is calculated at each moment in order to ensure that the total power input is equal to 825 W for the geometry considered. On all external surfaces of the mesh except the symmetry plane, an exchange coefficient of 10 Wm −2 ( • C) −1 is applied with an outside temperature of 20 • C. Metallurgical transformations are not included in the simulations, since the 316 L stainless steel has only one (austenitic) phase.
The temperature distributions obtained with the three meshes during the deposit of the first weld pass are very close, as shown in Figure 6.

Mechanical Simulations
A thermoelastoplastic constitutive law with isotropic hardening is used to describe the behavior of the AISI316L stainless steel. Figure 7 illustrates the variations of Young's modulus and the thermal strain versus temperature and provides the hardening curves (stress versus strain) at various typical temperatures. The duration of the welding and cooling simulation is taken as 1000 s, in order for the cooling to be complete. The calculation based on the nodal approach is performed without any stabilization (the coefficient α of Equation (14) is taken to be nil); further, the nodal thermal strains are obtained simply from the nodal temperatures, following the first option discussed in Section 2.3 above.
Let us first compare the CPU times and computer resources needed for each calculation. Such a comparison is provided in Table 3. All the calculations were performed on the same computer, with the same accuracy required for convergence of the equilibrium iterations, and the same number of time steps. The only difference in the methods of solution used was that for the classical FEM with Q1P0 elements and the nodal approach, a BFGS method [7] was used to solve the global nonlinear problem at each time step, whereas for the classical FEM with P1P1 elements, since the BFGS method did not permit reaching convergence at all time steps, a Newton-Raphson procedure with full calculation and inversion of the tangent matrix was used. The CPU time required by the simulation based on the standard FEM with Q1P0 elements is the smallest of all. However, such a hexahedral mesh cannot be obtained from current-day automatic mesh generators. It was not difficult to create it manually for the very simple structure considered here, but doing the job for a complex welded structure would be arduous and extremely time-consuming. With regard to the other two simulations, the CPU time required by the nodal approach with meshes B and C is shorter than that necessary with the classical FEM with P1P1 elements, in spite of the increased bandwidth of the left-hand-side matrix. The increase of the matrix bandwidth within the nodal approach is compensated on the one hand by an easier convergence of the global iterations, and on the other hand by a much smaller number of integration points where the constitutive equations must be solved.
Concerning the computer resources, with the nodal approach, since all results (stresses, strains, internal variables) are stored at nodes instead of Gauss points, the disk space required is much smaller than with the standard FEM. However, the situation is opposite for the RAM needed, obviously because of the increased matrix bandwidth required by the nodal approach.
Two features may be compared at the end of the welding process: distortions and residual stresses. Distortions are very small in the case considered, because of the small size of the weld bead compared to that of the plate; for instance, in the classical FE calculation based on Q1P0 elements, the maximum vertical displacement in the plate is 0.104 mm, to be compared to the largest dimension of the plate, 150 mm. Therefore, we concentrate on the comparison of residual stresses. This is a crucial aspect in practice, since accurate stress calculations are required for such applications as future cracking predictions, and difficulties may be anticipated due to the very localized nature of the thermal strain fields resulting from welding. Figure 8 displays a comparison of the distributions of the longitudinal residual stress σ yy obtained in the five simulations, after the deposit of the first bead (see Figure 4 for the definition of x, y, z directions). Moreover, Figure 9 provides a more detailed view of the distributions of both the transverse and longitudinal residual stresses along a vertical line. The line chosen, located on the plane of symmetry, is shown in blue at the top of Figure 9.   Figure 8 shows that the results of the five simulations are globally similar. However, a more careful look at the various stress distributions reveals that they are more or less smooth, depending on the simulation. This is confirmed by inspection of Figure 9, which permits better appreciating the differences. The superiority (in the context of the classical FEM based on Gaussian integration) of brick elements, such as Q1P0 elements, has been established in numerous numerical studies, and for various types of problems, especially in the context of welding; on this point, see notably [31][32][33]. Therefore, the solution obtained with Q1P0 elements may be adopted as a reference. The "reference" residual stresses thus calculated by the classical FEM with Q1P0 elements (Figure 8a) are the smoothest of all. The distributions obtained with the standard FEM with P1P1 elements (Figure 8b,c) are slightly less regular, especially that corresponding to the cruder mesh (B), but still quite acceptable. On the other hand, the nodal approach generates somewhat more chaotic stress distributions (Figure 8d,e), especially when the cruder mesh (B) is employed-note also that the red curves in Figure 9, corresponding to this case, accordingly present the largest oscillations of all. This is somewhat of a problem.
Several remedies may be tested to try and resolve this specific issue of the nodal approach: 1. The first consists of ascribing a nonzero (though small) value to the parameter α of Equation (14), since it is supposed to have a "stabilizing" virtue [25,26]. Unfortunately, numerical experience shows that this is insufficient to erase stress oscillations; 2.
Since in the simulations just presented, the nodal thermal strains were simply evaluated from the nodal temperatures, another possible remedy consists of using the other two options presented in Section 2.3 for the calculation of these thermal strains. Figure 10 compares the distributions of the residual mean stress (pressure) obtained with the three different options, on the cruder mesh B. The comparison clearly shows that the stress distribution obtained is slightly smoother when the thermal strain at a given node is obtained from some average of the temperature over all the subvolumes (option 2) or elements (option 3) containing this node, instead of simply the nodal temperature there (option 1). Option 2 (using an average of the temperature over the subvolumes containing the node considered), which warrants a smoother stress distribution while being logically consistent with the SCNI approach of Chen et al. [19][20][21], is therefore recommended; 3. Finally, one may note that the stress distributions provided by the standard FEM are basically obtained from average values of the stress components in the elements (before these values are transported to the nodes for the plot), rather than from truly discrete values at Gauss points; this means that with the classical FEM, more confidence is put in volume values than point values.
(It is well-known, for instance, that for quadratic subintegrated elements, stress values at Gauss points are much less reliable than average values in elements.) Therefore, it seems logical for consistency reasons, when using the nodal approach, to first average the stresses in the elements before plotting their distribution in the same way as with the standard FEM. (The suggestion was first made by Puso et al. [25].) Figure 11 compares the distribution of the mean stress (pressure) obtained in this way to the reference distribution obtained with the classical FEM with Q1P0 elements. The slight smoothing resulting from the suggested averaging of stresses over elements may be seen to lead to some stress distribution of quality comparable to that of the reference.

Conclusions
A nodal-integration-based finite element method was developed for the numerical solution of thermomechanical problems and implemented into the SYSWELD TM software [27]. The approach uses linear triangular (in 2D) or tetrahadral (in 3D) meshes, thus permitting the use of standard present-day tools for automatic meshing. As a consequence of the placement of integration points at the nodes instead of the Gauss points of the elements, it avoids volumetric locking phenomena currently encountered with the standard finite element method, for elastoplastic or viscoplastic constitutive equations involving plastic incompressibility.
Nodal subdomains, as required by Chen's et al. [19][20][21] SCNI approach to nodal integration, were defined for both triangular and tetrahedral meshes, as a natural byproduct of the discretization of the domain of study into finite elements. It was shown that with this definition of nodal subdomains, the nodal strains evaluated through the SCNI approach are identical to those calculated using Dohrmann's et al. [12] nodal averaging technique. A natural weak formulation of the identity of the strain fields in the nodal approach and the standard FEM was also shown to ultimately lead to the same definition of nodal strains. In addition, three different methods were proposed for the calculation of the nodal thermal strains required by thermomechanical analyses.
Several calculations of residual stresses induced by some welding process were performed, using both the standard FEM and the nodal approach. The results obtained with the nodal approach are globally consistent with those obtained using the classical FEM with mixed displacement/pressure elements (Q1P0 hexahedra and P1P1 tetrahedra). With the nodal approach, some spatial fluctuations of the residual stresses may be obtained when using a coarse mesh. These oscillations may be reduced and virtually eliminated by refining the mesh and/or refining the method of calculation of nodal thermal strains and/or using a "smoothing procedure" of the stress components, similar to that used within the classical FEM when evaluating the average stress components in the elements from their values at Gauss points.
The possible use of linear triangular or tetrahedral meshes, as permitted by the nodal approach without the need for precautions for the prevention of locking phenomena, is a major advantage in the context of numerical simulations of the mechanical consequences induced by welding processes involving complex geometries and/or welding paths. In spite of the increased bandwidth of the left-hand side matrix, the CPU time required by the nodal approach is less than that necessary with the classical FEM with P1P1 elements. Additionally, convergence of the global "equilibrium iterations" is generally easier, which permits the use of more economical quasi-Newton methods, instead of a full Newton-Raphson method requiring calculation and inversion of a new tangent matrix at each iteration. Finally, with the nodal approach, time-integration of the constitutive equations is performed at the nodes instead of the Gauss points; since the ratio of the total number of nodes over the total number of Gauss points is generally considerably smaller than unity (of the order of 1/5 for a typical large tetrahedral mesh), the disk space necessary to store the results (stresses, strains, internal variables) at each time-step is appreciably reduced, which permits envisaging the simulation of larger welded structures.

Conflicts of Interest:
The authors declare no conflict of interest.