Performing an Indirect Coupled Numerical Simulation for Capacitor Discharge Welding of Aluminium Components

: Capacitor discharge welding (CDW) for projection welding provides very high current pulses in extremely short welding times. This requires a quick follow up behaviour of the electrodes during the softening of the projection. The possibilities of experimental process investigations are strongly limited because of the covered contact zone and short process times. The Finite Element Method (FEM) allows highly resoluted analyses in time and space and is therefore a suitable tool for process characterization and optimization. To utilize this mean of optimization, an indirect multiphysical numerical model has been developed in A NSYS M ECHANICAL APDL. This model couples the physical environments of thermal–electric with structural analysis. It can master the complexity of large deformations, short current rise times and high temperature gradients. A typical ring projection has been chosen as the joining task. The selected aluminium alloys are EN-AW-6082 (ring projection) and EN-AW-5083 (sheet metal). This paper presents the investigated material data, the model design and the methodology for an indirect coupling of the thermal–electric with the structural physic. The electrical contact resistance is adapted to the measured voltage in the experiment. The limits of the model in A NSYS M ECHANICAL APDL are due to large mesh deformation and decreasing element stiffness. Further modelling possibilities, which can handle the limits, are described.


Introduction
Resistance projection welding is a versatile welding process. It is characterized by various geometries (e.g., ring projection, segment projection and long projection), machine characteristics (electrical source and machine construction) and materials (steel, aluminium, mixed joints). Resistance projection welding for aluminium components is hardly used. The reason is very high welding currents are needed to generate enough heat within materials with high thermal and electrical conductivities. That can be realized by capacitor discharge welding (CDW). The investigation of the process with experimental methods is strongly limited. In most cases, it is only possible to measure process parameters at discrete locations. Furthermore, methods of destructive test, such as strength tests or structural analysis of the joints are used to investigate their quality. As a result, many of the current investigations analyze the influence of welding parameters on the joined components [1][2][3]. Besides these correlations between the process parameters and the characteristics of the joint a more detailed understanding of the process is not possible with the mentioned methods. There is no way to investigate the temperature distribution, voltage drops at the interfaces of the components or even the deformation of the projections. To overcome these barriers for a more detailed understanding of the process, the FEM provides a suitable tool to investigate inner processes as well as the hardly measurable quantities such as deformation, stresses, temperatures in the joining zone and current densities. It allows the local and temporal resolution of the whole welding process. In the scientific publications [4][5][6][7][8], projection welding by capacitor discharge welding is investigated. In [4,5] a mathematical model for determining the shear stress and temperature distribution is presented. However, no temperature-dependent deformations are considered. In [6][7][8] an iterative model is developed to determine the weld parameters of a wave profile. The width of the connection is also investigated. Therefore, the FEM might be used for process optimization as well. Scalable parameters of the model such as geometrical conditions, input parameters like welding current and electrode forces, and boundary conditions provide an easy and fast way to investigate the influences of the named parameters on the process without expensive and time-consuming experiments. Projection welding by capacitor discharge produces a high and temperature-dependent deformation, which causes convergence problems in direct coupled numerical models. For these reasons, an indirect coupled FEM Simulation regarding structural, thermal and electrical phenomena is established in this paper. ANSYS MECHANICAL APDL was used as the software environment. Practical, as well as theoretical, aspects in modelling the complex phenomena are described. Methods of modelling a process are shown, firstly examined by experiments only by adapting the calculation of contact resistance. Thereby common models of contact resistance are used and we take advantage of the already included numerical parameters.

Capacitor Discharge Welding (CDW)
CDW belongs to the resistance welding processes and is mainly used for projection applications. Figure 1 shows, as an example, the welding current I and electrode force F as a function of time t for one of the conducted experiments. With the starting time t 0 (t = 0 ms) the current flow starts. At t p (t = 3.6 ms) the maximum current I of 113 kA is reached. At t h (t = 7.1 ms) the current decreased to 50% of its maximum, which is defined as the welding time [9]. Additionally Figure 1 shows four characteristic process phases, which are typical for steel alloys and were elaborated in [10,11] In comparison to spot welding with the medium-frequency inverter technique, projection welding with capacitor discharge has some complexities concerning the characteristic process phases:
High deformations of the projection 3.
Activation of the contact surface by metal vapor

Geometry and Material Data
The technical drawing of the projection component is shown in Figure 2. The projection component is made of aluminium alloy EN-AW-6082. The thermophysical properties of the aluminium alloys were determined by a material simulation using the software JMATPRO. The thermophysical properties of EN-AW-6082 are shown in Table A1 and the flow stress in Table A4. The geometric dimensions of the projection were optimized within a research project conducted also by the authors at Technische Universität Dresden [12]. The experimental data result in a recommendation of the width of the contact area of 0.5 mm and a projection angle of 120°. For the sheet component, the aluminium alloy EN-AW-5083 is used. The thermophysical properties of EN-AW-6082 are shown in Table A2 and the flow stress in Table A5. A common resistance welding electrode type C20 according to [13] is used. This electrode is made of CuCr1Zr. The thermophysical properties of CuCr1Zr are represented in Table A3 and the flow stress in [14].

Application Specific Modelling
The investigated 2D axisymmetric finite element model is shown in Figure 3. The quality of the mesh affects the solution of the simulation. Therefore, a geometry specific mesh is generated by analysing the convergence of the displacement. As a critical area, where the results have a big effect on the whole simulation, the penetration of the projection into the sheet is evaluated. The size of the elements was reduced until the penetration results did not change with further refinement of the mesh. This analysis was done at room temperature. The optimized element size in the contact area for this joining task is 0.0178 mm. With further distance from the contact area, the element size is growing. From the numerical point of view, the edges of the projections contacting surface are singularities. The settings of the contact elements have to be chosen carefully to keep the difficulties under control coming with singularities. Therefore, the contact pressure of the joining zone was investigated at room temperature. The elaborated settings (key options) is shown in Table 1. Adapting the model to the conditions at the beginning of the process enables a stable simulation until temperatures of about 600 K. However, with rising temperatures, and resulting from this the lowering of element stiffnesses as well as thermal expansion, the simulation terminates due to unconverged loadsteps of the structural physic. At the upper electrode the force F(t) and the current I(t) is applied as load from the measured data ( Figure 1). As a boundary condition the displacement u y is blocked at the lower electrode and the electrical potential is 0 V. The uniform temperature T uni is 298 K. There are three contact areas, which are made of contact and target elements: Projection (Contact) -Sheet (Target) 3.
Sheet (Contact) -Electrode (Target) The behaviour of all contact surfaces is standard (not bonded, no roughness, no separation).

Implementation of Electrical Contact Resistance
There are many theories for calculating the ECR in resistance welding [15]. There is no validated theory for projection welding by capacitor discharge for aluminium components. The contact theory in [14,16] describes the ECR according to Equation (1). This theory has been validated for aluminium spot welding by the medium frequency inverter technique [17,18].
The contact theory in [19] describes the ECR according to Equation (2). This theory has been validated for welding of steel and brass projections on steel sheets with a direct current.
The discussion of the contact theory in [19] is done by means of calculations with κ set to 0.5. SONG inserts a parameter for the film resistance ρ f ilm . A good agreement with experiments is achieved. The pressure and temperature dependency of the ECR can also be described by the polynomial in Equation (3) with a varying exponent z. The author states the need of the variation in order to describe the temperature dependency of the pressure influence on the ECR. Transferring this insight to the application of Equation (2) in a coupled FE-model the variation of the exponent κ seems to be appropriate for modelling all phenomena taking place at the interfaces.
Despite the different theoretical principles behind the theories, the effects of the exponent κ are similar. It balances the pressure against the temperature dependency of the ECR. It relates to the surface condition of the interfaces [14]. With progress in the simulation time, changes in the surfaces conditions are expectable from experimental investigations [10]. Therefore, a time-varying exponent corresponds to the phenomena in the real process. Namely, the disruption of surface film and flattening of asperities. For using the ECR as element parameter in the ANSYS MECHANICAL simulation, the parameters ECC and TCC are necessary. ECC is calculated by Equation (4). TCC is calculated by Equation (5) with the LORENZ Number L of Equation (6) and the current temperature T of the element [14,16].

Indirect Coupling of a Numerical Simulation
There are two methods of coupling possible to model the interactions between structural, thermal and electrical phenomena. First, a direct coupling condenses the governing equations in one equation system. This is realized by the coupled field elements. By neglecting of the magnetic effects, the degrees of freedom are: the displacement u, the electric potential Φ and the temperature field T [20].
Secondly, an indirect coupling uses more than one equation system while the coupling is executed by transferring results of one equation system to the other and vice versa. Two systems of equations are established with Equations (8) and (9) which will be solved individually. In this case, thermal-electrical effects are coupled direct while the coupling with structural phenomena is established indirect. This results in a separated modelling of the two different physics, thermal-electrical and structural. Investigations of welded joints show big deformations of the projection near the joining zone [20].
Modelling big deformations with finite elements leads to disorted meshes, where bad qualities of the results or even unconverged loadsteps might result from [21]. ANSYS MECHANICAL provides the methods of nonlinear adaptive meshing and rezoning to overcome those difficulties. For two-dimensional simulations, both methods are restricted to the solid structure elements PLANE182 and PLANE183 as well as the contact elements TARGE169, TARGE170, CONTA172 and CONTA174. However, to model thermal-electrical as well as structural effects other kinds of elements are needed (e.g., PLANE223). The separation of the physics provides the possibility to apply the rezoning method provided by the software at least for the structural physic. With the requirement of equally meshed physics used in an indirect coupled simulation, a way to transfer results of the thermal-electrical physic onto the new mesh still has to be found. The indirect coupling offers the following advantages: 1.
Minimization of the computing effort 2.
Adaptation of the solver to the requirements of each system of equations 3. Dynamic calculation of the contact resistance depending on the physic conditions without tabular precalculations (material data, temperature distribution and contact pressure) 4.
Rezoning in structural physic (No rezoning for Coupled Field Elements) Figure 4 shows the schematic process of the developed indirect simulation model. This is realized by a do loop. With the start of the simulation, the input data (F(t), I(t)) were imported. The input data were provided by the measurements of the experiments as shown in Figure 1. After this, the two physical environments (structure and thermal-electric), short physics, were defined. At the beginning of the first iteration loop, the structure physic was read and solved. As a result, the contact pressure σ cont (normal stress in contact plane) of all contact zones was exported. This enables the first calculation of the electrical contact resistance (ECR) at room temperature (see Section 3.3). After calculating the ECR, the solution of the structure physic was saved and the thermal-electrical one was loaded. This was followed by an update of the geometry due to the temperature-dependent deformation and the solution. As a result, the temperatures of the nodes are exported and the solution of the thermal-electrical physic is saved. The next loop iteration follows by applying the temperature field from the thermal-electrical physic as a boundary condition to the structure physic. The rapid rise of the temperature leads to thermal expansion, softening of the elements and therefore to bigger deformations compared to room temperature. This results in a new contact situation, so that the ECR was updated with the new temperature field and the new contact pressures. These iteration steps ∆t were repeated until the current simulation time t n reached the end time t end of the simulation. Afterwards, the post-processing started. The physical environments were established by declaring the geometry of the model, meshing of this geometry, application of boundary conditions and solver settings. As mentioned above, the meshes of both physics have to be identical regarding the coordinates of the nodes and elements. Therefore, meshing has to be executed in the structural physic only. Afterwards, the established thermal-electrical physic built up onto the existing mesh by declaring different attributes to the elements. Both physics can be saved and loaded during the simulation as described in Figure 4. Within an indirect coupled simulation, special care has to be taken to the transfer of results between the physics:

1.
To simulate the thermal-electrical phenomena with the current state of deformation, the calculated displacement in the structural physic is transferred by the UPGEOM command as a change of geometry.

2.
Coupling in the contrary direction is implemented by applying the calculated temperatures as body loads by the LDREAD command. The temperatures are transferred to each node. Therefore the current state of thermal strain can be involved in the structural solution. 3.
The coupling of effects within the contact zones was established by using the electrical contact conductivity ECC and the thermal contact conductivity TCC. Their values were calculated elementwise by Equations (4) or (5) and were therefore location-dependent. Before solving the thermal-electric physic starts, ECC and TCC were handed over to the contact elements by the RMODIF"19 and RMODIF"14 command respectively. Both parameters were declared as table arrays using the primary variable x-coordinate (2D) to guarantee every element gets their corresponding value.
The alternating calculation of the physics allows updating of the solver setting depending of the other physics than the one being solved. e.g., assuming large thermal strain with steeply rising temperatures the loadsteps and substeps of the structural solution might be chosen smaller when large changes of the temperature were detected in the thermal-electric solution.

Experimental Determination of the Electrical Voltage
In order to adapt the numerical factors of the contact theory, the electrical voltage drop was measured during the welding of the projection component (Figure 5a). Similar method for resistance spot welding of steel alloys was investigated by [22]. The voltage was measured at three positions (electrode-electrode U e,e ; projection-sheet U p,s ; sheet-electrode U s,e ). From these, the voltage drop at the electrode-projection U e,p interface was calculated.  The electrical voltage was measured with 1 MHz and a resolution of 16 bit (measuring range from 0 V to 10 V). Figure 5b shows the standard deviation of all (7) measured electrical voltages and the mean value. The sum of all individual voltages gives the total voltage between the electrodes U e,e . The measured voltage curves include inductivities. These were corrected according to [23]. The electrical voltage U p,s is lower than the others because the contact area is the smallest. Therefore, the contact pressure increases and the contact resistance decreases compared to the other contact zones (see Equation (2)). A strong increase of the individual voltage can be seen in the first 0.5 ms. Afterwards, the increase becomes flatter. The characteristic curve of U p,s is different to U e,p and U s,e .

Comparison of Experimental and Simulative Results
All simulations were executed until their termination due to convergence problems arising from large distortion of the mesh, low stiffnesses of elements with high temperatures and high contact pressure. The range of values for the numerical parameters d and κ used in the Equation (1) as well as ξ and κ in Equation (2) are declared by the authors SONG and WANG. It was not possible to calculate the measured voltages across the interfaces of the experiments within this range. The simulation results shown here were determined according to SONG (Equation (2) with ρ f ilm = 0). SONG recommends the numerical parameters ξ from 0.1 to 10 and κ from 0.5 to 1. In this range, the simulated electrical voltages were strongly outside the measurement. For this reason, the parameters were adapted to the measured voltages to reproduce the experiment. The new parameters are shown in Tables A6 and A7. In simulation 1, 2 and 3 the value of κ was set to 1.
Uniformly defined values lead to a distribution of voltage drops, which is not consistent with the experiments. This is shown in Figure 6a. The three contact zones differ from each other regarding materials as well as the surface treatment and roughness. This implies that the numerical parameters should also be treated separately. This can be seen in Figure 6b-d. This results in a distribution of the simulated voltage drops at each contact as measured in the experiments. Figure 5b shows the characteristic voltage curves of the contacts U p,s , U e,p and U s,e are different from each other. The contacts are modelled with smooth surfaces. Therefore, it is not possible to include specific effects in the contact zone. One of these effects is the different behaviour of the asperities at the different contact surfaces, that is due to their different material characteristics. Furthermore, the behaviour of the surface films cannot be modelled, due to neglecting the specific film resistance ρ f ilm . To model the change of these effects, a time variable exponent κ is introduced. The simulation results in Figure 6b,c show that a constant set of numerical parameters does not lead to the measured voltages over process time. Either the simulated voltage, or their gradient does not match the experimental results. For bigger values of the factor ξ in Equation (2), the strong increase of voltage drops at the beginning of the process can be modelled. It results in a moment of time where the calculated voltages exceed the measured values while maintaining their high gradient. After this point of time, the values of the simulation were constantly to big. Choosing low values results in not reaching the measured voltages at all. With κ = 1, the voltages are smaller then the experimental values. If κ < 1, the voltages become bigger. Therefore, the values of the exponent are changed as shown in the Table A6. Figure 6d shows the effect of the final configuration of κ. With the simulation, four values were found which can model the experimental measured voltages in a satisfying manner. The simulation terminated at a time of 1.67 ms and a maximum temperature of 671 K located in the joining zone. The convergence problems were due to thermal expansion and low element stiffnesses resulting from high temperature. This is accompanied by large deformations.

Discussion
With the elaboration of the FEM model in this work the most suitable settings where found regarding the mesh quality, calculation methods, boundary conditions and contact settings. Additionally assuming correct material data one can make the assumption that the numerical errors in the results are neglectable for this type of investigation. Because of this, the only way to reproduce the measured voltages is to adapt the contact resistance. For the reason of adaption the authors of contact the resistance models discussed in this work already included parameters to be changed in their equations. Taking advantage of this, it is possible to reach very good agreement of the FEM simulation with the real process regarding the electrical phenomena. By calculating the thermal conductance from the electrical conductance by the commonly used WIEDEMANN-FRANZ-Law the quality of the electrical results applies for the thermal phenomena too.
The results show big influence of the numerical parameters on the whole simulation. Thus they are critical for the quality of the simulation results and further investigations. Figure 7a It was already stated that the adaption of the numerical parameters has to be done according to the characteristics of the contacting surfaces. With change of this characteristics e.g. by using different materials, surface roughness or different surface treatments it is necessary to adapt the numerical parameters recursively. Thus the quality of the simulation results is not only dependent from the numerical modelling methods but also highly influenced by the calibration of the contact resistance parameters.
All in all, a separated treatment of the contact resistance calculation at each contact, as well as the implementation of time dependent, numerical parameters allows the reproduction of the experiments concerning the thermal-electric and structural aspects. By interpreting the adaption of the numerical parameters, future investigations might show further insights into contact resistance effects not implemented by physical parameters in Equation (2).

Conclusions
Despite the discontinuous geometry of the projection an indirect coupled model, allowing stable and efficient simulation of the CDW-process up to 671 K was built. With a reasonable adaption of the numerical parameters to the experimental measured voltages the simulation time was extended from 0.57 ms to 1.67 ms. Furthermore, the indirect coupling allows the dynamic calculation of the contact resistance and therefore a suitable way of determining the numerical parameters. A set of parameters was found, which yields to voltage drops at the contacts as measured in experiments. Therefore, it is possible to consider effects which can't be seen using constant factors, such as the evolution of surface conditions during the process. For the first time, this provides the possibility to adapt the simulation model to the complexity (large deformations, short current rise times and high temperature gradients) of the CDW. The presented model cannot describe the CDW process completely due to convergence problems. The convergence problems appear because of the high deformation at high temperatures. The finite elements deform heavily on the edge of the projection area. Additionally, the contact settings have been determined at room temperature, and therefore might be inadequate for high temperatures. Furthermore, rezoning was evaluated as a method to possibly overcome convergence issues within the structural physic [24].

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

Abbreviations
The following abbreviations are used in this manuscript:  Table A1. Thermophysical properties of EN-AW-6082.