Mechanics-Based Failure Model of Tube Hydro-Bulging Test

Tube hydroforming has been widely applied by the automotive sector to produce hollow parts. As a popular tube hydroforming test method, tube hydro-bulging needs an analytical failure model to analyze the formability of tubular materials. In the present work, a failure prediction model has been developed to predict the bulging height limit (BHL) of the hydro-bulging test. The model utilized Hill’s orthogonal anisotropic model to characterize the tube material, a geometry model to characterize the non-loading path and the M-K model to predict failure. Defects in multiple directions were taken into consideration. The developed model was applied on two tubes of different materials as case studies to verify its validity. It is shown that the developed model is capable of predicting the forming limit or determining the imperfection factor of tubular materials.

1 Introduction

Hydroforming technology has been widely adopted by the automotive sector in recent decades due to the capability of forming complex geometries from lightweight materials, whilst avoiding joining processes. Tube hydroforming (THF) uses tubes as the raw material and applies internal pressure and/or axial compression to form hollow parts, including tubular parts, irregular cross-sections, or multi-way tubes.

Forming limit test methods of THF is essential to the prediction of failure. Nakajima test is the standard way to determine the forming limit (International Standard Organization, 2021) for sheet materials. However, this standard does not apply to tubular materials. One of the most popular methods to solve this problem is the tube hydro-bulging test (or hydraulic/tube bulge test). Numerous research studies have developed experimental devices to perform such tests (Fuchizawa et al., 1993; Sokolowski et al., 2000; Fllice et al., 2001; Aue-u-lan, 2007), and the hydro-bulging test has been effectively used to predict the forming limit (Zhu et al., 2020).

However, an analytical derivation of the forming limit prediction for the hydro-bulging test has not been developed yet, which is the primary aim of the article. There are many existing analytical models for forming limit prediction, including Swift’s diffuse necking model (Swift, 1952), Hill’s localized necking model (Hill, 1948), and the M-K model (Marciniak and Kuczyński, 1967). The original M-K model assumes that an imperfection area only exists across the width of the sheet. Hutchinson and Neale (1978) modified the M-K model by assuming that the imperfection area direction is arbitrary. The modified M-K model was adopted in the present work to study the potential necking behavior in all directions.

The M-K model takes the loading path of a point as input. In the context of the hydro-bulging test, the point is the pole point, as it undergoes the largest deformation. The loading path of the point can be obtained from stress and strain analysis, which depends on the geometry model of the tube. The key difference between the published geometries is the assumption of the bulging zone’s profile, such as the circular arc (Boudeau and Malécot, 2012), cosine-like function (Strano and Altan, 2004), and elliptical curve (Hwang and Lin, 2007; He et al., 2014a). He et al. (2014a) proposed two models with elliptical curve assumption. The first assumed that the filleted corner was negligible, and the bulging tube can be approximated by only one elliptic arc. The second took the filleted corner into account, and the bulging profile was approximated by the combination of two circular arcs and an elliptic arc. The latter geometry model has better accuracy and was adopted in the present study.

After determining the geometry model, stress and strain analysis should be conducted at the pole point to obtain the loading path. As one of the most important conditions to be analyzed, the end-conditions have four different types, namely, free-ends (Imaninejad, Subhash and Loukus, 2004), closed-ends (Fuchizawa, Narazaki and Yuki, 1993), fixed-ends (Sokolowski et al., 2000; Imaninejad, Subhash and Loukus, 2004; Hwang and Lin, 2007; Hwang and Wang, 2009; Boudeau and Malécot, 2012; He et al., 2014a), and forced-ends (Imaninejad, Subhash and Loukus, 2004; Kuwabara et al., 2005), which can be applied in tube hydro-bulging test. The fixed-ends condition was mostly used in previous studies due to the simplest mathematical expression. Thus, it was adopted in the present study.

Hill’s orthogonal anisotropic model (Hill, 1950) was used to describe the plastic mechanical behavior of the tubular material. For the convenience of model derivation, a tube-friendly version of the model was derived in advance.

Combining Hill’s orthogonal anisotropic model, hydro-bulging geometry model, and M-K model, the present study developed a mechanics-based hydro-bulging test failure model to enable the prediction of bulging height limit (BHL) in tube hydro-bulging tests. Two case studies were reviewed to show the prediction capability of the model. The case studies illustrated how to fit initial imperfection factors $f_0$ from experimental results, found a way to simplify the model, and verified the validity of the model by comparing the required pressure evolution and pole point thickness prediction with experiments.

2 Mechanics-Based Failure Model

2.1 Hill’s Orthogonal Anisotropic Model for Tubular Materials

The original form of Hill’s orthogonal anisotropic model is too general to be applied in the specific file of the tube hydro-bulging test. In this section, a specific form for tubular materials was derived, as a footstone for further strain and stress analysis.

2.1.1 Derivation of Equivalent Stress

The general form of the Hill yield function (Hill 1950) is as follows:

${\left({\sigma }_{22}-{\sigma }_{33}\right)}^2+G{\left({\sigma }_{33}-{\sigma }_{11}\right)}^2+H{\left({\sigma }_{11}-{\sigma }_{22}\right)}^2+2L{\sigma }_{23}^2+2M{\sigma }_{31}^2+2N{\sigma }_{12}^2=1(1)$

where $F,G,H,L,M,N$ are anisotropy constants and determined experimentally. By using the coordinates/subscripts $z\theta r$ to replace the numerals and the elimination of shear stress due to the material undergoing plane stress conditions, the yield criterion can be derived as follows:

$f=F{\left({\sigma }_{\theta }-{\sigma }_{\mathrm{r}}\right)}^2+G{\left({\sigma }_{\mathrm{r}}-{\sigma }_{\mathrm{z}}\right)}^2+H{\left({\sigma }_{\mathrm{z}}-{\sigma }_{\theta }\right)}^2+2N{\sigma }_{\mathrm{z}\theta }^2=1(2)$

Since the material anisotropy is assumed to be orthogonal, once the uniaxial normal yield stress is achieved, the yield criterion can be expressed as follows:

$\left(G+H\right){\left({\sigma }_z^y\right)}^2=(F+H){\left({\sigma }_{\theta }^y\right)}^2=(F+G){\left({\sigma }_{\mathrm{r}}^y\right)}^2=1(3)$

where ${\sigma }_{\theta }^y,{\sigma }_{\mathrm{z}}^y,{\sigma }_{\mathrm{r}}^y$ are the uniaxial yield stresses on the axes of anisotropy. By solving for the constants $F,G,H$ from Eq. 3, we derive the following:

$F=\frac{1}{2}\left[\frac{1}{{\left({\sigma }_{\theta }^y\right)}^2}+\frac{1}{{\left({\sigma }_{\mathrm{r}}^y\right)}^2}-\frac{1}{{\left({\sigma }_z^y\right)}^2}\right](4)$

$G=\frac{1}{2}\left[\frac{1}{{\left({\sigma }_{\mathrm{r}}^y\right)}^2}+\frac{1}{{\left({\sigma }_z^y\right)}^2}-\frac{1}{{\left({\sigma }_{\theta }^y\right)}^2}\right](5)$

$H=\frac{1}{2}\left[\frac{1}{{\left({\sigma }_z^y\right)}^2}+\frac{1}{{\left({\sigma }_{\theta }^y\right)}^2}-\frac{1}{{\left({\sigma }_{\mathrm{r}}^y\right)}^2}\right](6)$

Subsequently, the associated flow rule is applied:

$d{\epsilon }_{ij}=d\lambda \frac{\partial f}{\partial {\sigma }_{ij}}\to d\lambda =d{\epsilon }_{ij}/\frac{\partial f}{\partial {\sigma }_{ij}}(7)$

where the ratio between $d{\epsilon }_z$, $d{\epsilon }_{\theta }$, $d{\epsilon }_r$, and $d{\epsilon }_{z\theta }$ can be obtained as follows:

$\frac{d{\epsilon }_z}{2\left(G+H\right){\sigma }_{\mathrm{z}}-2H{\sigma }_{\theta }-2G{\sigma }_{\mathrm{r}}}=\frac{d{\epsilon }_{\theta }}{2\left(F+H\right){\sigma }_{\theta }-2H{\sigma }_z-2F{\sigma }_{\mathrm{r}}}=\frac{d{\epsilon }_r}{2\left(F+G\right){\sigma }_{\mathrm{r}}-2G{\sigma }_z-2F{\sigma }_{\theta }}=\frac{d{\epsilon }_{z\theta }}{2N{\sigma }_{z\theta }}(8)$

Note that the last item is not $\frac{d{\epsilon }_{z\theta }}{4N{\sigma }_{z\theta }}$, because the term of $2N{\sigma }_{\mathrm{z}\theta }^2$ is divided into $N{\sigma }_{\mathrm{z}\theta }^2+N{\sigma }_{\theta z}^2$ when calculating the partial derivative $\frac{\partial f}{\partial {\sigma }_{ij}}$.

The wall thickness at the pole point of the tube is low ($\frac{t_0}{R_0}<\frac{1}{10}$), and thus the radial stress components can be ignored (${\sigma }_r=0$), resulting in the state of plane stress (Zhu et al., 2020). For plane stress ${\sigma }_{\text{r}}=0$, Eq. 8 is simplified as follows:

$\frac{d{\epsilon }_z}{2\left(G+H\right){\sigma }_{\mathrm{z}}-2H{\sigma }_{\theta }}=\frac{d{\epsilon }_{\theta }}{2\left(F+H\right){\sigma }_{\theta }-2H{\sigma }_z}=\frac{d{\epsilon }_r}{-2G{\sigma }_z-2F{\sigma }_{\theta }}=\frac{d{\epsilon }_{z\theta }}{2N{\sigma }_{z\theta }}(9)$

by defining the ratio $r_z=\frac{d{\epsilon }_{\theta }}{d{\epsilon }_r}$, $r_{\theta }=\frac{d{\epsilon }_z}{d{\epsilon }_r}$, and $r_{z\theta }=\frac{d{\epsilon }_{45}}{d{\epsilon }_r}$ (note that ${\epsilon }_{z\theta }$ is shear strain while ${\epsilon }_{45}$ is normal strain) under uniaxial tensile stress. ${\sigma }_{\text{z}},$ ${\sigma }_{\theta }$, or ${\sigma }_{45}$ can be obtained by uniaxial tensile tests with directions shown in Figure 1.

FIGURE 1

FIGURE 1. Directions of the anisotropy coefficients.

Note that in the uniaxial tensile tests used to calculate $r_z$, ${\sigma }_{\theta }$ equals to zero; similarly, in the uniaxial tensile tests used to calculate $r_{\theta }$, ${\sigma }_z$ equals to zero. Thus, simplified expressions for $r_z$ and $r_{\theta }$ can be obtained from Eq. 9:

$r_z=\frac{d{\epsilon }_{\theta }}{d{\epsilon }_r}=\frac{H}{G};r_{\theta }=\frac{d{\epsilon }_z}{d{\epsilon }_r}=\frac{H}{F}(10)$

For $r_{z\theta }$, the expression is more complex because the principal axes are not coincident with the axes of anisotropy. Additional transformation is needed to express $r_{z\theta }$ by $F,G,H,N$:

$\left[\begin{array}{c}{\sigma }_z\\ {\sigma }_{\theta }\\ {\sigma }_{z\theta }\end{array}\right]=\mathit{S}\left[\begin{array}{c}{\sigma }_{45}\\ 0\\ 0\end{array}\right](11)$

$\left[\begin{array}{c}d{\epsilon }_z\\ d{\epsilon }_{\theta }\\ d{\epsilon }_{z\theta }\end{array}\right]=\mathit{S}\left[\begin{array}{c}d{\epsilon }_{45}\\ -d{\epsilon }_{45}-d{\epsilon }_r\\ 0\end{array}\right](12)$

where

$\mathit{S}=\left[\begin{array}{ccc}{\text{cos}}^2\left(45°\right)& {\text{sin}}^2\left(45°\right)& 2\text{sin}\left(45°\right)\text{cos}\left(45°\right)\\ {\text{sin}}^2\left(45°\right)& {\text{cos}}^2\left(45°\right)& -2\text{sin}(45°)\text{cos}\left(45°\right)\\ -\text{sin}\left(45°\right)\text{cos}\left(45°\right)& \text{sin}\left(45°\right)\text{cos}\left(45°\right)& {\text{cos}}^2\left(45°\right)-{\text{sin}}^2\left(45°\right)\end{array}\right]$

is the transformation matrix.

By substituting ${\sigma }_z$, ${\sigma }_{\theta }$, ${\sigma }_{z\theta }$, $d{\epsilon }_{z\theta }$ obtained from Eqs 11, 12 into Eq. 9, $r_{z\theta }$ can be expressed as follows:

$r_{z\theta }=\frac{d{\epsilon }_{45}}{d{\epsilon }_r}=\frac{\text{N}}{\text{G}+\text{F}}-\frac{1}{2}(13)$

By eliminating ${\sigma }_{\theta }^y$ and ${\sigma }_r^y$ in Eqs 4, 5, 6, 10 and 13, $F,G,H,N$ can be expressed by using anisotropy coefficients and ${\sigma }_z^y$:

$F=\frac{r_z}{r_{\theta }\left(1+r_z\right){\left({\sigma }_z^y\right)}^2}(14)$

$G=\frac{1}{\left(1+r_z\right){\left({\sigma }_z^y\right)}^2}(15)$

$H=\frac{r_z}{\left(1+r_z\right){\left({\sigma }_z^y\right)}^2}(16)$

$N=\frac{\left(2r_{z\theta }+1\right)\left(r_z+r_{\theta }\right)}{2r_{\theta }\left(1+r_z\right){\left({\sigma }_z^y\right)}^2}(17)$

By defining ${\overline{\sigma }}^2=\frac{3r_{\theta }\left(1+r_z\right)}{2\left(r_z{r}_{\theta }+r_{\theta }+r_z\right)}{\left({\sigma }_z^y\right)}^2$, the anisotropy constants and the equivalent stress can be expressed as follows (Hwang and Lin, 2006; Hwang and Wang, 2009):

$F=\frac{3r_z}{2\left(r_z{r}_{\theta }+r_{\theta }+r_z\right){\overline{\sigma }}^2}(18)$

$G=\frac{3r_{\theta }}{2\left(r_z{r}_{\theta }+r_{\theta }+r_z\right){\overline{\sigma }}^2}(19)$

$H=\frac{3r_z{r}_{\theta }}{2\left(r_z{r}_{\theta }+r_{\theta }+r_z\right){\overline{\sigma }}^2}(20)$

$N=\frac{3\left(2r_{z\theta }+1\right)\left(r_{\theta }+r_z\right)}{4\left(r_z{r}_{\theta }+r_{\theta }+r_z\right){\overline{\sigma }}^2}(21)$

$\overline{\sigma }=\sqrt{\frac{3}{2}}{\left[\frac{\left(1+\left(1/r_z\right)\right){\sigma }_z^2-2{\sigma }_z{\sigma }_{\theta }+\left(1+\left(1/r_{\theta }\right)\right){\sigma }_{\theta }^2+\left(\left(1/r_{\theta }\right)+\left(1/r_z\right)\right)\left(2r_{z\theta }+1\right){\sigma }_{z\theta }^2}{\left(1/r_{\theta }\right)+1+\left(1/r_z\right)}\right]}^{1/2}(22)$

Consequently, the associate flow rule Eq. 9 can be expressed by $r_{\theta }$, $r_z$, and $r_{z\theta }$:

$\frac{d{\epsilon }_z}{r_{\theta }\left({\sigma }_z+r_z{\sigma }_z-r_z{\sigma }_{\theta }\right)}=\frac{d{\epsilon }_{\theta }}{r_z\left({\sigma }_{\theta }-r_{\theta }{\sigma }_z+r_{\theta }{\sigma }_{\theta }\right)}=\frac{d{\epsilon }_{z\theta }}{{\sigma }_{z\theta }\left(2r_{z\theta }+1\right)\left(r_z+r_{\theta }\right)}(23)$

2.1.2 Derivation of Equivalent Strain Increment

The equivalent stress can also be expressed in matrix notation (Mohr et al., 2010):

$\overline{\sigma }=\sqrt{\left(\mathbf{p\sigma }\right)\cdot \mathbf{\sigma }}(24)$

where $\mathbf{\sigma }$ is the vector form of stress components:

$\mathit{\sigma }=\left[\begin{array}{c}{\sigma }_z\\ {\sigma }_{\theta }\\ {\sigma }_{z\theta }\end{array}\right]$

$\mathbf{p}$ is the factor matrix:

$\mathbf{p}=\frac{3}{2[\left(1/r_{\theta }\right)+1+\left(1/r_z\right)]}\left[\begin{array}{ccc}1+\left(1/r_z\right)& -1& 0\\ -1& 1+\left(1/r_{\theta }\right)& 0\\ 0& 0& \left[\left(1/r_{\theta }\right)+\left(1/r_z\right)\right]\left(2r_{z\theta }+1\right)\end{array}\right]$

By applying the associate flow rule under matrix notation, the following is derived:

$\mathit{d\epsilon }=d\lambda \frac{\text{d}\overline{\sigma }}{\partial \mathit{\sigma }}=d\lambda \left[\frac{d\left(\left(\mathbf{p\sigma }\right)\cdot \mathbf{\sigma }\right)}{2\overline{\sigma }d\mathbf{\sigma }}\right]=d\lambda \left(\frac{\mathbf{p\sigma }}{\overline{\sigma }}\right)(25)$

where $\mathit{d\epsilon }$ is the vector form of strain increment components:

$\mathit{d\epsilon }=\left[\begin{array}{c}d{\epsilon }_z\\ d{\epsilon }_{\theta }\\ 2d{\epsilon }_{z\theta }\end{array}\right]$

By applying work conjugation $\mathit{\sigma }\cdot \mathit{d\epsilon }=\overline{\sigma }\cdot d\overline{\epsilon }$ , the equivalent strain increment $d\overline{\epsilon }$ can be expressed as follows:

$d\overline{\epsilon }=\frac{\mathit{\sigma }\cdot \mathit{d\epsilon }}{\overline{\sigma }}=\frac{\mathit{\sigma }\cdot \left[d\lambda \left(\frac{\mathbf{p\sigma }}{\overline{\sigma }}\right)\right]}{\overline{\sigma }}=\frac{\mathit{\sigma }\cdot \mathbf{p\sigma }}{{\overline{\sigma }}^2}d\lambda =d\lambda \left(\frac{{\overline{\sigma }}^2}{{\overline{\sigma }}^2}\right)=d\lambda (26)$

By Substituting $d\mathbf{\lambda }=d\overline{\mathbf{\epsilon }}$ into Eq. 25

$\mathit{d\epsilon }=d\overline{\epsilon }\left(\frac{\mathbf{p\sigma }}{\overline{\sigma }}\right)$

$\mathit{\sigma }=\frac{\overline{\sigma }}{d\overline{\epsilon }}\left({\mathbf{p}}^{-1}\mathit{d\epsilon }\right)(27)$

By Applying Work Conjugation $\mathit{\sigma }\cdot \mathit{d\epsilon }=\overline{\mathit{\sigma }}\cdot \mathit{d}\overline{\epsilon }$ Again

$\mathit{\sigma }\cdot \mathit{d\epsilon }=\frac{\overline{\sigma }}{d\overline{\epsilon }}\left[\left({\mathbf{p}}^{-1}\mathit{d\epsilon }\right)\cdot \mathit{d\epsilon }\right]=\overline{\sigma }\cdot d\overline{\epsilon }$

$d\overline{\epsilon }=\sqrt{\left({\mathbf{p}}^{-1}\mathit{d\epsilon }\right)\cdot \mathit{d\epsilon }}(28)$

where

${\mathbf{p}}^{-\mathbf{1}}=\frac{2\left(r_{\theta }+r_z+r_{\theta }{r}_z\right)}{3}\left[\begin{array}{ccc}\frac{r_{\theta }+1}{r_{\theta }+r_{\theta }{r}_z{+}_{\theta }^2\mathrm{r}}& \frac{1}{r_{\theta }+r_z+1}& 0\\ \frac{1}{r_{\theta }+r_z+1}& \frac{r_z+1}{r_z+r_{\theta }{r}_z{+}_z^2\mathrm{r}}& 0\\ 0& 0& \frac{1}{r_{\theta }+r_z+2\left(r_{\theta }+r_z\right)r_{z\theta }}\end{array}\right]$

Thus, the equivalent strain increment can be expressed as follows:

$d\overline{\epsilon }=\sqrt{\frac{2\left(r_{\theta }+r_z+r_{\theta }{r}_z\right)}{3}}\sqrt{\frac{1}{\left(r_{\theta }+r_z+1\right)}\left(2d{\epsilon }_{z}d{\epsilon }_{\theta }+\frac{\left(r_{\theta }+1\right)d{\epsilon }_z}{r_{\theta }}+\frac{\left(r_z+1\right)d{\epsilon }_{\theta }}{r_z}\right)+\frac{4d{\epsilon }_{z\theta }}{\left(2r_{z\theta }+1\right)\left(r_{\theta }+r_z\right)}}(29)$

2.2 Geometry Model of Hydro-Bulging Tube

The geometry model was applied to describe the plastic deformation of the bulging profile. Once the bulging profile can be expressed in mathematical ways, the strain can be calculated, and thus the loading path for the M-K model can be established.

The geometry model adopted in the present study is shown in Figure 2. The failure prediction model begins from the bulging height at the pole point $h_i=0$. In the following iteration, $h_i=h_{i-1}+dh$, where $i$ indicates the step number. The expression of the circular and the elliptical arc in the first quadrant in step $i$ is expressed as follows:

${\left[r-\left(R_0+R_d\right)\right]}^2+{\left(z-\frac{L}{2}\right)}^2=R_d^2\left(z<=\frac{L}{2},r<=r_A\right)(30)$

$\frac{z^2}{a_i^2}+\frac{r^2}{b_i^2}=1\left(0<z<z_A,r>=r_A\right)(31)$

where $L$ is the total length of the tube, $R_0$ is the initial radius of the tube, $R_d$ is the corner radius, and $a_i,b_i$ are the length of the major and minor axes of the elliptic arc. The elliptic arc passes through the pole point $P\left(r_{P,i},0\right)$, which gives the explicit expression of $b_i$:

$b_i=r_{P,i}=R_0+h_i(32)$

FIGURE 2

FIGURE 2. Geometry of the die-related model.

The profile equations can be written in the form of $r=F\left(z\right)$:

$r=\left(R_0+R_d\right)-\sqrt{R_d^2-{\left(z-\frac{L}{2}\right)}^2}(33)$

$r=\sqrt{b_i^2-\frac{z^2{b}_i^2}{a_i^2}}(34)$

These functions pass through the intersection point A $\left(r_{A,i},z_{A,i}\right)$:

$r_{A,i}=R_0+R_d-\sqrt{R_d^2-{\left(z_{A,i}-\frac{L}{2}\right)}^2}(35)$

$r_{A,i}=b_i\sqrt{1-\frac{z_{A,i}^2}{a_i^2}}(36)$

The derivative of the two profile functions at point A $\left(r_{A,i},z_{A,i}\right)$ is continuous:

$\frac{2z_{A,i}-L}{2\sqrt{R_d^2-{\left(z_{A,i}-\frac{L}{2}\right)}^2}}=-\frac{b_i{z}_{A,i}}{a_i\sqrt{a_i^2-z_{A,i}^2}}(37)$

Once bulging height and the geometry of the initial tube and die are determined, $r_{A,i},z_{A,i},a_i$ can be solved numerically from Eqs 35–37. This system of equations can be simplified as one equation with only one unknown $z_{A,i}$:

$\frac{z_{A,i}}{{\text{A}}_i}-\frac{L}{2{\text{A}}_i}=-\frac{R_0+h_i}{z_{A,i}\sqrt{B_i}\sqrt{B_i-1}}(38)$

where $A_i=\sqrt{R_d^2-{\left(\frac{L}{2}-z_{A,i}\right)}^2}$ and $B_i=\frac{{\left(R_0+h_i\right)}^2}{{\left(R_0+R_d-{\text{A}}_i\right)}^2-{\left(R_0+h_i\right)}^2}$.

Substituting the function of the elliptic arc Eq. 36 into ${\rho }_z=\frac{{\left(1+r{\prime }^2\right)}^{3/2}}{\left|r″\right|}$, the curvature radius in the axial direction at the pole point $P\left(r_{P,i},0\right)$ can be obtained:

${\rho }_{z,i}=\frac{a_i^2}{b_i}(39)$

The curvature radius in hoop direction at the same point is as follows:

${\rho }_{\theta ,i}=b_i(40)$

2.3 Strain and Stress Analysis at the Pole Point

The pole point of the tube-hydro-bugling test undergoes the largest plastic deformation and thus neck first. The loading path of the pole point can be derived based on the geometry model. Figure 3 shows the geometry and stress of the infinitesimal element at the hydro-bulging tube’s pole. By analyzing this infinitesimal element, the expression of stress and strain components can be derived (Hwang and Lin, 2007; He et al., 2014a; He et al., 2014b; Zhu et al., 2020).

FIGURE 3

FIGURE 3. Infinitesimal element at the hydro-bulging tube’s pole.

Radial strain and hoop strain components on the pole point can be written as follows:

${\epsilon }_{\theta ,i}=\mathit{ln}\left(\frac{{\rho }_{\theta ,i}-\frac{t_{p,i}}{2}}{R_0-\frac{t_0}{2}}\right)(41)$

${\epsilon }_{t,i}=\mathit{ln}\left(\frac{t_{p,i}}{t_0}\right)(42)$

The expression of axial strain can be calculated through volume constancy:

${\epsilon }_{z,i}=\mathrm{-}\left({\epsilon }_{\theta ,i}+{\epsilon }_{t,i}\right)=ln\left(\frac{t_0\left(R_0\mathrm{-}\frac{t_0}{2}\right)}{t_{p,i}\left({\rho }_{\theta ,i}\mathrm{-}\frac{t_{p,i}}{2}\right)}\right)(43)$

where $t_0,t_{p,i}$ are the thickness at the undeformed and deformed stage, respectively.

The force equilibrium equation in the radial ($r$) direction of the element is as follows:

$p_i\left({\rho }_{z,i}-t_{p,i}\right)\left({\rho }_{\theta ,i}-t_{p,i}\right)d\phi d\theta =2{\sigma }_{\theta ,i}d\theta \left({\rho }_{z,i}-\frac{t_{p,i}}{2}\right)t_{p,i}\sin \frac{d\phi }{2}+2{\sigma }_{z,i}d\phi \left({\rho }_{\theta ,i}-\frac{t_{p,i}}{2}\right)t_{p,i}\sin \frac{d\theta }{2}(44)$

where $p_i$ is the internal pressure during the process, and ${\sigma }_{\theta ,\text{i}},{\sigma }_{z,\text{i}}$ are stress components in the hoop and axial direction. The equilibrium equation can be simplified as follows:

$\frac{p_i}{t_{p,i}}=\frac{{\sigma }_{\theta ,i}\left({\rho }_{z,i}-\frac{t_{p,i}}{2}\right)}{\left({\rho }_{z,i}-t_{p,i}\right)\left({\rho }_{\theta ,i}-t_{p,i}\right)}+\frac{{\sigma }_{z,i}\left({\rho }_{\theta ,i}-\frac{t_{p,i}}{2}\right)}{\left({\rho }_{z,i}-t_{p,i}\right)\left({\rho }_{\theta ,i}-t_{p,i}\right)}(45)$

By applying fixed-ends boundary condition on the axial direction force equilibrium equation, the following is derived:

${\sigma }_{z,i}\pi \left(2{\rho }_{\theta ,i}-t_{p,i}\right)t_{p,i}=p_i\pi {\left({\rho }_{\theta ,i}-t_{p,i}\right)}^2(46)$

where $R_0$ and $t_0$ are the initial radius and thickness of the tube. From Eq. 46, the expression of ${\sigma }_z$ can be written as follows:

${\sigma }_{z,i}=\frac{p_i{\left({\rho }_{\theta ,i}-t_{p,i}\right)}^2}{t_{p,i}\left(2{\rho }_{\theta ,i}-t_{p,i}\right)}(47)$

By substituting Eq. 47 into Eq. 45, the hoop stress components can be expressed as follows:

${\sigma }_{\theta ,i}=\frac{p_i\left({\rho }_{\theta ,i}-t_{p,i}\right)\left(2{\rho }_{z,i}-{\rho }_{\theta ,i}-t_{p,i}\right)}{t_{p,i}\left(2{\rho }_{z,i}-t_{p,i}\right)}(48)$

By substituting Eqs 47, 48 into Eq. 23, the associated flow rule can be expressed as follows:

$\frac{d{\epsilon }_{z,i}}{d{\epsilon }_{\theta ,i}}=\frac{-\left(d{\epsilon }_{\theta ,i}+d{\epsilon }_{r,i}\right)}{d{\epsilon }_{\theta ,i}}=\frac{r_{\theta ,i}\left(-2{\rho }_{z,i}{\rho }_{\theta ,i}+2{\rho }_{z,i}{t}_{p,i}+{\rho }_{\theta ,i}{t}_{p,i}-2r_{z,i}{\rho }_{\theta ,i}^2-t_{p,i}^2+2r_{z,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}{r_{z,i}\left(2{\rho }_{z,i}{t}_{p,i}-4{\rho }_{z,i}{\rho }_{\theta ,i}+{\rho }_{\theta ,i}{t}_{p,i}+2r_{\theta ,i}{\rho }_{\theta ,i}^2+2{\rho }_{\theta ,i}^2-t_{p,i}^2-2r_{\theta ,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}(49)$

By converting the differential equation into difference form, the associated flow rule can be expressed as follows:

$\frac{{\epsilon }_{r,i}-{\epsilon }_{r,i-1}}{{\epsilon }_{\theta ,i}-{\epsilon }_{\theta ,i-1}}=\frac{r_{\theta }\left(2{\rho }_{z,i}{\rho }_{\theta ,i}-2{\rho }_{z,i}{t}_{p,i}-{\rho }_{\theta ,i}{t}_{p,i}+2r_{z,i}{\rho }_{\theta ,i}^2+t_{p,i}^2-2r_{z,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}{r_{z,i}\left(2{\rho }_{z,i}{t}_{p,i}-4{\rho }_{z,i}{\rho }_{\theta ,i}+{\rho }_{\theta ,i}{t}_{p,i}+2r_{\theta ,i}{\rho }_{\theta ,i}^2+2{\rho }_{\theta ,i}^2-t_{p,i}^2-2r_{\theta ,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}-1(50)$

By substituting Eqs 41, 42 into Eq. 50, a non-linear equation with only one unknown, $t_{p,i}$, can be obtained as follows:

$\frac{ln\left(\frac{t_{p,i}}{t_0}\right)-{\epsilon }_{r,i-1}}{ln\left(\frac{{\rho }_{\theta ,i}-\frac{t_{p,i}}{2}}{R_0-\frac{t_0}{2}}\right)-{\epsilon }_{\theta ,i-1}}=\frac{r_{\theta }\left(2{\rho }_{z,i}{\rho }_{\theta ,i}-2{\rho }_{z,i}{t}_{p,i}-{\rho }_{\theta ,i}{t}_{p,i}+2r_{z,i}{\rho }_{\theta ,i}^2+t_{p,i}^2-2r_{z,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}{r_z\left(2{\rho }_{z,i}{t}_{p,i}-4{\rho }_{z,i}{\rho }_{\theta ,i}+{\rho }_{\theta ,i}{t}_{p,i}+2r_{\theta ,i}{\rho }_{\theta ,i}^2+2{\rho }_{\theta ,i}^2-t_{p,i}^2-2r_{\theta ,i}{\rho }_{z,i}{\rho }_{\theta ,i}\right)}-1(51)$

Once $t_{p,i}$ is solved, all the strain components can be calculated using Eqs41–43. The strain component increment can be calculated by subtracting the total strain components in the previous step. Then, the equivalent strain increment can be obtained by referring to Eq. 29:

$d{\overline{\epsilon }}_i=\sqrt{\frac{2\left(r_{\theta }+r_z+r_{\theta }{r}_z\right)}{3}}\sqrt{\frac{1}{\left(r_{\theta }+r_z+1\right)}\left(2d{\epsilon }_{z,i}d{\epsilon }_{\theta ,i}+\frac{\left(r_{\theta }+1\right)d{\epsilon }_{z,i}}{r_{\theta }}+\frac{\left(r_z+1\right)d{\epsilon }_{\theta ,i}}{r_z}\right)+\frac{4d{\epsilon }_{z\theta ,i}}{\left(2r_{z\theta }+1\right)\left(r_{\theta }+r_z\right)}}(52)$

The Total Equivalent Strain

${\overline{\epsilon }}_i={\overline{\epsilon }}_{i-1}+d{\overline{\epsilon }}_i(53)$

The equivalent stress in Zone $a$ can be expressed by referring to Eq. 22:

${\overline{\sigma }}_i=\sqrt{\frac{3}{2}}{\left[\frac{\left(1+\left(1/r_{\theta }\right)\right){\sigma }_{\theta ,i}^2-2{\sigma }_{\theta ,i}{\sigma }_{z,i}+\left(1+\left(1/r_z\right)\right){\sigma }_{z,i}^2}{\left(1/r_{\theta }\right)+1+\left(1/r_z\right)}\right]}^{\frac{1}{2}}(54)$

The flow stress curve of the tubular material was expressed by the following:

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

where ${\sigma }_0$ is the initial yield stress, $K$ is the strength coefficient, and $n$ is the strain hardening exponent. By substituting Eqs 47, 48, 52, 54 into Eq. 55, the expression of $p_i$ can be obtained as follows:

$p_i=\frac{2\left({\sigma }_0+K{\overline{\epsilon }}_i^n\right)\sqrt{\frac{1}{r_{\theta }}+\frac{1}{r_z}+1}}{\sqrt{6}\sqrt{\frac{\left(\frac{1}{r_z}+1\right){\left({\rho }_{\theta }-t_{p,i}\right)}^4}{t_{p,i}^2{\left(2{\rho }_{\theta }-t_{p,i}\right)}^2}+\frac{\left(\frac{1}{r_{\theta }}+1\right){\left({\rho }_{\theta }-t_{p,i}\right)}^2{\left({\rho }_{\theta }-2{\rho }_z+t_{p,i}\right)}^2}{t_{p,i}^2{\left(2{\rho }_z-t_{p,i}\right)}^2}+\frac{2{\left({\rho }_{\theta }-t_{p,i}\right)}^3\left({\rho }_{\theta }-2{\rho }_z+t_{p,i}\right)}{t_{p,i}^2\left(2{\rho }_z-t_{p,i}\right)\left(2{\rho }_{\theta }-t_{p,i}\right)}}}(56)$

All the stress components and equivalent stress can be obtained by back substitution of $p_i$. The collection of ${\epsilon }_{\theta ,i}$ and ${\epsilon }_{z,i}$ under different step $i$ formed the loading path that the M-K model needs.

2.4 M-K Model

After acquiring the loading path at the pole point, the necking prediction can be started. Hutchinson and Neale (1978) gave a modified M-K model, which is shown in Figure 4.

FIGURE 4

FIGURE 4. Schematic diagram of the M-K model.

The necking speed in Zone $b$ is greater than in Zone $a$. The fracture criterion can be expressed as follows (Barata da Rocha et al., 1985; Graf and Hosford, 1990):

$C_i=\frac{d{\epsilon }_{rb,i}}{d{\epsilon }_{ra,i}}>10(57)$

By discretising, the following is derived:

$C_i=\frac{{\epsilon }_{rb,i}-{\epsilon }_{rb,i-1}}{{\epsilon }_{ra,i}-{\epsilon }_{ra,i-1}}>10(58)$

Once the fracture criterion is fulfilled, the material is deemed to be necking.

2.4.1 The Imperfection Factor

The initial imperfection factor is defined as follows:

$f_0=\frac{t_{b_0}}{t_{a_0}}(59)$

where $a$, $b$ are the subscripts indicating Zone $a$ and $b$, and $t_{b_0}$, $t_{a_0}$ are the initial thickness of these zones. $f_0$ is a material constant to describe the defect in the material and can be obtained by regression. The imperfection factor is defined as follows:

$f_i=\frac{t_{b,i}}{t_{a,i}}(60)$

Subsequently the relationship between $f_{\text{i}}$ and $f_0$ can be derived as follows:

${\epsilon }_{ra,i}=\ln \left(\frac{t_{a,i}}{t_{a_0}}\right)\Rightarrow t_{a,i}=t_{a_0}\exp \left({\epsilon }_{ra,i}\right)(61)$

${\epsilon }_{rb,i}=\ln \left(\frac{t_{b,i}}{t_{b_0}}\right)\Rightarrow t_{b,i}=t_{b_0}\exp \left({\epsilon }_{rb,i}\right)(62)$

$f_i=\frac{t_{b,i}}{t_a,i}=f_0\exp \left({\epsilon }_{rb,i}-{\epsilon }_{ra,i}\right)(63)$

2.4.2 Strain and Stress Analysis in Zone a

Zone $a$ is applied by external loading. Thus, its strain and stress state are the same as the macroscopic geometric model in the previous section:

${\epsilon }_{\theta a,i}={\epsilon }_{\theta ,i},{\epsilon }_{za,i}={\epsilon }_{z,i},{\epsilon }_{ra,i}={\epsilon }_{r,i}$

${\sigma }_{\theta a,i}={\sigma }_{\theta ,i},{\sigma }_{za,i}={\sigma }_{z,i}$

In order to calculate the strain and stress state in Zone $b$, the compatibility and equilibrium relations are applied in the $ntr$ coordinates. Thus, it is necessary to transform Zone$a$’s strain and stress components from $z\theta r$ coordinates to $ntr$ coordinates. Matrix notations are used to describe the transformation (Ganjiani and Assempour, 2008).

${\mathit{\sigma }}_{a,i}^{ntr}={\mathbf{T}}_{\mathbf{\text{i}}}{\mathit{\sigma }}_{a,i}^{z\theta r}{\mathbf{T}}_i^T(64)$

${\mathit{\epsilon }}_{a,i}^{ntr}={\mathbf{T}}_{\mathbf{\text{i}}}{\mathit{\epsilon }}_{a,i}^{z\theta r}{\mathbf{T}}_i^T(65)$

where

${\mathit{\sigma }}_{a,i}^{z\theta r}=\left[\begin{array}{ccc}{\sigma }_{za,i}& 0& 0\\ 0& {\sigma }_{\theta a,i}& 0\\ 0& 0& 0\end{array}\right]\mathit{ }{\mathit{\epsilon }}_{a,i}^{z\theta r}=\left[\begin{array}{ccc}{\epsilon }_{za,i}& 0& 0\\ 0& {\epsilon }_{\theta a,i}& 0\\ 0& 0& {\epsilon }_{ra,i}\end{array}\right]$

${\mathit{\sigma }}_{a,i}^{ntr}=\left[\begin{array}{ccc}{\sigma }_{na,i}& {\sigma }_{nta,i}& 0\\ {\sigma }_{tna,i}& {\sigma }_{ta,i}& 0\\ 0& 0& 0\end{array}\right]\mathit{ }{\mathit{\epsilon }}_{a,i}^{ntr}=\left[\begin{array}{ccc}{\epsilon }_{na,i}& {\epsilon }_{nta,i}& 0\\ {\epsilon }_{tna,i}& {\epsilon }_{ta,i}& 0\\ 0& 0& {\epsilon }_{ra,i}\end{array}\right]$

${\mathbf{T}}_{\mathbf{\text{i}}}=\left[\begin{array}{ccc}cos\left({\phi }_i\right)& sin\left({\phi }_i\right)& 0\\ -sin\left({\phi }_i\right)& cos\left({\phi }_i\right)& 0\\ 0& 0& 1\end{array}\right]$

${\phi }_i$ is the angle between axes $n$ and $z$, which changes during the loading process.

As shown in Figure 5, ${\phi }_i$ can be calculated from the geometry relation:

$\text{tan}\left({\phi }_i\right)=\frac{l_{z,i}}{l_{\theta ,i}}=\frac{l_{z0}}{l_{\theta 0}}exp\left(\frac{{\epsilon }_{za,i}}{{\epsilon }_{\theta a,i}}\right)=\text{tan}\left({\phi }_0\right)exp\left(\frac{{\epsilon }_{za,i}}{{\epsilon }_{\theta a,i}}\right)(66)$

FIGURE 5

FIGURE 5. Top views of the M-K model in (A) undeformed stage and (B) deformed stage.

Note that the width of Zone $b$ is neglected.

2.4.3 Strain and Stress Analysis in Zone b

The stress and strain states in Zone $b$ are calculated from Zone $a$ through compatibility of strain and force equilibrium:

${\epsilon }_{tb,i}={\epsilon }_{ta,i}(67)$

$d{\epsilon }_{tb,i}={\epsilon }_{tb,i}-{\epsilon }_{tb,i-1}(68)$

${\sigma }_{nb,i}=\frac{{\sigma }_{na,i}}{f_i}=\frac{{\sigma }_{na,i}}{f_0\text{exp}\left({\epsilon }_{rb,i}-{\epsilon }_{ra,i}\right)}=F_1\left({\epsilon }_{rb,i}\right)(69)$

${\sigma }_{ntb,i}=\frac{{\sigma }_{nta,i}}{f_i}=\frac{{\sigma }_{nta,i}}{f_0\text{exp}\left({\epsilon }_{rb,i}-{\epsilon }_{ra,i}\right)}=F_2\left({\epsilon }_{rb,i}\right)(70)$

where the capital $F$ represents the functional relation with unknowns in its bracket. In addition to the stress components, all the strain and strain increment components can also be expressed in the form of $F\left({\epsilon }_{rb,i}\right)$:

${\epsilon }_{nb,i}=-\left({\epsilon }_{tb,i}+{\epsilon }_{rb,i}\right)=F_3\left({\epsilon }_{rb,i}\right)(71)$

$d{\epsilon }_{nb,i}=F_3\left({\epsilon }_{rb,i}\right)-{\epsilon }_{nb,i-1}=F_4\left({\epsilon }_{rb,i}\right)(72)$

$d{\epsilon }_{rb,i}={\epsilon }_{rb,i}-{\epsilon }_{rb,i-1}=F_5\left({\epsilon }_{rb,i}\right)(73)$

Eq. 23 in Zone $b$ can be transformed as follows:

$\begin{array}{l}\frac{d{\epsilon }_{nb,i}+d{\epsilon }_{tb,i}-2\text{d}{\epsilon }_{ntb,i}\text{sin}\left(2{\phi }_i\right)+d{\epsilon }_{nb,i}\text{cos}\left(2{\phi }_i\right)-d{\epsilon }_{tb,i}\text{cos}\left(2{\phi }_i\right)}{r_{\theta }\left[{\text{σ}}_{\text{nb},\text{i}}+{\text{σ}}_{\text{tb},\text{i}}+\left({\text{σ}}_{\text{nb},\text{i}}-{\text{σ}}_{\text{tb},\text{i}}+2{\text{r}}_{\text{z}}{\text{σ}}_{\text{nb},\text{i}}-2{\text{r}}_{\text{z}}{\text{σ}}_{\text{tb},\text{i}}\right)\text{cos}\left(2{\text{φ}}_{\text{i}}\right)-\left(2{\text{σ}}_{\text{ntb},\text{i}}+4{\text{r}}_{\text{z}}{\text{σ}}_{\text{ntb},\text{i}}\right)\text{sin}\left(2{\text{φ}}_{\text{i}}\right)\right]}\\ =\frac{d{\epsilon }_{nb,i}+\text{d}{\epsilon }_{tb,i}+2d{\epsilon }_{ntb,i}\text{sin}\left(2{\phi }_i\right)-d{\epsilon }_{nb,i}\text{cos}\left(2{\phi }_i\right)+\text{d}{\epsilon }_{tb,i}\text{cos}\left(2{\phi }_i\right)}{r_z\left[{\sigma }_{nb,i}+{\sigma }_{tb,i}+\left(-{\sigma }_{nb,i}+{\sigma }_{tb,i}-2r_{\theta }{\sigma }_{nb,i}+2r_{\theta }{\sigma }_{tb,i}\right)\text{cos}\left(2{\phi }_i\right)+\left(2{\sigma }_{ntb,i}+4r_{\theta }{\sigma }_{ntb,i}\right)\text{sin}\left(2{\phi }_i\right)\right]}\\ =\frac{d{\epsilon }_{nb,i}\sin \left(2{\phi }_i\right)-\text{d}{\epsilon }_{tb,i}\sin \left(2{\phi }_i\right)+2\text{d}{\epsilon }_{ntb,i}\cos \left(2{\phi }_i\right)}{\left(2r_{z\theta }+1\right)\left(r_{\theta }+r_z\right)\left(2{\sigma }_{ntb,i}\cos \left(2{\phi }_i\right)+{\sigma }_{nb,i}\sin \left(2{\phi }_i\right)-{\sigma }_{tb,i}\sin \left(2{\phi }_i\right)\right)}\end{array}(74)$

There are only two equivalence relations and three unknowns (${\sigma }_{tb,i},d{\epsilon }_{ntb,i},{\epsilon }_{rb,i}$) in Eq. 74. Thus, $\text{d}{\epsilon }_{ntb,i}$ and ${\sigma }_{tb,i}$ cannot be solved explicitly and have to be expressed by ${\epsilon }_{rb,i}$:

$d{\epsilon }_{ntb,i}=F_6\left({\epsilon }_{rb,i}\right)(75)$

${\sigma }_{tb,i}=F_7\left({\epsilon }_{rb,i}\right)(76)$

The equivalent stress and strain increment can be calculated on non-principal axes of anisotropy to make use of the constitutional relationship. The transformed expression of equivalent stress and strain increment components on the axes $ntr$ is derived. The stress and strain components in $ntr$ coordinates are as follows:

${\mathit{\sigma }}^{ntr}={\mathbf{T}}_{\text{i}}{\mathit{\sigma }}^{z\theta r}{\mathbf{T}}_i^T\to {\mathit{\sigma }}^{z\theta r}={\mathbf{T}}_i^T{\mathit{\sigma }}^{ntr}{\mathbf{T}}_{\mathbf{\text{i}}}(77)$

$\mathit{d}{\mathit{\epsilon }}^{ntr}={\mathbf{T}}_{\text{i}}\mathit{d}{\mathit{\epsilon }}^{z\theta r}{\mathbf{T}}_i^T\to \mathit{d}{\mathit{\epsilon }}^{z\theta r}={\mathbf{T}}_i^T\mathit{d}{\mathit{\epsilon }}^{ntr}{\mathbf{T}}_{\mathbf{\text{i}}}(78)$

where the strain and stress are in matrix form:

${\mathit{\sigma }}_{b,i}^{z\theta r}=\left[\begin{array}{ccc}{\sigma }_{zb,i}& {\sigma }_{z\theta b,i}& 0\\ {\sigma }_{\theta zb,i}& {\sigma }_{\theta b,i}& 0\\ 0& 0& 0\end{array}\right]\mathit{d}{\mathit{\epsilon }}_{b,i}^{z\theta r}=\left[\begin{array}{ccc}d{\epsilon }_{zb,i}& d{\epsilon }_{z\theta b,i}& 0\\ d{\epsilon }_{\theta zb,i}& d{\epsilon }_{\theta b,i}& 0\\ 0& 0& d{\epsilon }_{rb,i}\end{array}\right]$

${\mathit{\sigma }}_{b,i}^{ntr}=\left[\begin{array}{ccc}{\sigma }_{nb,i}& {\sigma }_{ntb,i}& 0\\ {\sigma }_{tnb,i}& {\sigma }_{tb,i}& 0\\ 0& 0& 0\end{array}\right]\mathit{d}{\mathit{\epsilon }}_{b,i}^{ntr}=\left[\begin{array}{ccc}d{\epsilon }_{nb,i}& d{\epsilon }_{ntb,i}& 0\\ d{\epsilon }_{tnb,i}& d{\epsilon }_{tb,i}& 0\\ 0& 0& d{\epsilon }_{rb,i}\end{array}\right]$

${\mathbf{T}}_{\mathbf{\text{i}}}=\left[\begin{array}{ccc}\mathit{cos}\left({\phi }_i\right)& \mathit{sin}\left({\phi }_i\right)& 0\\ -\mathit{sin}\left({\phi }_i\right)& \mathit{cos}\left({\phi }_i\right)& 0\\ 0& 0& 1\end{array}\right]$

The solutions are shown in vector form:

${\mathit{\sigma }}_{b,i}=\left[\begin{array}{c}{\sigma }_{zb,i}\\ {\sigma }_{\theta b,i}\\ {\sigma }_{z\theta b,i}\end{array}\right]=\left[\begin{array}{ccc}\text{cos}{\left({\phi }_i\right)}^2& \text{sin}{\left({\phi }_i\right)}^2& -\text{sin}\left(2{\phi }_i\right)\\ \text{sin}{\left({\phi }_i\right)}^2& \text{cos}{\left({\phi }_i\right)}^2& \text{sin}\left(2{\phi }_i\right)\\ \frac{\text{sin}\left(2{\phi }_i\right)}{2}& -\frac{\text{sin}\left(2{\phi }_i\right)}{2}& \text{cos}\left(2{\phi }_i\right)\end{array}\right]\left[\begin{array}{c}{\sigma }_{nb,i}\\ {\sigma }_{tb,i}\\ {\sigma }_{ntb,i}\end{array}\right]={\mathit{Q}}_{\mathit{i}}{\mathit{\sigma }}_{b,i}^{\prime }(79)$

$\mathit{d}{\mathit{\epsilon }}_{b,i}=\left[\begin{array}{c}d{\epsilon }_{zb,i}\\ d{\epsilon }_{\theta b,i}\\ 2d{\epsilon }_{z\theta b,i}\end{array}\right]=\left[\begin{array}{ccc}\text{cos}{\left({\phi }_i\right)}^2& \text{sin}{\left({\phi }_i\right)}^2& -\text{sin}\left(2{\phi }_i\right)\\ \text{sin}{\left({\phi }_i\right)}^2& \text{cos}{\left({\phi }_i\right)}^2& \text{sin}\left(2{\phi }_i\right)\\ \frac{\text{sin}\left(2{\phi }_i\right)}{2}& -\frac{\text{sin}\left(2{\phi }_i\right)}{2}& \text{cos}\left(2{\phi }_i\right)\end{array}\right]\left[\begin{array}{c}d{\epsilon }_{nb,i}\\ d{\epsilon }_{tb,i}\\ 2d{\epsilon }_{ntb,i}\end{array}\right]={\mathit{Q}}_{\mathit{i}}\mathit{d}{\mathit{\epsilon }}_{b,i}^{\prime }(80)$

Note that ${\sigma }_{b,i}^{\prime }$ and $\mathit{d}{\mathit{\epsilon }}_{b,i}^{\prime }$ can be written in the form of $F\left({\epsilon }_{rb,i}\right)$:

${\mathit{\sigma }}_{b,i}^{\prime }=\left[\begin{array}{c}{F}_1\left({\epsilon }_{rb,i}\right)\\ F_7\left({\epsilon }_{rb,i}\right)\\ F_2\left({\epsilon }_{rb,i}\right)\end{array}\right]\mathit{ d}{\mathit{\epsilon }}_{b,i}^{\prime }=\left[\begin{array}{c}{F}_4\left({\epsilon }_{rb,i}\right)\\ {\epsilon }_{tb,i}-{\epsilon }_{tb,i-1}\\ 2F_6\left({\epsilon }_{rb,i}\right)\end{array}\right]$

Thus, by substituting Eqs 79, 80 into Eqs 24, 28, respectively, the equivalent stress and strain increment can be expressed in $ntr$ axes:

${\overline{\sigma }}_{b,i}=\sqrt{\left(\mathbf{p}{\mathit{Q}}_{\mathit{i}}\mathit{\sigma }{\prime }_{b,i}\right)\cdot {\mathit{Q}}_{\mathit{i}}\mathit{\sigma }{\prime }_{b,i}}=F_8\left({\epsilon }_{rb,i}\right)(81)$

$d{\overline{\epsilon }}_{b,i}=\sqrt{\left({\mathbf{p}}^{-1}{\mathit{Q}}_{\mathit{i}}\mathit{d\epsilon }{\prime }_{b,i}\right)·{\mathit{Q}}_{\mathit{i}}\mathit{d\epsilon }{\prime }_{b,i}}=F_9\left({\epsilon }_{rb,i}\right)(82)$

The total strain is as follows:

${\overline{\epsilon }}_{b,i}={\overline{\epsilon }}_{b,i-1}+d{\overline{\epsilon }}_{b,i}=F_{10}\left({\epsilon }_{rb,i}\right)(83)$

The constitutional relationship in Eq. 55 can be transformed as follows:

$F_8\left({\epsilon }_{rb,i}\right)={\sigma }_0+K{F}_{10}{\left({\epsilon }_{rb,i}\right)}^n(84)$

It is clear that Eq. 84 is a non-linear equation with ${\epsilon }_{rb,i}$ as the only unknown. Solving ${\epsilon }_{rb,i}$ and substituting the result into all of the $F$ functions, the stress and strain state in zone $b$ can be determined.

2.5 Numerical Process

All the equations necessary for the failure prediction model have been derived in the previous sections and the procedure of utilizing the model is demonstrated in the flow chart of Figure 6. The core of the model requires the solution of three non-linear equations, Eq 38, 51, 84. The equations are simplified by eliminating to only one unknown that can be solved numerically.

FIGURE 6

FIGURE 6. Flow chart of the failure prediction model.

3 Case Studies

3.1 Case Description

In the case studies, the derived model was applied on two tubes made from annealed C26800 zinc copper and AISI 1215 carbon steel, the material properties and geometry parameters of which are shown in Table 1 and Table 2 (Hwang and Wang, 2009). $r_{z\theta }$ is not a commonly used parameter as it is challenging to experimentally determine. In the present model, $r_{z\theta }$ is evaluated as the average of $r_{\theta }$ and $r_z$ as it is found to negligibly affect the results, which will be explained in the next section.

TABLE 1

TABLE 1. Mechanical properties of C26800 and AISI 1215.

TABLE 2

TABLE 2. Geometry parameters of the tube.

3.2 Results and Discussion

Hwang and Wang (2009) conducted hydro-bulging tests on C26800 and AISI 1215. They provided the evolutions of bulging height, inner pressure, and pole point thickness. The last bulging height recorded for each test was taken as the BHL. Thus, the BHL of C26800 and AISI 1215 are 11.1 and 5.7 mm, respectively. The failures in the experiments are developed along the axial direction. Therefore, Figure 7 used the derived model to predict the BHLs for both tubes under different $f_0$ and fixed ${\phi }_0=90°$. It is found that when $f_0$ are 0.932 and 0.982, respectively, the BHL predictions are corresponding to the experimental result. The factor can be used in other studies that are associated with the M-K model, such as a post-FE failure prediction module (Gao et al., 2017).

FIGURE 7

FIGURE 7. Effect of the initial imperfection factor ${\mathit{f}}_0$ on the BHL.

After determining $f_0$ of each tube, Figure 8 studied the influence of ${\phi }_0$ to explain why the failure appeared along the axial direction. In Figure 8, both curves of C26800 and AISI 1215 undergo an increase, a plateau, and a decrease. Both curves reach their minimum at ${\phi }_0=90°$, which suggests that the failure should appear in the axial direction first and agrees with the experimental result. In fact, in experiments on many other materials, crack formations are also along the tube axis (Mori et al., 2007; Zhu et al., 2020). Thus, when utilizing the model, it is reasonable to cease the iteration of ${\phi }_0$ and set it to a constant value. Once ${\phi }_0=90°$, all the terms including $r_{z\theta }$ are eliminated in the model. Consequently, the value of $r_{z\theta }$ is not an essential factor and the use of an average of $r_{\theta }$ and $r_z$ in the previous section is a reasonable assumption.

FIGURE 8

FIGURE 8. BHL under changing.${\mathit{\phi }}_0$

Figure 9 and Figure 10 predict the required pressure and pole point thickness evolution during the predictions under given $f_0$ and ${\phi }_0$ and then compare the data with experimental results. In Figure 9, AISI 1215 shows a faster growth and a higher BHL than C26800. In Figure 10, both tubes undergo a steady pole point thickness decrease. The average relative error for the required pressure and the pole point thickness are 3.81% and 5.75%, respectively. Thus, it is reasonable to say that the model can be used to reflect the process of hydro-bulging precisely.

FIGURE 9

FIGURE 9. Required pressure evolution during the hydro-bulging process.

FIGURE 10

FIGURE 10. Pole point thickness evolution during the hydro-bulging process.

4 Conclusion

A failure model to predict the bulging limit of the tube hydro-bulging test is needed to evaluate the formability of a tubular material without conducting real tests. In the present study, a failure prediction model for the tube hydro-bulging test was developed by the combination of Hill’s orthogonal anisotropic model, the geometry model, and the M-K model. The main conclusions can be summarized as follows:

• Given $f_0$, the model can predict the BHL of the tube hydro-bulging test. On the contrary, the model can also be used to fit $f_0$ once the BHL is known. Figure 7 is an example of the mapping relation between the BHL and $f_0$.

• By assuming that the necking can only appear along the axial direction, the model can be simplified by stopping the iteration of ${\phi }_0$ and setting ${\phi }_0=90°$.

• In the case studies, the predictions of the required pressure and the pole point thickness evolution demonstrate marginal errors compared with the experimental results, which are 3.81% and 5.75%, respectively. This verified the validity of the model.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author Contributions

BD is the main developer of the model and the first author of the article. HL contributed to the discussion part and figures. SD contributed to the establishment of the model. KW contributed to the establishment of the model. XL is the corresponding author of the article. He guided BD in the model derivation and article writing. DP contributed to the checking of the model and its improvement.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Aue-u-lan, Y. (2007). Hydroforming of Tubular Materials at Various Temperatures Dissertation/doctoral Thesis. Columbus (OH): The Ohio State University.

Google Scholar

Barata da Rocha, A., Barlat, F., and Jalinier, J. M. (1985). Prediction of the Forming Limit Diagrams of Anisotropic Sheets in Linear and Non-linear Loading. Mater. Sci. Eng. 68, 151–164. doi:10.1016/0025-5416(85)90404-5

CrossRef Full Text | Google Scholar

Boudeau, N., and Malécot, P. (2012). A Simplified Analytical Model for Post-processing Experimental Results from Tube Bulging Test: Theory, Experimentations, Simulations. Int. J. Mech. Sci. 65, 1–11. doi:10.1016/j.ijmecsci.2012.08.002

CrossRef Full Text | Google Scholar

Fllice, L., Fratini, L., and Micari, F. (2001). A Simple Experiment to Characterize Material Formability in Tube Hydroforming. CIRP Ann. 50, 181–184. doi:10.1016/s0007-8506(07)62100-3

CrossRef Full Text | Google Scholar

Fuchizawa, S., Narazaki, M., and Yuki, H. (1993). “Bulge Test for Determining Stress-Strain Characteristics of Thin Tubes,” in Advanced Technology of Plasticity, 1, 488–493.

Google Scholar

Ganjiani, M., and Assempour, A. (2008). Implementation of a Robust Algorithm for Prediction of Forming Limit Diagrams. J. Materi Eng Perform 17, 1–6. doi:10.1007/s11665-007-9121-4

CrossRef Full Text | Google Scholar

Gao, H., El Fakir, O., Wang, L., Politis, D. J., and Li, Z. (2017). Forming Limit Prediction for Hot Stamping Processes Featuring Non-isothermal and Complex Loading Conditions. Int. J. Mech. Sci. 131-132, 792–810. doi:10.1016/j.ijmecsci.2017.07.043

CrossRef Full Text | Google Scholar

Graf, A., and Hosford, W. F. (1990). Calculations of Forming Limit Diagrams. Metallurgical Transactions A 21, 87–94. doi:10.1007/BF02656427

CrossRef Full Text | Google Scholar

He, Z., Yuan, S., Lin, Y., Wang, X., and Hu, W. (2014a). Analytical Model for Tube Hydro-Bulging Test, Part I: Models for Stress Components and Bulging Zone Profile. Int. J. Mech. Sci. 87, 297–306. doi:10.1016/j.ijmecsci.2014.05.009

CrossRef Full Text | Google Scholar

He, Z., Yuan, S., Lin, Y., Wang, X., and Hu, W. (2014b). Analytical Model for Tube Hydro-Bulging Tests, Part II: Linear Model for Pole Thickness and its Application. Int. J. Mech. Sci. 87, 307–315. doi:10.1016/j.ijmecsci.2014.05.010

CrossRef Full Text | Google Scholar

Hill, R. (1948). A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 193, 281–297. doi:10.1098/rspa.1948.0045

CrossRef Full Text | Google Scholar

Hill, R. (1950). The Mathematical Theory of Plasticity. Oxford: Clarendon.

Google Scholar

Hutchinson, J. W., and Neale, K. W. (1978). “Sheet Necking-II. Time-independent Behavior,” in Mechanics of Sheet Metal Forming: Material Behavior and Deformation Analysis. Editors D. P. Koistinen, and N. Wang (Boston, MA: Springer US), 127–153. doi:10.1007/978-1-4613-2880-3_6

CrossRef Full Text | Google Scholar

Hwang, Y.-M., and Lin, Y.-K. (2006). Analysis of tube bulge forming in an open die considering anisotropic effects of the tubular material. Int. J. Mach. Tools Manuf. 46, 1921–1928. doi:10.1016/j.ijmachtools.2006.01.025

CrossRef Full Text | Google Scholar

Hwang, Y.-M., and Lin, Y.-K. (2007). Evaluation of Flow Stresses of Tubular Materials Considering Anisotropic Effects by Hydraulic Bulge Tests. J. Eng. Mater. Technol. 129, 414–421. doi:10.1115/1.2744406

CrossRef Full Text | Google Scholar

Hwang, Y. M., and Wang, C. W. (2009). Flow Stress Evaluation of Zinc Copper and Carbon Steel Tubes by Hydraulic Bulge Tests Considering Their Anisotropy. J. Mater. Process. Technol. 209, 4423–4428. doi:10.1016/j.jmatprotec.2008.10.033

CrossRef Full Text | Google Scholar

Imaninejad, M., Subhash, G., and Loukus, A. (2004). Influence of End-Conditions during Tube Hydroforming of Aluminum Extrusions. Int. J. Mech. Sci. 46, 1195–1212. doi:10.1016/j.ijmecsci.2004.08.001

CrossRef Full Text | Google Scholar

International Standard Organization (2021). Metallic Materials — Determination of Forming-Limit Curves for Sheet and Strip — Part 2: Determination of Forming-Limit Curves in the Laboratory.

Google Scholar

Kuwabara, T., Yoshida, K., Narihara, K., and Takahashi, S. (2005). Anisotropic Plastic Deformation of Extruded Aluminum Alloy Tube under Axial Forces and Internal Pressure. Int. J. Plasticity 21, 101–117. doi:10.1016/j.ijplas.2004.04.006

CrossRef Full Text | Google Scholar

Marciniak, Z., and Kuczyński, K. (1967). Limit Strains in the Processes of Stretch-Forming Sheet Metal. Int. J. Mech. Sci. 9, 609–620. doi:10.1016/0020-7403(67)90066-5

CrossRef Full Text | Google Scholar

Mohr, D., Dunand, M., and Kim, K.-H. (2010). Evaluation of Associated and Non-associated Quadratic Plasticity Models for Advanced High Strength Steel Sheets under Multi-Axial Loading. Int. J. Plasticity 26, 939–956. doi:10.1016/j.ijplas.2009.11.006

CrossRef Full Text | Google Scholar

Mori, K., Maeno, T., and Maki, S. (2007). Mechanism of Improvement of Formability in Pulsating Hydroforming of Tubes. Int. J. Mach. Tools Manuf. 47, 978–984. doi:10.1016/j.ijmachtools.2006.07.006

CrossRef Full Text | Google Scholar

Sokolowski, T., Gerke, K., Ahmetoglu, M., and Altan, T. (2000). Evaluation of Tube Formability and Material Characteristics: Hydraulic Bulge Testing of Tubes. J. Mater. Process. Technol. 98, 34–40. doi:10.1016/s0924-0136(99)00303-9

CrossRef Full Text | Google Scholar

Strano, M., and Altan, T. (2004). An Inverse Energy Approach to Determine the Flow Stress of Tubular Materials for Hydroforming Applications. J. Mater. Process. Technol. 146, 92–96. doi:10.1016/j.jmatprotec.2003.07.016

CrossRef Full Text | Google Scholar

Swift, H. W. (1952). Plastic Instability under Plane Stress. J. Mech. Phys. Solids 1, 1–18. doi:10.1016/0022-5096(52)90002-1

CrossRef Full Text | Google Scholar

Zhu, H., He, Z., Lin, Y., Zheng, K., Fan, X., and Yuan, S. (2020). The Development of a Novel Forming Limit Diagram under Nonlinear Loading Paths in Tube Hydroforming. Int. J. Mech. Sci. 172, 105392. doi:10.1016/j.ijmecsci.2019.105392

CrossRef Full Text | Google Scholar