CONTROL OF DAMAGE IN FLANGES HYDROFORMING

C. Labergère*, A. Lejeune† and J.C. Gélin†

*Applied Mechanic Laboratory

ENSMM/CNRS/University Franche-Comté

24 Chemin de l'épitaphe 25000 Besançon (France)

*e-mail: carl.labergere@univ-fcomte.fr

e-mail : arnaud.lejeuene@univ-fcomte.fr

e-mail : jean-claude.gelin@univ-fcomte.fr

**Key words:** Hydroforming, control, meshless method, damage.

**Abstract.**

In the flanges hydroforming processes, the research of the way to control loading in pressure and blankholder load are often long and tiresome. One proposes a numerical control methodology based on an approximation per pieces of the objective and constraints functions compared to the control variables. This method allows the construction step by step ways of controlled loading with respect to the required objectives

A meshless method to build the objective function approximation in order to guarantee a good representation of the real objective function.

An application is proposed that consists to minimize the development of the ductile damage during a stamping operation with pressure to reach a rather high deep drawing. One carried out then the coupling between the ductile damage model and the finite element code and for reasons of simplicity, one works with isotropic damage assumptions.
# 1 INTRODUCTION

The flanges or tubes hydroforming processes are now increasingly used in automotive industry for manufacturing lightweight components. These processes permit to obtain complex structural parts comparatively to classical stamping ones and are used as example in automotive production for standard vehicles or prototype ones, and are often associated to the use of new incoming materials.

Among different hydroforming processes, one classically distinguishes two distinct technologies. In tube hydroforming technology, a tube generally with a circular cross section, is first shaped by bending and then clamped in a tooling system, filled by a fluid and then pressured. It results that the tube takes the inner shape of tools that allows producing complex seamless parts. The flange hydroforming technology consists to take benefit from hydraulic pressure to assist flange deep drawing by replacing the die cavity by a cavity filled by an hydraulic fluid. In such processes, the pressure in the die cavity decreases the sliding between the flange and the punch and reduce the friction between the flange and the die cavity under the blankholder, leading to the increase of Limit Drawing Ratio comparatively to classical deep drawing processes. The flange hydroforming process permit to get complex components in a single shape that cannot be produced by classical deep drawing processes. It also permits to get small size details due to the fact that pressure applies the flange on the dies.

In order to proceed with flange or tube hydroforming in the process window to avoid failure, bursting and buckling [5], one have proposed strategies to determine the limiting pressure in the die cavity in order to prevent risks of failure. Control strategy based on coupling FE simulation of flange hydroforming and optimization has been proposed that lead to a method that consists to build incrementally the pressure and blankholder force.

# 1 DUCTILE DAMAGE

In that part, one are interested to perform the implementation of a ductile damage made in the finite element code. For reasons of simplicity, one works with isotropic damage and isotropic hardening, the ductile damage model that has been chosen in the Lemaitre one [6]. Such the model induces few changes in the definition of mechanical equilibrium according to the principle of equivalence strain. In the case of the isotropic damage, the effective stress vector is connected to the stress vector by the relation:

$$ \tilde{\sigma} = \frac{\sigma}{1-d} \qquad (1) $$

where d is the damage variable.

Thus in the case of an orthotropic material obeying the generalized von Mises plasticity criterium, the behaviour of damaged material is expressed under the following form :
$$f(\sigma, \varepsilon) = \frac{\sqrt{\frac{3}{2} \boldsymbol{\sigma} : \mathbf{P} : \boldsymbol{\sigma}}}{(1-d)} - \sigma_0(\varepsilon^{-p}) = \frac{\Phi}{(1-d)} - \sigma_0(\varepsilon^{-p}) \quad (2)$$

where $\mathbf{P}$ is the orthotropic Hill tensor and $\sigma$ is the stress vector, $\Phi$ stands for the equivalent effective stress.

The generalized normality rule gives the flow law and the evolution of the internal variables as:

$$\begin{aligned} \dot{\varepsilon}^p &= \dot{\lambda} \frac{\partial f}{\partial \boldsymbol{\sigma}} = \frac{3}{2} \frac{\dot{\lambda}}{(1-d)} \frac{\mathbf{P} : \boldsymbol{\sigma}}{\Phi} = \frac{\dot{\Omega}}{(1-d)} \mathbf{P} : \boldsymbol{\sigma} \\ \dot{\Omega} &= \frac{3}{2} \frac{\dot{\lambda}}{\Phi} \end{aligned} \quad (3)$$

where $\dot{\lambda}$ is the plasticity multiplier and $\dot{\Omega}$ is the modified plasticity multiplier.

According to the equivalence strain principle, the equivalent plastic strain is calculated using the following realtionship:

$$\dot{\varepsilon}^p = \frac{\dot{\lambda}}{(1-d)} = \frac{2}{3} \frac{\dot{\Omega}}{(1-d)} \Phi \quad (4)$$

It results that stress vector and elastic strain tensor are related through the following expression that account damage effects:

$$\boldsymbol{\sigma} = (1-d)\mathbf{C}^e : \boldsymbol{\varepsilon}^e \quad (5)$$

If one derives relation (5) vs time, one thus obtain:

$$\dot{\boldsymbol{\sigma}} = -\dot{d}\mathbf{C}^e : \boldsymbol{\varepsilon}^e + (1-d)\mathbf{C}^e : \dot{\boldsymbol{\varepsilon}}^e = -\dot{d}\frac{\boldsymbol{\sigma}}{(1-d)} + (1-d)\mathbf{C}^e : \dot{\boldsymbol{\varepsilon}}^e \quad (6)$$

The integration constitutive equations consists in determining the stress state at time $t_{n+1}$ knowing stress state at time $t_n$. Several algorithms are possible, one have chosen an implicit integration algorithm. The preceding relation can be written in an incremental way at time $t_{n+1}$ as:

$$\Delta \boldsymbol{\sigma}_{n+1} = -\Delta d_{n+1} \frac{\boldsymbol{\sigma}_{n+1}}{(1-d_{n+1})} + (1-d_{n+1}) \mathbf{C}^e : \Delta \boldsymbol{\varepsilon}_{n+1} \quad (7)$$

For that, the strain increment at time $t_{n+1}$ is splitted into an elastic part and a plastic part so that one can write:

$$\Delta \boldsymbol{\varepsilon}_{n+1}^e = \Delta \boldsymbol{\varepsilon}_{n+1} - \Delta \boldsymbol{\varepsilon}_{n+1}^p \quad (8)$$
By combining relations (7) and (8), one obtain:

$$ \Delta\boldsymbol{\sigma}_{n+1} = -\Delta d_{n+1} \frac{\boldsymbol{\sigma}_{n+1}}{(1-d_{n+1})} + (1-d_{n+1})\mathbf{C}^e : (\Delta\boldsymbol{\varepsilon}_{n+1} - \Delta\boldsymbol{\varepsilon}_{n+1}^p) \quad (9) $$

If one replaces the plastic strain increment by his expression (3) according to the modified plasticity multiplier, one has:

$$ \boldsymbol{\sigma}_{n+1} - \boldsymbol{\sigma}_n = (d_n - d_{n+1}) \frac{\boldsymbol{\sigma}_{n+1}}{(1-d_{n+1})} + (1-d_{n+1}) \mathbf{C}^e : \left( \Delta\boldsymbol{\varepsilon}_{n+1} - \frac{\Delta\Omega}{(1-d_{n+1})} \mathbf{P} : \boldsymbol{\sigma}_{n+1} \right) \quad (10) $$

After transformation, one get the following expression for the stress tensor:

$$ \boldsymbol{\sigma}_{n+1} = \underbrace{\left[ \frac{(1-d_n)}{(1-d_{n+1})} \mathbf{I} + \Delta\Omega\mathbf{P} \right]}_{\Pi} = \underbrace{\boldsymbol{\sigma}_n + (1-d_{n+1})\mathbf{C}^e : \Delta\boldsymbol{\varepsilon}_{n+1}}_{\sigma_{n+1}^{\text{trial}}} \quad (11) $$

$$ \boldsymbol{\sigma}_{n+1} = \Pi^{-1} \boldsymbol{\sigma}_{n+1}^{\text{trial}} \quad (12) $$

where $\boldsymbol{\sigma}_{n+1}^{\text{trial}}$ is the elastic prediction of the stress vector, $\Pi$ is called the modified elastic operator.

## 1.1 Damage evolution

The damage evolution is chosen as a power function of the elastic energy restitution -Y, and is expressed in the following form [6]:

$$ \Delta d_{n+1} = d_{n+1} - d_n = \left( \frac{-Y}{S_0} \right)^{s_0} \Delta \bar{\varepsilon}_{n+1}^p \quad (13) $$

where $S_0$ and $s_0$ are the material parameters to be identified from experiments. In our case, one have chosen the following values for damage model parameters corresponding to a mild steel:

<table><caption>Table 1: Damage parameters used for a mild steel</caption><thead><tr><th>Material</th><th>S<sub>0</sub></th><th>s<sub>0</sub></th></tr></thead><tbody><tr><td>Steel</td><td>3.</td><td>1.</td></tr></tbody></table>

The elastic energy restitution expressed from effective stress $\Phi$ and triaxiality stress $\sigma_H$ is given by:
$$ -Y = \frac{\Phi^2}{2E(1-d_{n+1})^2} \left[ \frac{2}{3}(1+\nu) + 3(1-2\nu) \left( \frac{\sigma_H}{\Phi} \right)^2 \right] \quad (14) $$
with $$ \sigma_H = \frac{\text{tr}(\sigma_{n+1})}{3} $$

## 1.2 Determination of Stress tensor and damage variable

In a first time, one calculates the stress vector corresponding to an elastic prediction by using the following elastic stress prediction $\sigma_{n+1}^{\text{trial}} = \sigma_n + (1 - d_{n+1})C^e : \Delta\epsilon_{n+1}$. To get the quadratic Hill criteria, one gives a quadratic form of the plasticity criterion (2):

$$ F(\sigma_{n+1}, \bar{\epsilon}_{n+1}^p) = \frac{\Phi^2}{(1 - d_{n+1})^2} - \sigma_0^2 (\bar{\epsilon}_{n+1}^p) \quad (15) $$

The law hardening can be chosen as a Ludwick one, such as:

$$ \sigma_0(\bar{\epsilon}_{n+1}^p) = \sigma_y + K(\bar{\epsilon}_{n+1}^p)^n \quad (16) $$

If $\sigma_{n+1}^{\text{trial}}$ exceeds the plasticity limit, i.e. $F(\sigma_{n+1}^{\text{trial}}, \bar{\epsilon}_{n+1}^p) > 0$, one realize a correction by modifying equation (3). To get the stress tensor $\sigma_{n+1}$, it is necessary to calculate the increment of the modified plasticity multiplier corresponding to the damage variable $d_{n+1}$, for that one uses the plastic flow condition $F(\sigma_{n+1}, \bar{\epsilon}_{n+1}^p) = 0$ and relation (13). While replacing $\sigma_{n+1}$ by its expression resulting from equation (3), one obtain:

$$ \left\{
\begin{aligned}
G_1(\Delta\Omega, d_{n+1}) &= \frac{\Phi^2}{(1-d_{n+1})^2} - \left[ \sigma_y + K \left( \frac{2}{3} \frac{\Delta\Omega}{1-d_{n+1}} \Phi \right)^n \right]^2 = 0 \\
G_2(\Delta\Omega, d_{n+1}) &= d_{n+1} - d_n - \left( \frac{-Y}{S_0} \right)^{s_0} \left( \frac{2}{3} \frac{\Delta\Omega}{1-d_{n+1}} \Phi \right) = 0
\end{aligned}
\right.
\quad (17) $$
with $$ \Phi = \sqrt{\frac{3}{2} \sigma_{n+1}^{\text{trial}} \Pi^{-1} : \mathbf{P} : \Pi^{-1} \sigma_{n+1}^{\text{trial}}} $$ $$ -Y = \frac{\Phi^2}{2E(1-d_{n+1})^2} \left[ \frac{2}{3}(1+\nu) + \frac{(1-2\nu)}{3} \left( \frac{\text{tr}(\Pi^{-1}\sigma_{n+1}^{\text{trial}})}{\Phi} \right)^2 \right] $$
The resulting system with two non-linear equations in $\Delta\Omega$ and $d_{n+1}$ is solved by a Newton iterative method expressed as:

$$ \left[ \frac{\partial G_1(\Delta\Omega^i, d_{n+1}^i)}{\partial\Omega} \quad \frac{\partial G_1(\Delta\Omega^i, d_{n+1}^i)}{\partial d_{n+1}} \right] \cdot \left[ \frac{\delta\Delta\Omega^{i+1}}{\delta d_{n+1}^{i+1}} \right] = - \left[ \frac{G_1(\Delta\Omega^i, d_{n+1}^i)}{G_2(\Delta\Omega^i, d_{n+1}^i)} \right] \quad (18) $$

$$ \text{with} \quad \begin{cases} |\Delta\Omega^{i+1}| = \Delta\Omega^i + \delta\Delta\Omega^{i+1} \\ d_{n+1}^{i+1} = d_{n+1}^i + \delta d_{n+1}^{i+1} \end{cases} \quad \text{and} \quad \begin{cases} |\Delta\Omega^0| = 0 \\ d_{n+1}^0 = d_n \end{cases} $$

## 2 THE CONTROL ALGORITHM

This chapter develops a control strategy to determine the command laws to obtain an hydroformed part with better quality. This objective consists to include a control loop inside the incremental loop used in the finite element method. The developed method proposes a coupling between the response surface method and optimal control algorithm [3] and [4].

### 2.1 Problem formulation

In the case of hydromechanical deep drawing, one search to determine the evolution of the pressure cavity $p(\tau)$ and the blankholder force $F(\tau)$ (Figure 1), to minimize the development of the ductile damage, the problem can be formulated like a succession of optimization problems.

<figure><img src="image_1.png" /><figcaption>Figure 1: The hydromechanical deep drawing process with axisymmetric tools</figcaption></figure>
The control algorithm consists in determining along the successive step "j" of the simulation the values of the pressure $p^{j+1}$ and blankholder force $F^{j+1}$ that minimize the objective function with respect to the constraint inequality:

$$ \left\{ \begin{array}{l} \min_{(p^{j+1}, F^{j+1})} f^j(p^{j+1}, F^{j+1}) \\ \text{avec } g^j(p^{j+1}, F^{j+1}) \le \varepsilon_{\text{tol}} \end{array} \right. \qquad (19) $$

where $p^j$ and $F^j$ correspond respectively to the pressure and punch stroke for a loading factor equal to $\tau_j$, such as:

$$ \left\{ \begin{array}{l} p^j = p(\tau_j) \\ F^j = F(\tau_j) \end{array} \right. \qquad (20) $$

If one uses the objective function described in equation (21), the function that characterizes the damage for any spatial position in the component is given by:

$$ f^j(p^{j+1}, F^{j+1}) = \frac{1}{n} \sqrt{\left(\sum_{i=1}^{n} d_i^2(p^{j+1}, F^{j+1})\right)} \qquad (21) $$

where $d_i$ is the damage variable corresponding to the ith element.

The objective function can be expressed accordingly to the values $p^{j+1}$ and $\delta^{j+1}$ at any time of the calculation process. One thus will try to build an analytical approximation of the objective function $f^j(p^{j+1}, F^{j+1})$.

## 2.2 Meshless method and approximation of the objective function

As previously indicated, the proposed strategy is based on an approximation of the objective function through Response Surface Modelling. The response surface modelling uses the MLSA Method (Moving Least Square Approximation) [1,2].

One use a discretisation domain to compute the approximation of the cost function. In case of two parameters $\alpha$ and $\beta$ which respectively correspond to pressure and blankholder force (22), this domain corresponds to a square [0,1]x[0,1] where one realize the discretization by using regular intervals, see Figure 2. Each node represents the center of an element associated to a circle of radius $r$ representing its influence zone.

$$ \alpha = \frac{p^{j+1} - p_{\min}^{j+1}}{p_{\max}^{j+1} - p_{\min}^{j+1}}, \qquad \beta = \frac{F^{j+1} - F_{\min}^{j+1}}{F_{\max}^{j+1} - F_{\min}^{j+1}} \qquad (22) $$

with $p_{\min}^{j+1} < p^{j+1} < p_{\max}^{j+1}$, $F_{\min}^{j+1} < F^{j+1} < F_{\max}^{j+1}$
<figure><img src="image_2.png" /><figcaption>Figure 2: Meshless representation for 2 variables optimization problem</figcaption></figure>

The Moving Least Square Approximation (MLSA) is a technique which provides an approximation of a function with a set of weight functions, and a set of basis functions, avoiding the use of nodal points. At each point of the domain, a set of nodes $x_i$, $I=1,...,N$, is considered and the parameter associated with the approximation at node $I$ is denoted by $f_i$. The coefficients for the basis functions in this linear combination differ from point values and are computed using the weight functions and the nodal values of the function under consideration. For this, an equations set has to be solved at each point of the domain in order to get the coefficients for the basis functions.

The approximation $\tilde{f}(x)$ of the function $f(x)$ in the domain $\Omega$ is defined from a linear combination of basis functions $P(x)$ as :

$$ \tilde{f}(x) = \sum_{i=1}^{m} P_i(x) a_i(x) = P^T(x) a(x) \qquad (23) $$

where $m$ is the basis number of terms, $P(x)$ is the vector of the basis polynomial functions and $a(x)$ is the vector of the coefficients which are functions of the spatial coordinates $x$.

The influence zone of a node $I$ is defined by a weight function which can take the following form:

$$ w_I(x) = \begin{cases} 1 - \left( \frac{\|x - x_I\|}{r} \right)^2 & \text{if } \|x - x_I\| \le r \\ 0 & \text{else} \end{cases} \qquad (24) $$
The coefficients $\mathbf{a}(x)$ are obtained by minimizing the square of the difference between the local approximation $\tilde{f}(x)$ and the value $f_1$ of the function $f(x)$ corresponding to $x = x_1$. Therefore, it is necessary to minimize the following functional:

$$J(x) = \frac{1}{2} \sum_{i=1}^{n} w_i(x) \left[ \sum_{i=1}^{m} p_i(x_i) a_i(x) - f_1 \right]^2 \quad (25)$$

This equation can be rewritten in matrix form:

$$\mathbf{J}(x) = \frac{1}{2} \left( \left[ \mathbf{P} \right]_{n \times m} \mathbf{a}(x) - \mathbf{F} \right)^{\mathrm{T}} \left[ \mathbf{W}(x) \right] \left( \left[ \mathbf{P} \right]_{n \times m} \mathbf{a}(x) - \mathbf{F} \right) \quad (26)$$

where $\mathbf{F}^{\mathrm{T}} = \{F_1, F_2, \dots, F_n\}$, is the vector values function at the nodal points. $[\mathbf{P}]$ is a polynomial vector.

## 2.3 Presentation of the control algorithm

This algorithm is inspired by the algorithms used for buckle control of a actuator. Methodology employed allows to minimize an objective function at certain moments of the process. To present the algorithm, we will take again the development with two parameters. One will consider that the process duration is described thanks to the loading factor $\tau$ such as $\tau \in [0,1]$ and the algorithmic diagram corresponding to Figure 3.

* **1st stage :** One divides the loading factor $\tau$ into $N_{sub}$ subintervals, noted such as in the case of the example represented in Figure 3, ($N_{sub}=5$).

* **2nd stage :** Initially, one begin calculation in the first subinterval. One realize several simulations on the subinterval by imposing different values different of pressure $p^0$ and blankholder force $F^0$. One obtain thus various evaluations of the objective function.

* **3rd stage :** From the various evaluations of the objective function objective $f_i^j$, one build an analytical approximation $\tilde{f}^j$ of the objective function $f^j$ given from the equation (21).

* **4th stage :** One uses a genetic algorithm [7] to determine the optimum parameters optimum in order to obtain the global minimum of the analytical approximation of the objective function. One thus obtains the first value of the pressure curve $p^1_{opt}$.

* **5th stage :** One realize a last simulation corresponding to a duration $I_{sub}^0$ while imposing as values for the pressure $p^0$ and punch stroke $\delta^0$ the optimal parameters given in stage 5. Then the first step of control is realized and one has obtained the first values for the command laws.
## 6th stage :

One places now at the following step, therefore with the following subinterval $I_{sub}^1 = [0.2,0.4]$. One takes the same step by repeating stages 2 to 6.

<figure><img src="image_3.png" /><figcaption>Figure 3 : Schematic representation of the control algorithm over one normalized time interval</figcaption></figure>

This algorithm thus allows to build step by step the commands laws of pressure, to see Figure 4. $p_{opt}^j$ is the optimal value of pressure determined for a particular value of loading factor $\tau^j$, so that:

$$p_{opt}^j = (p_{max}^j - p_{min}^j)\alpha_{opt} + p_{min}^j \qquad (27)$$
<figure><img src="image_4.png" /><figcaption>Figure 4 : Build step by step the pressure curve</figcaption></figure>

## 3 DEEP DRAWING OF A CYLINDRICAL CUP WITH A COUNTER PRESSURE IN THE DIE CAVITY

The test reported in this paragraph consists with deep drawing of a circular blank with a diameter 330 mm and thickness 0,6 mm using a semi-spherical punch. The geometry of the tools is described **Figure 5** according to a cross-section.

<figure><img src="image_5.png" /><figcaption>Figure 5 : Geometry used for the deep drawing of a cylindrical cup with a counter pressure</figcaption></figure>
For this test, one uses a XES steel whose hardening law is given by the relation (16) and whose material parameters is indexed in the **Table 2** :

<table><caption>Table 2: material Parameters for XES</caption><thead><tr><th>Material</th><th>E (MPa)</th><th>ν</th><th>σ<sub>y</sub> (MPa)</th><th>K (MPa)</th><th>n</th></tr></thead><tbody><tr><td>XES</td><td>210000</td><td>0.3</td><td>169</td><td>366,55</td><td>0.406</td></tr></tbody></table>

Because of symmetry of the problem and supposing an isotropic material behaviour, only one quarter of the problem is discretized by finite elements. The flange is discretized by shell elements with three nodes, the complete mesh for the quarter of tube comprises 883 elements. In this sample, the loading factor τ is subdivides into 30 subinterval (N<sub>sub</sub>=30) and one divides space domain, see **Figure 2**, into 4 intervals.

We can see in **Figure 6** the evolution of the objective function $\tilde{f}^j$ at different loading factor $\tau_j$, this function only depends on pressure and blankholder force. One represents by red points, **Figure 6**, which one notes of A with I the couple solution obtained with the Genetic algorithm. So, the control algorithm permits to build step by step the command laws for pressure and blankholder force.

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6 : Evolution of the objective function $\tilde{f}^j$ at different loading factor $\tau_j$</figcaption></figure>

The results concerning the pressure profile and those concerning to the blankholder force profile are given in Figure 7.

<figure><img src="image_8.png" /><figcaption>Figure 7 : Pressure and blankholder load vs. Loading factor obtained from meshless approximation</figcaption></figure>

The damage variable contours for the deep drawn cup after application of proposed control algorithm are represented in Figure 8. One can observe that the use of the proposed control procedure really leads to obtain a better distribution of damage during tube hydroforming.
<figure><img src="image_9.png" /><figcaption>Figure 8 : Damage contours at the end of the hydroforming process corresponding to the meshless approximation</figcaption></figure>

## 4 CONCLUSION

A new strategy and the associated procedures are proposed for the optimal control of hydroforming processes. This strategy is based on the combined use of incremental optimization over the processing time divided into subintervals. So the algorithm that has been proposed consists first to build an approximation of the variables entering in the control functions and then to incrementally determine the process parameters over the time interval permitting to effectively build in an incremental way the command laws for the process variables. It is really a new approach in the sense that these process variables could be inserted in the command of the hydroforming machine tool controller system to get the required component. The obtained result proves the efficiency of the proposed approach to control damage contours and values at the end of process.

## REFERENCES

[1] G. Ayad, A. Lejeune, T. Barriere, J.C. Gelin, *An optimization strategy for the determination of material and process parameters to avoid segregation defects during metal injection*, 8th Int. Conf. In Numerical Methods in Forming processes, NUMIFORM 2004, Ohio, U.S.A. 13-17 June 2004, 2004.
[2] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl. *meshless methods :An overview and recent developments*. Comput. Methods Appl. Mech. Engrg. Vol. 139, pp 3-47, 1996.

[3] J.C. Gelin, C. Labergere, *Numerical design and optimal control for sheet forming and tube hydroforming processes*, 7th Int. Conf. In Numerical Methods in Forming processes, NUMIFORM 2001, Toyohashi, Japan, 18-20 June 2001, ed. by K.I. Mori, pp 897-902, 2001.

[4] J.C. Gelin, C. Labergere, *Application of optimal design and control strategies to the forming of thin walled metallic components*, J. Mater. Process. Technol., Vol.125-126, pp 565-572, 2002.

[5] M. Koç, T. Altan, *Prediction of forming limits and parameters in tube hydroforming process*, Mach. Tool Des. Res., Vol. 42, pp. 123-138, 2002.

[6] Lemaitre J., Chaboche J.L., *Mécanique des Matériaux Solides*, Dunod, 1985.

[7] Michalewicz Z., Nazhiyath G., *Genocop III, A Co-evolutionary Algorithm for Numerical Optimization Problems with Nonlinear Constraints*, Proceedings of the 2nd IEEE International Conference on Evolutionary Computation, Vol.2, pp.647-651, 1995.