Towards sensitive prediction of wrinkling instability in sheet metal forming by introducing evolution of triple nonlinearity: Tube forming

Journal Pre-proof

Towards sensitive prediction of wrinkling instability in sheet metal forming by introducing evolution of triple nonlinearity: tube forming Heng Li , Hao Ran Liu , Nan Liu , Hong Sun , Heng Yang , Bi Ying Liu PII: DOI: Reference:

S0020-7403(19)31198-1 https://doi.org/10.1016/j.ijmecsci.2019.105054 MS 105054

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

8 April 2019 8 July 2019 30 July 2019

Please cite this article as: Heng Li , Hao Ran Liu , Nan Liu , Hong Sun , Heng Yang , Bi Ying Liu , Towards sensitive prediction of wrinkling instability in sheet metal forming by introducing evolution of triple nonlinearity: tube forming, International Journal of Mechanical Sciences (2019), doi: https://doi.org/10.1016/j.ijmecsci.2019.105054

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights •

A hybrid integrated framework of wrinkling prediction in sheet forming by introducing evolution of triple nonlinearity.

•

Transferring features of buckling modes under less and more degrees of constraints in sheet forming operations.

•

Superposing of multiple geometric micro-imperfections (GMI) for sheet metal forming with less and more constraints.

•

Considering anisotropic/asymmetric properties (AAP) for wrinkling predication closely related to compressive deformation.

•

Coupling role of GMI and AAP on wrinkling prediction capability in sheet metal forming with less and more constraints.

Towards sensitive prediction of wrinkling instability in sheet metal forming by introducing evolution of triple nonlinearity: tube forming Heng Li *a, Hao Ran Liu a, Nan Liu b, Hong Sun a, Heng Yang a, Bi Ying Liu a a School of Materials Science and Engineering, Northwestern Polytechnical University, Xi 'an 710072 b School of Materials Science and Engineering, Xi‟An University of technology E-mail address: [email protected] (H. Li)

Abstract: Wrinkling instability is one of key defects restricting the formability of sheet metallic materials How to sensitively predict the wrinkling occurring is one of the most fundamental and challenging issues in sheet metal forming operations to ensure defect-free deformation. In the study, taking typical axial compression and draw bending of tubular materials as the cases, the buckling modes of sheet metal forming under less and more degrees of constraints are experimentally compared, then a hybrid integrated prediction framework of wrinkling in sheet metal forming is constructed by introducing evolution of triple nonlinearity in geometry, material and boundary conditions, viz., multiple geometric micro-imperfections (GMI), anisotropic/asymmetric properties (AAP) and explicit FE algorithm, respectively. The results show that transferred wrinkling modes may occur for less constrained process, while the multiple constraints make the wrinkling occur in a higher order buckling mode consuming more energy. Multiple GMI modes constructed by linear buckling analysis are superimposed into explicit FE models, viz., the first three buckling modes superposed for axial compression, while the higher order buckling mode analogous to actual wrinkling selected for draw bending. Three constitutive models, viz, Mises, normal anisotropic Hill‟48-r and stress invariant based anisotropic/asymmetric models (SIM), are comparatively implemented in the explicit FE models. Compared with the perfect model without introducing GMI, the GMI-implemented FE model provides more sensitive prediction of wrinkling for both forming cases. The SIM-implemented FE model predicts the instability more sensitively than other constitutive models in tension-compression asymmetry dominated bending case. By coupling the GMI and SIM into explicit FE models, the desired wrinkling prediction capability for both forming cases is further confirmed since evolution of triple nonlinearity is simultaneously considered. Keywords: Wrinkling instability; Sheet metal forming; Geometric micro-imperfection; Anisotropy; Asymmetry; Contact conditions; Tube compression; Tube bending

1. Introduction The urgent demands of heavy loading, low consumption and long life of high-end equipment in aviation, aerospace and automobile call for the precision forming of high-performance lightweight thin-walled sheet metal components. As one of the key lightweight components for “bleeding” transforming or loading carrying, tubular parts are widely used in the above industries with increasing needs for production of high-strength and thin-walled parts[1,2]. However, non-uniform tension-compression stress/strain states during forming operations easily induce the wrinkling instability for thin-walled metallic materials. It is known that the initiation and growth of wrinkles are interactively influenced by multiple factors including geometric parameters, material properties and boundary conditions, and sheet metal forming is generally a dynamic tri-nonlinear physical process in geometry, material and boundary conditions. Especially, the forming of thin-walled high-strength materials generally undergoes more complicated loading conditions, which involves higher order tri-nonlinearity and causes more confused mechanisms of wrinkling initiation and propagation. Thus, wrinkling instability has nowadays been being one of the key defects restricting formability of thin-walled sheet metallic materials, and better understanding the wrinkling mechanism in sheet metal forming and accurate prediction of the wrinkling is one of the most fundamental and challenging issues for achieving defect-free deformation processing. Much effort has been taken on the study of wrinkling mechanisms, predictions and controlling in sheet metal forming by using analytical, numerical and experimental methods. With regards to analytical approach, static equilibrium methods were early proposed based on the differential equilibrium equation involving small deflection consumption or non-linear large deflection theory[3-5]. The energy methods were also raised by comparing the internal energy of buckled parts and the external work to judge whether it was a stable condition and could be utilized to solve massive wrinkling problems of thin shells under simple boundary conditions[6-8] and complex forming process with multiple constraints[9]. To predict and control the wrinkling quantitatively in unsupported region for thin-walled shells with curved surface, a theoretical model on critical wrinkling stress was proposed by considering proper “reverse bulging effect” based on energy method[10]. Kong et al. proposed a theoretical flange wrinkling prediction method based on the plastic buckling theory and the energy method with the boundary conditions extracted from the FE simulation results[11]. An investigation into material wrinkling failure mechanics of conventional metal spinning and

the effects of process parameters and material properties are presented by Watson et al[12]. However, the analytical methods are of limited applications due to many simplifications and assumptions, which are difficult and inaccurate in dealing with complicated cases. As for numerical methods, much progress is gained in recent years for accurate wrinkling prediction and sheet metal forming optimization, including the implicit methods and explicit methods. Liu et al. proposed an imperfection-based perturbation method for plastic wrinkling prediction in tube bending under multi-die constraints. By providing an over prediction of wrinkling, a “lower bound” forming conditions are obtained under multi-die constraints[13]. The implicit analysis owns abilities to catch the bifurcation point along the loading path by checking the determinant of stiff matrix during iteration procedure[14-16]. Nevertheless, implicit method is hard to deal with the convergence problems when facing the complex boundary conditions[17,18]. Comparatively, explicit models are easy to overcome the convergence problems due to its solution algorithm and were widely used in complicated forming occasions[19,20]. However, the explicit algorithm based FE simulation generally underestimates the wrinkling risks especially for forming operations under multiple constraints due to the absence of the singularity check for tangent stiffness. Also, it was put that the triggering of the out-plane instability could depend on the relative small accumulating numerical error[7]. Thus, it is of significance to realize sensitive wrinkling prediction under multiple constraints. As for the experimental ways, it is an intuitionistic way to understand the occurrence of wrinkling and its evolution, and much effort has been put in the experimental studies[21,22]. The wrinkling modes of tube axial compression were classified into seven categories and their different energy absorbing abilities were studied [23]. Abebe et al. carried out that a dimpling-and-wrinkling-free technique for order preferences by similarity to ideal solution is numerically verified by comparison with results of a full model simulation and experimentally validated by its application to a manufactured product[24]. Cao et al. presented a contact die to offer normal pressure to the buckling tests, and they found that the one-side contact could delay the initiation of buckling[25]. Simonetto et al. presented dynamic detection of wrinkling in tube bending process and it was capable to give out dynamic feedback and fast adjustment of the forming parameters during the process[26]. By utilizing the digital image correlation (DIC) analysis, Zhou et al. studied the wrinkling behavior of origami crash boxes[27]. The experiments could provide an intuitionistic understanding of wrinkling instability and help improve the precision of analytical model and numerical model.

In order to explain the difference between the analytical or numerical results, the researches mainly focused on the initial imperfections and nonlinearities. As for initial imperfections, it comprises of the dimensional deviation, thickness unevenness, residual stress and loading imperfection. The initial imperfections were studied about both its detection and construction. Considering the geometric imperfection, Donnell et al. reduced the discrepancy of limit load between experimental and analytical values[28]. By using the eigenvalue analysis, the buckling modes from low order to high order were applied to simulate the onset and growth of wrinkling[29-32]. However, the linear buckling analysis comprises no plasticity of the material and a majority of the researches are concerning on the critical point of buckling, few are about the post-buckling deformation or ultimate forming shape. Initial imperfections were also explored using random field model[33-35], which was based on the statistical characteristics from measured data. And by measuring the initial imperfections experimentally, Wheeler and Pircher studied the structural behavior of thin-walled tubes[36], and Bisagni studied buckling and post-buckling behavior of composite cylindrical shells under axial compression[37]. Orifici and Bisagni used single perturbation load analysis (SPLA) to model the imperfection pattern[38], and Paquette and Kyriakides used analytical expression to model the initial imperfection for accurate prediction of limit load and limit strain, turning out that the analysis fitted the experiments well when the amplitude of imperfection was 0.1%[39]. The above researches give out many novel and useful methods to characterize the imperfections, and among the methods, eigenvalue analysis is employed in this article to construct the initial imperfection with both sole and superposed buckling modes because of its easy achievements and convenient applications. As for nonlinear properties, the researches about geometric nonlinearities are much more than that of material nonlinearities and boundary nonlinearities. However, the nonlinear material properties and boundary conditions still influence the wrinkling behavior a lot. The plastic paradox showed that the buckling load based on simple J2 flow theory exceeds realistic values, but simple J2 deformation theory is in good agreement[40]. The reason for these discrepancies were thought to be the simple J2 flow theory, which is unable to describe the strain-stress relations accurately under inhomogeneous and complicated deformations when buckling occurs, especially for the anisotropic and asymmetric materials. Aiming at catering for these demands, many complicated and accurate constitutive models are raised by scholars[41-43]. As for titanium material focused in this study, the multi-pass thermal-mechanical processing and hexagonal close-packed (HCP) crystal structure result in obvious anisotropy and asymmetry[44]. Especially, it is noted that the compressive

instability is closely related to the compressive stress/strain states, while most researches were conducted based on sole uniaxial tension stress-strain constitutive relations. Therefore, it is necessary to introduce the anisotropy and asymmetry behaviors and their evolution in wrinkling analysis and prediction. Tube forming is a typical kind of sheet metal forming processes such as tube axial compression and draw bending. Tube axial compression with simple tools is easy to acquire complicated tubular parts with double walls. However, due to the poor structural stability of thin-walled tubes and weak constraints on deformation zones, multiple instable defects could happen before reaching the ultimate strength of material such as wrinkling and collapse[45,46]. Furthermore, the perturbations of material properties, geometric parameters and process parameters all could lead to buckling mode changing of the instabilities with same compression conditions[47,48]. With respect to tube bending methods, draw bending is attracting more and more applications since it is possible to achieve more precision and efficiency in bending radius and bending angle under multiple tools[49]. However, it is a complex dynamic process with inhomogeneous tension-compression deformation, which may induce wrinkling instability[50]. Thus, to solve the wrinkling instability is the most fundamental and challenging issue for both of the tube forming processes. The process of tube axial compression is with less constraints, while the process of tube draw bending is with more constrains. Thus in this paper, taking typical axial compression and draw bending of tubular materials as the cases, the wrinkling features about buckling modes under less and more constraints are experimentally explored. Then a hybrid integrated prediction framework of wrinkling in sheet metal forming is constructed by introducing evolution of triple nonlinearity in geometry, material and boundary conditions. Finally, the prediction capability of the above framework is evaluated when initial GMI and AAP are individually or simultaneously implemented into explicit 3D-FE models of the above tube forming cases.

2. Wrinkling occurring characteristics under less and more constraints 2.1 Experimental facilities and procedure Tube axial compression is regarded as a forming process with less boundary constraints and draw bending is of more boundary constraints. As shown in Fig. 1, tube axial compression is conducted under an axial force P supplied by simple tools. During the initial stage, the tube deforms uniformly and the compressive load is rising stably. The wrinkling instability happens with the load going to fall and the maximum load achieved in the process is called the

limit load, which is an important index to measure the resistance of tubes to instability. As shown in Fig. 2, tube draw bending consists of multiple dies, including bend die, clamp die, pressure die, wiper die and flexible mandrel (including shank and balls). After clamped by the clamp die and bend die, the tube goes past the tangent point and rotates along the groove of bend die with designed bending radius. When the bending angle is reached, the mandrel is withdrawn from working position and then the dies are removed to unload the tube.

Fig. 1 Tube axial compression of (a) schematic diagram and (b) forming dies with less constraints

Fig. 2 Tube draw bending of (a) schematic diagram and (b) forming tools with more constraints

The test specimens for axial compression and draw bending are made of commercial pure CP3 titanium with the thickness and diameter of 1.07 mm and 76.2 mm, respectively. Regarding axial compression, three categories of tubes with different heights were used, i.e. 40 mm, 80 mm, and 120 mm. For further studying the influence of material properties on wrinkling behaviors of tube during axial compression, 6061-T4 aluminum tubes with the same geometrical specification to the above titanium alloy ones were also considered for comparison. Fig. 3 shows the plastic stress-strain curves of tension and compression of

titanium tubes and tension curves for 6061-T4 aluminum tube.

The yield stresses of CP3 Ti

are 372 MPa and 374 MPa for tension and compression, respectively. The tension yield stress of 6061 T4 Al is 157 MPa.

Fig. 3 Plastic stress-strain curves of tubular materials

Axial compression tests were carried out with an upper die and a lower die fixed with a 200KN universal testing machine. All tests are under displacement control with the speed value of 2mm/min, repeating three times for each specification. Considering the efficiency and cost, the first two tests were compressed with the height to gain at least one whole wrinkling wave, and the third test was compressed with the height of no more wrinkles occurring. The tube draw bending is performed on servo-hydraulic bender machine, which can detect the bending speed and relative pushing speed dynamically. The bending radius, bending angle, bending speed, and push assistant speed of pressure die is 152.4 mm, 90º, 0.8 rad/s and 121 mm/s, respectively. The length of pressure die, mandrel, wiper die, clamp die, diameter of mandrel, diameter of balls, thickness of balls and pitch of balls are 700 mm, 360 mm，350 mm, 305 mm, 73.75 mm, 73.75 mm, 22.4 mm and 25.2 mm, respectively. The lubricant is imposed on the wiper die, mandrel and flexible balls only.

2.2 Experimental results and discussion 2.2.1 Tube axial compression under less constraints Fig. 4 and Fig. 5 are the experimental results with two visual views of tube axial compression for titanium tubes and aluminum tubes, respectively. Fig.4 shows the experimental results of axial compression forming of Φ76.2×t1.07 CP3 titanium tubes with different height. For the results of 40mm height tubes, the forming quality of No.1 is the best with no obvious wrinkling while the other two both have severe wrinkling. For the results of

80mm height tubes, No.1 and No.2 have four wrinkling and No.3 has three. For the results of 120mm height tubes, No.1 and No.2 have three wrinkling and No.3 has four. Fig.5 shows the experimental results of axial compression forming of Φ76.2×t1.07 6061 aluminum tubes with different height. For the results of 40mm height tubes, the first two tubes have good forming quality while the third one wrinkled. For the results of 80mm height tubes, the forming quality of No.1 is the best with no obvious wrinkling but the other two both have severe wrinkling at bottoms. For the results of 120mm height tubes, No.1 and No.2 have three wrinkling and No.3 has four. It can be found that the appearing of wrinkling defects has some degree of randomness. The detailed information about the wrinkling mode and limit load is listed in Table 1 and Table 2. It is worth noting that transferred buckling modes may occur even for the same geometric parameters and material types, which was also mentioned by Al Galib and Limam[47]. This transferred performance is assumed to be caused by many uncertainties in microstructure distributions, GMI, contact and boundary conditions or loading imperfections, of which GMI are considered to be the main cause and assumed to be equivalent index representing the above uncertainties in this paper. The wrinkling mode is found to transfer from axisymmetric mode to diamond mode with increasing tube height. The reason may be the relatively weakened constraints from clamp ends and the accumulating effects of initial GMI with increasing tube height. With respect to effects of the material properties on wrinkling mode, the numbers of the axisymmetric and diamond modes are compared between the aluminum tubes and titanium tubes. There are totally 8 diamond modes occurring for titanium tubes and 6 diamond modes for aluminum tubes, indicating that titanium tubes tend to obtain diamond wrinkling mode.

Fig. 4 Experimental results of Φ76.2×t1.07 CP3 titanium tubes axial compression forming: (a) H=40 mm, (b) H=80 mm, (c) H=120 mm

Fig. 5 Experimental results of Φ76.2×t1.07 6061 aluminum tubes axial compression forming: (a) H=40 mm, (b) H=80 mm, (c) H=120 mm

Table 1 Comparison of wrinkling mode and critical force of Φ76.2×t1.07 tubes with different height

Material

Titanium tube

Height /mm

Wrinkling mode

Limit load /kN

40

four-edged diamond

87.2878

40

axisymmetric

94.1532

40

four-edged diamond

88.8705

80

four-edged diamond

92.723

80

four-edged diamond

88.8205

80

three-edged diamond

89.5262

120

three-edged diamond

91.9878

120

three-edged diamond

90.9326

120

four-edged diamond

92.1493

*n-edged diamond represents n circumferential wrinkling waves Table 2 Comparison of wrinkling mode and critical force of Φ76.2×t1.07 tubes with different height Material

Aluminum tube

Height /mm

Wrinkling mode

Limit load /kN

40

axisymmetric

42.1211

40

axisymmetric

42.3602

40

four-edged diamond

40.9455

80

four-edged diamond

41.8647

80

axisymmetric

42.1634

80

four-edged diamond

42.1609

120

four-edged diamond

41.6615

120

four-edged diamond

42.4038

120

four-edged diamond

42.382

2.2.2 Tube draw bending under more constraints Fig. 6 shows that, for multi-tool constrained bending with different tubular materials, geometric specifications and loading conditions, though the occurring positions of wrinkling are various, including front portion, straight portion, curved portion, whole portion, upper curved portion and lower curved portion, the wrinkling forms are all similar with regular waves at the intrados of bent tubes. Fig. 6 (a) is the titanium tubular part focused in this study.

Fig. 6 Locations where wrinkles occurred in draw bending process: (a) front portion, (b) straight portion, (c) curved portion, (d) the whole portion, (e) upper curved portion, (f) lower curved portion

2.3 Comparison of characteristics under less constraints and more constraints As mentioned above, for tube axial compression with less constraints, different buckling modes such as axisymmetric mode, diamond mode would occur in the same geometric parameters and material types. Regarding tube draw bending, though occurring in different portions, the wrinkling waveforms belong to the similar buckling mode which can be unified as sinusoidal wave along the tangent direction of the bending center. It can be found in literature[51] that filling the hollow tube with concrete, which could be seen as a soft constraint, can delay the local buckling and change the buckling mode from outwards to inwards. Similar to tube compression experiments, multiple possible instability modes were also mentioned in Kyriakides‟s research[45,46] for tube pure bending with less constraints, such as short wavelength rippling (wrinkling), sharp local kink and kink. So the potential buckling modes are supposed to be various for deformation process with less constraints, and unified for deformation process with more constraints. It was found that the boundary constraint can affectively influence the buckling mode for a plate under edge compression with lateral constraint[52]. If the binder force is rising, the half-wave number of the buckling mode is incrementally increased, also leading to the increasing of buckling load and decreasing of the wavelength. This understanding provides a theoretical guide for the GMI construction of the different sheet metal forming process with less and more constraints. As shown in Fig. 7, the wrinkling mode of thin-walled tube axial compression is determined by the lower order buckling mode (point A). Since the wrinkling is sensitive to the change of the geometric parameters, material properties and boundary conditions, the tiny perturbation of these

parameters will cause different buckling modes with close buckling load (point B and C). The different buckling modes correspond to different bifurcation points and different experienced post-buckling behaviors, and different wrinkling results. To ensure the loading path proceed along the bifurcation path, the low order superposed buckling modes are seeded in the primary path. As for tube draw bending process constrained by multiple tools, the wrinkling tendency is suppressed to some extent. So the bifurcation path corresponds to a high order buckling mode with high deformation energy and low wavelength, which is the „nth path‟ shown in Fig. 7.

Fig. 7 Schematic of loading paths for less and multiple constraints

2.4 Integrated hybrid methodology for sensitive wrinkling prediction under different degrees of constraints As show in Fig. 8, a framework of integrated hybrid methodology for wrinkling prediction is explored by introducing evolution of triple nonlinearity in geometry, material and

boundary

conditions,

viz.,

multiple

geometric

micro-imperfections

(GMI),

anisotropic/asymmetric properties (AAP) and explicit FE algorithm, respectively. To evaluate the capability of wrinkling prediction under different degrees of constraints, the above cases of large-diameter thin-walled titanium tubular parts are tested under less and multiple constraints, and comparisons between experimental and numerical results are carried out.

Large-diameter thin-walled aluminum and titanium tubular parts Plastic wrinkling characteristics under less and multiple constraints Influenced by

Stress state evolution Geometry evolution

Interaction

Material evolution

Contact evolution

GMI

Explicit FE algorithm

Asymmetry /anisotropy

Coupling model by considering triple nonlinearity Comparisons between experimental and numerical results Capability of wrinkling prediction for less and multiple constraints are elevated Fig. 8 Framework of integrated hybrid methodology for wrinkling prediction in sheet metal forming

3. Micro-imperfection related geometrical nonlinearity for wrinkling prediction 3.1 Appropriate methods for constructing geometric micro-imperfections In a mathematical sense, buckling is a bifurcation in the solution to equations of static equilibrium. At a certain point, under an increasing load, any further load is able to be sustained in one of two states of equilibrium. To introduce GMI, linear buckling analysis based on the capability of eigenvalue extraction is established with implicit algorithm. Critical load and buckling modes can be obtained from the following equation[53]:

( Ke  i K g )i  0

(5)

In Eq. (5), Ke is the elastic stiffness matrix regarding the base state, K g is the differential initial stress and load stiffness matrix due to the incremental loading pattern, is the eigenvalue and

i



is the ith buckling mode shape. To search for a singularity point of

the coefficient matrix, the value of

 is increasing gradually and when the equation has a

nontrivial solution, the buckling load and buckling mode is extracted. These buckling modes can be implemented into explicit FE models as superposed initial GMI with the form as: N

xi   Aii i 1

(6)

where Ai is the scaling factor for GMI and N is the order of the buckling mode. Thus, by setting the boundary conditions and the loading conditions, different buckling modes are obtained through the above analysis. As for tube compression with less constraints, the low buckling modes could be relatively accurate and useful for constructing GMI. Since three wrinkling modes occurred in the experiments, i.e. axisymmetric, three-edged diamond and four-edged diamond mode, the buckling modes of the first three orders are constructed and superposed for tube compression. The scale factor is set to be 0.01 of the thickness according to ABAQUS User‟s Manual[53]. As for draw bending with more constraints, the analysis doesn‟t take the normal/lateral constraint into consideration such as the wiper die and mandrel die in tube draw bending under multiple constraints. Therefore, a high order buckling mode is supposed to be constructed for complex bending with multiple constraints. In order to realize sensitive and accurate wrinkling prediction, the buckling mode should own analogous waveform with the experiment and the right wavelength and scale factor are filtered based on the minimum energy principle. The solution procedure is accordingly shown in Fig. 9. Start Establish buckling model Write the buckling modes in the global system as nodal data Less constraints Multiple low order superposed buckling modes

Multiple constraints

High order buckling mode with analogous waveform to experiments Minimum energy

Update the wave length of imperfection N

Imperfection λ1 Imperfection λ2 ……

Amplitude = 0.01 of the thickness

Select amplitude according to experiment or minimum energy

The appropriate construction for initial GMI End

Fig. 9 Solution procedure of imperfection selected method under different constraints

3.2 Geometric micro-imperfections construction 3.2.1 Tube axial compression under less constraints Regarding the tube axial compression, the boundary conditions of buckling analysis are fixed support at one side, and free axial degree at the other side, as shown in Fig. 12(a). The element type for the deformable tube is S4R and the mesh size is set as 1 mm × 1 mm. The material properties and geometric parameters are referred to Fig. 3. Hence the buckling modes from 1 to 10 are obtained and the results are given in Table 3. Then according to ABAQUS User‟s Manual, the scaling factor of imperfection is set as 0.01 of the wall thickness. The 3D FE model for axial compression is established with three parts, i.e. the thin-walled tube, upper die and lower die (see Fig. 10). The tube possesses the same element type and mesh setting with buckling analysis and the two dies are set as rigid parts with R3D4 element. The reference point of the upper die is set in the geometric center and is locked of all its degrees of freedom except its displacement along axial direction. Considering the efficiency and accuracy, the mass scaling factor is set to be 10000. Other parameters are the same as experiments. By taking the 40 mm height CP3 tube as a typical case, the different GMI corresponding to different buckling modes, i.e. the 1st, 2nd and 6th buckling mode, are seeded individually into the FE model. As shown in Fig. 11, different wrinkling modes are obtained, viz. the five-edged diamond mode, four-edged diamond mode and axisymmetric mode. Even though the difference is little from the beginning, it is effective to change the wrinkling mode. Hence, we can infer that the wrinkling mode is sensitively influenced by the GMI. If the tube naturally owns different GMI, the various wrinkling modes would occur with the same experimental conditions. To predict the potential wrinkling modes as completely as possible, the first three buckling modes are simultaneously seeded into the model through their linear combination. The simulation results for CP3 tube with height of 40 mm, 80 mm and 120mm, are given in Table 4. When tube height is 40 mm, it is worth noting that both axisymmetric and diamond mode occurred in the opposite half part of the tube. The results agree well with the experiments, which show both the axisymmetric and diamond wrinkling mode. As for 80 mm height tube, the GMI-implemented model only gives a four-edged wrinkling mode. However, the three-edged and four-edged diamond wrinkling mode occurred in the experiments. Similarly, regarding 120 mm height tube, the GMI-implemented model only shows a four-edged diamond wrinkling mode, but the experiments have shown two wrinkling modes, i.e. the

three-edged diamond mode and four-edged diamond mode. A possible explanation for this might be the magnifying difference of the introduced GMI corresponding to different buckling modes due to the increase of the tube height. This may result in the domination by one certain form of buckling mode and the weakening of the other form of buckling mode in GMI.

Fig. 10 FE model of CP3 tube axial compression forming Table 3 Buckling results of 40 mm height CP3 tube under axial compression Mode

1

2

3

4

5

Mode

6

7

8

9

10

Fig. 11 Simulation results of 40 mm height CP3 tube with sole GMI induced from buckling (a) mode 1, (b) mode 2, (c) mode 6

Table 4 Experimental and simulation results of tube axial compression H=40 mm H=80 mm

H=120 mm

GMI

Exp

3.2.2 Tube draw bending under more constraints In draw bending process, the tube is generally subjected to two combined forces, i.e. the axial compression from pressure die and the bending moment from bend die and clamp die. Thus in linear buckling analysis, the two loading conditions are set to achieve different buckling modes. For axial loading condition, the linear buckling analysis under axial load is the same as tube axial compression. For bending condition, the bending moment is coupled to a reference point of the clamped center. According to the above discussion, buckling mode under multiple constraints will cross over the lower order buckling mode to a high order one. Therefore, the buckling mode is chosen with a high order that owns an analogous waveform shape as the experiments, of which 14th mode is for axial load and the 4th mode is for bending load, as shown in Fig. 12 (b). Fig. 12 (a) shows that the boundary conditions of tube axial compression are symmetric as there is no more constrain surrounding the tube. Fig 12. (b) shows that non-axisymmetric boundary conditions are set to the tube bending process because the outside and inside of the tube are under different deformation process. As a result of the compressive stress distribution of the two loading conditions, the axial load gives an axisymmetric GMI and the bending moment gives a non-axisymmetric GMI since it is only subjected to a compressive load at intrados. The FE model includes a deformable S4R tube, and rigid R3D4 bend die, clamp die, wiper die, pressure die, and mandrel. The material properties are referred to Fig. 3. Thus, a 3D FE explicit model is established (see Fig. 13).

Fig. 12 The boundary conditions and the buckling results under more constraints: (a) axisymmetric GMI and (b) non-axisymmetric GMI

Fig. 13 Cutaway view of thin-walled CP3 tube bending FE model

The wavelength and scale factor for the GMI should be filtered based on Timoshenko‟s energy method after the determination of the buckling mode. The reasonable GMI should lead to the least energy during forming process. In the forming process, the sum of ALLSE and ALLPD is an indicator of strain energy, of which ALLSE is the recoverable elastic strain energy and ALLPD is the plastic dissipation energy in energy output variables of ABAQUS/explicit. A series of GMI with different wavelengths are constructed around the initial wavelength. The strain energy of axisymmetric and non-axisymmetric GMI with different wavelengths is shown in Fig. 14 (a) and Fig. 15 (a) respectively. As the time goes on, the strain energy increases quickly at the beginning and slow down appropriate at the point of 2s. It can be found that as the wavelength increases, the strain energy is going up. The energy consumption of non-axisymmetric GMI is generally higher than that of axisymmetric GMI. The strain energy at the ending time of forming process is extracted and plotted in wavelength-strain energy curves, as shown in Fig. 14 (b) and Fig. 15 (b). Then, through polynomial function fitting of relations between strain energy and wavelength, the appropriate wavelength for the geometric imperfection with the lowest strain energy is obtained. From the fitting curves of these two figures, the lowest energy of axisymmetric GMI condition is with

the wavelength of 20.33mm, while that of the non-axisymmetric GMI condition is with the wavelength of 12,37mm. As shown in Fig. 16, for axisymmetric GMI, the strain energy generally shows an increasing tendency as the scale factor increases and then the minimum energy cost is difficult to distinguish. Combined with the experimental results, when the scale factor is selected as 0.3, the prediction shows the best agreement with the experiment. As for non-axisymmetric GMI, the energy consuming falls from the scale factor 0.01 to 0.3 and then increases. So both the scale factors for two GMI are set to be 0.3. 3.00

(a)

(b) FE model result Fitting curve

2.9

wavelength wavelength wavelength wavelength wavelength wavelength

2.8

2.7

1.95

2.10

Strain Energy /kJ

Strain Energy /kJ

3.0

16.0121mm 19.0131mm 22.0105mm 24.0134mm 28.0148mm 31.0185mm

2.25

2.98

The lowest strain energy with the wavelengthcr=20.33mm 2.96

2.94

2.40

15

20

Time /s

25

30

Wavelength /mm

Fig. 14 Energy consumption with different wavelength for axisymmetric GMI: (a) strain energy-time curve, (b) strain energy-wave length fitting curve (b)

(a)

2.8

1.95

2.10

Time/s

2.25

2.40

Strain energy/kJ

Strain energy/kJ

wavelength 10.0071mm wavelength 11.0077mm wavelength 12.0063mm wavelength 13.0098mm wavelength 14.0153mm wavelength 16.0139mm wavelength 18.0065mm

3.0

1.80

FE model result Fitting curve

3.25

3.2

3.20

The lowest strain energy with wavelength cr=12.37 mm

3.15

3.10

10

12

14

16

Wavelength/mm

Fig. 15 Energy consumption with different wavelength for non-axisymmetric GMI: (a) strain energy-time curve, (b) strain energy-wave length fitting curve

18

3.04

(a)

(b)

Strain Energy /kJ

Strain Energy /kJ

3.00

2.96

2.92

3.30

3.15

2.88 0.0

0.1

0.2

0.3

Scale factor b

0.4

0.5

3.00

0.0

0.1

0.2

0.3

0.4

0.5

Scale factor b

Fig. 16 Strain energy-scale factor curve of different geometric imperfections: (a) axisymmetric GMI with wavelength 20.33mm, (b) non-axisymmetric GMI with wavelength 12.37mm

The contours of wall thickness for axisymmetric GMI model, non-axisymmetric GMI model and perfect geometry model are shown in Fig. 17. The perfect geometry model gives a smooth ridge line at intrados with homogeneous thickness distributions increasing from extrados to intrados. As for the two GMI-implemented models, the wall thickness changes unevenly at the intrados of the tube. As shown in Table 5, comparisons are made from the aspect of wrinkling position, wrinkling height and the wall thickness. If the wrinkling height is under 0.5 mm, the wrinkling angle is indicated by the position where the wall thickness increases most. The axisymmetric GMI model predicts three wrinkling waves at the intrados, but the wrinkling height is much lower than the experimental value and the predicted wrinkling position angle ranges from 8.7º to 36.6º, which is wider than the real wrinkling position that focuses between 22.0º and 33.6º. The non-axisymmetric GMI model does not give discernable wrinkles at intrados, but the wall thickness fluctuates from an initial angle of 7.2º to 27.7º for about eight times. Regarding the wall thickness, if the tube deforms homogeneously without wrinkling, the wall thickening at intrados is small and uniform. As seen from the results, the GMI-implemented models endure more severe wall thickening than perfect geometry model and the values are closer to experiment. Generally speaking, compared with the perfect geometry mode, the GMI-implemented models show more sensitive prediction for the wrinkles or wrinkling tendency, but the models are not reliable on the wrinkling height and wrinkling position.

Fig. 17 Thickness distributions for different models (Mises yield criterion): (a) perfect geometry model, (b) axisymmetric GMI, (c) non-axisymmetric GMI and (d) experiment result. Table 5 Comparisons of three waves with different geometric imperfection First wave Model

Second wave

Angle (°)

Height

Thickness

(mm)

(a)

8.6

(b)

8.7

Height

Thickness

(mm)

Angle (°)

(mm)

Unobvious

1.25

13.5

0.74

1.36

11.6

(c) Exp

Third wave Height

Thickness

(mm)

Angle (°)

(mm)

(mm)

Unobvious

1.29

20.1

Unobvious

1.22

0.63

1.49

36.6

0.21

1.26

2.48

1.40

Unobvious (Wall thickness fluctuates from 7.2º to 27.7º) 22.0

2.30

1.56

28.6

3.08

1.58

33.7

4. Anisotropy/asymmetry related material nonlinearity for wrinkling prediction 4.1 Construction for constitutive models Sensitive wrinkling prediction not only requires the design for appropriate GMI but also precise characterization of plastic deformation. Especially for hexagonal close-packed metals that are sensitive to the loading directions, common isotropic yield functions are not able to

depict its distorted behavior well such as Mises and Tresca yield criterion. Hence in this section, the isotropic model Mises, the normal anisotropic model Hill‟48-r and the discontinuous anisotropic/asymmetric stress invariable based model (SIM)[54] are applied and compared with experimental results to study the effects of the AAP related constitutive features. The Hill‟48 yield criterion is a pioneering anisotropic plastic yield criterion firstly raised by Hill as an extension of the von Mises isotropic yield criterion. In this section, only the normal anisotropy is taken into account for titanium tube. The SIM is based on a discontinuous framework and aimed at capturing the irregular evolution of yield loci under various plastic strain levels, especially the large deformation. And the newly proposed Yoon yield function[41] is able to accurately describe both the asymmetry and anisotropy in plastic deformation. As for pressure insensitive titanium alloy, the yield function can be expressed as follows: f ( )  ( J 2 ( Σ A ))3/2  J 3 ( Σ B )  ( Y )3  0

(7)

where: J2 (Σ A ) 

2 2 2 2 2 2 1 A 1 Σ : Σ A   Σ11A    Σ 22A    Σ 33A     Σ12A    Σ13A    Σ 23A    2 2

2 2 2 J3 ( Σ B )  det( Σ B )  Σ11B Σ22B Σ33B  2 Σ12B Σ23B Σ31B   Σ11B  Σ23B   Σ22B  Σ13B   Σ33B  Σ12B    

(8) (9)

In the formulation, J 2 (Σ A ) and J 3 (Σ B ) are the second and third stress invariant for stress tensor Σ A and Σ B . Σ A and Σ B are transformed from the stress tensor σ by two fourth-order linear transformation tensors:

Σ A  LA   ; Σ B  LB  

(10)

Linear transformation tensors can be expressed as follows: a3 / 3 a2 / 3 0 (a3  a2 ) / 3  a / 3 (a1  a3 ) / 3 a1 / 3 0 3   a2 / 3 a1 / 3 (a1  a2 ) / 3 0 LA =  0 0 0 a4   0 0 0 0  0 0 0 0 

0 0 0 0 a5 0

0 0  0  0 0  a6 

(11)

b3 / 3 b2 / 3 0 (b3  b2 ) / 3  b / 3 (b1  b3 ) / 3 b1 / 3 0 3   b2 / 3 b1 / 3 (b1  b2 ) / 3 0 LB =  0 0 0 b4   0 0 0 0  0 0 0 0 

0 0 0 0 b5 0

0 0  0  0 0  b6 

(12)

With respect to the in-plane stress state, eight material coefficients are required to be calibrated, i.e. a1, a2, a3, a4, b1, b2, b3 and b4. Combined with experimental results, visco-plastic self-consistent (VPSC) model and genetic algorithm, the material coefficients are obtained, as shown in Table 6. The schematic of yield loci with different plastic stain for Hill‟48-r model and SIM are plotted in Fig. 18. yy(MPa)

（ b）

yy(MPa)

(a)

800

30%

800 30%

20%

20%

400

400

-800

-400

10%

400

800

10%

xx(MPa)

-800

-400

800

xx(MPa)

5%

5%

-400

400

-400

2.5%

2.5%

0.2%

0% -800

-800

Solid dot: Exp+VPSC Temperature: 25℃

Solid dot: EXP+VPSC Temperature: 25℃

Fig. 18 Calibrated yield loci of CP3 titanium tube constructed with experimental and VPSC results under plane stress for (a) Hill‟48-r model and (b) SIM Table 6 Material coefficients in yield function for CP3 titanium tube under various levels of plastic strain



p

a1

0 1.287

0.03 1.335

0.05 1.382

0.1 1.435

0.2 1.457

0.3 1.462

a2

1.609

1.623

1.592

1.535

1.474

1.476

a3

1.944

1.930

1.903

1.816

1.714

1.678

b1

-0.460

-0.427

-0.291

-0.988

-0.974

-1.031

b2

1.701

1.635

2.094

-1.069

-1.295

-1.342

b3

0.585

0.618

0.180

-1.081

-1.579

-1.645

a4

1.784

1.797

1.795

1.775

1.708

1.657

b4

0.434

0.460

-0.360

-0.782

1.250

-1.350

4.2 Tube axial compression under less constraints The wrinkling prediction results of different constitutive models are listed in Table 7. Regarding 40 mm height tube, both Hill‟48-r and SIM model show axisymmetric wrinkling mode. With respect to 80 mm height tube, the Hill‟48-r model shows a three-edged wrinkling mode and SIM model shows a four-edged wrinkling mode, which is same with the tubes of 120 mm height. Compared with the GMI-implemented model, it indicates that the AAP is not so sensitive as GMI to change the wrinkling mode under less constraints.

Model

Table 7 Simulation results of tube axial compression (without GMI) H=40 mm H=80 mm H=120 mm

Hill‟48-r

SIM

4.3 Tube draw bending under more constraints Fig. 19 shows the wall thickness distributions of Hill‟48-r model and SIM. Compared with the Mises model shown in Fig. 17，Hill‟48-r model shows the wrinkling tendency with fluctuating wall thickness distribution at the intrados and SIM model shows three discernable wrinkling waves in the initial bending angle. Listed in Table 8 are the detailed information of wrinkling position, wrinkling wave height, and wall thickness for Hill‟48-r model, SIM model and the experimental results. Since the Hill‟48-r model shows no discernable wrinkling waves, the related information is not given in Table 8. As for the wrinkling position, Hill‟48-r model gives wrinkling tendency in whole bending area, and the wrinkles predicted by SIM are occurring below 15º. However, the experimental wrinkles are concentrating in the range from 22.01º to 33.66 º. Regarding the wrinkling height, only the SIM gives countable wrinkling height, but the experimental value is larger than that of SIM. Further, the wall thickness predicted by SIM is larger than that by Mises model and closer to the experiments. Overall, although the wrinkling angle range is lower, the SIM gives out

discernable wrinkling wave, the right wrinkling number, and the closer wall thickness. The results verify the effectiveness of the SIM with consideration of both the asymmetry and anisotropy. Despite the less efficiency of Hill‟48-r model compared to SIM about discernable wrinkling wave, it still shows more severe wrinkling tendency compared to Mises model. Compared to axial compression with less constraints, the added anisotropy and asymmetry are more sensitive for multiple constraints than less constraints. The reason may be the intricate stress/strain conditions, involving non-uniform tension and compression in draw bending under multiple constraints.

(a) (b) Fig. 19 Thickness distributions for different constitutive models (without GMI): (a)Hill‟48-r, (b)SIM Table 8 Comparison of three wrinkle waves between different constitutive models and experiment First wave Model SIM

Third wave

Angle

Height

Thickness

Angle

Height

Thickness

Angle

Height

Thickness

(°)

(mm)

(mm)

(°)

(mm)

(mm)

(°)

(mm)

(mm)

5.5

1.74

1.42

8.7

1.97

1.41

11.8

2.43

1.31

2.48

1.40

Hill‟48-r Exp.

Second wave

Unobvious (Wall thickness fluctuates from 0º to 90º) 22.0

2.30

1.56

28.6

3.08

1.58

33.7

5. Coupling of GMI and anisotropy/asymmetry under different boundary conditions 5.1 Tube axial compression under less constraints As for tube axial compression, the lower order superposed GMI are implemented into the coupling models as well as the Hill‟48-r normal anisotropy or anisotropic/asymmetric SIM.

By considering the accuracy and efficiency, only the 40 mm and 80 mm tubes with two typical wrinkling modes, i.e. axisymmetric and diamond mode, are taken into consideration. Fig. 20 illustrates simulation results of the 40 mm height tubes for the coupling model (a), representing the Hill‟48-r and GMI, and the coupling model (b), representing the SIM and GMI. It can be seen that the coupling model (a) shows both diamond mode and axisymmetric mode, but the left diamond mode part is weakened when compared to the GMI-implemented model in Table 4. However, coupling model (b) is totally showing axisymmetric mode. With respect to 80 mm height CP3 tube, the results are shown in Fig. 21. Two types of coupling models all show a four-edged diamond wrinkling mode. Broadly speaking, since tube compression is of simple loading conditions and the stress state is simple before buckling or wrinkling, the anisotropy/asymmetry seems not to be sensitive to predict the wrinkling mode and may also weaken the effects of GMI. To the author‟s view, introduced GMI is enough for wrinkling prediction under less constraints and simple stress/strain state.

(a) Hill‟48-r +GMI (b) SIM+ GMI Fig. 20 Simulation results for thickness distribution of 40 mm height CP3 tube: (a) coupling model of Hill‟48-r and GMI, (b) coupling model of SIM and GMI

(a) Hill‟48-r + GMI (b) SIM+ GMI Fig. 21 Simulation results for thickness distribution of 80 mm height CP3 tube: (a) coupling model of Hill‟48-r and GMI, (b) coupling model of SIM and GMI

5.2 Tube draw bending under more constraints As for tube draw bending, axisymmetric GMI and non-axisymmetric GMI are individually implemented into the coupling models as well as Hill‟48-r normal anisotropy or anisotropic/asymmetric SIM. Fig. 22 illustrates simulation results of the bent tubes for the coupling model (a), representing the Hill‟48-r and axisymmetric GMI, and the coupling model (b), representing the Hill‟48-r and non-axisymmetric GMI, coupling model (c), representing the SIM and axisymmetric GMI, and coupling model (d), representing the SIM and non-axisymmetric GMI. The coupling model (a) and (b) do not give discernable wrinkling waves, but reveal the wrinkling tendency at the intrados of tubes (uneven wall thickness distribution). The coupling model (c) and (d) both give the discernable wrinkling waves in the front bending region with heavier wall thickness changing. To make comprehensive analysis for the coupling effects, the local amplification figures of all the numerical results mentioned and the cutting wrinkled titanium tube are illustrated in Fig. 23 and the detailed information of the wrinkling characteristics are given in Table 9 with three indices, i.e. wrinkling range, maximum wrinkling height and wrinkling number. As for the wrinkling range, the experimental wrinkling waves are occurring from 22.0º to 33.7º and the span of the wrinkling angle is about 10 degrees. Compared to this, all models show a relatively small span of wrinkling range except Hill‟48-r model and coupling model of Hill‟48-r and axisymmetric GMI, but the wrinkles all initiate before 10º of the bending angle. Regarding the maximum wrinkling height, only four FE models (marked using greyhound) show the availabilities by predicting discernable wrinkling waves, i.e. axisymmetric GMI model, SIM, the coupling model of SIM and axisymmetric GMI and the coupling model of SIM and non-axisymmetric GMI. Among the four models, coupling model of non-axisymmetric GMI and SIM gives the most reasonable predictions. Regarding the wrinkling number, the four models with discernable wrinkling waves are all right about the wrinkling numbers except coupling model of non-axisymmetric GMI and SIM. Regarding different initial GMI, the non-axisymmetric GMI is better at improving the wrinkling angle compared to axisymmetric GMI. And regarding the constitutive models, the SIM considering the anisotropy and asymmetry can predict the discernable wrinkling height better. Taken together, coupling model of non-axisymmetric GMI and SIM is of the best fitness to the experiment. More intuitionistic results can be seen in Fig. 23 about the wall thickness fluctuating and the ultimate geometric shape at tube intrados and the part enclosed by a red circle is the wrinkling position.

Conclusions can be made through the above comparisons: 1) as for draw bending process with more constraints, the proper constitutive relations with anisotropy and asymmetry are more sensitive to wrinkling prediction than the GMI, but the axisymmetric and non-axisymmetric GMI could help improve its precision. The reason may be that the enhanced multiple constraints limit the evolution of initial GMI and the intricate stress/ strain state of inhomogeneous deformation; 2) the non-axisymmetric GMI induced by bending moment shows better fitness with experiments than the axisymmetric GMI induced by axial compression. That is because bending moment is the main load that tube endures during bending, but not the axial compression.

(a) Hill‟48-r + axisymmetric GMI

(b) Hill‟48-r + non-axisymmetric GMI

(c) SIM + axisymmetric GMI (d) SIM + non-axisymmetric GMI Fig. 22 Thickness distribution of different coupling models

Mises

Mises axisymmetric

Mises non-axisymmetric

Hill‟48-r

Hill‟48-r + axisymmetric GMI

GMI

GMI

Hill‟48-r SIM SIM SIM Experiment non-axisymmetric axisymmetric non-axisymmetric GMI GMI GMI Fig. 23 Local amplification of thickness distribution of different prediction methods Table 9 Three indices for wrinkling waves with different FE models and experiment Model Wrinkling range (º) Max wrinkling height (mm) Wrinkling number Mises ( 8.6, 20.1) 0 Mises + axis-GMI ( 8.8, 36.6) 0.7 (too low) Mises + non-axis-GMI ( 7.2, 27.7) 0 Hill‟48-r ( 0, 90) 0 Hill‟48-r + axis-GMI ( 0, 90) 0 Hill‟48-r + non-axis-GMI ( 6.8, 21.5) 0 SIM ( 3.5, 11.0) 2.4 SIM + axis-GMI ( 4.2, 14.9) 1.9 SIM + non-axis-GMI ( 6.2, 20.1) 2.9 Experiment ( 22.0, 33.7) 3.1 *non-axis-GMI represents non-axisymmetric GMI, axis-GMI represents axisymmetric GMI

0 3 0 0 0 0 3 3 4 3

6. Conclusion Taking tube axial compression and tube draw bending as two cases, this study focuses on sensitive wrinkling prediction in sheet metal forming by introducing inevitable triple nonlinearity in geometry, material and contact into 3D-FE model. Through coupling the GMI related geometrical nonlinearity, anisotropy/asymmetry related material nonlinearity and contact nonlinearity, an integrated hybrid framework is established based on explicit FE platform for realizing sensitive wrinkling prediction of sheet metal forming under less and more boundary constraints, then the capability for wrinkling prediction of the above framework is thoroughly evaluated. The main contents and results are summarized as followed: (1) Through experimental exploration in tube axial compression and draw bending, the wrinkling characteristics under less constraints and multiple constraints are articulated. As for axial compression, the wrinkling mode would transfer from axisymmetric mode to

diamond one with increasing tube height for each tubular material. Especially various low order wrinkling modes occur even with the same parameters, which is assumed to be caused by initial GMI. While, regarding draw bending, the closely imposed constraints are found to change the instability mode of tube from a lower order buckling mode to a higher order one, which consumes more energy, and under multiple constraints, the wrinkling waveforms are revealed to be the same regardless of various occurring segments. (2) Based on the linear buckling analysis, multiple GMI modes are constructed and superimposed into explicit FE models corresponding to the wrinkling characteristics under different degrees of constraints, viz., the first three superposed buckling modes are simultaneously coupled into FE model for less constraints, one high order buckling mode with analogous wrinkling waveform is selected for forming operation with multiple constraints. Regarding the axial compression, the sensitivity of sole GMI and three superposed GMI for wrinkling is studied respectively. Regarding draw bending, axisymmetric and non-axisymmetric GMI are implemented and are proved to be sensitive to the wrinkling waves or the wrinkling tendency. (3) Three constitutive models, viz, Mises, normal anisotropic Hill‟48-r and stress invariant based anisotropic/asymmetric models (SIM), are comparatively implemented in explicit FE models without initial GMI. As for tube axial compression, AAP models are less sensitive to the wrinkling prediction compared to GMI, which may be attributed to the simple stress/strain state under less constraints. With respect to multiple constraints, the Hiil48-r model reflects the wrinkling tendency in the whole bending region and SIM gives out discernable wrinkling waves, which shows the effectiveness for wrinkling prediction. (4) By coupling the GMI and SIM into explicit FE models, the desired wrinkling prediction capability for both forming cases is further confirmed since evolution of triple nonlinearity is simultaneously considered. For axial compression under less constraints, coupling effects are not obvious since its simple constraint and loading condition, whereas for tube draw bending under multiple constraints, the coupling model of non-axisymmetric GMI and SIM is of the best agreement with experiments. The non-axisymmetric GMI improves the prediction accuracy of wrinkling range better and AAP improves the wrinkling height more effective.

Acknowledgments The authors gratefully acknowledge the financial support from the National Science Fund for Excellent Young Scholars (51522509), the Commercial Aircraft Research and

Development Project of China (MJ-2016-G-64), and the National Natural Science Foundation of China (51275415) and (51605382).

Reference [1] Ahmetoglu M, Altan T. Tube hydroforming: state-of-the-art and future trends. J Mater Process Technol 2000;98:25-33. [2] Hashmi MSJ. Aspects of tube and pipe manufacturing processes: meter to nanometer diameter. J Mater Process Technol 2006;179:5-10. [3] Karman T Von, Tsien H-S. The buckling of thin cylindrical shells under axial compression. J Aeronaut Sci 1941;8:303-12. [4] Donnell LH, Akron O. A new theory for the buckling of thin cylinders under axial compression and bending. Trans ASME 1934;56:795-806. [5] Karman T Von, Tsien H-S. The buckling of spherical shells by external pressure. J Aeronaut Sci 1939;7:43-50. [6] Cao J, Boyce MC. A predictive tool for delaying wrinkling and tearing failures in sheet metal forming. J Eng Mater Technol 1997;119:354. [7] Wang X, Cao J. On the prediction of side-wall wrinkling in sheet metal forming processes. Int J Mech Sci 2000;42:2369-94. [8] Wang X, Cao J. An analytical prediction of flange wrinkling in sheet metal forming. J Manuf Process 2000;2:100-7. [9] Wang X, Cao J. Wrinkling limit in tube bending. J Eng Mater Technol 2001;123:430. [10] Chen Y Z, Liu W, Zhang Z C, et al. Analysis of wrinkling during sheet hydroforming of curved surface shell considering reverse bulging effect. Int J Mech Sci 2017;120:70-80. [11] Kong Q, Yu Z, Zhao Y, et al. Theoretical prediction of flange wrinkling in first-pass conventional spinning of hemispherical part. J Mater Process Technol 2017;246:56-68. [12] Watson M, Long H, Lu B. Investigation of wrinkling failure mechanics in metal spinning by Box-Behnken design of experiments using finite element method. Int J Adv Manuf Tech 2015;78(5-8):981-995. [13] Liu N, Yang H, Li H, et al. An imperfection-based perturbation method for plastic wrinkling prediction in tube bending under multi-die constraints. Int J Mech Sci 2015;98:178-194. [14] Riks E. An incremental approach to the solution of snapping and buckling problems. Int J Solids Struct 1979;15:529-51.

[15] Kim JB, Yang DY, Yoon JW, Barlat F. The effect of plastic anisotropy on compressive instability in sheet metal forming. Int J Plast 2000;16:649-76. [16] Kim JB, Yoon JW, Yang DY. Investigation into the wrinkling behaviour of thin sheets in the cylindrical cup deep drawing process using bifurcation theory. Int J Numer Methods Eng 2003;56:1673-705. [17] Neto D M, Oliveira M C, Santos A D, et al. Influence of boundary conditions on the prediction of springback and wrinkling in sheet metal forming. Int J Mech Sci 2017;122:244-254. [18] Najafabadi H M, Naeini H M, Safdarian R, et al. Effect of forming parameters on edge wrinkling in cold roll forming of wide profiles. Int J Adv Manuf Tech 2019;101(1-4):181-194. [19] Liu N, Yang H, Li H, Li ZJ. A hybrid method for accurate prediction of multiple instability modes in in-plane roll-bending of strip. J Mater Process Technol 2014;214:1173-89. [20] Liu N, Yang H, Li H, Yan S. Plastic wrinkling prediction in thin-walled part forming process: A review. Chinese J Aeronaut 2016;29:1-14. [21] Savadkoohian H, Arezoodar A F, Arezoo B. Analytical and experimental study of wrinkling

in

electromagnetic

tube

compression.

Int

J

Adv

Manuf

Tech

2017;93(1-4):901-914. [22] Neto D M, Oliveira M C, Dick R E, et al. Numerical and experimental analysis of wrinkling during the cup drawing of an AA5042 aluminium alloy. Int J Mater Form 2017;10(1):125-138. [23] Andrews KRF, England GL, Ghani E. Classification of the axial collapse of cylindrical tubes under quasi-static loading. Int J Mech Sci 1983;25:687-96. [24] Abebe M, Lee K, Kang B S. Surrogate-based multi-point forming process optimization for dimpling and wrinkling reduction. Int J Adv Manuf Tech 2016;85(1-4):391-403. [25] Cao J, Cheng SH, Wang HP, Wang CT. Buckling of sheet metals in contact with tool surfaces. CIRP Ann - Manuf Technol 2007;56:253-6. [26] Simonetto E, Ghiotti A, Bruschi S, Gemignani R. Dynamic detection of instability defects in tube rotary draw bending. Procedia Manuf 2017;10:319-28. [27] Zhou CH, Wang B, Luo HZ, Chen YW, Zeng QH, Zhu SY. Quasi-static axial compression of origami crash boxes. Int J Appl Mech 2017;09:1750066. [28] Donnell LH, Wan CC. Effect of imperfections on buckling of thin cylinders and columns under axial compression. J Appl Mech 1950; 17: 73-83.

[29] Yuan L, Kyriakides S. Plastic bifurcation buckling of lined pipe under bending. Eur J Mech A/Solids 2014;47:288-97. [30] Huang H, Han Q, Wei D. Buckling of FGM cylindrical shells subjected to pure bending load. Compos Struct 2011;93:2945-52. [31] Hilberink A, Gresnigt a. M, Sluys LJ. Liner wrinkling of lined pipe under compression: a numerical and experimental investigation. vol. 5. ASME; 2010. [32] Wong W, Pellegrino S. Wrinkled membranes III: numerical simulations. J Mech Mater Struct 2006;1:63-95. [33] Papadopoulos V, Stefanou G, Papadrakakis M. Buckling analysis of imperfect shells with stochastic non-Gaussian material and thickness properties. Int J Solids Struct 2009;46:2800-8. [34] Vryzidis I, Stefanou G, Papadopoulos V. Stochastic stability analysis of steel tubes with random initial imperfections. Finite Elem Anal Des 2013;77:31-9. [35] Chen G, Zhang H, Rasmussen KJR, Fan F. Modeling geometric imperfections for reticulated shell structures using random field theory. Eng Struct 2016;126:481-9. [36] Wheeler A, Pircher M. Measured imperfections in six thin-walled steel tubes. J Constr Steel Res 2003;59:1385-95. [37] Bisagni C. Numerical analysis and experimental correlation of composite shell buckling and post-buckling. Compos Part B Eng 2000;31:655-67. [38] Orifici AC, Bisagni C. Perturbation-based imperfection analysis for composite cylindrical shells buckling in compression. Compos Struct 2013;106:520-8. [39] Paquette JA, Kyriakides S. Plastic buckling of tubes under axial compression and internal pressure. Int J Mech Sci 2006;48:855-67. [40] Zhang LC, Yu TX, Wang R. Investigation of sheet metal forming by bending-Part II. Plastic wrinkling of circular sheets pressed by cylindrical punches. Int J Mech Sci 1989;31:301-8. [41] Yoon JW, Lou Y, Yoon J, Glazoff M V. Asymmetric yield function based on the stress invariants for pressure sensitive metals. Int J Plast 2014;56:184-202. [42] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. Int J Plast 2006;22:1171-94. [43] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc R Soc Lond A 1948; 193(1033): 281-297. [44] Li H, Zhang HQ, Yang H, Fu MW, Yang H. Anisotropic and asymmetrical yielding and its evolution in plastic deformation: Titanium tubular materials. Int J Plast 2017;90:177-211.

[45] Ju GT, Kyriakides S. Bifurcation and localization instabilities in cylindrical shells under bending-II. Predictions. Int J Solids Struct 1992;29:1143-71. [46] Kyriakides S, G.T. J. Bifurcation and localization instabilities in cylindrical shells under bending-I: Experiments. Int J Solids Struct 1992;29:1117-42. [47] Al Galib D, Limam A. Experimental and numerical investigation of static and dynamic axial crushing of circular aluminum tubes. Thin-Walled Struct 2004;42:1103-37. [48] Bardi FC, Kyriakides S. Plastic buckling of circular tubes under axial compression-part I: Experiments. Int J Mech Sci 2006;48:830-41. [49] Yang H, Li H, Zhang Z, Zhan M, Liu J, Li G. Advances and trends on tube bending forming technologies. Chinese J Aeronaut 2012;25:1-12. [50] Li H, Ma J, Liu BY, Gu RJ, Li GJ. An insight into neutral layer shifting in tube bending. Int J Mach Tools Manuf 2018;126:51-70. [51] Fam A, Qie FS, Rizkalla S. Concrete-filled steel tubes subjected to axial compression and lateral cyclic loads. J Struct Eng 2004;130:631-40. [52] Cao J, Boyce MC. Wrinkling behavior of rectangular plates under lateral constraint. Int J Solids Struct 1997;34:153-76. [53] ABAQUS User‟s Manual. Version 6.4. vol. 1080. 2003. [54] Li H, Hu X, Yang H, Li L. Anisotropic and asymmetrical yielding and its distorted evolution: modeling and applications. Int J Plast 2016;82:127-58.