Analytical approach to damage prediction in incremental sheet metal forming

Incremental sheet metal forming (ISF) is known to exhibit higher formability compared to conventional stamping. It is established that the mechanism of failure during ISF is by fracture occurring at higher effective strain than the local necking observed in traditional forming processes. The deformation limit in ISF is therefore estimated numerically using a suitable continuum damage model. However, simulation of incremental forming process including the damage model is often computationally expensive. In the present work, a simpler alternative combining finite element simulation and numerical solution outside finite element software is utilized. The finite element simulation is performed ignoring the damage model. The stress and strain history thus obtained in critical locations is utilized to estimate the damage evolution. The proposed method is useful when utilizing uncoupled continuum damage models. The proposed method is validated for typical cases by comparing it against the predictions from finite element method. Excellent correlation is observed between the proposed method and finite element simulation results.


Introduction
Incremental sheet metal forming (ISF) is a highly versatile technique aimed at dieless sheet metal forming. It finds relevance in many industries, viz. automotive, aerospace, biomedical sectors for manufacturing complex sheet metal components [1,2]. In ISF, the maximum strain region is dynamic due to the moving tool that postpone the localization of strains allowing larger overall deformation before complete fracture.
Failure during ISF is a combined effect of material behaviour and the stress states to which they are subjected to [3]. In ISF the deforming specimen endures a multiaxial state of stress. Depending upon the geometry of the component being formed the loading conditions and strain path (mode) vary at different locations within the material. Based on experiments and finite element (FE) simulations, it has been observed that the plastic deformation in ISF is limited to a narrow zone in the vicinity of tool-sheet contact area [4]. This can be understood by tracing the strain evolution in the contact region. Since the failure mechanism during ISF is by fracture rather than local necking, fracture forming limit curve (FFL) are used to represent deformation limits as opposed to the conventional forming limit curve (FLC) [5,6]. Although the overall formability of ISF is higher than the conventional forming, the process is limited by gradual IOP Publishing doi: 10.1088/1757-899X/1238/1/012024 2 thinning and subsequent fracture [7,8]. An investigation of the failed parts confirms that the cracks in the part manufactured by ISF were ductile in nature [9]. It can be shown that typical FFL can be fit as a straight line (in the range of uniaxial to equibiaxial tension deformation mode as applicable to FLC) with a constant slope in the traditional two dimensional principal strain space [7].
An accurate prediction of fracture during ISF is essential for successful process planning. From the mechanics of materials perspective (based on stress strain distribution), continuum damage models can be effectively employed to obtain a theoretical estimate for the critical fracture limits [10]. The constitutive description of damage modelling is broadly classified as coupled and uncoupled models [11]. The damage variable evolves and mutually influences the plasticity relations in the former approach. In the case of uncoupled models such as modified Mohr-Coulomb [12] or Hosford Coulomb [13], the damage evolution is independent of the stressstrain relation and the onset of failure is estimated when the accumulated damage variable reaches a critical value. The uncoupled damage models are relatively simpler to implement numerically and are increasingly gaining interest in successfully modelling fracture prediction in ISF [14]. Huang et al. [10] provided an estimation of approximate stress distribution in the ISF deformation zone using force equilibrium. Martins et al. [15] estimated the damage in the material using a fracture criterion based on hydrostatic stress. Nguyen et al. [16] utilized Oyane fracture criterion to estimate damage in ISF for magnesium. Malhotra et al. [17] used damage approach in FE analysis of a simple cone. It was observed that the damage accumulates the most in the regions of local bending in the sheet during deformation. Similar results were shown by the study conducted by Shamsari et al. [18] using Modified Mohr-Coulomb fracture criterion. The numerical simulations indicated that significant damage accumulates in the transition region between the deformed and contact zone.
Determining the location of the damage initiation is critical and may help to analyze key factors leading upto material fracture. Although a great amount of research is reported on the deformation mechanism of ISF, initiation and evolution of the material fracture is scarcely discussed. Implementation of damage in FE simulations involves meticulously written user subroutines, which utilizes heavy computational resources. As mentioned earlier, damage accumulation in uncoupled approach is computed independent of the deformation behaviour of the material. This presents an opportunity to reduce the overall computational load by separating the FE analysis from the damage calculations. In this context, the current work presents a preliminary framework where the FE simulations are used to determine the deformation history, thereby damage progression is estimated via a standalone analytical code using the FE output. This procedure is validated using selected ISF formed components.

Methodology
Commercially pure aluminium alloy (AA1050) sheet of 1.2 mm thickness is considered in this study. A series of mechanical tests were conducted to record the mechanical behaviour of the material till fracture. Material response under different loading conditions was achieved by varying the specimen geometry (details are discussed in section 3). All the mechanical tests were performed using a 100 kN capacity (Zwick-Roell 100) universal testing machine (UTM) at an initial strain rate of 0.005s −1 and the displacements were recorded by video extensometer.

Experiments on single point incremental forming (SPIF)
Components with three geometrical variations namely, (i) variable wall angle conical frustum (VWACF), (ii) truncated pyramid and (iii) a hybrid shape with five lobes were designed and developed by incremental forming in the present work ( Fig.1). SPIF experiments were conducted on CNC machining center ('Jyoti Huron KMill'). The sheet metal blank of dimensions (150 mm × 150 mm × 1.2 mm) was rigidly held by a customized fixture, such that the material flow is IOP Publishing doi:10.1088/1757-899X/1238/1/012024 3 constrained along the periphery to avoid any draw-in during forming. Tool of 8 mm diameter was moved across the sheet with a constant feedrate of 750 mm/min. 3D spiral toolpath conformal to the designed component geometry were developed using commercially available CAM software, 'Mastercam' (ver. 19.0.7874.0). Hydraulic oil was applied to the forming surface to reduce the contact friction during forming. In order to characterize the forming limits by fracture all the component(ϵ 2f ) at the vicinity of the crack were estimated using circle grid analysis. Thickness strains (ϵ 3f ) are estimated by measuring final thickness at the crack location. Details of the procedure can be found in work by Cristino et al. [19]. Laser etched samples with grid of 1.5mm (d) diameter circles were used, and the shape change of circle to elliptical indicated the nature of strain induced.
where d 2 is the length of minor axis of the ellipses after deformation. The strain measurements are obtained using 3Com USB camera. The major inplane strains ϵ 1f at fracture were obtained by imposing volume constancy of plastic deformation.

Finite Element Simulations of tension tests and ISF
Commercially available explicit finite element code ABAQUS ® was used to simulate the mechanical tests and SPIF. 3D solid elements C3D8R (8 node linear Brick with hourglass control) were chosen with minimum mesh size of 0.2 mm and 0.5 mm for mechanical tests and SPIF respectively. Mesh is refined in critical regions with a minimum of 3 elements along the thickness direction. The boundary conditions for UT, shear and notch tests and SPIF are depicted in Fig.2.
A user-defined VUMAT subroutine was implemented to simulate the mechanical behaviour. The plastic deformation was modelled assuming vonMises yield criterion and isotropic hardening. An uncoupled phenomenological damage criterion (Hosford-Coulomb) is implemented to simulate the damage evolution.

Proposed methodology for computation efficiency
It will be subsequently shown that the use of VUMAT increases the computation time significantly even for relatively simpler shapes. Unlike coupled damage models, the uncoupled damage models such as that being investigated in the present work does not interfere with the plastic constitutive behaviour. Instead, the damage evolution is independently evaluated for each element based on the deformation history. Taking advantage of the uncoupled formulation, a combined approach is proposed in the present work to reduce the computation time significantly. The proposed methodology involves finite element simulation without the use of any damage model. The stress and strain history of the elements are exported and the damage evolution is estimated outside the framework of finite element code. In the present work, we have used commercial MATLAB code for estimating damage. The damage model of interest is independently coded in MATLAB and using the strain and stress history, the damage variable is computed. The deformation limit when the damage variable reaches a critical value is considered as failure. The procedure is shown in the flow chart below (Fig.3). Since damage is evaluated independently, the in-built material models in commercial software that are numerically efficient can be used to simulate ISF. In addition to the above, different combinations of yield criterion, hardening and damage models can be evaluated easily. In the previous approach, a separate VUMAT has to be written for each combination and is time consuming. One of the factors that contribute to increased computation time in VUMAT is the damage evaluation of all the elements in the model. In practice, only select few regions are subjected to large strain. Therefore, it is unnecessary to evaluate the entire model. In the proposed method, a minimum number of elements from the highly strained region is chosen and only these selected elements are tracked for damage accumulation.

Results and Discussion
The uniaxial stress-strain flow behaviour is described using a combined hardening law (Eq.3), as depicted in the Fig.(4). A linear combination of Swift and Voce hardening rules are used to describe the pre and post necking deformation of the material till the onset of fracture [13,20].
Here K, ϵ o ,n,ϵ p and σ y denote strength coefficient, strain offset, strain hardening coefficient, equivalent plastic strain and yield stress respectively. A and B are fitting constants and Q is a post necking hardening parameter. Table 2 lists the mechanical properties of the material.

Calibration of damage model
The stress state at a material point can be described in terms of dimensionless ratios of stress invariant, namely stress triaxility (η) and Lode angle parameter (θ). In terms of principal stress components(σ 1 ,σ 2 ,σ 3 ), vonMises effective stress (σ), η andθ can be defined as follows :   In the present work, an uncoupled phenomenological fracture initiation model, Hosford Coulomb fracture criterion (HC) is used to predict the onset of ductile fracture in ISF. It is based on a localization criteria defined in terms of Hosford equivalent stress and normal stress on the plane of maximum shear, Eq(7) ( [13]).
HC model is a stress based criterion which is transformed from principal stress space, using Lode angle parameter dependent functions (f 1 ,f 2 ,f 3 ) as given below, The principal stress may be expressed in (η,θ,ε) space, as given : Using Eqs. (11)(12)(13), in the Eq.(7) the stress at the onset of fracture (σ f ) can be formulated as: Parameter a controls the Lode angle sensitivity of the fracture strain, c is called the friction parameter and b is the strain to onset of fracture under uniaxial tension. Furthermore, strain at onset of fracture (ε f ) is obtained by taking inverse of isotropic hardening law (Eq.3), which can be expressed as,ε f = f −1 (σ f ). Here, f signifies hardening function rewritten in terms ofσ f ."

Figure 6. Specimen geometries a) Notch Test (NT) b) Shear Test (SH) c) Uniaxial Tension (UT)
Strain data at the onset of fracture for three monotonic tests namely uniaxial tension (UT), shear test (ST) and notch test (NT) are used to calibrate the damage model. The loading history denoted by (η i ,θ i ) of each calibration test is evaluated with help of FEA simulation. It is assumed that the displacement at the onset of fracture (δ f ) coincides with the location of highest equivalent plastic strain within the specimen near fracture Fig(7). The damage progression (D) for each loading case is calculated according to the Eq. (15): D is damage indicator, which attains unity at the onset of fracture. The unknown parameters of HC damage model obtained by fitting fracture surface to the experimental fracture strains for UT,SH NT tests is given in Table 3. The above formulation is coded in MATLAB environment for calibration and offline prediction of fracture strains for a given forming condition.

Validation of numerical code
In order to test the validity of the code, the fracture strains obtained by FE simulations are compared with analytically estimated strain values for the mechanical tests. Fig.(9) shows the fracture strain distribution obtained by FE simulations and the analytical results for each mechanical test.  Figure 9. Fracture strains and evolution for damage estimated for a) Uniaxial tension test (UT) b) Notch test (NT) and c) Shear test (SH) using MATLAB code Furthermore, to test the applicability of the analytical method to predict the onset of failure strains in ISF, strain path data (η,θ and ϵ p ) were extracted for the critical elements from FE simulations and used as the inputs for the MATLAB code to estimate the fracture stains. It can be seen from Table 4 that the fracture strains under various loading conditions are well predicted by the coded damage theory in MATLAB.  Table 5, lists the CPU time recorded for FE analysis of ISF, when damage calculations are implemented with material subroutine (VUMAT), compared to analyses when damage is computed outside of FE environment. A workstation running on a windows platform, Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz 32 GB RAM is used to run all the FE simulations. A significant difference in computational time can be noticed in both cases. Fig.(11) depicts the inplane major and minor strain distribution, and the comparative placement of FFLs along with the evolution of strain path for the critical elements at fracture obtained by experiments, numerical simulations and analytical code. The experimental fracture forming limit is marked by interpolating a straight line from the values of inplane strains obtained by grid circle analysis. Similarly, fracture limit lines (FFLs) developed based on numerical prediction and MATLAB code represent the onset of fracture stains. The analytical and numerically estimated fracture strains are depicted by violet and yellow colored diamond marks respectively. Dashed lines marks the average strain path ( strain ratio, β), recorded while   Figure 11. Fracture forming limits forming three ISF parts i,e. 5 lobe, VWACF and pyramid. Furthermore, the damage evolution after onset is implemented by extrapolating the strain paths to predict the inplane strains, as function of average strain ratio corresponding to the loading, given as : here e 1 and e 2 are the major and minor inplane strains.

Conclusions
The present work demonstrated the analysis of incremental forming using uncoupled continuum damage model. The trend of fracture initiation predicted and experiments match well, though the simulation estimate were conservative. In order to reduce the computation time, we have proposed a hybrid strategy to separate the FE simulation without damage and later use the stress-strain history to analytically estimate the failure. It is demonstrated that the computation time reduced approximately by 50% when FE simulation is carried without damage model. The results were consistent with three different geometries analyzed for incremental forming. The deformation history of critical elements is analyzed using a MATLAB code after FE simulation. The predictions in the proposed approach is compared with that predicted using a user defined material model. It was shown that the results correlate very well. Therefore, the presented methodology can be utilized to significantly reduce the computation time. The added advantage is the flexibility of analyzing different combination of constitutive relations without detailed development of individual subroutines.

Acknowledgments
Authors would like to acknowledge Department of Heavy Industries, Government of India and Mahindra & Mahindra for funding the project (sanction number: 7(8)/2019-AEI (19310)).