Numerical simulation of self-piercing riveting process ( SRP ) using continuum damage mechanics modelling

The extended Bonora damage model was used to investigate joinability of materials in self-piercing riveting process. This updated model formulation accounts for void nucleation and growth process and shearcontrolled damage which is critical for shear fracture sensitive materials. Potential joint configurations with dissimilar materials have been investigated computationally. In particular the possible combination of DP600 steel, which is widely used in the automotive industry, with AL2024-T351, which is known to show shear fracture sensitivity, and oxygen-free pure copper, which is known to fail by void nucleation and growth, have been investigated. Preliminary numerical simulation results indicate that the damage modelling is capable to discriminate potential criticalities occurring in the SPR joining process opening the possibility for process parameters optimization and screening of candidate materials for optimum joint.


INTRODUCTION
elf-piercing riveting (SPR) is a cold-forming operation used to fasten two or more sheets of material by driving a semi-tubular rivet through the top sheet(s), piercing the bottom sheet, and spreading the rivet skirt under the guidance of a suitable die.As the process relies on a mechanical interlock rather than fusion, it can be used for a wide range of advanced materials that are dissimilar, coated, and hard to weld.Unlike the "traditional" riveting process, no pre-process, such as hole drilling and alignment, is necessary.The complex joint geometry and its three-dimensional nature combine to increase the difficulty of obtaining an overall system of governing equations for predicting the properties of the SPR joints.The effective way to analyze SPR joint during forming process is to perform numerical simulation.This has been addressed by several authors using different numerical techniques.Mori et al. [1] used finite element simulation to optimize the die shape provided the mechanical properties of the rivet.They found that the strength of the joint is influenced by both mechanical properties of the sheet metals and the ratio of the lower sheet thickness and the overall thickness.Atzeni et al. [2] investigated the possibility to predict the strength of a SRP joint simulating both process and shear test finding a reasonable agreement with experiments in terms of deformed geometry and force-displacement response.More recently, He et al. [3] simulated SPR S process in AL 5754 using LS-DYNA showing the possibility to reproduce with good accuracy the load vs displacement response during the process.However, finite element analysis can be particularly useful in providing a better understanding of the potential mechanisms and conditions that may result in a poor joining.In fact, several defects may occur during SRP.These include: the penetration through the lower sheet, the necking of the lower sheet, and the separation of sheets.In the penetration through the lower sheet, the leg of the rivet goes out from the bottom surface of the lower.This is found to be a driver for corrosion in service.The necking of the lower sheet is a thinning phenomenon, similar to the penetration defect.In the separation of the sheets, the spread of the leg of the rivet is too small to join the sheets, and the riveted sheets are easily separated [5].A sketch of these defects is shown in Fig. 2. The SPR technique is mainly applied to joining of two sheets, although application of this technique to multi-layer sheets has also been investigated [6].The possibility to achieve a joint without defects defines the "joinability" of candidate materials.At present, for given materials (i.e.steel and aluminum), this ability seems to depends on the ratio between the upper and the lower sheets thickness.For steel-aluminum it was found that penetration, necking and separation are caused by small total thickness, small thickness of the lower sheet and large total thickness, respectively [5].The possibility to predict the occurrence of such defects by means of numerical simulation relies mainly in the ability to model plastic deformation in the materials and the capability to predict the occurrence of fracture phenomena.This requires an accurate description of material flow curve and a failure model.In SPR of joinable materials, the upper sheet material initially flows around the rivet with a progressive reduction of the ligament eventually followed by fracture.In the literature, this process has been simulated using arbitrary separation criteria based on geometrical consideration.Cacko [7] presented a review of selected material separation criteria available in commercial MSC software applied for the SPR process simulation.His conclusions confirm that geometrical separation criteria, although very simple to use, relays on experiments for calibration.This makes simulation suitable for reproducing a known process but incapable of any prediction for different scenarios in which process parameters or material characteristics are changed or unknown.Alternatively, the conditions for separation, occurring in the material under large plastic deformation as such in SPR process can be accurately predicted using continuum damage mechanics (CDM).In this field, there is an extensive literature on modelling ductile damage in metals and alloys.Among available approaches to model ductile damage and its effects, Bonora [8], in the framework of CDM introduced by Lemaitre [9], proposed a damage model formulation in which a general non-linear law of evolution for damage that was demonstrated to be suitable to predict ductile fracture in different classes of metals and alloys under variable geometry and loading conditions [10][11][12][13][14][15].Very recently, Testa et al. [16] extended the original Bonora damage model in order to account for the effect of the third invariant of the deviatoric stress tensor on the fracture strain of shear fracture sensitive materials.This new feature is particularly suited to investigate fracture in deformation processes largely dominated by shear such as those occurring in the self-piercing process.In this work, the extended Bonora damage model (XBDM) was used to investigate the joinability of materials in the SPR process.
To this purpose, an extensive numerical investigation was performed analyzing different combinations of material sheets including also materials that are known to show susceptibility to shear controlled fracture.Results seem to indicate that the model is capable to predict the occurrence of damage during the process and in particular before the rivet flaring stage.The possibility to predict the poor joinability of materials is also shown.

DAMAGE MODELLING
he Bonora damage model (BDM) is derived according to the thermodynamics framework of continuum damage mechanics (CDM) initially introduced by Lemaitre [17].In this context, from the definition of the potential of dissipation, the evolution law for the damage, which is a state variable, can be derived.Under the assumption of isotropic damage, D is a scalar defined as the ratio between the damaged and the nominal net resisting area of the material reference volume element (RVE), D ranges from 0, for the material in the undamaged state, to cr D at rupture where the material has zero load carrying capability.Testa et al. [16] proposed the following expression for the damage dissipation potential, Here, Y is the damage associated variable, which represents the strain energy density release rate, derived from the state potential [18], where, E is the Young's modulus and R  is the function that accounts for stress triaxiality effects, where  is the Poisson's ratio and  is the stress triaxiality factor defined as the ratio of the means stress m  and equivalent von Mises stress  , The first term on the right hand side is the original expression of damage dissipation potential proposed by Bonora [8].This accounts mainly for the dissipation associated with damage state ascribable to nucleation and growth of microvoids (NAG) and it is led by stress triaxiality.Here, p is the "active plastic strain" defined as the equivalent plastic strain that accumulates under tensile state of stress (T≥0), is the Heaviside function that is equal to 1 when the stress triaxiality is positive and 0 otherwise.The unilateral condition for damage states that under compressive state of stress, damage does not accumulate and its effects are temporarily restored ( ).The second term on the right hand side of eqn.( 2) accounts for dissipation associated with void sheeting under shear dominated stress state.Here,  is the parameters, bounded between 0 and 1, that accounts for the for the influence of the third invariant of the deviatoric stress tensor J 3 on material ductility, which is defined as a function of the Lode parameter L as, The damage evolution law is obtained from the generalized normality rule where   is the plastic multiplier equal to the equivalent plastic strain rate scaled by damage effect,   For constant  and T (>0) deformation process, the damage rate equation can be integrated to obtain the following expression for damage evolution, Here, , D cr ,  th ,  f ,  f , , k are the material damage parameters:  th is the plastic strain threshold under uniaxial state of stress at which damage process is initiated, f is the failure strain under constant stress triaxiality equal to 1/3,  is the damage exponent that defines the shape of damage evolution law for NAG,  and k are the damage exponents for shear controlled damage contribution, and f is the critical strain for pure shear.These parameters can be determined performing selected experimental tests.In particular,  th and  f can be determined by fitting on the stress triaxiality vs failure strain plot of fracture data obtained with round notched bar in tension for which w=0, using the expression for the failure locus Similarly,  f and k can be determined on the same plot fitting fracture data in the negative stress triaxiality regime with the following expression of the fracture locus, Other damage parameters, which define the shape of the evolution of different damage contribution with plastic strain, can be determined performing measurement of loss of stiffness as discussed in [19].In a first approximation, =0.3 is a suitable value for several classes of materials.More difficult is the determination of .Barsoum and Faleskog [8] found that at low levels of stress triaxiality, failure seems to be governed by a simple criterion based on a critical measure of plastic shear deformation.Therefore, following the Occam's razor principle, =0 is assumed.

Geometry
he self-piercing riveting process was simulated with FEM using XBDM model.All finite element simulations have been performed with MSC MARC v2017.1.In this FEM code, the BDM is ready available, however the XBDM was implemented using user subroutines.Here, a self-pierce H4 C-SKR rivet was selected to simulate joining of T different combinations of hard to soft metal sheets.The reference self-piercing riveting configuration is shown in Fig. 4.
Here, the holder and the punch are simulated as rigid bodies while deformable-deformable contact is considered between the rivet, the upper, and the lower metal sheets.Since the problem is axisymmetric, all simulations have been performed using 2D four node axisymmetric elements.This element type has four Gauss points, bilinear shape functions and it is particularly suited for large strain and contact problems.Analyses have been carried out using large displacement, finite strain and Lagrangian updating.Since elements in the metal sheets undergo large plastic deformation, automatic remeshing was used to avoid severe element distortion and analysis interruption due to the impossibility to reach convergence.Remeshing was performed using advanced front quadrilateral prescribing the average (0.04 mm) and minimum element size (0.01 mm) for the new created elements.In order to investigate the modelling capabilities to predict different material joinability, the combinations of materials given in Tab. 1, have been investigated.The upper and lower sheet thicknesses are 27 mm and 15 mm, respectively.These are unchanged during the analyses.

Materials
Materials considered in this study are DP600, AL 2024-T351 and OFHC 99.98% pure copper.The DP600 is a dual phase steel consisting of a ferrite matrix containing a hard second phase, usually islands of martensite, widely used in the automotive industry for different parts and safety cage components (B-pillar, floor panel tunnel, engine cradle, front subframe package tray, shotgun, seat).Compared with high strength low alloy steels (HSLA), DP steel exhibits higher initial work hardening rate, higher ultimate tensile strength, and higher ultimate stress over yield stress ratio than a similar yield strength HSLA.AL 2024-T351, usually provided as plate or sheet, is ideal for everyday applications where a high strength to weight ratio is necessary.The material offers a fair level of workability as well as high level of corrosion resistance.This alloy is used in a variety of applications such as fuselage structurals, wing tension members, shear webs and ribs, and structural areas where stiffness, fatigue performance, and good strength are required.This material has been widely characterized.It susceptibility to shear controlled fracture was investigated in [20] while Testa et al. [16] identified the XBDM model parameters.OFHC 99.8% Cu has been characterized extensively for application in the dynamic deformation range.[21][22][23][24][25][26][27].Bonora et al. [28] determined the BDM model parameters for OFHC CU in the fully annealed state.This material was selected because fracture process is controlled by NAG and does not show shear sensitivity.Flow curves for these material at room temperature and low strain rate (10 -3 /s) is given in Fig. 5.

DP600 damage parameters identification
Damage parameters for AL 2024 and OFHC were identified elsewhere.For DP600, the experimental data given in [29] where used for the identification of damage parameters.Authors performed tensile tests on several specimen geometries to obtain different combinations of stress triaxiality and Lode parameter.Among all available experimental data, those relative to geometries for which w=0, i.e. axisymmetric round notched bar (RNB), have been selected to determine the damage parameters for DT.Damage parameters for D, would requires experimental data in the regime for T<0.
Unfortunately, no such data where available.Therefore, k could not be fitted and was assumed equal to 0.1 which sounds reasonable for the few data reported in the near zero stress triaxiality range.The choice of this value was verified simulating with FEM some of the flat notched samples tested by Basaran and Weichert [29] that in Fig. 6 showed a failure strain in the intermediate stress triaxiality range.In Fig. 6 the comparison of predicted failure locus for shear and stress triaxiality controlled stress states with available experimental data is plotted.Here, the red line is the failure locus given by eqn.(12) where there is no contribution of the third invariant of deviatoric stress (w=0) while the blue line is the failure locus predicted for shear controlled fracture (=1) without the contribution of damage due to stress triaxiality under the assumption of plane stress condition ( 3 =0).
Light symbols are relative to specimen geometries in which both T and  contribute.In Fig. 6 the traced failure loci have to be interpreted as follow.For the stress triaxiality range where only a solution is drawn, this represents the locus where fracture strain is expected to occur.In the range of stress triaxiality where the two solutions superimpose, fracture strain will be determined by the concurrent action of stress triaxiality and Lode parameter.Therefore, the expected fracture strain in this range is much lower than the values retuned by the two limit solutions.In Tab. 2, the summary of damage parameters for DP600 is given.

RESULTS AND DISCUSSION
n Fig. 7 the calculated global response in terms of applied punch force vs vertical displacement for different combination of materials sheets in the joint is given.These results are in a general good agreement with the response measured in experiments.For instance, for the DP600-OFHC combination, all the essential phases of the process can be recognized: 1) clamping, 2) piercing, 3) rivet flaring and 4) compression.When less joinable materials are used, some of the phases no longer occur.For instance, for AL2024-AL2420 configuration, the shear sensitivity of the material causes the aluminum sheet to fracture when in contact with the rivet.In this case full penetration occurs without rivet flaring.On the opposite side, if the material in contact with the rivet is ductile but with lower strength, piercing occurs during clamping and the compression phase starts well before the rivet has begun to flare.In Fig. 8 the comparison of the calculated damage contour map for the different material combination for a punch vertical displacement of 1.25 mm is shown.Here, it can be noted that for shear sensitive materials (AL-AL), damage controlled by Lode parameter accumulates almost immediately resulting in the piercing at the contact without appreciable deformation of the metal sheets.Opposite situation is predicted for the OFHC-DP600 configuration.Here the ductility of the upper sheet combined with that of DP steel, result in much larger deformation in the die with a limited indentation of the rivet.This indentation is larger for the DP600-AL configuration.In fact, the limited deformability of the AL restraints that of the upper steel plate promoting the rivet penetration at this stage.
Increasing the punch displacement (=2.26mm), the rivet is predicted to penetrate without any relevant deformation in the AL-AL configuration, while in the DP600-OFHC the lower sheet has touched the die and the complete penetration of the upper sheet has not occurred yet Fig. 9.It is interesting to note the difference in the deformed shape of the upper sheet material that flow into the rivet for the OFHC-DP600 and DP600-OFCH.In the latter case, shear bands (damage contours) initiated from the rivet edge and developing along the upper metal sheet are clearly visible while in OFHC-DP600 configuration, damage is much more localized along the contact surface between the copper and the rivet.Interestingly, in the case of DP600-AL, damage start to accumulate on the back surface of the AL sheet.This damage state is mainly controlled by stress triaxiality and it is due to the plastic strain accumulation due to the local bending.
Increasing further the punch displacement, Fig. 10 and Fig. 11, fracture is predicted to occur in AL-AL on both metal sheets with almost no deformation of the rivet.Similar fracture is predicted to occur for the DP600-AL configuration.
The main difference is that in the AL-AL damage is driven by shear in the DP600-AL manly by NAG.Complete clutching is predicted to occur for the DP600-OFHC, although extensive damage bands along the ligament of the lower metal sheet is predicted to develop.This model formulation allows accounting also for shear-controlled fracture in materials that suffer shear fracture sensitivity.SRP process, with several combinations of materials, has been investigated to demonstrate that the model is capable to predict different piercing characteristics that could not be predicted with simple separation criteria.Numerical simulation results show that it is possible to anticipate the behavior of the material during the joining process allowing to tune process parameters and to screen materials in a virtual environment.Of particular interest is the case of AL2024-T351.This material showed to result in uncomplete joint when coupled with another AL or DP steel sheet but for completely different mechanisms.Further experimental investigation may provide indications to support model validation.

Figure 1 :
Figure 1: Schematic representation of SRP fastening process: a) the semi-tubular rivet is place in the holder; b) applied load forces the rivet to deform with the metal sheets that eventually pierce, c) usually, only the first sheet is pierced and while the rivet lock-in the second metal sheet, d) the load is finally released [4].

Figure 5 :
Figure 5: True stress-strain curves for materials used in the simulation.

Figure 6 :
Figure 6: Comparison of predicted failure locus for shear and NAG dominated regimes and experimental data for DP600.

Figure 7 :
Figure 7: Comparison of punch force vs displacement curves for different combinations of materials.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figure 8: Comparison of deformed state for different material combinations at the end of the clamping stage, punch displacement: 1.25 mm.

Table 1 :
Summary of investigated SPR configurations.

Table 2 :
Summary of damage model parameters for DP600.