Prediction of fracture limits of IN625 alloy by coupling ductile damage models and anisotropic yielding functions with the classical Marciniak and Kuczyński model

Abstract

The fracture forming limit diagram (FFLD) is gaining special attention in high strength materials where the necking tendency rarely occurs during sheet metal forming processes. In the present work, the classical Marciniak and Kuczyński (MK) model has been modified by coupling it with different ductile damage models (Cockcroft and Latham, Brozzo, Oyane, Ko, Oh, Rice and Tracey, McClintock and Clift) and anisotropic yielding functions (Hill 1948 and Barlat 1989) to predict the fracture limits of Inconel 625 (IN625) alloy at different temperatures. Firstly, uniaxial tensile testing has been conducted for the determination of important mechanical properties. Consequently, stretch forming experiments have been performed to analyze the forming limits of a material. It has been found that the safe and fracture forming limits of the material increased by approximately 17.26% and 22.22%, respectively, on increasing the temperature from 300 to 673 K. From the comparative analysis of different combinations of ductile damage models and yielding functions, the Cockcroft and Latham (C-L) damage model in combination with the Barlat 1989 yielding function helped in best predicting the theoretical FFLD as it displayed the least average root mean square error (RMSE) of 0.033. The other ductile damage models used for predicting the theoretical fracture limits displayed large error; hence, they should not be considered while designing a critical component in the manufacturing industry using IN625 alloy.

Introduction

Inconel 625 (IN625) alloy is a precipitation-hardening alloy with major chemical contents of Ni-Cr-Fe along with some secondary precipitation in the metal matrix of Ni-Al/Ti/Nb. These precipitates in the metal matrix help achieve the excellent blend of different mechanical properties even at higher processing temperatures. Due to the excellent performance of IN625 alloy in extreme conditions, it has been used immensely in aerospace and nuclear applications. ¹ The sheet metal forming process is one of the major processes used in the manufacturing industry for making several complex components. The workability or the maximum deformation during any sheet metal forming process is often affected and limited by a failure occurring in the sheet metals. The necking phenomenon usually occurs due to the sheet’s localized thinning during the excessive deformation that occurs while forming complex shape components. The necking phenomenon also leads to the loss of the load bearing capacity of sheet metal.

The forming of a high strength material such as IN625 alloy is difficult at room temperature conditions. ² In the past, several forming techniques, such as hydro-forming,^3,4 electro-magnetic forming, ⁵ incremental forming,^6,7 hot and warm forming, ⁸ deep drawing, ⁹ etc., have been used to improve the formability of the material. Hydro-forming and electro-magnetic forming techniques help in improving the quality and formability by many folds. Components with highly complex shapes can be easily produced using these forming methods, but they have some disadvantages, such as high equipment cost and poor production rate. Hot forming is usually carried out at a temperature higher than the recrystallization temperature of the material. It is sometimes less preferred to form components of complex shapes, as instability in the material might arise due to the recrystallization, recovery, and grain growth. However, warm forming helps in improving the safe and fracture forming limits of the material by increasing the ductility of the material. ¹⁰ It is also simple and economic compared to other advanced forming techniques.

The forming limit diagram (FLD) usually describes the onset or the occurrence of localized necking. Goodwin ¹¹ and Keeler ¹² first established an experimental way of determining the FLD. The FLD is used extensively to evaluate the formability of different sheet metals. However, the determination of an experimental FLD involves extensive resource consumption and also requires an expert researcher to draw a trustworthy limiting curve. ¹³ Thus, various researchers in the past have proposed different theoretical models to determine the limiting curve in order to reduce rigorous experimental processes. By considering that sheet metals possess homogenous properties, the first such models for the theoretical prediction of localized necking and diffused necking were proposed by Hill and Swift respectively. However, several discrepancies have been found in the forming limits between the experimental results and these proposed models for the case of having the ratio of minor to major strains greater than zero. Further, the Marciniak and Kuczyński (MK) ¹⁴ proposed one of the most commonly used theoretical forming limit prediction model by having an assumption that the failure occurred due to the presence of non-homogeneity (structural and geometric) present in the sheet metals. The classically proposed MK model was later improved by Hutchinson and Neale^15,16 incorporating the influence of groove rotation angle. It is noteworthy to mention that the accuracy of the necking limit estimated by the MK model can be influenced by the choice of yield models, hardening equations, initial geometric imperfection factor, and strain rate sensitivity index etc. as demonstrated in the previous literature.^17–19 The MK model is mainly based on thickness variation and necking theory for FLD prediction. The MK model is unable to predict failure when there is a negligible change in the thickness of the considered metal, especially in the negative or low stress triaxiality region. Marciniak et al. ²⁰ modified the proposed MK model by using the Hill 1948 yield criterion and found very credible forming limits for aluminum, steel, and copper. Banabic et al. ²¹ used five different methods while predicting the forming limits, namely, Swift’s diffused necking, Hill’s localized necking, the finite element approach, the modified maximum force criterion, and the MK model for comparing the predictive efficiency of these methods for AA5182 metallic sheet. The normal anisotropy of the material also influences the forming limits of the material. Ávila and Vieira ²² used von-Mises, Hill 1948, Hill 1979, Logan and Hosford and Hill 1993 yielding functions while predicting the forming limits using the MK model and they also found that the criterion used highly affects the predicted forming limits of the material.

Recently, several challenges have arisen in forming complex-shaped components due to the sudden appearance of fractures without any prior hint of substantial necking. ²³ On the basis of the opening mode of the crack tip, the fracture limits can be classified into the shear fracture forming and tensile fracture forming limit diagram (FFLD). The continuum damage (ductile) models in combination with the yield criterion used help in the prediction of the theoretical FFLD. The Bao-Wierzbicki (BW) ²⁴ model is one of the most popular damage models used for the evaluation of the fracture locus. Three different empirical relations have been proposed in the BW model to estimate the fracture limits in the whole stress triaxiality (−0.33 < η < 0.67) range. Different researchers in the past have used several damage models in the estimation of the FFLD. Isik et al. ²⁵ in their study over AA1050-H111 introduced normal anisotropy into the damage models for FFLD prediction using the Nakazima, bulge, and Incremental sheet metal forming experiments. Park at al. ²⁶ combined the Hill 1948 and Barlat 1991 yielding functions into the Lou-Huh ductile damage model for the prediction of an accurate FFLD for stainless steel. Pandre et al. ²⁷ used seven different ductile fracture criteria to calibrate the fracture limits of DP590 steel. They found that among all the models used, BW damage criteria can be used in the future for evaluating the fracture limits of DP590 steel. Takuda et al. ²⁸ used five different proposed criteria for investigating the failure of cylindrical deep-drawn cups using A5052, A1100, and A2024-T4 alloys. They further found these criteria to be accurate in predicting the HDR after implementing the ductile damage models. Prasad et al. ²⁹ also compared the abilities of six models to predict the fracture limits of the solution-treated IN718 alloy for aerospace applications. They found that Oh’s model accurately followed the experimentally evaluated fracture limits of the material with the least average absolute error (AAE) of approximately 10%. Wu et al. ³⁰ used six different damage models to study the failure of material using the hydroforming technique. The Rice-Tracey model was found to have results approximately the same as those obtained by the experimental approach. These ductile damage models have also been used to evaluate the fracture for different processes such as biaxial stretching, bore expansion, hydroforming, tube bending, and deep drawing. These fracture models can predict fracture initiation with high accuracy, as they perform evaluations by considering the deformation energy. However, these damage models have certain limitations in the calculation of material constants which restricts their application in different manufacturing processes. Further, it is also evident from the literature that very limited work has been done for the prediction of fracture limits using the damage models (considering anisotropy) for the warm forming of Inconel alloys.

In the present work, different combinations of ductile damage models (Cockcroft and Latham, Brozzo, Oyane, Ko, Oh, Rice and Tracey, McClintock and Clift) and yield criteria (Hill 1948 and Barlat 1989) have been incorporated into the classical M-K model for predicting the fracture limits of IN625 alloy at forming temperatures of 300 and 673 K. The stretch forming experiments have been performed to find experimental fracture limits for validating the theoretical results. A comparative analysis has also been performed to find the best combination of yield criteria and damage models for predicting the fracture limits using a modified M-K model which could be used in the manufacturing industry in order to avoid costly and time-consuming experimental procedures for finding the limiting strains and forming complex shapes.

Theoretical framework of modified Marciniak and Kuczyński (m-MK) Model

The fracture prediction in the present work has been carried out with the help of a modified Marciniak and Kuczyński (m-MK) model. ²³ The theoretical framework to obtain the fracture limits of the material using the proposed m-MK model has been described step by step below in Figure 1. The non-linear equations used in this theoretical model have been solved using the Newton-Raphson solver in the MATLAB software.

Figure 1.

Layout of the procedure to be followed for calculating the fracture limits using the proposed m-MK model. ²³

The classical MK model has been used in abundance to find the forming limits of a given material. When applying the classical MK model, inhomogeneity has been assumed in the material in terms of the groove. The groove region (region “B”) is the critical region; it plays a major role in the initiation of fracture, while region “A” is the safe part of the sheet metal, as shown in Figure 2. The ductile fracture criteria have been used for characterizing the condition of fractures in the stretch formed specimens. The principle strains in region A at the time of the fulfillment of the fracture conditions are considered to be the limiting strains of the material.

Figure 2.

Geometric imperfection assumed while solving MK model.

The classical MK model has been designed by assuming the plane stress condition in regions A and B. Hence, equations (1) and (2) hold true (for the indices 1, 2, and 3 in these equations, please refer to Figure 4).

$${\mathrm{\sigma }}_{3i}={\mathrm{\sigma }}_{i3}=0;\mathrm{i}=1,2,3$$

(1)

$${\epsilon }_{3i}={\epsilon }_{i3}=0;\mathrm{i}=1,2,3$$

(2)

While undergoing the deformation, the region A has been assumed to have gone under the proportional loading which further resulted in no tangential strains ( ${\epsilon }_{12}={\epsilon }_{21}=0$ ) and stresses ( ${\mathrm{\sigma }}_{12}={\mathrm{\sigma }}_{21}=0$ ) in the plastic orthotropic frame. The strains for each individual strain path have been updated in small increments and the limiting strain has been achieved on the initiation of fracture. The term [t·t + Δt] corresponds to the different time interval during the theoretical analysis. At the moment, t = 0, the principle strain ( ${\epsilon }_{1\mathrm{A}},{\epsilon }_{2\mathrm{A}},{\epsilon }_{1\mathrm{B}},\mathrm{and}{\epsilon }_{2\mathrm{B}}$ ) and the equivalent strains ( ${\overline{\epsilon }}_{\mathrm{A}}\mathrm{and}{\overline{\epsilon }}_{\mathrm{B}}$ ) are assumed to be zero. Further, the unknown quantities need to be evaluated at the time interval t + Δt.

The inhomogeneity factor or the initial thickness ratio f₀ is defined as f₀ = s_B0/s_A0; where 0 < f₀ < 1. Further, s_A0 and s_B0 are defined as the initial thickness in regions A and B respectively. Additionally, to have simplicity in solving the mathematical formulations, the deformation in region A is assumed to have occurred along linear strain paths. The slope of each path is defined as the principle stress ratio ( ${\mathrm{\beta }}_{\mathrm{A}}=\raisebox{1ex}{${\mathrm{\sigma }}_{2\mathrm{A}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\sigma }}_{1\mathrm{A}}$}\right.$ ), which lies between 0 and 1.

The computation of the fracture strains starts with a small increment in the principle strain ( $\mathrm{\Delta }{\epsilon }_{1A}$ ) of region A. The increment should be small, in the range of 10⁻⁴–10⁻³, in order to have a high accuracy. The equivalent stress in a load state can be expressed as $\overline{\mathrm{\sigma }}={\mathrm{\sigma }}_1.\mathrm{F}(\beta )$ where, ${\mathrm{\sigma }}_1>0$ . The function F is related to the adopted yield criteria. The Hill 1948 and Barlat 1989 yield criteria have been used in the present study to incorporate the anisotropic effect of sheet metal which is further elaborated below.

The Hill 1948 ³¹ yield criterion can be expressed using the mathematical function shown below in equation (3).

$${\overline{\mathrm{\sigma }}}_{\mathrm{Hill}}^2={\mathrm{F}\mathrm{\sigma }}_{\mathrm{yy}}^2+\mathrm{G}{\mathrm{\sigma }}_{\mathrm{xx}}^2+\mathrm{H}({\mathrm{\sigma }}_{\mathrm{xx}}-{\mathrm{\sigma }}_{\mathrm{yy}}{)}^2+2{\mathrm{N}\mathrm{\sigma }}_{\mathrm{xy}}^2$$

(3)

Here, F, G, H, and N are the material constants which describe the anisotropy of the material. Stoughton ³² suggested neglecting the component for shear stress (σ_xy) in order to reduce the complexity of the yield function. Further, in a plane of principle stress, the stress components ${\mathrm{\sigma }}_{\mathrm{xx}}$ and ${\mathrm{\sigma }}_{\mathrm{yy}}$ can be rewritten as ${\mathrm{\sigma }}_1$ and ${\mathrm{\sigma }}_2$ respectively. The associative flow law represented by equation (4) can be used for representing the strain ratio ( $\alpha =\raisebox{1ex}{${\epsilon }_2$}\!\left/ \!\raisebox{-1ex}{${\epsilon }_1$}\right.,\mathrm{where}{\epsilon }_1>{\epsilon }_2$ ) as a function of the stress ratio ( $\beta =\raisebox{1ex}{${\mathrm{\sigma }}_2$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\sigma }}_1$}\right.,\mathrm{where}{\mathrm{\sigma }}_1>{\mathrm{\sigma }}_2$ ) as per equation (5).

$$\mathrm{d}{\epsilon }_{\mathrm{ij}}=\mathrm{d}\mathrm{\lambda }\frac{\partial \varphi }{\partial {\mathrm{\sigma }}_{\mathrm{ij}}}$$

(4)

$$\alpha =\frac{(\mathrm{F}+\mathrm{H})\mathrm{\beta }+\mathrm{H}}{\mathrm{G}+\mathrm{H}-\mathrm{H}\mathrm{\beta }}$$

(5)

The parameter ϕ, derived from the Hill 1948 yielding function can be expressed as

$$\varphi =\frac{{\overline{\mathrm{\sigma }}}_{\mathrm{Hill}}}{{\mathrm{\sigma }}_1}=[(\mathrm{G}+\mathrm{H})+(\mathrm{F}+\mathrm{H}){\mathrm{\beta }}^2-2\mathrm{H}\mathrm{\beta }{]}^{0.5}$$

(6)

By assuming the plane stress condition ( ${\mathrm{\sigma }}_3=0$ ), the plastic work equation ( $\overline{\mathrm{\sigma }}d\overline{\epsilon }=\sum _{i=1}^3{\mathrm{\sigma }}_{i}d{\epsilon }_i$ ) can be used to derive the effective plastic strain ( $\overline{\epsilon }$ ) as shown below in equation (7).

$$\overline{\mathrm{\epsilon }}=({\mathrm{\epsilon }}_1(1+\alpha \mathrm{\beta }))/\varphi$$

(7)

The material constants shown in equation (3) can be evaluated using two different approaches, namely, the r-based and σ-based approach. In the present work, the r-based approach has been used. Hence, the anisotropic coefficients in the three directions, namely, the RD (Rolling direction) or 0°, DD (Diagonal direction w.r.t. RD) or 45°, and ND (Normal direction w.r.t. to RD) or 90°, have been evaluated. The formulation for the evaluation of all the constants is shown below.

$$\mathrm{G}\left(\mathrm{r}\right)=\left(\frac{1}{1+{\mathrm{r}}_0}\right)$$

(8)

$$\mathrm{H}\left(\mathrm{r}\right)=\left(\frac{{\mathrm{r}}_0}{1+{\mathrm{r}}_0}\right)$$

(9)

$$\mathrm{F}\left(\mathrm{r}\right)=\left(\frac{{\mathrm{r}}_0}{{\mathrm{r}}_{90}(1+{\mathrm{r}}_0)}\right)$$

(10)

$$\mathrm{N}\left(\mathrm{r}\right)=\left(\frac{({\mathrm{r}}_0+{\mathrm{r}}_{90})(1+2{\mathrm{r}}_{45})}{2{\mathrm{r}}_{90}(1+{\mathrm{r}}_0)}\right)$$

(11)

The stress triaxiality (η) using the Hill 1948 yield criterion can be evaluated using equation (12).

$$\eta =\frac{1}{3}\mathrm{x}\frac{1+\beta }{\sqrt{\frac{r_{90}+r_0{\beta }^2+r_0{r}_{90}{(1-\beta )}^2}{r_{90}(1+r_0)}}}$$

(12)

The Barlat 1989 ³³ yield criterion can be mathematically expressed using equation (13).

$${\overline{\sigma }}_{Barlat}={\left\{\frac{a}{2}\left[{\left({\sigma }_1-{\sigma }_3\right)}^M+{(-h)}^M{\left({\sigma }_2-{\sigma }_3\right)}^M\right]+\left(1-\frac{a}{2}\right){\left({\sigma }_1-{\sigma }_2-2{\sigma }_3\right)}^M\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.}$$

(13)

Here, σ₁, σ₂, and σ₃ are the principle stresses. M is an integer exponent with 6 as its value for the BCC material and 8 for the FCC material. IN625 is the FCC material; hence, the value of M is selected to be 8. The material constants a, c, and h can be calculated using equation (14)–(16).

$$\mathrm{a}=2(1-\sqrt{\frac{{\mathrm{r}}_0{\mathrm{r}}_{90}}{(1+{\mathrm{r}}_0)(1+{\mathrm{r}}_{90})}})$$

(14)

$$\mathrm{h}=\sqrt{\frac{{\mathrm{r}}_0(1+{\mathrm{r}}_{90})}{{\mathrm{r}}_{90}(1+{\mathrm{r}}_0)}}$$

(15)

$$\mathrm{c}=(2-\mathrm{a})=\sqrt{\frac{{\mathrm{r}}_0{\mathrm{r}}_{90}}{(1+{\mathrm{r}}_0)(1+{\mathrm{r}}_{90})}}$$

(16)

Under the plane stress condition, equation (13) can be transformed into equation (17).

$$\begin{array}{c}\multicolumn{1}{c}{{\overline{\mathrm{\sigma }}}_{\mathrm{Barlat}}=\{\frac{a}{2}[{\left({\mathrm{\sigma }}_1\right)}^M+(-h{)}^M{\left({\mathrm{\sigma }}_2\right)}^M]}\\ \multicolumn{1}{c}{+(1-\frac{a}{2})({\mathrm{\sigma }}_1-{\mathrm{\sigma }}_2^M{\}}^{1/M}}\end{array}$$

(17)

The parameter ξ, which is the ratio of equivalent to the major principle stress for the Barlat 1989 yielding function can be evaluated using equation (18).

$$\varphi =\frac{{\overline{\mathrm{\sigma }}}_{\mathrm{Barlat}}}{{\mathrm{\sigma }}_1}={\left[\frac{a}{2}+\frac{a}{2}{(-h\beta )}^M+\left(1-\frac{a}{2}\right){\left(1-\beta \right)}^M\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.}$$

(18)

Further, the strain ratio (α) using the Barlat 1989 yielding function can be evaluated as

$$\alpha =\frac{{\epsilon }_2}{{\epsilon }_1}=\frac{\frac{a}{2}{(-h)}^M{\beta }^{M-1}-(1-\frac{a}{2}){(1-\beta )}^{M-1}}{\frac{a}{2}+(1-\frac{a}{2}){(1-\beta )}^{M-1}}$$

(19)

The stress triaxiality (η) using the Barlat 1989 yield criterion can be evaluated using equation (20).

$$\mathrm{\eta }=\frac{0.33(1+\beta )}{{\left[\frac{1}{2}(a+a{\left(-h\beta \right)}^M+(2-a){(1-h\beta )}^M)\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$M$}\right.}}$$

(20)

The incompressibility constraint (dε₃ = −(dε₁ + dε₂) has been followed for the determination of the fracture limit. The non-zero components for the evaluation of the strains are defined using the flow rule $(\mathrm{d}{\epsilon }_{\mathrm{i}}=\mathrm{d}\overline{\epsilon }\frac{\partial \overline{\mathrm{\sigma }}}{\partial {\mathrm{\sigma }}_{\mathrm{i}}};\mathrm{i}=1,2)$ . The principle strain ratio can be used for obtaining the increment in the second principle strain which is given below by equation (21).

$$\mathrm{\Delta }{\overline{\epsilon }}_{2A}={\alpha }_A.\mathrm{\Delta }{\epsilon }_{1A}$$

(21)

The increment in the equivalent strain can be evaluated using equation (22) in which $\gamma$ is the ratio of equivalent to the major plastic strain.

$$\mathrm{\Delta }{\overline{\epsilon }}_A={\gamma }_A.\mathrm{\Delta }{\epsilon }_{1A}$$

(22)

At this particular stage, the strains in the area of the non-defective region have been updated using equations (23) and (24).

$${\epsilon }_{\mathrm{iA}}^{\mathrm{t}+\mathrm{\Delta }\mathrm{t}}={\epsilon }_{\mathrm{iA}}^{\mathrm{t}}+\mathrm{\Delta }{\epsilon }_{\mathrm{iA}};\mathrm{i}=1,2,\mathrm{and}3$$

(23)

$${\overline{\epsilon }}_{\mathrm{A}}^{\mathrm{t}+\mathrm{\Delta }\mathrm{t}}={\overline{\epsilon }}_{\mathrm{A}}^{\mathrm{t}}+\mathrm{\Delta }{\overline{\epsilon }}_{\mathrm{A}}$$

(24)

All the components of strain have been calculated up until now. Further, the components corresponding to the groove for t + Δt also need to be calibrated. Additionally, the link between regions A and B can be expressed as shown below.

• Continuity in the groove direction of strain; Δε_2B = Δε_2A.

• From the interface of both of the regions, equilibrium can be achieved; σ_1A·s_A = σ_1B·s_B, where, s_A and s_B are the current thicknesses of regions A and B respectively.

The continuity equation can be used for calculating the ratio of principle strains in the groove region.

$${\alpha }_B=\raisebox{1ex}{$\mathrm{\Delta }{\epsilon }_{2B}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{\Delta }{\epsilon }_{1B}$}\right.={\alpha }_A.\raisebox{1ex}{$\mathrm{\Delta }{\epsilon }_{2A}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{\Delta }{\epsilon }_{1B}$}\right.$$

(25)

Further, the relation of principle strains in region B can be written as

$$\frac{\partial {\overline{\mathrm{\sigma }}}_B}{\partial {\mathrm{\sigma }}_{1B}}\left(\mathrm{\Delta }{\epsilon }_{2B}\right)-\frac{\partial {\overline{\mathrm{\sigma }}}_B}{\partial {\mathrm{\sigma }}_{2B}}\left(\mathrm{\Delta }{\epsilon }_{1B}\right)=0$$

(26)

The principle stress ratio in the case of the Hill 1948 and Barlat 1989 yield criteria for region B can be calibrated using equations (5) and (19) respectively. Further, the unknowns, namely, α_B, β_2B, and Δε_1B can be calculated using equations (25) and (26) and the stress ratio relations for the respective yield criteria. Thus, the strains in the groove region get updated at this step.

$${\epsilon }_{\mathrm{iB}}^{\mathrm{t}+\mathrm{\Delta }\mathrm{t}}={\epsilon }_{\mathrm{iB}}^{\mathrm{t}}+\mathrm{\Delta }{\epsilon }_{\mathrm{iB}};\mathrm{i}=1,2,\mathrm{and}3$$

(27)

$${\overline{\epsilon }}_{\mathrm{B}}^{\mathrm{t}+\mathrm{\Delta }\mathrm{t}}={\overline{\epsilon }}_{\mathrm{B}}^{\mathrm{t}}+\mathrm{\Delta }{\overline{\epsilon }}_{\mathrm{B}}$$

(28)

To perform the iterative process to update the strains, the force compatibility equation $({\varnothing }_B{\overline{\mathrm{\sigma }}}_A={\varnothing }_A{\overline{\mathrm{\sigma }}}_{B}f)$ can be used to achieve the condition of equilibrium. The inhomogeneity factor (f) varies according to the relation (f = f₀ exp(ε_3B − ε_3A)). For predicting the hardening behavior of metal, the Swift equation (equation (29)) has been considered in the present study. It has the benefit of considering the combined effect of temperature and strain rate while predicting the flow stress.

$$\overline{\mathrm{\sigma }}=K({\epsilon }_0+\overline{\epsilon }{)}^n$$

(29)

Where, $\overline{\mathrm{\sigma }}$ and $\overline{\epsilon }$ are the equivalent stress and plastic strain respectively. K and n are the strength coefficient and strain hardening exponent respectively. The force compatibility equation can be rewritten as shown below in equation (30) by replacing the effective stress term using equation (29).

$${\varnothing }_B{({\epsilon }_0+\overline{\epsilon })}_A^n=f{\varnothing }_A{({\epsilon }_0+\overline{\epsilon })}_B^n$$

(30)

On achieving the parameters associated with the fracture of region B, the iterations must be stopped. The principle strains in region A at the current step are known as the limiting strains. The FFLD has been plotted using equation (31), which is shown below.

$$\begin{array}{c}\multicolumn{1}{c}{\underset{0}{\overset{{\overline{\epsilon }}_{\mathrm{B}}}{\int }}\mathrm{F}(\mathrm{process}\mathrm{parameters})\mathrm{d}{\overline{\epsilon }}_{\mathrm{B}}}\\ \multicolumn{1}{c}{\ge \mathrm{Critical}\mathrm{Damage}\mathrm{Value}(\mathrm{CDV})}\end{array}$$

(31)

The ductile damage models listed below are used in abundance for influencing the damage history of stress and strain on the occurrence of the fracture of sheet metals. The absolute value of material constants should be considered in the case of each damage model for easy mathematical computations.

• Cockcroft and Latham (C-L) ³⁴ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}{\mathrm{\sigma }}_{\max }\mathrm{d}\overline{\epsilon }=C_1$$

(32)

• Brozzo et al. (Bo) ³⁵ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}\frac{2}{3}{(1-\frac{{\mathrm{\sigma }}_h}{{\mathrm{\sigma }}_{\max }})}^{-1}\mathrm{d}\overline{\epsilon }=C_2$$

(33)

• Oyane et al. (Oy) ³⁶ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}(\frac{{\mathrm{\sigma }}_h}{\overline{\mathrm{\sigma }}}+C_3)\mathrm{d}\overline{\epsilon }=C_4$$

(34)

• Clift et al. (Cl) ³⁷ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}\overline{\mathrm{\sigma }}\mathrm{d}\overline{\epsilon }=C_5$$

(35)

• Rice and Tracey (RT) ³⁸ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}0.283\exp \left(\frac{3{\mathrm{\sigma }}_h}{2\overline{\mathrm{\sigma }}}\right)\mathrm{d}\overline{\epsilon }=C_6$$

(36)

• McClintock (MT) ³⁹ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}\frac{{\mathrm{\sigma }}_h}{\overline{\mathrm{\sigma }}}\mathrm{d}\overline{\epsilon }=C_7$$

(37)

• Oh et al. ⁴⁰ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}\frac{{\mathrm{\sigma }}_{\max }}{\overline{\mathrm{\sigma }}}\mathrm{d}\overline{\epsilon }=C_8$$

(38)

• Ko et al. ⁴¹ :

$$\underset{0}{\overset{{\overline{\epsilon }}_f}{\int }}\frac{{\mathrm{\sigma }}_{\max }}{\overline{\mathrm{\sigma }}}(1+\frac{3{\mathrm{\sigma }}_h}{\overline{\mathrm{\sigma }}})\mathrm{d}\overline{\epsilon }=C_9$$

(39)

In equation (32)–(39), ${\overline{\epsilon }}_f$ represent the equivalent plastic strain at fracture. $\overline{\mathrm{\sigma }}$ , ${\mathrm{\sigma }}_h$ , and ${\mathrm{\sigma }}_{\max }$ represent the equivalent stress, hydrostatic pressure, and maximum principle stress respectively. The C_i (i = 1–9) are the material constants or the critical damage values. On further rearranging the damage models mentioned above, a common form of integral equation shown by equation (40) has been obtained. The condition of sheet metal fracture is satisfied; in case, $I_i$ (i = 1–9) tends to unity.

$$I_i=\frac{1}{C_i}\underset{0}{\overset{{\overline{\epsilon }}_{\mathrm{f}}}{\int }}(\mathrm{damage}\mathrm{model}\mathrm{function})\mathrm{d}\overline{\epsilon }$$

(40)

Material and method

A 1 mm thick IN625 alloy sheet has been used in the present work for performing the uniaxial and stretch forming experiments at 300 and 673 K temperature. The chemical composition for the as received IN625 alloy by weight percentage has been determined according to the procedure mentioned in the ASTM E3047 ⁴² standard and is shown below in Table 1.

Table 1.

Chemical content of IN625.

Element	Ni	Cr	Nb	Mo	Ti	Al	Fe
wt.%	61.495	21.739	3.271	9.479	0.166	0.067	3.304
Element	C	Si	Mn	P	S	Co	Cu
wt.%	0.021	0.101	0.123	0.001	0.002	0.080	–

The uniaxial tensile testing has been performed using the dog bone shaped specimen, designed using the ASTM E8/E8M ⁴³ standard. The tensile test has been conducted using the specimens cut in three different directions, namely, 0° (parallel to the rolling direction of the sheet), 45° from the rolling direction, and 90° (perpendicular to the rolling direction of the sheet). The testing has been done at 300 and 673 K temperature. The testing has been done at a slow strain rate of 0.0001 s⁻¹.

A 40 ton hydraulic press has been used for to perform the experimental stretch forming process at 300 and 673 K. It has an inductive heating coil for high-temperature testing, and a K-type contact thermocouple has been used for finding the exact temperature over the whole forming setup. Figure 3 represents the schematic diagram of the stretch forming setup and design of the specimens used for finding the limiting strains in the present work. It consists of three primary components, namely, a hemispherical punch with the diameter of 50 mm, die, and blank holder. To restrict the easy flow of sheet metal, a circular draw bead with the diameter of 84 mm has been designed over the circular die. The specimen designs for stretch forming were selected based upon the ASTM E2218-15 ⁴⁴ standard and the method proposed by Sajun Prasad et al. ⁴⁵ The molybdenum disulfide based Molykote lubricant has been used in this study. The laser etching process has been used to mark the grids with the diameters of 2.5 mm over specimens to capture the major and minor strains across the deformed specimens accurately. The deformed grids over all the specimens were measured using an optical microscope with image analyzing software to obtain major and minor strains. The punch speed and the blank holding pressure (BHP) have been maintained to be 5 mm/min and 20 MPa respectively.

Figure 3.

Schematic representation of stretch forming setup and design of all the specimens.

In the present work, six types of specimens have been used to find the complete forming limit diagram of IN625 alloy at forming temperatures of 300 and 673 K. The Design 1, 2, and 3 specimens are of rectangular shape and have dimensions of 100 mm × 100 mm, 100 mm × 80 mm, and 100 mm × 60 mm, respectively. The 100 mm × 100 mm square-shaped specimens were used to induce biaxial tensile strain in the material, and they were lubricated using Molykote lubricant to impose deformation close to the equi-biaxial strain path. Rectangular specimens of the same length but with the narrower widths of 80 and 60 mm were deformed to change the strain path during deformation. However, failure in the sheet metals was consistently observed at the draw-bead position during the testing of rectangular specimens of 40-mm width. Hence, the Hasek specimens of Designs 4, 5, and 6, shown in Figure 3, were used. In these specimens, the material could flow into the die cavity from the lateral side, as the width was lower than the draw-bead diameter. However, the sheet metal gets locked in the draw bead and stretched in the longitudinal direction. The quantity of material drawing into the die cavity depended on the width of the specimen, and this induced the compressive strain and hence the strain path. ⁴⁵

Results and discussion

The engineering stress-strain curve of IN625 alloy is shown in Figure 4. The basic material properties at 300 and 673 K are shown in Table 2. On increasing the temperature from 300 to 673 K, the total elongation increased by 15.83%, while the yield and ultimate stress decreased by 19.74% and 14.81%, respectively. It is worth mentioning that serrated flow stress behavior has been reported at 673 K. This serrated behavior is popularly described by the Portevin–Le Chatelier effect or dynamic strain aging (DSA). According to the classification of Rodriguez, ⁴⁶ serrations type A, B, and C are generally observed in metal. The oscillation type of the flow stress behavior is found from 673 K, which indicates the B-type of serrated flow behavior. The fluctuation of flow stress is because of rapid dynamic recrystallization, where the former process ends before the latter one starts. Similar observations were reported in previous work done by Lin et al. ⁴⁷

Figure 4.

Engineering stress-strain curve of IN625 alloy.

Table 2.

Material properties of IN625 alloy.

Temperature (K)	σyield (MPa)	σultimate (MPa)	Total elongation (%)	Anisotropic coefficient $(r_{\theta }=\frac{{\epsilon }_{\mathrm{width}}}{{\epsilon }_{\mathrm{thickness}}})$
				r 0	r 45	r 90
300	785.21	972.14	48.62	0.85	0.96	0.74
673	630.48	828.62	56.32	1.08	0.66	0.70

The material constants for Hill 1948 and Barlat 1989 yield criteria have been estimated and are shown in Tables 3 and 4, respectively.

Table 3.

Hill 1948 yield criterion material constants.

Temperature	H	G	F	N
300K	0.614015	0.54010	0.45989	1.694937
673K	0.739487	0.47918	0.52082	1.415115

Table 4.

Barlat 1989 yield criterion material constants.

Temperature	a	c	h
300K	0.887581	1.112419	1.036294
673K	1.025793	0.974207	0.881727

The true major and minor strains for the safe and failed regions have been plotted in the principle strain graph using six different designs of specimens. The specimens were stretch formed until the fracture of the specimens. It was of particular interest that a sub-size hemispherical punch of ϕ50 mm was used in the present work instead of the standard punch of ϕ101.4 mm proposed by Hecker. ⁴⁸ Hence, the bending strain induced on the outer surface of the sheet metal while wrapped around the smaller punch during deformation was calculated using the procedure mentioned by Morchhale et al. ⁴⁹ and Prasad et al. ⁵⁰ The approximately V-shaped black line separating the safe and fractured points is called the forming limit diagram (FLD). The FLDs obtained at temperatures of 300 and 673 K after the inclusion of the bending strain are shown in Figure 5(a) and (b), respectively. The whole FLD space can be distributed into the Tension-Compression (T-C) region, Tension-Tension (T-T) region, and plane-strain (P-S, minor strain = 0) region. The point of intersection between the experimental FLD and the P-S axis is known as the FLD₀. This point can also be used for representing the formability of the material. At 300 and 673 K, the FLD₀ has been found to be 0.386 and 0.458, respectively. The FLD₀ displayed an improvement of approximately 17.26% with an increase in the temperature from 300 to 673 K. The lubricated condition has been used during the experimental stretch forming process; hence in the T-T region, the strain ratios (α) obtained are 0.72 and 0.63 at 300 and 673 K, respectively. On reducing the width of the specimens, the strain path, marked by differently colored long arrows, shifted toward the T-C region. The D-4 specimen exhibited deformation near the P-S region at both of the temperatures.

Figure 5.

FLD at (a) 300 K and (b) 673 K.

The experimental fracture forming limit diagram (FFLD) has also been evaluated using the same stretch formed specimens at 300 and 673 K temperature. Silva et al. ⁵¹ reported large errors in the direct measurement of the major fracture strains (ε_1f) because of the inhomogeneous deformation around the crack region. Thus, the volume consistency equation has been followed for evaluating the fracture strains of the specimens. ⁵² No significant stretching of the material in the lateral direction across the necked region has been observed after the appearance of necking in the specimens at both of the temperatures. Rather, the thinning of the specimens took place because of the excessive strain localization. Thus, the minor fracture strain (ε_2f) has been assumed to be equal to the necking strain (ε_2n). Further, the fractured specimens were cut in the perpendicular direction of the fracture line as shown in the representative diagram (Figure 6). The thickness (t_f1 and t_f2) across the maximum thinned region has been measured and the minimum of the two thicknesses has been considered as the fracture thickness of the specimen. The true strain along the thickness direction has been measured as (ε_3f = ln(t_f/t₀)), where t_f = min(t_f1, t_f2) and t₀ = initial thickness of the sheet. Thus, the major fracture strain has been estimated using the equation (ε_1f = −ε_2f − ε_3f). For each specimen, several fracture points have been collected and represented in the plain stain space. A curve passing below each cluster of fracture data points has been drawn and named as experimental fracture forming limit diagram (FFLD).

Figure 6.

Thickness measurement through the thickness direction of the fractured dome.

The FFLDs of IN625 alloy at 300 and 673 K have been estimated and are shown in Figure 7(a) and (b) respectively. The estimated fracture points have been fitted with the help of a straight line in the principle fracture strain locus also known as the FFLD. Table 5 represents the maximum safe forming major strains for IN625 alloy. The fracture limits of the material improved by approximately 22.22% as the forming temperature increased from 300 to 673 K.

Figure 7.

FFLD at (a) 300 K and (b) 673 K.

Table 5.

Safe forming and fracture major strains obtained from the specimens deformed along different deformation paths.

Curve	Temperature	D-1	D-2	D-3	D-4	D-5	D-6
FLD	300K	0.395	0.392	0.391	0.386	0.405	0.453
	673K	0.511	0.491	0.473	0.458	0.492	0.516
FFLD	300K	0.430	0.444	0.451	0.453	0.461	0.467
	673K	0.532	0.538	0.541	0.545	0.557	0.569

The theoretical prediction of the FFLD has been done using a modified MK model. The curve predicted using the MK model is highly sensitive to the assumed inhomogeneity factor (f₀). A proper theoretical FFLD can be obtained only if an appropriate value of f₀ has been chosen before running the whole modified MK model algorithm. A detailed representative figure for the prediction of the FFLD using a modified MK model by introducing a combination of the Cockcroft and Latham damage model and Hill 1948 criterion into the MK model at 300 K is shown in Figure 8. The FFLD found using the modified MK model with f₀ to be 0.992 followed the experimental fracture locus with the highest accuracy. A similar trial and error process needs to be performed to find the appropriate f₀ value while implementing the modified MK model using different combination of damage models and yielding functions to find the accurate and best fitted FFLD. The material constants (Table 5) for different ductile damage models shown in equations (32)–(39) have been calibrated at both 300 and 673 K using the Hill 1948 and Barlat 1989 yielding functions in order to find the fracture strains using the stretch forming process. The finalized f₀ value after applying the trial and error approach for all of the cases is shown in Table 6.

Figure 8.

Representative figure to find appropriate f₀ value using the trial and error approach.

Table 6.

Material constants for different damage models.

Temperature	Yield criteria	C 1	C 2	C 3	C 4	C 5	C 6	C 7	C 8	C 9
300 K	Hill1948	308.6	0.42	0.87	1.28	198.6	0.39	0.21	0.23	0.62
	Barlat 1989	262.4	0.58	1.34	1.66	209.4	0.24	0.39	0.39	0.94
673 K	Hill 1948	322.5	0.49	0.62	0.98	230.4	0.47	0.32	0.36	0.87
	Barlat 1989	287.2	0.71	0.99	1.15	242.5	0.35	0.48	0.57	1.27

Figure 9 shows the theoretical fracture limits predicted with the help of the modified MK model using different ductile damage models and yielding functions at 300 and 673 K. These fracture loci have been predicted using the f₀ value shown in Table 7. Notably, none of the utilized models is able to accurately predict the fracture limits (on comparing it with the experimental results) in any of the graphs shown below. Different models seem to predict the fracture limits accurately in the different regions of the FFLD. Thus, to confirm the damage model having best prediction ability of fracture strains while applying the modified MK model, the root mean square error (RMSE) has been calculated using equation (41) in each case.

$$\mathrm{RMSE}=\sqrt{\frac{(\sum ^{\}{({\epsilon }_{1e}-{\epsilon }_{1p})}^2})}{N}}$$

(41)

Figure 9.

Prediction of theoretical FFLD using modified MK model in case of (a) 300 K, Hill 1948, (b) 300 K, Barlat 1989, (c) 673 K, Hill 1948, and (d) 673 K, Barlat 1989.

Table 7.

Finalized f₀ value for all the considered combination of ductile damage models and anisotropic yield criterion.

Temperature	Yield criteria	C-L	Bo	Oy	Cl	RT	MT	Oh	Ko
300 K	Hill1948	0.992	0.987	0.997	0.989	0.990	0.984	0.994	0.993
	Barlat 1989	0.989	0.995	0.986	0.992	0.990	0.998	0.994	0.989
673 K	Hill 1948	0.988	0.991	0.999	0.987	0.994	0.992	0.990	0.997
	Barlat 1989	0.996	0.994	0.989	0.990	0.990	0.991	0.988	0.992

Figure 10 shows the RMSE obtained using the modified MK model by incorporating different ductile damage models and using the Hill 1948 and Barlat 1989 yielding functions at 300 and 673 K. At 300 K, the Cockcroft and Latham (C-L) model displayed the least RMSE while predicting the fracture limits using both of the yielding functions. However, at 673 K, the C-L model had the best prediction only in the case of the Barlat 1989 criterion. On considering the overall average of RMSE, the C-L model (RMSE = 0.062) best predicted the results and hence should be preferred when predicting the theoretical fracture limits. The Oyane (Oy, RMSE = 0.134) and Oh (RMSE = 0.137) models also predicted results with high accuracy, though not as high as the C-L model, and can also be used for predicting the FFLD. Further, the Brozzo (Bo, RMSE = 0.183) model displayed the highest overall average RMSE; hence, it should be least preferred to predict the fracture limits. The combination of the C-L damage model and the Barlat 1989 yield criterion has been found to predict the theoretical fracture limits of IN625 alloy with the highest accuracy using the modified MK model.

Figure 10.

RMSE obtained in case of different considered cases.

Conclusions

In the present study, the experimental forming and fracture limits have been evaluated using the stretch forming process at 300 and 673 K using six different design specimens. On increasing the temperature from 300 to 673 K, the forming and fracture limits of the material improved by approximately 17.26% and 22.22%, respectively. The classical Marciniak and Kuczyński (MK) model has been modified by introducing eight different ductile damage models into it for the prediction of fracture limits using the Hill 1948 and Barlat 1989 yielding functions. This study helped in finding the best combination of the ductile damage model and yield criterion for the theoretical prediction of the FFLD using the modified MK model. Among the 16 different combinations of ductile damage models and yielding functions, the Cockcroft and Latham (C-L) damage model along with the Barlat 1989 yielding function best predicted the theoretical FFLD with the smallest average root mean square error (RMSE) of 0.033. Some of the other models, namely, the Rice-Tracey, Brozzo, McClintock, and Ko models, displayed an average RMSE beyond the acceptance limit, and hence, they should be the least preferred for the prediction of fracture limits while using the modified MK model.

Future work involves the further improvement of the modified MK model by introducing the critical thickness imperfection effect, which will be related to the equivalent plastic strain, initial grain size, and surface roughness, as these factors play an important role in the critical sheet metal forming process. Also, effect factors such as the initial groove angle, sheet thickness, and strain hardening coefficient, strain rate sensitivity, constitutive model, advanced yield criteria, and inhomogeneity factor over the predicted fracture limit diagram using the modified MK model will also be communicated in the future.