1 Introduction

One of the most effective ways of joining metals is to use an electric arc. Among the methods of welding metals, one of the best processes in terms of quality is the method of electric welding in shielding gases uses a non-consumable tungsten electrode (GTAW). Vendam et al. [1] and Weman [2] indicated in their papers that GTAW method allows one obtained the high quality welds.

The availability of equipment used in the GTAW method is the reason for its popularity. The most important factor for this method in the heating process is the arc efficiency. Stenbacka [3] analysed values of the arc efficiency coefficient published in year 1955-2011. This analysis showed that the efficiency decreases when the distance between the electrode and the heating material increases. To maintain the high efficiency of this process, the distance should be constant. This can be ensured by automating the process. Therefore, the trend of automating production processes allows the GTAW method to be widely used.

The Gas Tungsten Arc Welding process is of a much higher efficiency coefficient [3] than eg. laser heating process, which Li [4] presented in his work. There are a numerous papers on the GTAW methods and its modifications to improve the quality and productivity of welding or surface modification by remelting. Ciechacki and Szykowny [5] used the GTAW method for surface remelting of cooper ductile cast iron. This method may be used not only for joining but also in the heating of steel elements. Models that describe the classic GTAW technology with remelting can be inappropriate for the modelling of heating without fusion zone using a non-consumable tungsten electrode. In addition, the GTAW method is based on the use of a shielding gas (argon), which increases the area of its application, among others, to the point heating or hardening process without surface de-generation (oxidation, etc.). Due to such a non-classical application of the GTAW method, it is necessary to determine the parameters of the procedure/analytical method, which are helpful in the selection of technological parameters of the process. We assumed that one of the ways to determine these parameters is to carry out numerical calculations / simulations. In this work, we present the results of such calculations and their comparison with the results obtained on the research stand.

Another important problem is to determine the parameters of the boundary conditions that simulate the GTAW electrode and the cooling conditions. There are studies in which the topic of using inverse analysis to determine the parameters of boundary conditions for this process (with remelting) is discussed. There is also a publication in which the author’s determine parameters of simulation with using the artificial intelligence method. Tafarroj and Kolahan [6] used artificial neural networks and the regression model in modelling the heat source parameters for the GTAW process and also compared the precision of both models. Whereas Wróbel and Kulawik [7] used artificial neural networks to determine parameters of a hybrid heat source. However, artificial neural networks require immense training, testing and validation sets [8]. In this example there is not enough experimental data for GTAW methods without remelting.

We think that the examples analysed in the paper cannot be referred to the results of other authors for the GTAW method with remelting. Therefore, the results obtained from numerical simulations were compared with the results of the experiment.

The proper selection of heat source parameters during GTAW process requires verification of arc voltage, arc current and arc efficiency coefficient. It is also necessary to check that the results obtained from the applied method for one set of process parameters will be suitable for example for the large current change. In the literature, relationships between simulation parameters, e.g. arc voltage and arc current, are often observed. The paper compares these literature relationships with the results obtained from experimental measurements. Because the voltage, current, arc efficiency and radius of the heat source are mutually linked, their proper adoption must be the result of certain assumptions from an experiment or literature. The paper focuses on the comparison of different radius of the heat source, arc voltage assumed from the literature and experimental dependencies and one value of the arc efficiency coefficient. In the studies conducted up to now, there are no comprehensive solutions for the conversion of GTAW heating conditions into parameters of a simulation heat source.

2 Details of experimental procedure

In order to verify numerical calculations, experimental studies were carried out. A research stand consisting of four main blocks was built: mechanical part, control system, heat source and measuring system (Figs. 1, and 2).

Fig. 1
figure 1

The experimental stand structure: 1 - operational panel; 2 - programmable logic controller; 3 - power supply; 4, 5 - stepper motor drivers; 6, 7 - stepper motors; 8 - current source; 9 - compressed argonium tank with pressure reduction unit; 10 - welding grip with tungsten electrode; 11 - measuring head of laser pirometer; 12 - measuring transducer/amplifier of laser pirometer; 13 - analysed sample

Fig. 2
figure 2

Scheme of the analysed task: 10 - heat source; 11 - measuring system; 13 - analysed sample

The heat source was the inverter welding machine PIROTEC SIC 301/1 (element 8, Fig. 1). To measure the welding current, a digital clamp meter UNI-T UT203 was used. In addition, the voltage between the current output of the welding machine and the common wire was measured using a UNI-T UT139B digital multimeter. The welder has a high-frequency arc ignition function, which allows one maintain a predetermined distance between the electrode and the surface of the sample - no need to arc the arc by rubbing the electrode against the surface of the sample. The heat supply was started automatically. This was done by the PLC controller bypassing the switch located on the welding gun.

Steel specimens (Fig. 2) in the form of 0.245 m long shafts were considered and placed in grips to allow a given rotation and feed movement. To homogenise the surface of the sample, it was cleaned with 800 grit sandpaper. The geometric centre of the heating area of the sample (A) lies on the surface of the sample, opposite the top of the infusible electrode. This point moves along a helix line on the surface of the sample. Its movement is described by means of axial velocity (Vp) and rotational speed of the sample (ωp). By modifying the speed in ωp, the peripheral speed of point A is changed. The pitch of the helix is dependent on the linear speed Vp (Fig. 2).

The measuring point (M) is located at the other end of the diameter exiting the heating point of the sample. The temperature measurement system consisted of a RAYTEK MI3 pyrometer equipped with a measuring head MI31002MSF1 (element 11, 12, Fig. 1). Basic parameters of the measuring head are presented in Table 1.

Table 1 Basic parameters of the measuring head MI31002MSF1

The control system of the experimental stand was based on the PLC SIEMENS S7-1200 controller (element 2, Fig. 1) equipped with HMI operational panel TP600 BASIC SIEMENS (element 1, Fig. 1). The control system has some features which are helpful during experiments: high resistance to electrical disturbances, easy parametrization by human machine interface, built-in ready to use functions for numerical control of drives, build in analog to digital converter. The PLC controller has high resistance to electrical disturbances, which was very important during conducted experiments. Basic parameters of experimental research:

  • specimen material (rolled steel bar): C45,

  • specimen diameter: 0.012 m,

  • specimen length: 0.245 m,

  • distance of the initial heat source from the face of the sample: 0.04 m,

  • length of the heated area by the heat source: 0.14 m,

  • specimen rotation speed: ωp = 2.12 rot/s (0.08 m/s),

  • axial speed of the sample: Va = 0.001 m/s,

  • sharpening angle of the GTAW electrode: 30°,

  • arc length: g = 0.001 m,

  • focal length: f = 200 mm,

  • the pitch of screw line : p = 0.00047 m,

  • heat source parameters: current 30A and 50A,

  • GTAW shielding gas: argon,

3 Finite element simulation

On the basis of the analysis presented in the author’s earlier paper [9], the method of modelling a moving heat source was selected. Due to the accuracy of calculations, heat sources were taken into account in the form of a moving boundary condition along the perimeter of the analysed steel element. Whereas, the movement along the axis was assumed using the drift velocity. To modelling of thermal phenomena, a differential equation of heat conduction in Euler coordinates (Petrov-Galerkin formulation) with a convection term in the following form was used:

$$ \nabla\cdot\left( \lambda \nabla T \right)-\rho c \frac{\partial T}{\partial t} - \rho c \nabla T \cdot \textbf{V} =0 $$
(1)

The solution of equation (1) has been supplemented with the initial condition and boundary conditions in the form:

  • initial condition

    $$ T(x_{\alpha}, t_{0})=T_{0}(x_{\alpha}) $$
    (2)
  • Neumann boundary condition

    $$ -\lambda\frac{\partial T}{\partial n}(x_{\alpha},t)= - \lambda\nabla T \cdot\textbf{n} = q^{*}(x_{\alpha},t) $$
    (3)
  • Newton-Robin boundary condition

    $$ -\lambda \frac{\partial T}{\partial n}(x_{\alpha},t)=\alpha_{\infty}(x_{\alpha},t)(T(x_{\alpha},t)-T_{\infty}(x_{\alpha},t)) $$
    (4)

To solve the above mathematical model, the finite element method, cubic elements with linear interpolation for each of the directions were used. Finally, Eq. 1 was obtained in the following form [10, 11]:

$$ \begin{array}{@{}rcl@{}} &&\left( \upbeta\left( K_{ij}+B_{ij}^{N}\right)+M_{ij}\right)T_{j}^{s+1}=\\ &&\left( M_{ij}-\left( 1-\upbeta\right)\left( K_{ij}+B_{ij}^{N}\right)\right)T_{j}^{\left( s\right)}+\\ &&+\upbeta\left( B_{ij}^{N} T_{j}^{\infty\left( s+1\right)}-B_{ij}^{Q}q_{j}^{*\left( s+1\right)}+Q_{ij}q_{j}^{V\left( s+1\right)}\right)+\\ &&+\left( 1-\upbeta\right)\left( B_{ij}^{N}T_{j}^{\infty \left( s\right)}-B_{ij}^{Q}q_{j}^{*\left( s\right)}+Q_{ij}q_{j}^{V\left( s\right)}\right) \end{array} $$
(5)

The matrix in the Eq. 5 are defined as follows:

$$ \begin{array}{@{}rcl@{}} &&\displaystyle K_{ij}^{e}=\int\limits_{{\Omega}^{e}}\lambda w_{i},_{\alpha}{\Phi}_{j},_{\alpha} d{\Omega}+{\rho c}\int\limits_{{\Omega}^{e}}w_{i}{\Phi}_{j},_{\alpha} V_{k\alpha}{\Phi}_{k} d{\Omega},\\ &&\displaystyle M_{ij}^{e}=\frac{\rho c}{\Delta t}\int\limits_{{\Omega}^{e}}w_{i} {\Phi}_{j}d{\Omega},\\ &&\displaystyle Q_{ij}^{e}=\int\limits_{{\Omega}^{e}}w_{i}{\Phi}_{j}d{\Omega},\\ &&\displaystyle B_{ij}=\int\limits_{{{\Gamma}^{e}_{Q}}}w_{i}{\Phi}_{j}d{\Gamma}+\int\limits_{{{\Gamma}^{e}_{N}}}\alpha_{\infty} w_{i}{\Phi}_{j}d{\Gamma}= B_{ij}^{Q}+B_{ij}^{N}, \end{array} $$
(6)

The task was solved with Euler’s backward scheme assuming β = 1. The approximation and weighting functions were assumed for each direction according to equations:

$$ \begin{array}{@{}rcl@{}} &&\displaystyle \mathop{w}\limits_{i}(\xi)=\mathop{\Phi}\limits_{i}(\xi)+{\mathop{w}\limits_{i}}^{*}(\xi)\\ &&\displaystyle \mathop{\Phi}\limits_{i}(\xi)=\frac{1}{2}(1+\mathop{\xi}\limits_{i}\xi)\\ &&\displaystyle {\mathop{w}\limits_{i}}^{*}(\xi)=\mathop{\xi}\limits_{i}\mathop{\alpha}\limits_{i}P(\xi)\\ &&\displaystyle \mathop{\alpha}\limits_{i}=coth\left( \mathop{Pe}\limits_{i}|_{x}\right)-\frac{1}{\mathop{Pe}\limits_{i}|_{x}}, \mathop{Pe}\limits_{i}|_{x}=\frac{\rho c}{\lambda}\mathop{h}\limits_{i}v_{x} \end{array} $$
(7)

Finally, a complex weight function is obtained:

$$ \begin{array}{@{}rcl@{}} \mathop{w}\limits_{i}(\xi,\eta,\zeta)&=&(\mathop{\Phi}\limits_{i}(\xi)+{\mathop{w}\limits_{i}}^{*}(\xi)) (\mathop{\Phi}\limits_{i}(\eta)\\ &&+{\mathop{w}\limits_{i}}^{*}(\eta)) (\mathop{\Phi}\limits_{i}(\zeta)+{\mathop{w}\limits_{i}}^{*}(\zeta)) \end{array} $$
(8)

The Neumann boundary condition has been supplemented with the following equation describing its distribution in the source area of operation. Teixeira el al. [12] presented a surface heat source model for modeling welding processes. While, Weman [2] determined the dependence on arc voltage for the GTAW method:

$$ \begin{array}{@{}rcl@{}} &&\displaystyle q=\frac{Q}{2\pi R^{2}}exp\left( -\frac{(x-x_{0})^{2}+(z-z_{0})^{2}}{2R^{2}} \right)\\ &&Q=UI\xi\\ &&U=10+0.04I \end{array} $$
(9)

Equation 9 describes point superficial heating with a Gaussian distribution. In this equation the voltage, current and arc efficiency of the heat source are taking into account. The Fig. 3 shows the power distribution for the three tested diameters of the heat source.

Fig. 3
figure 3

Mathematical model of the distribution of the II type condition depending on the adopted radius and current

The cooling process was realized by a third type of condition, in which the heat transfer coefficient was assumed as for air. Li et al. [13] presented the distribution of the coefficient for the considered temperature range:

$$ \alpha_{air}=\left\lbrace \begin{array}{ll} 0.0668 \times T & T_{0}< T < 773 K \\ 0.231 \times T -82.1 & T \geq 773 K \end{array}\right. $$
(10)

Material properties from Eq. 5 are assumed to be temperature dependent according to Fig. 4.

Fig. 4
figure 4

Material properties adopted in the computations

4 Examples of calculation

Task conditions were defined as the closest to the conditions of the conducted experiment. In the simulation, a shaft with a diameter of 0.012m and a length of 0.245m was analysed. The Figs. 5 and 7 and Table 2 shows the divide of finite element mesh into cubic elements. Based on the assumed values of the heat source radius, width of the heat affected zone and the values of temperature gradients in the area of heat source operation, a concentrated mesh around the perimeter of the elements in segment 2 was adopted (Fig. 5, Table 2). The grid pitch length below 0.0005m was obtained. In the other segments, the distance between the nodes was assumed at the level from 0.008 to 0.01m. The length of the segment preceding the heat source’s impact area has been adopted so that the heat source can interact freely with this area without introducing errors using e.g. the Dirichlet condition or third type condition on the front of the element. Because the heat source is moved around the perimeter with the rate ωp = 0.08m/s, symmetrical arrangement of finite elements and an equal density of elements at the surface were applied. In order to take into account the movement of the heat source in a spiral path, using the heat transfer equation with the convection term, it was assumed that the axial speed is a heating rate of Va = 0.001m/s. This ensured the pitch of the heat source paths p = 0.00047m. By using the heating rate, it was possible to maintain the grid density only in the selected area (Fig. 5, Table 2).

Fig. 5
figure 5

Segmentation of the area due to the finite element size

Table 2 Segmentation of the area due to the finite element size

To ensure the movement of the heat source no larger than its diameter (20 positions around the perimeter of the shaft), the time was discretized in accordance with Table 3. The time step density did not increase the temperature or a change its characteristics. The adopted time step also does not cause temperature fluctuations at the place of measurement by the sensor.

Table 3 Segmentation of the area due to the finite element size

The presented paper focuses on the analysis of the influence of parameters of boundary conditions on the temperature distribution and the correlation of the simulation results with the results of the experiment. Two cases of current intensity (30A, 50A) and voltage for the heat source model were considered. The arc efficiency was assumed at the level of 0.9 [3]. For arc current modelling, both literature relations between arc voltage and arc current (9c) and measurements from experimental research were used.

During GTAW heating analysis, the heat affected zone (HAZ) and fusion zone (FZ) are most often determined. Therefore, there are empirical equations determining the size of these zones depending on: current intensity, arc length and traveling speed [14, 15]. However, from these relationships, it is impossible to deduce the size of the heat source in the boundary condition, especially since there is no fusion zone in the studied case. The radius of the heat source, especially in the considered Gaussian distribution limited to 2R, is smaller than the HAZ together with the FZ. Moreover, the radius cannot be too large, because the area of phase transformations will not occur. Yan et al. [16] showed that the maximum diameter of the direct effect of the GTAW for electrode length distance greater than 2mm does not exceed twice the arc length. For smaller values of the arc length, this proportion increases. Therefore, 3 cases of the heat source radius have been taken into account, i.e. the ratio of arc width to arc length at values: 2:1, 2.5:1, 1.5:1; less than 50% of the HAZ size. Taking into account the combination of the above parameters, 12 calculation cases were considered (Table 4).

Table 4 Summary of parameters of numerical simulations carried out

5 Results and discussion

A significant influence on the temperature distribution has the appropriate discretion of the area at the point of heat source application. The discretization of the boundary condition of the adopted finite element mesh was compared. The theoretical fields of heat source distribution were compared with the numerical distribution. The difference in approximate Gaussian distribution for particular variants of radius did not exceed (in the cross section of the heat source) 4%. The power distribution of the heat source was approximated by the diameter, respectively R = 0.001m - 11 nodes, R = 0.00125m - 14 nodes, R = 0.00075m - 9 nodes. The analytical function (3D space) of the heat source was integrated using the trapezoidal method, the same was done for numerical discretization. The difference between the integral of ideal heat source distribution and numerical distribution was not greater than 0.6%. The Fig. 6 shows that the poor discretization of the power of the heat source cannot be an explanation for the difference between the maximum temperature values achieved.

Fig. 6
figure 6

Discretion of the II type boundary condition

The analysis of the results obtained from numerical simulations was performed in the rod cross-section (Fig. 7). The Fig. 8 shows the values of temperature which were extracted in points located on the diameter connecting measuring point M (Fig. 2) with the centre point of the heat source A (Fig. 2). Based on the analysis of the results along the diameter D, (Fig. 8) it can be concluded that the assumed discretization after the depth of the element does not significantly affect the value of the temperature at the place of measurement during the experiment (D ∈< 0.0115 − 0.012 > m). If the measuring point were moved to a heating location, the higher density mesh should be considered - large derivative values in the area D ∈ < 0 − 0.001 > m. The Fig. 8 shows that the temperature decrease significantly when the steel rod makes half a turn. However, it is still higher than in the axis of the steel rod.

Fig. 7
figure 7

Discretization of the assumed area with the measurement line and temperature values for time step 5000

Fig. 8
figure 8

Temperature distribution along the diameter of the steel rod passing through the maximum temperature point

In the performed experiment and numerical simulation it was assumed that due to GTAW heating only HAZ will be created, whereas there will be no remelting area. There is no experimental data on the size of HAZ during GTAW without remelting. The color of the material during heating in comparison to the color standards indicated that the maximum temperature in the considered cases was not less than 200K below the solidus temperature (this applies to 50A current). This was also confirmed by the achieved remelting for slightly higher power of the heat source.

Analysing the obtained results, it can be seen that for smaller power values (30A current), reducing the radius value increases the maximum temperature value (Figs. 8, and 9). At the places of measurement, the slope of the heating curve also changes significantly in the area of the heat source (Figs. 91011 and 12). The use of the same power for a smaller radius must result in a local temperature increase due to the increased power density. The best matching of the experimental and simulation results was obtained when the voltage value measured during the experiment was taken into account (Figs. 9, and 10). For this case of arc current intensity, the best of the analysed radius is 0.001 m. The use of the empirical formula (6c) to determine the voltage (U = 11.2V,U = 12V) leads to obtaining much higher temperatures than in the case of values from the experimental measurement (U = 9V). This applies to the assumption that we have the same efficiency factor for both cases. The consequence of this is to obtain a higher maximum temperature during heating. This approach, in order to achieve similar convergence as in the case of the measured voltage (9b), arc efficiency coefficient of 0.72 should be used, not 0.9 as it was assumed. The value of this coefficient is within the ranges assumed in the literature [3]. Assuming the values of voltages determined from the empirical formula leads to incorrect conclusions regarding the efficiency of the heating process.

Fig. 9
figure 9

Time-temperature distribution for I = 30A - experiment and numerical simulations

Fig. 10
figure 10

Quantile-Quantile plot for the results presented in Fig. 9 in the temperature rage 500 − 1000K

Fig. 11
figure 11

Time-temperature distribution for I = 50A - experiment and numerical simulations

Fig. 12
figure 12

Quantile-Quantile plot for the results presented in Fig. 11 in the temperature rage 500 − 1200K

6 Conclusions

The results of GTAW heating simulation using a superficial heat source with higher current parameters do not allow for determining the direction of the heat source parameters selection (Figs. 11, and 12). The application of measurement data for voltage leads to better compatibility of the experiment and simulation in the range of heating rate of a steel element (only for 30A). However, the maximum temperature from the experiment differs significantly from the simulation results. Accurate determination of the radius of the heat source is more important for lower arc voltages.

In this paper we successfully matched the heating parameters only for a smaller value of the current intensity. The use of the same parameters for higher power values leads to significant errors, especially in the area of temperature stabilization during heating process. This difference equal about 62K for the considered case (Fig. 11). Figures 9 and 11 also shows significant differences during the cooling process of the sample. However, the authors of the paper analysed only heating process and its parameters. Therefore, Quantile-Quantile plot (Figs. 10, and 12) were performed only for a specific temperature range during heating process. The equations and ranges for the heat source parameters (e.g. arc efficiency: 0.36 − 0.9) given in the literature are so wide that their use in a particular process unfortunately requires additional experimental research.

It should be noted that the data from the experiment was obtained assuming a constant emissivity value for the measurement system. Therefore, the values of temperature should be verified by, for example, a thermocouples system.