Contents lists available at SciVerse ScienceDirect

Computers and Structures

journal homepage: www.elsevier.com/locate/compstruc

Computers & Structures

Finite element modelling and optimisation of net-shape metal forming processes with uncertainties

H. Ou a,b'*, P. Wanga, B. Luc, H. Longd

a Department of Mechanical, Materials and Manufacturing Engineering, University of Nottingham, Nottingham NG7 2RD, UK b State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing 400044, China c Department of Plasticity Forming Engineering, Shanghai Jiao Tong University, Shanghai 200030, China d Department of Engineering and Computing Sciences, Durham University, Durham DH1 3LE, UK

ARTICLE INFO ABSTRACT

The importance of the stochastic nature of forging operations, especially for net-shape accuracy, has been widely recognised. However it is a under-researched area in terms of effective methods and computational tools for easy industrial implementation and applications. In this research, Monte Carlo simulation, response surface method and most probable point analysis are used to quantify probabilistic characteristics of the shape and dimensional errors in net-shape metal forming processes. A two-step optimisation approach is developed to minimise systematic errors using a direct compensation method for die shape modification and to reduce random variations through a control variable method. Two industrial based case studies, i.e. forging of a 2D aerofoil component and forward extrusion of a cylinder are conducted with good results obtained for much improved accuracy. This approach can be easily applied to general metal forming processes for both 2D and 3D cases where achievable accuracy is an essential criterion.

© 2011 Elsevier Ltd. All rights reserved.

Article history: Received 22 October 2009 Accepted 10 October 2011 Available online 5 November 2011

Keywords: Net-shape forming Aerofoil blade Extrusion

Stochastic modelling Optimisation Reliability-based analysis

1. Introduction

Metal forming processes such as forging and extrusion are competitive routes for manufacturing structural products due to the advantages of enhanced mechanical properties and microstructures as well as substantial benefits in reduced production cost. The quality and accuracy of metal forming processes are dependent upon the material properties, die and tool behaviour as well as process conditions [1-3]. The shape and dimensional accuracy specification of the formed components is an essential requirement in net-shape metal forming process design and production together with other criteria such as required material microstructure and mechanical properties, minimal material loss. Due to the variability of material flow behaviour, interfacial properties, as well as uncertainties in process and operational conditions, the shape and dimensions of formed components inevitably exhibit considerable variations and scatter. To achieve repeatable shape and dimensional accuracies of net-shape formed components, trial and error based iterations are normally used in practical metal forming process design and validation. This often results in prolonged lead time and increased cost.

* Corresponding author at: Department of Mechanical, Materials and Manufacturing Engineering, University of Nottingham, Nottingham NG7 2RD, UK. Tel.: +44 (0)115 8467391; fax: +44 (0)115 9513800.

E-mail address: h.ou@nottingham.ac.uk (H. Ou).

0045-7949/$ - see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2011.10.014

Finite element (FE) method is now widely used to simulate forging and other metal forming processes for a wide range of industrial applications [4]. In the past decade FE based modelling and optimisation methods have been developed for improved material deformation and material properties and reduced forming forces using gradient based methods, evolutionary algorithms and response surface based approximations [5-9]. In metal forming related shape optimisation, gradient based shape optimisation methods were developed for preform design [10-12]. Preform design optimisation is obtained by varying the geometry of the billet. On the other hand, die shape optimisation for net-shape specification is aimed towards the specified shape and dimensions by small perturbation of the die shape so that the shape and dimensional errors of the formed component due to, for example, die elastic deformation and thermal distortions are minimised. This approach has been used for net-shape forging operation and control of springback in sheet metal forming [13-16].

However, the aforementioned FE based simulation and optimisation methods are deterministic without the consideration of the variability of material properties such as material flow stresses and the uncertainties of process conditions including temperature, stroke, friction and initial set-up. In real metal forming production, these uncertainties can cause increased possibility of die failure and especially defects of formed components and scatter in shape and dimensions. Therefore the goal is to achieve minimised process variability and maximised shape and dimen-

sional accuracy of the formed components. To quantify and simulate process uncertainties, the stochastic finite element method (SFEM) may be used [17]. SFEM incorporates uncertainties by defining the deterministic and the uncertainty terms separately in the formulation of global governing matrix equations. To solve the system of random algebraic equations, different SFEMs such as the Taylor series expansion based perturbation method and spectral stochastic finite element method (SSFEM) may be applied to such stochastic problems but this often requires additional substantial computing resource [18,19]. Currently there is no SFEM based software available that can be used for stochastic metal forming simulation. Instead a decoupled approach for FE simulation and stochastic computation may be developed for stochastic process modelling using such as Monte Carlo simulation (MCS) for both statistical characterisation and quantification of stochastic propagation [20,21]. However, it becomes computationally demanding when stochastic based modelling and optimisation are taken into account. To deal with optimisation of uncertainties, it is important to develop efficient and robust techniques for uncertainty characterisation and stochastic modelling using reliability-based optimisation (RBO) and robust design methods (RDM) [22].

Little research has been reported in quantifying and incorporating uncertainties in metal forming simulation and optimisation. Grandhi et al. first quantified the effect of affecting parameters on the uncertainties of forging using Design of Experiment (DOE) approach in order to control the geometric variations and to optimise preform design using a stochastic response surface (SRSM) and a reliability-based optimisation (RBO) method [23-25]. Similar research in prediction of stochastic responses in cold forging and stamping processes were also reported [26,27]. Recently a non-intrusive stochastic Galerkin (NISG) method was proposed involving piecewise continuous interpolation of the probability distribution function over the support space with deterministic function evaluations at the element integration points [28]. This method provides accurate and efficient stochastic response estimation of different metal forming processes and may be used as an attractive alternative to intrusive SFEM techniques but it suffers the curse of dimensionality for high dimensional uncertainty problems. In sheet metal forming, the quality and final shape of the formed parts are influenced by the uncertainty of material and interfacial properties, blank holder force and other process conditions. Different methods including Monte Carlo based metamodels [29], adaptive Monte Carlo simulations (AMCS) [30] and polynomial chaos expansion [31] were reported to quantify the probabilistic characteristics of uncertainty propagation in sheet metal forming processes. In net-shape metal forming processes, most deterministic methods were employed to minimise systematic errors [13-16]. Although many stochastic approaches were developed to quantify probabilistic characteristics and uncertainty propagations of metal forming processes [23-31], little research has been reported on FE based stochastic optimisation for net-shape metal forming requirements.

The aim of this paper is to develop a generic approach for quantifying probabilistic characteristics and optimising the forming process that both minimises the systematic errors and reduces the random variations. Section 2 outlines forming error specification, quantification of probabilistic characteristics using Monte Carlo simulation (MCS), response surface method (RSM) and most probable point (MPP) analysis and a two-step optimisation approach. Sections 3 and 4 present FE modelling and stochastic optimisation results of two case studies including a 2D aerofoil forging and a forward extrusion process. Section 5 summarises a few concluding remarks of this research.

2. Net-shape forming optimisation with uncertainties

2.1. Specification of shape and dimensional errors

The shape and dimensional errors of net-shape formed parts may be classified by systematic deviations and random variations. In the case of simple upsetting of a cylinder as shown in Fig. 1(a), systematic deviations of the upset cylinder may be attributed to a specific trend in metal forming such as the die elastic deflection, whilst random variations are largely due to the variability of material and interfacial properties and fluctuations of operational conditions. The shape and dimensional error of a formed part is defined as the discrepancy between the actual and the nominal surface profiles. The calculation of the shape and dimensional error at a particular point, i, may be given by

Ad(x,-) = S(Xi)~ S0(Xi) (1)

where 5(xi) and d0(xi) are the actual and nominal profiles. A tolerance measure A is often used to ensure that necessary shape and dimensional requirements are satisfied for a specific metal forming operation, i.e.

|Ad(Xi)| = |d(Xi)- d0(Xi)| < A

A limit state function g(xi), as an implicit function, can be written as follows

g(Xi) = A - |d(Xi) - d0(Xi)|

The probability of failure P(g(xi) 6 0) may be given in Fig. 2. In consideration of uncertainties, the shape and dimensional error Ad(Xi) is a random response of the final shape of the formed component d(Xi), an implicit function of factors such as material deformation, temperature and other process conditions, and can be obtained from FE modelling. Its statistical mean and standard deviation is denoted by i and axi, respectively. Assuming there are N random responses Sj(xi), the shape and dimensional error Adj(xi) can be seen as the deviation from the desired shape S0(xi) and the mean square deviation (MSD) can be defined as follows

MSD = NE (d(xi) - d°(xi))

and further derived as

MSD = ^ E { [Si(x>) - ftj + - dVi)] } j=1

= + (lXi - <5°(x*)) 2 where lx = 1Ejij) and al = (dKxO - lx.)2.

From Eq. (4b), it is clear that MSD is made up of two parts: the deviation of the mean from the target dimension ix - d0(x,) which

upper die

upper die

initial die shapes

(a) Systematic error in upsetting (b) Systematic error minimisation

Fig. 1. Definition of shape and dimensional error.

Fig. 2. Probability of failure due to stochastic variations.

is a measure of the systematic error of the formed part and the shape and dimensional error variance o"2 from the mean, an indication of random variations. Therefore the objective is to minimise the systematic error and at the same time to reduce random variations. This can be achieved in two steps: (1) adjust the mean of the shape and dimensional error, to fall on the target d0(xi), i.e.

- d0(xi) = 0 as illustrated in Fig. 1(b), (2) reduce the random variations that causes the dispersion of the shape and dimensional error about its mean as shown in Fig. 3.

2.2. Quantification of probabilistic characteristics

Different methods including the first- (FORM) and second-order reliability (SORM) methods, stochastic expansions, response surface method (RSM), Monte Carlo simulation (MCS) and stochastic finite element method (SFEM) may be used to quantify probabilistic characteristics and to carry out reliability design with varying degrees of precision, efficiency and robustness for different stochastic problems [19,20]. In metal forming processes, FE based modelling and optimisation involves large material deformation, complex workpiece and tool interaction and thermo-mechanical coupled analysis. Computational efficiency and easy implementation are two important considerations. In this research, Monte Carlo simulation (MCS), response surface method (RSM) and most probable point (MPP) analysis are used to quantifying probabilistic characterisation of net-shape forming processes. They are summarised in the following subsections.

05 0.45 0 i 0.35 0.3 0 25 0,2 0 15 01 0 05

LSL îigggjmean i /f№i // 'W / / I \ \ / I \ / M 1 \ //ii\ / / i i \ / / i \ \ / i \ / / i i \ / / 1 I \ USL / \ i \ i \ / Initial \ / M*) \ i \ i \ i \

: / / i \ \ / I i 1 \ / i \ \ f 1 1 \ \ ' f -PBîji^^- \ ^ / sm \ / i \ i \ / \ / \ / \ V \

2.2.1. Monte Carlo simulations (MCS)

MCS is a class of computational algorithms that rely on repeated random sampling to compute results. It follows a particular procedure by: (1) defining uncertainties of the problem and quantifying the probabilistic characteristics of the random variables in terms of their probability density functions (PDFs) and cumulative distribution functions (CDFs); (2) evaluating the problem deterministically for each set of realisations of all random variables; (3) extracting probabilistic information from results to calculate the required statistics: the histogram, the PDF and corresponding CDF, and finally to obtain the probability of reaching the design target in consideration of various performance criteria [17].

As the first step of MCS, random variables of material properties and process conditions in a metal forming process may be generated using the inverse transform method [17,19,32]. In this method, the pseudorandom number ui uniformly generated over (0 6 ui 6 1) is equated to the CDF of the random variable FX(xi). Thus the realisations of the random variables can be calculated as follows

xi = Fx\ui)

where F^ (xi) is the inverse cumulative function. If the random variable X is normally distributed, i.e. N(iX, aX), the random variable may be obtained by

Xi = 1X + Ox U 1(u)

where U-1(ui) is the inverse of the CDF of the standard normal variable. After the generation of random variables, deterministic metal forming simulations are conducted and the probabilistic characteristics of the system response in forms of PDF and CDF are derived. The accuracy and efficiency of MCS are dependent upon the number of simulations and the following equation may be used as an efficiency estimator [17] ™ p > - ffiffi

where COV(pf) is the coefficient of variation, pf is the estimated probability of failure and N is the number of simulations. Especially

Fig. 3. Two-step stochastic optimisation for net-shape forming operations.

Fig. 4. Flow chart of the stochastic optimisation.

for small failure probability problems, a large number of simulations are required with a relatively small value of COV(pf). In metal forming simulation, however, it may be impractical to run a large number of simulations using direct MCS. Methods including importance sampling method [17,32] and adaptive Monte Carlo (AMC) simulation [30] may be used to improve both simulation efficiency and accuracy. In this research, different numbers of simulations are conducted in the two case studies and the results are used to compare with RSM and MPP methods.

2.2.2. Response surface method (RSM)

RSM is also used in reliability analysis of metal forming processes [8,26,27,29,33]. In RSM, system probabilistic response is represented by linear or quadratic approximation. This is sufficient for most metal forming processes due to relatively small variability of material and interfacial properties and process conditions. The RSM procedure involves the following steps: (1) selecting a number of the most influential random variables for system response;

Table 1

Mechanical and material properties of the workpiece and dies.

(2) evaluating the system response using deterministic analysis for all the sets of values of the selected random variables; (3) constructing linear or quadratic approximation to represent the system response by regression analysis with data collected from the deterministic analyses; (4) once the approximate close-form representation is obtained the system probabilistic characteristics are computed and the probability of failure is evaluated using such as MCS.

In net-shape metal forming, the shape and dimensional errors Adi(xi) may be represented as a polynomial function of the random variables in the following form,

Adi(xi ) = Adi(x1, x2,..., xp)

= bi,0 +J2 bi.«x>< +J2 kp+xxl + •••; (i = 1; 2, ..., n) (8)

k=1 k=1

where constants bi0, bitk.....bip+k,... are all m x 1 vectors and xk is a

m x 1 vector denotes the selected m sets of random values for random variable xk. n denotes the number of check points of the forged component and p the number of chosen random variables for probabilistic analysis. Coupled terms of different random variables are not involved if they are assumed independent of each other. The

regression analysis is to find out the unknowns (bi0,bik.....bip+k,...)

through the obtained m results. For each result, there is an error e. m sets of results collected from m deterministic FE analyses give m equations as follows in the matrix form,

Ad = Xb + e (9)

Note bold X in the matrix form of the above equation contains all terms in Eq. (8): constant term (an m x 1 vector with unit values), random variable terms and squared random variable terms and so on. Least square method may be used to minimise the norm of the vector of error e, i.e. to minimise || Ad - Xb||. This requires XT(Ad - Xb) = 0 or XTXb = XTDd. Therefore,

b = (XTX) -1Xt Dd (10)

However, as the normal matrix XTX is often ill conditioned, instead of solving Eq. (10) directly a QR (orthogonal and triangular) decomposition in MATLAB is used for the solution of unknowns

Young's modulus E (GPa) Poisson's ratio m Density q (kg/m3) Specific heat Cp (J/kg °C) Conductivity k (W/m °C) Thermal expansion (1/°C)

Ti6Al4V 50 0.3 4400 700 15.0 1.0 x 10"5

A199.5 69 0.3 2707 896 204.0 8.418 x 10~5

Dies 206 0.3 7900 500 48.1 1.2 x 10~5

Fig. 6. Flow stresses of Ti-6Al-4V.

Fig. 7. Location of check points for thickness measurement.

(bi0,bik.....bUp+k,...). Then similar to MCS, MATLAB is used to

generate random values for the chosen uncertain parameters

(x1,x2.....xp) and substitute them into Eq. (8) to complete the RSM

based stochastic simulations of the metal forming process. These results are used to obtain the system probabilistic characteristics: the histogram, the PDF and its corresponding CDF and to estimate the probability of failure.

Table 2

Case study 1: selected sets of random parameters for RSM.

Temperature T (°C) Friction F Stroke D (x103m)

1000 0.290 1.10

1000 0.290 1.00

1000 0.110 1.10

1000 0.110 1.00

900 0.290 1.10

900 0.290 1.00

900 0.110 1.10

900 0.110 1.00

950 0.200 1.05

950 0.290 1.10

950 0.290 1.00

950 0.110 1.10

950 0.110 1.00

Temperature T (°C) Friction F Stroke D (x103m)

1000 0.200 1.10

1000 0.200 1.05

1000 0.200 1.00

900 0.200 1.10

900 0.200 1.05

900 0.200 1.00

950 0.200 1.10

950 0.200 1.00

1000 0.290 1.05

950 0.290 1.05

900 0.290 1.05

1000 0.110 1.05

950 0.110 1.05

900 0.110 1.05

Fig. 8. Distributions of von Mises stress; strain and temperature.

Fig. 9. Forging force and displacement results.

To construct the approximate equations using the standard three level (low-central-high) approach, the total number of simulations to evaluate the performance function is 3p. Whilst this approach may be adequate to extract necessary system response data, it may not be enough to provide sufficient information for estimation of such as the effect of tail distributions. Therefore, the Latin Hypercube Sampling (LHS) as a stratified sampling technique may be used for the selection of simulation points within the random space [19,34]. In LHS, the distribution for each random variable is divided into n intervals and one random value is selected from each interval corresponding to its probability density. By repeating this procedure for all p random variables, n pairs of p random variables combined through a non-overlapping manner are used for multivariate sampling. This approach ensures that the full range of the probability distribution is sampled and hence results in relatively small variance in the system response. In addi-

tion, optimal LHS procedure may be used for achieving well-stratified sampling [19,35].

2.2.3. Most probable point (MPP) analysis

With suitable approximation such as RSM, the system probabilistic response may be estimated using the first-order Taylor series expansion at the mean value point as given by

ASi(Xi ) = Ddi (ix)

-A dASi ,

@Ad-(Xk x 1

1 A P d2 ASi

+ 2 h h @XkdXj

xk x lXk Xj x 1

where i is the mean value of xk. Truncating the series to linear terms, the first-order approximate mean and standard deviation of Adi(xi) may be obtained as follows

1as ASi ( 1x) = ASi ( 1X1 ; 1X2 ; ••• ; 1x„ ) t

OAS ffi

(12a) (12b)

In Eq. (12b), the random variables are assumed to be uncorrelated. This approach is often referred as mean value first-order second-moment (MVFORM) method. Using the MVFORM approach, direct correlation may be obtained between random variables and system response of the shape and dimensional error of formed parts. However, the linearisation of the system probabilistic response may result in considerable inaccuracy especially for low probability of failure. Improved results may be obtained by using the most probable point (MPP) also called Hasofer-Lind (H-L) method [17,19,36]. In the MPP analysis, the original random variables xi are transformed into normalised ones xj with zero mean and unit standard deviation as given by

(c) Quadratic approximation using RSM Fig. 10. Comparison of thickness (at the 5th check point) from MCS and RSM (m).

Table 3

Case study 1: stochastic optimisation for average thickness error d

(unit: x!0~3 m).

Check point i

Desired thickness d0 Initial thickness if Thickness, 1st iteration i(1) Thickness, 2nd iteration i(2) Thickness, 3rd iteration i(3

Standard deviation, 3rd iteration ct,(3)

0.96494

1.0909

1.1995

1.2535

1.2667

1.2247

1.1412

1.0516

0.95162

0.86397

0.79441

1.0894 1.2321 1.3244 1.4335 1.4678 1.3973 1.2819 1.2024 1.1020 1.0026 0.90261

0.89961

1.0323

1.1803

1.2528

1.2785

1.2441

1.1472

1.0381

0.89912

0.78771

0.67822

0.94920

1.0730

1.1814

1.2463

1.2605

1.2199

1.1552

1.0735

0.94086

0.86135

0.79789

0.95532

1.0949

1.1868

1.2477

1.2675

1.2270

1.1461

1.0617

0.97192

0.85362

0.79754

0.020981 0.023407 0.018563 0.014960 0.013772 0.014702 0.015310 0.017823 0.021621 0.020494 0.021756

0.15138

0.05298

0.01287

0.00937

■Initiaidie —»—1st die compensation optimized die --2nddie compensation

Fig. 11. Thickness error results of die shape compensation iterations.

Therefore the original limit state function as given in Eq. (3) is transformed from the original coordinate system into the transformed coordinate system, where the minimum distance also referred to as the reliability index b from the origin of the transformed coordinate system to the limit state function represents the most probable point (MPP) of failure. Hence this approach is also called MMP analysis and commonly used in reliability analysis. The reliability index b may be obtained by minimising b(xj) subject to the con-

straint of the limit state function g(xi) = 0. Using the method of Lagrange multipliers and linearisation of the limit state function [19], the reliability index b can be calculated by [17]

g(xr) -s » ^ i

Mx?) i

is the ith partial derivative of the original coordinates

(x£) and the asterisk means that it is evaluated at the most probable point xj*. The design point in the transformed coordinate system is given by

xi = -a iß (i = 1,2,..., n) where ai =-

is the direction cosine along the

coordinate axes x£. Using Eq. (13), the most probable point in the original coordinate system may be obtained by

x* = ix, - Oiax.b (i = 1, 2,..., n) (16)

By defining the appropriate limit state function of Eq. (3) and assuming an initial design point at the mean values of the random variables, an iterative procedure may be used to compute the most probable point (MMP) of failure x* [17,19].

2.3. Two-step stochastic optimisation

The aim of metal forming modelling and optimisation for net-shape accuracy involving probabilistic characteristics is to minimise the systematic errors and to reduce the random variations or control the probability of failure. As given in Eq. (4), this may be achieved by first minimising the mean values of the shape

and dimension measures, i.e. (jxx, - d0(xi^ close to zero and then reducing the dispersion of the dimensions about the mean values, i.e. ax2 as illustrated in Fig. 3.

In the first step, the systematic response (blue dashed) obtained from RSM based on the initial die design due to such as die elasticity and thermal distortion is minimised by using a direct die shape compensation method [13,14]. As a number of check points are used to measure the shape and dimensional errors along the cross-section of a formed part, it may be difficult to ensure that the mean values at all check points are all minimised to their targets simultaneously. Instead, the average error of the mean values is used to measure the specified tolerance A. Hence, the following objective function is used to measure the overall error as

1 £ (I- - d0)2 < A

In the second step, a control variable method is used to reduce the variance, o^ [32]. For one dimensional case, let Z be an unbiased estimator of lZ, to be obtained from a simulation run. A random variable C is called a control variable for Z if it is correlated with Z and its expectation, lC, is known. The control variable is used to construct an unbiased estimator of lZ with a variance smaller than that of Z

Zk = Z x k(C x 1c)

where k is a scalar parameter, called the linear control variable. The variance of Zk is given by

Var(Z>) = Var(Z) x 2kCo v(Z, C) + k2 Var(C) (19)

Consequently, the value k* that minimises Var(Zk) is k* = Cov(Z, C)/Var(C) (20)

and the minimum variance is

Var(Zk) = (1 x pjc)Var(Z)

where Cov(Z, C) is the covariance of two random variables Z and C, pZC is the correlation coefficient of Z and C and they have the relationship as follows

Cov(Z, C) E((Z x iz)(C x ic))

Eqs. (18)-(22) can be extended to the case of multiple control variables. It is noted from Eq. (21) that the larger |qZC| is, the greater is the variance reduction. In metal forming processes, the shape and dimensional errors are closely correlated with material and interfacial properties as well as process conditions. By selecting the most influential parameter or parameters and defining the amount of reduction of the variance of each parameter that can be achieved, the overall variance reduction of the system response, i.e. the shape and dimensional errors of formed part, can be obtained using the

(e) Thickness after variance reduction

Fig. 12. Comparison of thickness results in two-step optimisation (m).

control variable method. Therefore, the two-step optimisation procedure can be described as follows:

(1) Defining ranges for the random variables and selecting sets of random parameters.

(2) Finding the mean values and standard deviations of the shape and dimensional errors at the cross-section check points or the overall error using MCS, RSM or MMP methods.

(3) Calculating the shape and dimensional errors and checking them whether within the tolerance.

(4) If not, modifying the die shapes using the direct die shape compensation method and repeating steps (2) to (3).

(5) If yes, calculating the probability of failure and checking whether it meets the reliability target, if not, goes back to step (1) to change the ranges of the random variables to achieve the variance reduction. The flow chart is shown in Fig. 4.

3. Case study 1: forging of 2D aerofoil section

3.1. FE modelling

In this case study, hot forging of Ti-6Al-4V into 2D aerofoil sections is simulated using ABAQUS. The meshed models of the work-piece and two dies using CPE4T elements are shown in Fig. 5. The FE analysis is defined into four steps: forging, removing upper die, removing lower die and cooling. The time in each step is 0.2, 0.01,

0.01 and 2200 s, respectively. The workpiece is defined to be elastic-plastic and the forging dies are defined to be elastic. Young's modulus, thermal conductivity and specific heat of the workpiece and dies are shown in Table 1. The flow stresses of the workpiece as a function of strain, strain rate and temperature are shown in Fig. 6. The thickness errors of the forged aerofoil section are measured by the difference between the specified and actual thickness at each check point of the cross-section as shown in Fig. 7.

To quantify the probabilistic characteristics of the thickness error of the forged aerofoil section, three random variables, i.e. the initial temperature (T) of workpiece, the friction coefficient (F) and the stroke length (D) are selected. All three random variables are assumed to be normally (Gaussian) distributed and independent of each other. Their mean values are iF = 0.2, iT = 950 °C and iD =1.05 x 10-3 m and their standard deviations are oF = 0.03, oT = 16.7 °C and oD = 1.67 x 10-5 m, respectively. For Monte Carlo simulation (MCS), MATLAB function randn(3000,3) is used to generate a 3000 x 3 matrix containing three columns of independent data to represent the random variables. For response surface method (RSM), 27 (m = 33) sets of random parameters are used as listed in Table 2.

3.2. Material deformation

Fig. 8 shows the von Mises stress, the equivalent plastic strain and the temperature distributions relating to stochastic mean val-

Fig. 13. Thickness results (iAd + 3oAdi) at all check points between USL and LSL.

123456789 10 11 Initial die ■ Optimised die ■ optimised die & variance reduction

Fig. 14. Changes of standard deviation (oAdi) at all check points.

workpiece

extrusion die

Fig. 15. FE mesh of the workpiece and extrusion die.

ues using the initial die design. Similar results are observed with the optimised die design because to achieve minimised systematic error only small amount of die shape modification is necessary. Fig. 9 shows the forging force and displacement relationship during forging, in which the two solid lines denote the cases of using the initial die and the optimised die shapes with mean values of the random variables and the two dashed lines denote the cases of the optimised die with two limit sets of the random variables, i.e. one with the highest temperature and largest stroke length which results in the thinnest aerofoil blade, the other with the lowest temperature and the shortest stroke length which leads to the thickest blade.

3.3. Probabilistic characterisation

3.3.1. Monte Carlo simulation (MCS)

Using MCS, after running 3000 FE analyses deterministically, the probabilistic characteristics of the thickness errors are extracted. The thickness of the aerofoil cross-section at the 5th check point as defined in Fig. 7 is shown in Fig. 10(a). It is noted that the results of thickness are also normally distributed. Its statistical

Fig. 16. Flow stresses of Al199.5.

Fig. 18. Extrusion force history results.

Fig. 17. Distributions of (a) von Mises (MPa) and (b) plastic strain.

mean and standard deviation are also obtained and shown in the figure.

3.3.2. Response surface method (RSM)

After running 27 (m = 33) FE analyses deterministically for each set of the selected random variables listed in Table 2, RSM approximate representations of the thickness errors of the forged aerofoil section are derived. Using the linear terms in Eq. (8) the thickness response at the 5th check point can be approximately represented as

<5(xs) - 2.5764 x 10-3 + 2.4253 x 10-5F - 2.5572 x 10-7T

- 0.82908D (23)

Its mean and standard deviation are calculated to be 15 - 1.4678 x 10-3 m and o5 - 1.45 x 10-5 m, respectively. When taking all the linear and squared terms in Eq. (8), regression analysis gives the quadratic approximation at the 5th check point as

<5(xs) - 1.9879 x 10-3 + 7.2401 x 10-5F + 9.7674 x 10-7T

- 0.82908D - 1.2037 x 10-4F2 - 6.4867 x 10-10T2 (24)

For comparison purpose, the same 3000 realisations of the three random variables generated in Section 3.3.1 are substituted into Eqs. (23) and (24) and 3000 response results are calculated. Their histograms and fitted PDF, with the response mean and standard deviation at 5th check point are shown in Fig. 10(b) and (c). It shows that the thickness responses are also normally distributed. It is noted that there is no significant difference between the linear and quadratic approximations as the coefficient of the squared stroke length (D2) term equals to zero. To check the correlation between the thickness response and random variables, the correlation coefficients at the 5th check point are calculated to be qd5 F = 0.0502,

Ps5,t = -0.2941 and pd5D = -0.9535, respectively. This suggests that the stroke length has the biggest effect on the aerofoil thickness, whilst the effect of friction coefficient is negligible.

3.3.3. Most probable point (MPP) analysis

Using the MPP analysis, reliability analysis is carried out using the linear approximation of the aerofoil thickness given in Eq. (23). With a specified tolerance of A - 5.0 x 10-5 m at all check points, the limit state function at the 5th check point can be given as

g(xs) - 0.05 x 10-3 - (1.3097 x 10-3 + 2.4253 x 10-5F

- 2.5572 x 10-7T - 0.82908D) (25)

Table 4

Case study 2: selected sets of random parameters for RSM.

Friction Flow Heat transfer Friction Flow Heat transfer

F stress coefficient H F stress coefficient H

scaling (kW/m °C) scaling (kW/m °C)

factor w factor w

0.0968 1.0069 950.5 0.0961 1.0213 987.1

0.1020 1.0121 860.0 0.0842 1.0033 828.0

0.0845 1.0149 841.1 0.1052 0.9955 1019.5

0.1109 1.0219 866.6 0.1071 0.9847 893.9

0.1191 0.9891 798.7 0.0866 1.0064 811.3

0.0634 1.0175 917.7 0.1025 0.9975 766.6

0.0806 0.9990 1046.2 0.1123 1.0050 931.4

0.0903 0.9929 712.6 0.0991 0.9711 1120.5

0.0922 0.9915 925.3 0.1139 0.9835 978.6

0.0891 1.0008 778.7 0.0747 0.9936 875.4

0.1003 1.0279 846.3 0.1299 0.9803 966.4

0.1168 0.9730 880.6 0.0941 1.0007 941.5

0.1036 1.0082 1006.1 0.1225 1.0104 910.0

0.1092 0.9867 899.6

Fig. 19. Comparison of the overall error from MCS and RSM (mm).

A b value of -10.41 is obtained and this suggests that the most probable point is in the failure region and the forging design reliability is close to zero, i.e. P(g(x5) >0) « 0. Therefore with the current forging design, none of the products will meet the specified tolerance. Therefore, the two-step stochastic optimisation is used to improve the net-shape quality of the forged aerofoil components.

3.4. Stochastic forging optimisation

3.4.1. Minimisation of systematic error (mean)

Using the linear RSM approximation, 106 sets of random values of (F, T,D) are generated using MATLAB to obtain the thickness means ix. in Eq. (4b) at all check points as given in Table 3 with the 5th and 9th check point results shown in bold font. After three iterations using the direct die shape compensation method, the thickness means (6th column) are close to the target thickness d0 (2nd column) as given in the table. The last row shows the results of the objective function the square root of the average error of the thickness mean after three iterations. The thickness results and their variations during of die shape compensation iterations are shown in Fig. 11. The histograms and fitted PDF at the 5th check point are shown in Fig. 12(a)-(d). Using the MPP analysis, the reliability of the optimised die design at the 5th check point is P(g(x5) > 0) = 0.9998, which satisfies the 3r quality level. However, with a b value of 1.37, the reliability at the 9th check point is P(g(x9) > 0) = 0.9147. It indicates that 8.5% of forged aerofoil sections will be out of the tolerance specification at the 9th check point and further optimisation is needed so as to reach 3r quality level at all check points.

3.4.2. Reduction of random variation (variance)

The check point with maximum thickness error, i.e. the 9th check point as shown in Table 3 and Fig. 11, is chosen for variance reduction and in doing so variance reduction is also expected over the other check points. At the 9th check point, its mean thickness error is iAdg = 2.0 x 10-5 m and its standard deviation is aAs9 = 2.16 x 10-5 m. Therefore (lAdg + 3aAS^ = 8.5 x 10-5 m is larger than its specified tolerance, i.e. A = 5.0 x 10-5 m. To satisfy

the thickness error tolerance (lAdg + 3aAS^ = 5.0 x 10 5 m, the

standard deviation has to be within rAdg < 1.0 x 10-5 m. Hence the required standard deviation of the stroke length, rD1, can be calculated to be 3rD1 = 2.0 x 10-5 m using the relationship of the variance of the linear approximation. Using the control variable method, the scalar parameter 1* in Eq. (20) is calculated to be 1* = 0.7445. This requires the initial setting range of the stroke length D to be changed from D = (iD ± 3rD) = (1.05 x 10-3 ± 5.0 x 10-5^ m to (1.05 x 10-3 ± 2.0 x 10-5) m.

Another iteration of optimisation computation is carried out by using RSM with a new range of stroke length D. Fig. 13 shows that the thickness at all eleven check points have satisfied

(lAd, ± 3rAdj^ between upper and lower limits, i.e. LSL <

(lAd, ± 3aASj^ < USL. The histograms and fitted PDF after variance reduction of the 5 th check point is shown in Fig. 12(e). The changes of the standard deviation using the initial die and optimised die shapes without and with variance reduction are shown in Fig. 14. Using the MPP analysis, the reliability of the design is calculated to be P(g(xi) > 0) = 0.9999 and therefore satisfies the specified 3r quality level.

4. Case study 2: forward extrusion

4.1. FE modelling

In this case study, an axis-symmetric forward extrusion of part at room temperature is evaluated as shown in Fig. 15. To quantify the dimensional changes in the whole extrusion cycle, the FE analysis is carried out in four steps: extrusion (0.2 s), punch unloading (0.01 s), die removal (0.01 s) and cooling (1200 s). The workpiece material is Al199.5. The material properties of both the workpiece and the extrusion die are given in Table 1. The flow stresses as a function of strain, strain rate and temperature are shown in Fig. 16. During extrusion, the die elastic deformation and thermal induced deformation of the extruded component are major factors that determine the geometrical accuracy. The uncertainties of material flow stress (S), friction coefficient (F) and heat transfer

-0.00002 0 0.00002 0.00004 O.OOOOS 0.0000S 0.0001

Error (m)

Fig. 20. Geometric errors in die shape compensation iterations.

coefficients (H) are considered. All three random variables are assumed normally distributed and independent of each other. The mean values of friction coefficient and heat transfer coefficient are iF = 0.1 and iH = 9 kW/m °C and their standard deviations are oF =0.015 and oH = 0.91 kW/m °C, respectively. As a function of strain, strain rate and temperature, the material flow stresses can't be defined as a single value. Therefore a scaling factor w that obeys normal distribution N(iw, ow) is used to represent the variation of flow stresses (S), S(e, e, T) ■ (iw + ow x N(0,1)), where iw and ow are the mean and standard deviation of w, i.e. iw = 1, ow = 0.03. 8000 x 3 samples for the three random values are generated using MATLAB for Monte Carlo simulation (MCS).

4.2. Material deformation

Fig. 17 shows the von Mises stresses and equivalent plastic strain distributions at the mean values of the random variables by using the initial die design. Fig. 18 compares the force history results between the initial and optimised die shapes at the mean values of the random variables. It is observed that similar material deformation, stress and strain distributions are obtained by the variation of the random variables and in comparison with the cases between initial and optimised die shapes. This is because the changes of process conditions and die shapes for deterministic optimisation are relatively small in terms of stress, strain and tem-

Fig. 21. Comparison of overall error results in two-step optimisation (mm).

perature distributions. However, such changes can be significant to the final shape and dimensions of extruded parts.

4.3. Probabilistic characterisation

4.3.1. Monte Carlo simulation (MCS)

After 8000 deterministic FE analyses, the probabilistic characteristics of the geometrical error of the extruded component are obtained as shown in Fig. 19(a). It is noted that the results of the shape error are also normally distributed. Its statistical mean and standard deviation are given in the figure.

4.3.2. Response surface method (RSM)

Using the Latin Hypercube Sampling (LHS) method, 27 (m = 33) sets of parameters are generated for RSM based computation as given in Table 4. After running 27 FE analyses deterministically, the geometric errors of the extruded part are obtained for regression analysis. The linear approximation of the overall geometric error Ad of the extruded part is obtained as

Ad = 0.0159 + 0.371 F + 0.0057w - 5.285 x 10-6H (26)

The overall error mean and standard deviation are calculated to be 1Ad = 0.0587 mm and rA~ = 0.0058 mm. The quadratic approximation of the geometric error is represented as

Ad = -5.47 + 4.481 F + 9.9761w + 0.0006H - 5.045w • F

- 0.0007 w • H + 0.0006F • H - 4.41 w2 + 2.32F2

- 5.7 x 10-10H2 (27)

The histograms, the fitted PDF with mean and standard deviation are obtained by using the same 8000 realisations of the three random variables as shown in Fig. 19(b) and (c). Similar to the aerofoil forging case, the response shape errors are also normally distributed. Concerning the correlation between the geometric error and the random variables, the correlation coefficient are obtained to be pAiw = 0.0138, pAiF = 0.9824 and pAiH = -0.1007, respectively. This shows that the uncertainty of the friction coefficient is a major factor to affect the geometrical accuracy of the extruded part while the effect of the variation of material properties is negligible.

4.3.3. Most probable point (MPP) analysis

With a specified overall tolerance of A = 0.03 mm, the limit state function may be constructed based on the linear RSM approximation as follows

g(Ad) = 0.0092 - 0.371 F - 0.0057w + 5.285 x 10-6H (28)

Using the MPP analysis iterative procedure, the b value is found to be -2.434 suggesting that the MPP is at the failure region. Thus with the current extrusion design, the shape error of almost all products will be out of the specified tolerance. To improve the reliability of the extrusion design, the two-step stochastic optimisation is employed to reduce its overall geometric error and associated random variations.

4.4. Stochastic extrusion optimisation

4.4.1. Minimisation of systematic error (mean)

Similar to the blade forging case, the direct die shape compensation method is used to reduce the overall geometric error of the extruded part from original iAd = 0.0636 mm to 0.0175 mm after three die shape optimisation iterations. The distributed shape errors of the extruded part in each compensation iteration are presented in Fig. 20. The histograms and fitted PDF of the overall error are shown in Fig. 21(a)-(d). With the optimised die geome-

try, the reliability of the geometric accuracy of the extrusion process is obtained to be > 0^ = 0.992 (b = 2.43) using

the MPP analysis. This presents a significant improvement of the reliability of the extrusion process as compared to that of the initial design without die shape optimisation. However, it is still not sufficient to meet the 3r quality level due to geometric error variations.

4.4.2. Reduction of random variation (variance)

After die shape compensation iterations, the mean error is 0.0175 mm and its standard deviation is 0.0056 mm. So (lAd + 3rAd) = 0.0343 mm is larger than the specified tolerance Ad = 0.03 mm and variance reduction is necessary for improved quality of the extruded part. As shown in Section 4.3.2, the overall error of the extruded part is strongly correlated with the friction coefficient. Therefore friction coefficient is used as the control variable for variance reduction. To satisfy the tolerance, i.e. (lAd + 3rAd) = 0.030 mm, i.e. aAd 6 0.0041 mm. Using the control variable method, the scalar parameter 1* in Eq. (20) is calculated to be 1* = 0.4. This requires the initial setting range of the friction F to be changed from F = iiF ± 3rF) = (0.1 ± 0.045) to (1.0 ± 0.03) .By using the MPP analysis, a b value of 3.3 is obtained and this corresponds to a reliability of P(g('Ad) > 0^ = 0.9995 of the optimised extrusion process, which satisfies the 3r quality level. The histogram and fitted PDF of the overall error after the variance reduction is shown in Fig. 21(e) as compared to the results with and without die shape optimisation.

5. Conclusions

This study presents three stochastically based approaches, i.e. Monte Carlo simulation (MCS), response surface method (RSM) and most probable point (MPP) analysis, for quantifying the probabilistic characteristics of dimensional errors in net-shape metal forming processes. A two-step optimisation approach is developed first to minimise the systematic geometric errors using the direct die shape compensation method and then to achieve variance reduction through a control variable method.

Two industrial based metal forming cases, i.e. forging of 2D aerofoil shape and forward extrusion are presented. Although MCS can be easily implemented to quantify probabilistic characteristics of metal forming processes, it can be computationally demanding due to the large number of simulations required especially for 3D metal forming simulations. Alternatively, RSM can be employed using either standard low-central-high approach or the Latin Hypercube Sampling (LHS) with much less computing time. The results obtained from this research show a good degree of convergence. The reliability of the design or the probability of failure can be effectively estimated using the MPP analysis based on RSM approximation. In the two-step stochastic optimisation, the results show that the mean value of shape error on target is achieved using the die shape compensation, while the variance reduction can be obtained using the control variable method by reducing the variability of the most influencing factor. In the 2D aerofoil forging case, the stroke length is most sensitive to the aerofoil thickness, whilst in the cold extrusion case, the friction is the most sensitive to the shape error. Therefore this suggests that different control variables may be chosen for effective variance reduction of the shape and dimensional errors for different metal forming processes. This study also demonstrates that the developed stochastic optimisation approach can be easily employed to complex 3D metal forming processes with the consideration of uncertainties.

Acknowledgement

The authors would like to thank the UK Engineering and Physical Science Research Council (EPSRC) for funding of the present research (EP/C004140/2). The first author acknowledges the Visiting Fellowship Fund of the State Key Laboratory of Mechanical Transmission (2010).

References

[1] Wright DC, Smith DJ. Forging of blades for gas turbines. Mater Sci Technol 1996;2:742-7.

[2] Steffens K, Wilhelm H. Next engine generation: materials, surface technology, manufacturing processes. MTU aero engine report; 2003, p. 1-17.

[3] Voigtlander O, Gunther G. Isothermal precision forging of aero-engine compressor blade made in titanium alloy. Metallurgia 1983;50:322-6.

[4] Hartley P, Pillinger I. Numerical simulation of the forging process. Comput Methods Appl Mech Eng 2006;195:6676-90.

[5] Fourment L, Chenot J-L. Optimal design for non-steady-state metal forming processes: I. Shape optimization method. Int J Numer Methods Eng 1996;39:33-50.

[6] Chung JS, Hwang SM. Process optimal design in forging by genetic algorithm. J Manuf Sci Eng 2002;124:397-408.

[7] Antonio CAC, Dourado NM. Metal-forming process optimisation by inverse evolutionary search. J Mater Process Technol 2002;121:403-13.

[8] Wei DL, Cui ZS, Chen J. Optimization and tolerance prediction of sheet metal forming process using response surface model. Comput Mater Sci 2008;42:228-33.

[9] Breitkopf P, Naceur H, Rassineux A, Villon P. Moving least squares response surface approximation: formulation and metal forming applications. Comput Struct 2005;83:1411-28.

[10] Zhao G, Wright E, Grandhi RV. Preform die shape design in metal forming using an optimization method. Int J Numer Methods Eng 1997;40:1213-30.

[11] Zhao G, Huff R, Hutter A, Grandhi RV. Sensitivity analysis based preform die shape design using finite element method. J Mater Eng Perform 1997;6:303-10.

[12] Srikanth A, Zabaras N. Shape optimization and preform design in metal forming processes. Comput Methods Appl Mech Eng 2000;190:1859-901.

[13] Ou H, Lan J, Armstrong CG, Price MA, Walloe SJ, Ward MJ. Reduction in post forging errors for aerofoil forging using finite element simulation and optimisation. Model Simul Mater Sci Eng 2006;14:179-93.

[14] Lu B, Ou H, Armstrong CG, Rennie A. 3D die shape optimisation for net-shape forging of aerofoil blades. Mater Des 2009;30:2490-500.

[15] Karafillis AP, Boyce MC. Tooling and binder design for sheet metal forming processes compensating springback error. Int J Mach Tools Manuf 1996;34:503-26.

[16] Cheng HS, Cao J, Xia ZC. An accelerated springback compensation method. Int J Mech Sci 2007;49:267-79.

[17] Haldar A, Mahadevan S. Reliability assessment using stochastic finite element analysis. New York: John Wiley & Sons Inc.; 2000.

[18] Kleiber M, Hien TD. The stochastic finite element method: basic perturbation technique and computer implementation. New York: John Wiley & Sons Inc.; 1992.

[19] Choi S, Grandhi RV, Canfield RA. Reliability-based structural design. London: Springer-Verlag Ltd.; 2007.

[20] Schueller GI. Computational stochastic mechanics - recent advances. Comput Struct 2001;79:2225-34.

[21] Stefanou G. The stochastic finite element method: past, present and future. Comput Methods Appl Mech Eng 2009;198:1031-51.

[22] Schueller GI, Jensen HA. Computational methods in optimization considering uncertainties - an overview. Comput Methods Appl Mech Eng 2008;198:2-13.

[23] Bulur BK, Grandhi RV. Geometric deviations in forging and cooling operations due to process uncertainties. J Mater Process Technol 2004;152:204-14.

[24] Repalle J, Grandhi RV. Optimum design of forging process parameters and preform shape under uncertainties. In: Proceedings of NUMIFORM conference, Columbus, Ohio; 2004. p. 2032-7.

[25] Repalle J, Grandhi RV. Reliability-based preform shape design in forging. Commun Numer Methods Eng 2005;21:607-17.

[26] Gantar G, Kuzman K. Optimization of stamping processes aiming at maximal process stability. J Mater Process Technol 2005;167:237-43.

[27] Krusic V, Masera S, Pristovsek A, Rodic T. Adjustment of stochastic responses of typical cold forging systems. J Mater Process Technol 2009;209: 4983-93.

[28] Acharjee S, Zabaras N. A non-intrusive stochastic Galerkin approach for modelling uncertainty propagation in deformation processes. Comput Struct 2007;85:244-54.

[29] Jansson T, Nilsson L, Moshfegh R. Reliability analysis of a sheet metal forming process using Monte Carlo analysis and metamodels. J Mater Process Technol 2008;202:225-68.

[30] Kleiber M, Rojek J, Stocki R. Reliability assessment for sheet metal forming operations. Comput Methods Appl Mech Eng 2002;191:4511-32.

[31] Wei DL, Cui ZS, Chen J. Uncertainty quantification using polynomial chaos expansion with points of monomial cubature rules. Comput Struct 2008;86:2102-8.

[32] Rubinstein RY, Kroese DP. Simulation and the Monte Carlo method. Hoboken, New Jersey: John Wiley & Sons, Inc.; 2008.

[33] Kleiber M, Knabel J, Rojek J. Response surface method for probabilistic assessment metal forming failures. Sheet metal forming operations. Int J Numer Methods Eng 2004;60:51-67.

[34] Wyss GD, Jorgensen KH. A user's guide to LHS: Sandia's Latin Hypercube sampling software (SAND98-0210), Sandia National Lab, New Mexico; 1998.

[35] Park JS. Optimal Latin-Hypercube designs for computer experiments. J Stat Plan Infer 1994;39:95-111.

[36] Hasofer AM, Lind NC. Exact and invariant second-moment code format. J Eng Mech Div, ASCE 1974;100:111-21.