Numerical prediction of Joule heating effect in electric hot incremental sheet forming

Since the deformation region involves the interaction of electric-thermal-force coupling in electric hot incremental sheet forming, the numerical simulation of the forming process is unusually difficult. Currently, the thermal-force coupling method is adopted to simulate approximately the whole forming process, and the Joule heating effect is often ignored. Therefore, the numerical simulation of the Joule heating effect is especially significant for the prediction accuracy of forming process. In this paper, a novel numerical simulation method considering electric-thermal-force parameters was proposed to update the thermal-force condition of the forming region instantly. Meanwhile, the model of contact thermal conductance was established, combining geometrical and electric-thermal parameters. Then a high-precision finite element model was obtained to predict the Joule heating effect of the forming region. In addition to this, the impact of thermal superposition on forming temperature was further analyzed, and a modified model of contact thermal conductance was established in electric hot incremental sheet forming.


Introduction
Incremental sheet forming (ISF) is viewed as laminated manufacturing technology and is often used to fabricate small-batch products. The manufacturing cost is effectively reduced since the corresponding die could be simplified or removed [1]. Meanwhile, the technology is used to fabricate materials with good formability at room temperature, such as steel and aluminum alloys. With the development of science and technology, some high-performance materials, including magnesium and titanium alloys, are gradually applied in various sectors. However, this material has poor formability at room temperature and good ductility at elevated temperatures. Therefore, electric hot incremental forming (EHIF) is proposed to fabricate parts with magnesium and titanium alloys, for example, Ao et al. [2], Saidi et al. [3], Bao et al. [4], Ajay [5], Adams and Jeswiet [6], Xu et al. [7], and Li et al. [8].
During electric hot incremental forming, some features, such as forming temperature, forming accuracy, and material properties, change with forming process parameters in time and space. Therefore, the numerical method can visually analyze the variation of characteristic variables. More and more commercial software with higher computational efficiency is adopted to simulate the forming process of EHIF. Fan et al. [9] adopted the hexahedral element of SOLID 7 of the MSC. Marc software to simulate the forming process of EHIF based on the coupling model of electricity, displacement, and temperature. However, this module could not simulate the complex deformation process, and then the first machining path is only calculated. Afterward, Honarpisheh et al. [10] simplified Joule heating effect as heat flow to simulate the EHIF process of Ti-6Al-4 V, namely, that the coupling effect for three fields was converted to the coupling effect of two fields. However, the interaction of electric-thermal-force coupling is involved in EHIF, and the conventional simulation cannot real-time describe the process of electric-thermal-force coupling caused by the Joule heating effect.
To solve the above problem, a novel numerical simulation, namely, a direct simulation method of electric-thermalforce coupling, was proposed to describe the Joule heating effect of EHIF and its flow path shown in Fig. 1. Thermal interaction boundary conditions, combining geometrical and electric-thermal parameters, were established to predict the change rule caused by the Joule heating effect. In addition to this, the effect of thermal superposition on forming temperature was further analyzed, and the direct simulation method with thermal superposition was obtained in EHIF.

Electric-thermal boundary models
In the process of EHIF, the heat of forming region is calculated based on Joule's law, and it is written as where dQ is the instantaneous heat of the forming region due to the fact that dt is tiny. I is the value of currents, R is the equivalent resistance of the forming region, and dt is the active time of currents. The room temperature is viewed as 25 °C and it is viewed as the reference temperature (T 0 ). Based on the reference temperature, the temperature difference (ΔT) is obtained according to the law of temperature rise of solid metal, and the corresponding model is written as where c and ρ are separately the specific heat and density of materials and ΔT is the difference between the forming temperature (T) and T 0 . Combining Eqs. (1) and (2), the computational model of ΔT is further written as According to the above analysis, dV/dt and R are two key factors determining the value of ΔT when I, c, and ρ have been determined. Figure 2 shows the sketch of the contact region between the tool and the sheet, and the volume of the forming region is obtained by multiplying the deformation zone and the machining path in a short time. The geometric profile of the forming region is viewed as a cube since the active time of currents is extremely short. According to the study of Min et al. [11], the dV/dt is written as where a and b are two constants and separately obtained Eqs. (5) and (6), in which r T is the radius of forming tool.
Because of the characteristics of local contact, the equivalent resistance value is mainly the sum of the material resistance and the contact resistance in a short time. The contact resistance between the tool and the sheet largely depends on the electric properties of materials and the temperature of forming regions. The contact region would not change significantly with the time scale change due to a uniform feed rate and constant forming pressure. Therefore, the effective area of material resistance is described through the moving site of tools in a bit of time and shown in Fig. 3. Based on the resistance calculation formula and Eq. (4), the resistance value of materials is obtained according to the following expression: where R s is the resistance value of contact regions and ρ s (T) is the resistivity of material as a function of temperature, and the two parameters have been calculated according to the study of Fan and Gao [9]. Meanwhile, the contact resistance between the tool and the sheet is calculated according to the following model [9]: where R j is the contact resistance, R j (T r ) is the contact resistivity at room temperature, H(T) is the hardness of material as a function of temperature, and H(T r ) is the hardness of material at room temperature. Combining Eqs. (7) and (8), the equivalent resistance of forming region is written as Figure 4 shows the thermal interaction principle in the numerical simulation of EHIF, in which T T and T S are separately the temperatures of tools and sheets and K j is the contact thermal conductance between the tool and the sheet. According to the heat balance relationship, the heat relation of numerical simulation is written as where S j is the contact area between the tool and the sheet in numerical simulation, and it is further obtained according to Eq. (11).
In the simulating process, the model of K j is written in FORTRAN, in which some parameters, such as I, t 0 , θ 1 , R, dV/dt, T T , and T S , are imported, and the surface heat flux between the tool and the sheet is calculated to simulate the Joule heating effect of the forming region.

Representation of friction model
The Coulomb friction model is usually adopted to the interaction between the tool and the sheet in the numerical simulation of EHIF. The contact surface of tools is set as a rigid surface. The sheet with a deformable body is viewed as a slave surface. Figure 5 shows the normal position of master-slave surfaces, which is used to judge the contact relationship between the tool and the sheet, in which ɡ N and ε N are separately the normal gap of master-slave surfaces and the penalty stiffness, and the judging rule is established based on Eq. (12). The contact between the tool and the sheet is determined when ɡ N is greater than zero, and then the corresponding normal stress can be calculated by multiplying ɡ N and ε N . On the contrary, the contact stress is zero.
For isotropic materials, the tangential stress, namely, the friction stress between the tool and the sheet, is calculated according to the Coulomb friction model, and it is written as where τ sl is the scalar of tangential stress and μ is the friction coefficient between the tool and the sheet. ɡ T is the tangential displacement of master-slave surfaces, and it is calculated according to Eq. (15) where t 1 , t 2 , and ̇i i (Fig. 6) are separately time nodes and the tangential velocity of contact points. τ f is the tangential stress between the tool and the sheet.

Materials and simulation setting
In EHIF, the cone with Ti-6Al-4 V was fabricated to verify the accuracy of numerical simulation, and the corresponding geometric profile is shown in Fig. 7. The test platform (Fig. 8), including forming equipment, forming devices, and thermal imagers, was adopted to fabricate the part designed and collect the deformation region's forming temperature. The thermal imager consists of two cameras, such as lowand high-temperature cameras. The low-temperature camera with a temperature range of 0 to 200 °C is adopted to collect the temperature image of the forming region, and the hightemperature camera is used to collect the max temperature of the forming region.
In this work, the dynamic-temp-disp-explicit module was used to simulate the forming process. The solid element (C3D8RT) with eight nodes was used to simulate the deformation process of materials, and Johnson-Cook was adopted to calculate the equivalent stress varying with temperature based on the study of Honarpisheh et al. [10]. Meanwhile, thermal-force properties of Ti-6Al-4 V are separately shown in Tables 1 and 2. Figure 9 shows the finite element model of EHIF, and the holder region of sheets is divided by large-scale mesh (5-mm global size), and the forming region is refined by small-scale mesh (1.5-mm approximate size). The position of the sheet is between the binder sheet and the holder sheet. Meanwhile, binder and holder plates with 1-mm thickness are separately set as discrete rigid bodies (5-mm global size). In addition to this, the tool is set as a rigid analytic body, and the initial position of the tool is the initial point of the forming path. The thermal exchange area is mainly the sheet-air surface, which is set to describe the heat loss of the forming region. The same process parameters were adopted in the experiment and simulation and shown in Table 3, and the current value could be calculated according to Eqs. (3), (4), and (9). The ABAQUS-VUINTER with Eqs. (12), (13) and (14) was used to simulate the interaction of electric-thermal-force of EHIF.

Effects of thermal superposition
According to the characteristic of layered processing, a small portion of the forming region is influenced by the continuous action of tools. Therefore, if the interaction period between the tool and the contact region is short, the region will be affected by the thermal superposition effect, which is shown in Fig. 10. The tool moves from point C to point D, and the overlapping zone between two points is influenced through the thermal superposition effect due to a short cycle. Meanwhile, since the step size is far less than the contact area, the effect of the step size on the thermal superposition area can be ignored. Therefore, the thermal superposition area is approximately equal to the contact area between the tool and the sheet. According to the above analysis, the forming temperature increases with the increase of processing time, which is shown in Fig. 11. Meanwhile, the formability of materials can be influenced due to high-temperature oxidation.
Since the thermal loss of contact region is mainly the thermal exchange between the sheet and the air, the following assumptions need to be proposed before analyzing the heat superposition effect, as follows:  2. The thermal exchange is far greater than the thermal radiation, and then the effect of thermal radiation on thermal loss is ignored; 3. The analysis of the thermal superposition effect is based on room temperature; 4. The contact area between the tool and the sheet is content.
According to Fig. 10, the heat of point C can be expressed as follows: In unit time, the thermal exchange (Q air ) between the sheet and the air can be written as follows: where W A is the convective coefficient of air, and it is set as 28 × 10 −6 W·mm −2 , and S A is the thermal exchange area between the sheet and the atmosphere.
Combining Eqs. (16) and (17), the thermal loss of unit time is calculated in the cycle period of N, and it is written as Fig. 8 The test platform of EHIF where Q air-N is the thermal loss of the N-th unit time in the cycle period of N and Q air-N-1 is the thermal loss of the N-1'th unit time. According to the above analysis, the total thermal loss is obtained in the cycle period of N based on the superposition principle, and it is written as where Q air-T is the total thermal loss in the cycle period of N and Q air-i is the thermal loss of unit time in the cycle period of N. Figure 12 shows the analysis of the thermal superposition of EHIF, in which the thermal loss increases with the increase of the cycle period, and the superposed temperature gradually decreases. Meanwhile, the temperature fluctuation is gentle when the cycle period is greater than 20 s, and the error of forming temperature is less than 5%. Then the thermal superposition effect is ignored in this cycle. On the contrary, the thermal superposition effect needs to be considered.
According to the above analysis, the relationship between the total heat and the cycle period is established by the Fourier transform method when the cycle period is less than or equal to 20 s, and the corresponding formula is written as where Q D-i is the total heat after thermal superposition in the ith cycle period and i is a cycle period whose value is less than or equal to 20 s. Therefore, Eq. (12) is further modified considering the thermal superposition effect, and it is written as (20)    Figure 13 shows the comparative analysis between two conditions, considering and non-considering superposed heat. The difference between considering superposed heat and the designed temperature is slight, and the corresponding error is less than 5%. However, the instantaneous forming temperature without considering superposed heat increases gradually with the increase of forming time, and the maximum difference reaches 60 °C. Therefore, the forming temperature considering superposed heat is more uniform, and then the temperature accuracy is further ensured in EHIF. Figure 14 shows the comparative analysis between simulating and experimental temperature distribution in the initial forming stage, in which the error between the simulating result and the actual temperature is less than 5%. Meanwhile, the temperature distribution of the two cases is similar, and the high-temperature point, which always appears in the area after the contact between the tool and the sheet, has a hysteresis.

Prediction of Joule heat effect
The maximum temperature of four-quadrant points is further analyzed to verify further the accuracy of simulated results, shown in Fig. 15. The maximum temperature of the first four layers is collected to examine the simulating accuracy. The average difference between the simulation and the reality is less than 5% in each machining layer. Although the difference between the actual temperature and the simulating temperature is obtained, the average error rates are both less than 5%. Therefore, the difference between the actual temperature and the simulating temperature is not significant and it can be ignored for the simulation method. Meanwhile, the maximum temperature between four-quadrant points has a stable fluctuation in the simulation of EHIF, which also conforms to the actual tendency. Figure 16 shows the temperature distribution of forming part with the thermal superposition effect. According to Fig. 16a, the simulating temperature of forming region presents a symmetrical distribution, and the highest temperature appears at the endpoint. The bottom region unformed would obtain a higher temperature due to heat conduction. Meanwhile, the bottom temperature is higher than that of the sidewall, consistent with the actual result from Fig. 16b. The color at the bottom of the part is the most significant according to Fig. 16b, which cannot show that each area of the bottom reaches the maximum temperature. The phenomenon is obtained through the working principle, which is stated in Sect. 2.3, of the thermal imager. Meanwhile, the heat of the region of the max temperature is quickly passed to other regions at the bottom, and the time of this transfer process is very short, and the state of the transfer process is not easy to capture. Therefore, the difference between Fig. 16a, b is obtained. However, the max temperature of the end point is similar for the simulating result and the actual result, and the difference rate is about 1.2% and is less than 5%. Meanwhile, the gradient temperature distribution characteristics of the forming region is also consistent for the simulating result and the actual result. According to the above analysis, the simulation method proposed can accurately describe the Joule heating effect in EHIF.
Finally, the axial forming force between simulation and reality is further analyzed according to Fig. 17, in which the change tendency between simulated and actual results is similar. Change characteristics of axial forming force are accurately predicted based on the numerical simulation method in the initial and final stages. In addition to this, the coincidence degree between the prediction force and the actual force is higher. Therefore, the numerical simulation method proposed can accurately simulate the electric-thermal-force interaction of EHIF. Fig. 12 The analysis of thermal superposition of EHIF Fig. 13 The comparative analysis between considering and non-considering superposed heat

Conclusions
Aiming at the complex numerical simulation problem of EHIF, the mapping relationship between the geometric profile of deformation regions and the forming process parameters is analyzed from the perspective of the macro deformation mechanism, and the corresponding calculation model is established. Meanwhile, the equivalent resistance of contact regions is also established based on material resistance and contact resistance models. On this basis, an analytical model of contact thermal conductance, considering electrical-thermal parameters, forming process parameters, and thermal superposition effect, is established. In addition to this, the ABAQUS-VUINTER subroutine with the above models is adopted to update simulated results of the electric-thermal-force interaction in EHIF, the accuracy of which is verified through the corresponding physical experiment. Therefore, the finite element simulation technology for electric-thermal-force coupling is proposed to predict more accurately the forming process of EHIF.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.