Numerical simulation of high-frequency induction welding in longitudinal welded tubes

In the present paper the high-frequency induction welding process is studied numerically. The mathematical model comprises a harmonic vector potential formulation of the Maxwell equations and a quasi-static, convection dominated heat equation coupled through the Joule heat term and nonlinear constitutive relations. Its main novelties are a new analytic approach which permits to compute a spatially varying feed velocity depending on the angle of the Vee-opening and additional spring-back eﬀects. Moreover, a numerical stabilization approach for the ﬁnite element discretization allows to consider realistic weld-line speeds and thus a fairly comprehensive three-dimensional simulation of the tube welding process.


Introduction
High-frequency induction welding is widely used, especially in the production of superior quality oil and gas pipes and structural tubes.A steel strip is cold-formed into a tubular shape in a continuous roll forming mill.The strip edges are electromagnetically heated and joined mechanically by pushing the strip edges against each other to form the longitudinally welded tube (Fig. 1).
The welded joint as seen in the transverse cross-section of a welded tube (Fig. 2), is a very narrow zone compared to the tube diameter.The strip edges are heated to almost melting temperature and are pushed against each other in the viscoplastic state to form the welded joint where crystallographic texture and microstructural changes appear.It can be divided into three zones namely base metal, thermomechanically affected zone (TMAZ) and heat affected zone (HAZ).The transition zone between the base metal and HAZ is the TMAZ.In the HAZ, solid-solid phase transitions occur because of rapid temperature changes.These microstructural changes lead to change of mechanical properties in the HAZ.The characteristics of the HAZ depend heavily on how the tube and pipe manufacturer set up the welding line and how the welding machine is being used.The electromagnetic heating of the tube is analogous to transformer theory.The coil is the primary current source, the strip is where the current is induced and the impeder acts as a magnetic core.The entire setup, coil current and frequency determine the amount of induced current in the tube.
High-frequency alternating current is supplied to the induction coil.This induces eddy currents in the strip under the coil.The induced current can follow the principal paths indicated in Fig. 3 to complete the circuit.Along the strip edges, it can flow downstream from the coil towards the welding point or away from the coil in the upstream direction.At any strip cross-section, the current can follow a path either along the outer circumference or the inner circumference.The goal of high-frequency induction welding is to maximize the current density in the strip edge downstream, towards the welding point.
The impeder has a high relative magnetic permeability and guides the magnetic flux inside the formed strip.Thereby, the amount of current in the inner circumference of the strip to complete the circuit is reduced.Hence, the relative positioning of the strip, induction coil and impeder is very important to obtain an efficient heating process.Another factor which can affect the current concentration in the strip edges is the proximity of the strip edges to each other.The geometric shape of the opening between the strip edges is usually a Vee shape.Sometimes the Vee shape is distorted by spring-back due to the mechanical forming of the strip.This also affects the current distribution.
Further important process parameters are the coil current and welder frequency.The coil current is directly related to the power needed to obtain the desired temperature at the welding point.This would be in accordance with the required welding-line feed velocity for a given combination of tube diameter and wall thickness.Welder frequency is mainly based on what is required by the overall welding process, taking into account the tube/pipe product range, the available space for the impeder and the specified production rate.The choice of frequency is also influenced by considerations regarding the characteristics of the HAZ.
For better understanding of the complex interactions between the above parameters, numerical simulation is an indispensable tool.
One of the earliest studies was by Scott [12], where the two-dimensional (2D) heat equation is numerically solved and analytical formula is obtained for the power density in the tube in terms of some of the operating parameters.Parametric study on the Vee-angle, spring-back, weld speed and frequency with relative temperature distribution has been carried out in 2D for example, by Asperheim et al. [1].A simplified three-dimensional (3D) numerical study of the HF welding was done by Lessmann et al. [9].Jürgens et al. [8] also proposed a 3D FEM model but this had convergency problems with iterative solution process for high welding-feed velocities.One of the most recent 3D simulation approaches has been proposed by Nikanorov et al. [10].All of these papers have in common that they have to make severe assumptions with respect to geometry or process data like weldingline speeds, which is not used in common practice.
In this paper we present the first comprehensive simulation approach for high-frequency induction welding in 3D.Its main novelties are a new analytic expression for the spacedependent velocity of tubes accounting for arbitrary Vee-angle and spring-back.In addition, we utilize a stabilization strategy which allows us to consider realistic welding-line speeds.
The outline of the paper is as follows: In Sect. 2 we describe the mathematical model based on a harmonic vector potential formulation of the Maxwell equations and a nonlinear heat equation.We derive an analytic characterization for a spatially varying tube velocity depending on the Vee-opening and spring-back and its influence on the system of equations.In Sect.3, the approach to the numerical solution is discussed.Results of the simulations are presented in Sect. 4.

Mathematical model of induction welding
The high-frequency induction welding of steel tubes occurs in two stages: electromagnetic heating followed by mechanical joining of the tube edges.It is a multi-physics problem because of the electromagnetic, thermal, mechanical and metallurgical phenomena taking place during the process.Here, we focus only on the electromagnetic heating which is characterized by a coupling of the Maxwell equations with the heat equation.The welding of the tube is a part of a continuously moving production line.Hence, we will consider a quasi-stationary model.

Geometry
Figure 4 shows the different components considered for the electromagnetic heating of the tube.Steel strip and tube is the workpiece, denoted as t .Until the welding point, it is the formed steel strip.After the welding point, it is assumed that the strip edges join to form the tube.The induction coil c is the non-contact electromagnetic heat source which heats the workpiece.Impeder i is a cylindrical solid structure made of an electrically non-conductive material with high magnetic permeability and acts as a magnetic flux concentrator.It is placed inside the formed strip, very close to the inner circumference of the strip.Its length extends until the welding point and ensures the concentration of current on the tube edges and prevents the current from closing the loop on the inside surface strip.The hold-all domain is defined as = t ∪ i ∪ c ∪ a , where a is the surrounding air.
The strip geometry needs to be described such that it represents the forming of the strip during induction welding.The shape of the strip opening before the tube edges join at the welding point is an important factor that influences the welding.Ideally, it should appear as a "Vee" with a very narrow Vee-angle as shown in Figs. 4 and 5 (left).But effects like the spring-back effect can cause a distortion of the Vee shape, see Fig. 5 (right).A parametrized geometry description for t will be given in Sect.2.3.It is needed as a basis for the numerical simulations, moreover we will see that the deviation from a cylindrical shape leads to spatially varying feed velocity.

Electromagnetic induction heating problem
The problem domain is shown in Fig. 4. The problem is defined in a time-harmonic regime where electric and magnetic fields are coupled.Alternating current is supplied in the coil, Figure 5 Variations of tube opening c , which generates a magnetic field around the coil.The time-varying magnetic field induces an electric field creating eddy currents in the electrically conductive steel tube, t .The impeder, i , is used to control the path of the current in the tube.It is made of non-conducting material with high magnetic permeability.The surrounding air, a , is also considered in the model.
We describe high-frequency welding by the time-harmonic Maxwell equations.To this end, we assume that all the fields behave periodically with respect to time with the same frequency ω and all these fields are complex-valued.In eddy current problems, the displacement current is neglected.Hence, the harmonic Maxwell equations are given by Here, E [V/m] is the electric field, B [Vs/m 2 ] magnetic flux density, H [A/m] magnetic field, and J [A/m 2 ] the current density.The system is completed by using the constitutive law.We assume Here, μ [H/m] is the magnetic permeability which is material dependent and may be tensorial in nature.In general, for ferromagnetic materials such as steel, the magnetic permeability will also depend on the magnetic field intensity in a nonlinear way, which would make the harmonic approach not feasible.In order to still consider the equation in frequency domain and to simplify the numerical computation, we introduce a space dependent effective magnetic permeability.This effective value is adapted in a fixed point iteration using the time average of the magnetic permeability, see, e.g., [3,11].
To account for a conductor moving with a relative velocity v, Ohm's law has to be rewritten as Hence, in addition to the induced field resulting from the current within the coil, in principle there is also a contribution resulting from the movement of the tube.However, a dimensional analysis carried out in Sect.2.4 below will confirm that the contribution due to movement of the tube can be neglected for typical weld line speeds.
In view of (1b) we introduce the magnetic vector potential A such that Since A is not uniquely defined by this relation we impose in addition the Coulomb gauge div A = 0. Using (1a) and ( 4) we define the scalar potential φ by Together with Ohm's law (3) we obtain the following expression for the total current density Here, the first term represents the eddy currents and the last one the impressed source current, which lives only in the induction coil c and is given by Now, we utilize (1c) to obtain the vector potential formulation of the Maxwell equations (1a)-(1c), In Sect.2.4 we will show that the velocity term v × curl A can actually be neglected for typical weld-line speeds.Therefore, the final partial differential equation for the timeharmonic Maxwell's equation in vector potential formulation is given by We assume that the tangential component of the vector potential vanishes on the boundary ∂ , i.e., we define The source current J ext as defined in (7) can be precomputed analytically.To this end we assume that the coil is a toroid with rectangular cross-section, where r o and r i are the outer and inner radius, and h is the height.Then assuming continuity for the source current, we demand div(σ ∇φ) = 0.In cylindrical coordinates, we obtain where C 1 = U 0 /(2π) for a given voltage.For a given source current in the coil, C 1 is computed from In cartesian coordinates the source current is then given as Next we consider the thermal part of the problem.The aim of the simulation is to find the temperature distribution in the tube.Hence, only the strip and tube is considered in the thermal part of the simulation.The heat distribution in the tube is determined from the quasi-stationary convection-diffusion equation for heat transport, Here, θ is the temperature The term v • ∇θ represents the transport of the temperature field with v [ms -1 ] denoting the velocity of the tube.With this modification of the heat equation, the movement of the tube is taken into account, and we will see in Sect.2.4 that in contrast to the Maxwell system here the velocity term cannot be neglected.
We have to complement the heat equation with boundary conditions on ∂ t .We distinguish three different parts of the tube boundary, two faces and the lateral surface.At the strip inlet face the temperature is assumed to be at room temperature.Hence a Dirichlet boundary condition of 25 • C is assigned here.On the tube outlet face we assume a homogeneous Neumann condition.On the lateral surface a Robin boundary condition is prescribed, i.e., -λ∇θ Here, h denotes the heat transfer coefficient between the tube and air, and θ 0 is the surrounding temperature.

Tube parametrization and spatially varying velocity
As we have seen already, the strip is bent until it finally becomes a tube at the welding point, see Fig. 5.To parametrize the domain, we assume that the centerline of the tube coincides with the y-axis, such that cross-sections with y < 0 correspond to circular crosssections with (outer) tube radius r t .We assume that half the size of the tube opening is measured by a prescribed function w(y), such that w(y) = 0 for y ≤ 0. It allows to describe arbitrary strip openings including spring back (see Fig. 5, right).We assume that for y > 0 the bent strip is already of circular shape with a radius r(y) and (half ) an opening angle φ(y), both depending on the width of the strip opening w(y) (see Fig. 6).
A solution can only be found numerically.To this end, we fix N > 0 and define equidistant nodes y i = i y max N , i = 0, . . ., N .Note that φ 0 = φ(y 0 ) = 0. Then for 1 ≤ i ≤ N we compute φ i = φ(y i ) by applying Newton's method to solve where we choose the result at the previous node as the initial value for the next one.Then, we apply cubic spline interpolation to obtain a smooth function φ(y), y ∈ [0, y max ] and use (14a) to obtain r(y) (see Fig. 7), i.e.,
Now we are in a position to parameterize the tube domain using cylindrical coordinates with respect to the y-axis.Note that in view of (14a) for a fixed cross-section y ∈ [0, y max ], r(y) is the outer radius and (πφ(y))r(y) the arc length of the open tube arc.Hence, for a tube of wall thickness h, any point in the tube at cross-section y can be described as Here, ρ ∈ [0, 1] defines a radial location between outer (ρ = 1) and inner tube radius (ρ = 0), and α ∈ [0, π] is the angular position between opening edge (α = 0) and tube bottom (α = π ).This parameterization allows for an easy computation of the locally varying velocity.After the welding point is reached the velocity has only a component in the feeding direction y, while before the velocity varies locally with non-vanishing radial and angular components.To obtain the correct temperatures, especially close to the strip edge, it is crucial to use the correct locally varying velocity for the simulation.To this end we consider the trajectory of a point in a tube cross-section moving with constant weld-speed v, i.e., with velocity and Figure 8 shows the resulting local velocity vectors for selected points on the strip opening for a spring back tube opening.Instead of a constant velocity solely in y-direction one can see that now the velocity follows the contour of the opening.

Non-dimensionalization of governing equations
In this subsection we analyze the influence of tube movement on the solutions of the Maxwell system and the heat equation, respectively.Consider first the time-harmonic Maxwell equations in vector potential formulation, Substituting variables for non-dimensionalization, i.e., x = Lx and the magnetic potential vector to be A = A 0 Ã, we arrive at The term i Ã indicates the induced current and v Lω × curl Ã denotes the induced current due to movement.Consider a length scale L to be the penetration depth.The numerical value of the factor v Lω when calculated for high-frequency welding of steel (L = 1 × 10 -3 m and ω = 100 × 10 3 s -1 ) is of the order of 10 -2 .Hence, the current induced due to the movement can be neglected for HFI tube welding simulation.
Fourier's heat equation with convection is given by equation (12), Introducing the non-dimensional constants for length as x = L * x and temperature θ = θ , the equation ( 21) can be rewritten by, Here, the Péclet number is given by the coefficient, If the Péclet number is greater than 1, then this will cause a significant oscillation in the temperature.Indeed, for typical welding-line speeds and materials data for carbon steels it takes values ranging in the interval [12,26].Hence, the influence of the heat transport term during the electromagnetic heating simulation cannot be neglected.

System parameters and material data
Low-carbon steel is one of the most commonly used materials for manufacturing welded steel tubes.In the present work, non-linear material properties of steel strip have been considered.The electrical conductivity and thermal conductivity are shown in Fig. 9 which is for low-carbon steel [4].For the volumetric heat capacity and the definition of magnetization we use parametrized curves as described in [2].Note that these include a number of parameters which have been identified for induction heating processes in a costly process.They are confidential and cannot be disclosed by EFD Induction A.s., but similar data can be found in literature, cf., e.g., [6].
For the volumetric heat capacity we consider with constants ρC p0 , ρC Pi , E, σ , T c , and τ .The formula for magnetic flux density is given by

Discretization and solution strategy
The electromagnetic heating is governed by Maxwell's equation in the time-harmonic regime ( 9) and the quasi-stationary heat equation with transport (12).To reduce the computational effort we only consider one half of the domain, cf.Fig. 10, hence additional symmetry boundary conditions have to be added at this artifical boundary.The computational mesh is generated by the grid generator TetGen [13] and the coupled system is solved using the C++ library pdelib, a finite element toolbox developed at WIAS.
The heat equation is discretized by linear nodal finite elements.The heat in the tube is generated by the effect of Joule heating from the eddy currents induced.This heat is distributed by diffusion.The tube itself is not stationary but moving, which also causes convective transport of heat.The tubes in the welding-line move at very high speeds somewhere in the range of 40 m/min to 200 m/min.As we have remarked earlier, in this range the Péclet number is greater than 1, it means that the process is convection dominated.When using the Galerkin Finite Element method, as Péclet number increases, the solution has a sharp gradient and fails to capture the non-linear solution accurately [5].To stabilize the solution, we utilize the Streamline Upwind Petrov Galerkin (SUPG) method, which adds a residual based stabilization term to the finite element formulation [7].It uses a stabilization parameter which is a function of the local Péclet number, the local velocity and the mesh size and is also referred to as the intrinsic time scale.
In contrast to temperature, where continuity across element boundaries holds, for the magnetic vector potential A continuity across elements can only be expected in tangen- tial direction, in normal direction discontinuities might occur.To this end we follow [6] and use lowest order Nédélec elements for the discretization of the vector potential equation (9).The magnetization depends both on the temperature and the magnetic field (24a)-(24b).For fixed temperature this nonlinearity is resolved numerically based on an averaging approach.For details we refer again to [6].
Finally, the coupled system is iteratively decoupled and solved using a fixed point iteration as depicted in Fig. 11.

Results
For a better visibility most of the results are shown for tube 1, i.e., for the bigger tube diameter.Examples of temperature distribution in the welded strip are shown in Fig. 12.The strip edges are heated to very high temperatures because of Joule heating from eddy current concentration.In the simulations, a spatially dependent velocity is used as explained in Sect.2.3.Figure 13 shows the comparison of temperature distribution results in case of the model with non-spatially varying velocity (velocity only in longitudinal direction) to that of the model with spatially varying velocity.It is seen that the temperature profile follows the strip in the case where spatially varying velocity function is used.
Figure 14 shows the magnetic field lines entering the Vee-opening of the strip, passing through the impeder and then exiting.It must be noted that the number of field lines is not representative of the flux density or distribution of the magnetic flux.At the strip and tube boundaries with the Dirichlet condition, the magnetic field is contained in the welding station.The magnetic field lines around the coil on the underside of the strip are left out of the illustration.A simulation result of current density distribution in the strip edge is shown in Fig. 15.It shows current concentration both in the downstream and upstream directions.As explained in Sect. 1, the principal current paths are already well-known based on general understanding of magnetic field and current.With the developed model, it is possible to quantify the different current levels in the principal current paths (Fig. 15).For closer investigation of the current in different paths, models with four different domain boundaries corresponding to cross-sections 14 , 12 , 11 , 7 are considered (Fig. 16).In these models (tube radius, r t = 20 mm), the strip lengths beyond the impeder in the upstream direction are given in Table 1.The longitudinal current at different cross-sections ( 1 , 2 , . . ., 14 ) of the strip is calculated from equation (25).
In the table of Fig. 16, the current of equation ( 25) is expressed as a percentage of the external coil current (J ext ).Ideally, we would like to have a model where the current at the upstream domain boundary is almost nil.This is seen in the simulation result of the model with domain boundary at 14 .Because the model length is very long, the mesh quality maybe compromised.Hence, we try to choose a model with shorter strip length  2).As seen above, higher current density along the strip edge is observed both upstream and downstream of the coil.It reaches a peak downstream at the welding point whereas it decreases upstream.This is in accordance with the already known principal current paths as shown schematically in Fig. 3.The division of current path in the upstream and downstream directions is also reflected in consequent Joule heating of the strip edges in Fig. 17.The temperature drops midway of the coil because of the currents dividing in the two opposite directions and a considerably lower strip edge current density just upstream of coil center.The highest temperatures are reached close to the welding point.The previous results were analysis of the current density and temperatures along the length of the strip.The next set of results are analyzed at the cross-section of the strip.Figure 18 shows the temperature for different cross-sections in the strip.It is seen that the temperature profile in the strip edge has an hourglass shape.This is related to the current distribution in  the cross-section of the strip edges.The current distribution is influenced by the current penetration depth, the geometry of the strip edges and the position of the impeder.The temperature distribution is also a result of thermal conduction in the strip cross-section where power generation and heat transfer take place.The heat affected zone at the welding point is shown for two different tube dimensions but same operating frequency in Fig. 19.The hourglass shape is more pronounced in case of the thicker tube for same current frequency.At the welding point, the temperature of the strip edge at the outer circumference and inner circumference may differ depending on the closeness of the impeder to the inner circumference.When the distance between the impeder and the inner circumference of the strip is large, the magnetic field lines are spread out in this gap before concentrating into the impeder.This causes more heating   of the strip edges at the inner circumference.When the impeder is close enough to the inner circumference of the strip, most of the magnetic field lines are directed more vertically out of the gap between the two strip edges, with less heating of the corners at the inner circumference as the result.However there are physical limitations of how small the gap between the impeder and the inner circumference of the strip can be.At this point, a parameter called HAZ width is defined to quantify the size of the heat affected zone (HAZ).It is defined at the welding point cross-section of the strip as the region where the temperature is greater than 650°C.It is the width through mid-thickness as indicated in the Fig. 20.The simulation results for Tube 1 (Table 2) at a frequency of 250 kHz is taken as a reference case.The imposed current to the coil is tuned to get 1250°C at the welding point.The resulting HAZ width is taken to be of value 1.The imposed current for the reference case is also taken as 1.The role of the impeder which is essentially a ferrite core is to direct or concentrate the magnetic field.Or in other words, it has to function as a good conductor of magnetic flux.In absence of the impeder, the magnetic field will flow within a certain depth of the inner circumference of the strip.As a result, a high-frequency leakage currentwill flow in the inner circumference and generate Joule losses and unwanted heating of the inner surface of the strip.However in the presence of the impeder, most of the magnetic field flows within its magnetic material, reducing leakage current and unwanted heating.The impeder improves the efficiency of the process.The impeder's effect on the shape of the HAZ is seen in Fig. 21.It is also observed that the imposed current needs to be increased by 40% to achieve the desired temperature of 1250°C at the welding point for this specific simulation set-up.There are several operating parameters that affect the shape and size of the HAZ.The parameters of frequency and vee-angle were chosen to analyze the HAZ width.For the dimensions of Tube 1 (Table 2), as mentioned earlier the reference simulation case for a frequency of 250 kHz and 3°vee-angle is taken.The HAZ width is compared for frequencies: 150 kHz, 250 kHz, 350 kHz and vee-angles: 3°, 5°, 7°.From Fig. 22 it is observed that for any given vee-angle, the HAZ width reduces with higher frequencies.At any given frequency, the HAZ width significantly increases with a wider vee-angle.For example, at a frequency of 250 kHz the HAz width is 60% higher for a 7°vee-angle than that for a setup with 3°vee-angle.

Conclusions
A three-dimensional model has been developed for high-frequency induction welding.It is a non-linearly coupled system of Maxwell's electromagnetic equation and the heat equation.The system equations are discretized by FEM and solved using the fixed point iteration based on an incremental decoupling.The model includes the use of non-linear magnetic and thermal material properties.The model also considers the high welding-line velocities which accounts for the heating by advection and conduction in the strip.The results show temperature distribution in the strip edges that develops as expected from previous studies and visual observations of the process.First, the corners of the strip are heated and as the welding point approaches the characteristic hourglass shape of the heat affected zone is formed.The boundary conditions at the strip and tube ends are considered to be Dirichlet conditions as it represents the process correctly.The strip length used in the model also decides the amount of induced current that goes to the welding point.
For the thinner wall tube, the hourglass shape of heat affected zone is less pronounced.This is in line with what is expected from literature.The observations from the results are dimension specific and not necessarily the same for tubes of other dimensions.This new three-dimensional simulation tool will provide a basis for an optimization of the design of the welder, especially with respect to the dimensioning of induction coil, impeder and the configuration of these relative to the steel strip.Future work will include the study of the mechanics of the material squeeze-out when the strip edges are joined together after heating.

Figure 1 Figure 2 Figure 3
Figure 1 High Frequency induction welding of steel tube

Figure 4
Figure 4 Weld setup for electromagnetic heating

Figure 6
Figure 6 Parameterization of strip geometry

Figure 8
Figure 8 Local velocity to represent feed velocity and forming

Figure 9
Figure 9 Electrical and thermal conductivity of low-carbon steel

Figure 10
Figure 10 Weld setup for current controlled simulations

Figure 12 Figure 13
Figure 12 Temperature distribution in the tube

Figure 14 Figure 15
Figure 14 Path of magnetic field lines entering and exiting the Vee-opening of the strip

Figure 16
Figure 16 Current in Y-direction in different strip cross-sections (slices)

Figure 17 Figure 18 Figure 19
Figure 17 Temperature in the strip edge

Figure 21
Figure 21 Effect of impeder on HAZ

Table 1
Models with different domain boundaries