当前位置:首页 >> 工学 >>

3D FEM simulation of stationary metal forming processes with applications to slitting and rolling


Journal of Materials Processing Technology 148 (2004) 328–341

3D FEM simulation of stationary metal forming processes with applications to slitting and rolling
H.H. Wisselink? , J. Huétink
Netherlands Institute for Metals Research, University of Twente, P.O. Box 217, 7500 AE, The Netherlands Received 2 September 2002; accepted 12 February 2004

Abstract A computational method is presented for the 3D FEM simulation of stationary metal forming processes. The described method has been applied to model the slitting and shape rolling processes. Some results of these simulations are given, which show that the ALE formulation used is suitable for the simulation of such processes. In the ALE formulation the grid displacement is independent of the material displacement. Therefore, the displacement of the nodes has to be de?ned such that the movement of the free surfaces can be followed (mesh management) and the values of the state variables have to be calculated at the new nodal positions (transfer of state variables). Both topics are addressed in this paper. Also the possibility of modelling stationary crack propagation is treated. ? 2004 Elsevier B.V. All rights reserved.
Keywords: ALE method; Shape rolling; Slitting

1. Introduction Some different formulations, as updated Lagrangian (UL), Eulerian or arbitrary Lagrangian Eulerian (ALE) formulations, are used for FEM simulations of metal forming processes. As each of these formulations has its bene?ts and drawbacks, for each application the most suitable formulation can be chosen. In forming processes a UL formulation is often used. In such a formulation the mesh is moving with the material. With a UL formulation free surfaces are followed automatically and history dependent behaviour can easily be taken into account. A drawback of the UL formulation is that the shape of the elements can become distorted in the case of large deformations, leading to less accurate results or even to a premature end of the calculation. The UL formulation can be extended with a remeshing procedure, in which the old grid is replaced by a completely new grid, and the information on the old grid is transferred to this new grid. Remeshing is a ?exible method, as during the whole simulation the mesh can be adapted to best ?t the current situation. It can be used to avoid too much grid distortion and for mesh re?nement in critical areas. Impor? Corresponding author. E-mail addresses: h.h.wisselink@utwente.nl (H.H. Wisselink), j.huetink@utwente.nl (J. Hu? etink). URLs: http://www.nimr.nl, http://www.tm.ctw.utwente.nl.

tant issues in the use of remeshing are the automatic generation of the new ?nite element mesh and the accuracy of the transfer of information between the old and new meshes. Remeshing is widely used for the simulation of forming processes in two dimensions. Although there are some recent examples of the use of remeshing for 3D forming processes with hexahedral elements [1,2], the use of remeshing for three-dimensional problems is not yet well developed, due to a lack of robust mesh-generators for hexahedral elements. Using a UL formulation for a stationary process is less attractive if a large number of remeshings is required, due to localised deformations. In a Eulerian formulation the material ?ows through a frame of reference which is ?xed in space. Therefore, problems with grid distortion due to large deformations do not exist in such a formulation, but in general material boundaries are not equal to the element edges, so special procedures are required to follow free surfaces or boundaries between different materials, including contact [3,4]. To handle history dependent material properties convection must be taken into account. The arbitrary Lagrangian Eulerian formulation is a combination of the Lagrangian and the Eulerian descriptions. In such a formulation the mesh displacement is not necessarily equal to the material displacement (Lagrangian formulation) nor equal to zero (Eulerian formulation), but can be chosen independently from the material displacement. This formulation was ?rst applied in the ?nite element method to cal-

0924-0136/$ – see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jmatprotec.2004.02.036

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

329

culate ?uid-structure interactions [5,6] and later to 2D [7,8] and 3D [9–11] metal forming processes. The ALE formulation can be used to solve the problem with free surfaces in a Eulerian formulation or to avoid grid distortion in a Lagrangian formulation. As in the Eulerian formulation, convection of history dependent variables must be taken into account. A characteristic of the ALE formulation is that the topology of the mesh (number of elements and their connectivity) is constant during the entire simulation. This means that the initial domain as well as the ?nal domain must be described with the same topology, which limits the applicability of the ALE formulation. The transfer of information between the old mesh and the new mesh is in general more accurate for an ALE formulation, where convection techniques can be used, than for a remeshing step in the UL formulation. An ALE formulation is chosen as the most appropriate formulation for the calculation of the steady state of stationary metal forming processes, like rolling and slitting. Other formulations are more suitable for modelling the start and end of these processes. The steady state of a stationary process can be calculated by continuing a transient calculation until a steady state is reached. Stationary processes can also be seen as ?ow problems, when in- and out?ow boundaries are added. When the calculation is started with a volume suf?ciently close to the steady state volume, it is possible with the ALE method to avoid grid distortion and to describe the free surfaces correctly. This means that an estimation of the steady state geometry has to be made for the initial mesh (the exact shape of the free surfaces is not known beforehand). The speci?c ALE formulation used for the simulations is presented in Section 2. The transfer of state variables is addressed in Section 3 and the mesh management procedures are elaborated in Section 6. Next some attention is paid to the modelling of the contact between the material and the tools (Section 4), the possibility to incorporate ductile fracture (Section 5) and the solving of the FEM equations (Section 7). Finally in Section 8 some examples are shown of the application of the proposed simulation method to the slitting and rolling process. All simulations have been carried out with the ?nite element code DiekA, developed at the University of Twente.

ential point corresponds to one material point and one spatial point, but not always to the same points. The Lagrangian and Eulerian coordinates are special cases of the referential coordinates. The material velocity vm and the referential velocity vg are de?ned as ?x ˙ =x (1) vm = ?t X vg = ?x ?t
χ

=x

(2) (3)

vc = vm ? vg

The time derivatives, denoted by dot and star, are, respectively, for a constant material coordinate X and a constant referential coordinate χ. The convective velocity vc is the difference between the material velocity and the grid velocity. The material time derivative of an arbitrary quantity ξ depends on the coordinates in which it is written: d ? ξ(X, t) = ξ(X, t) ?t dt
X

˙ =ξ

(4) (5)

? d ξ(x, t) = ξ(x, t) + vm · ? ξ(x, t) dt ?t x d ? ξ(χ, t) = ξ(χ, t) dt ?t
χ

+ vc · ? ξ(χ, t) = ξ + vc · ? ξ(χ, t) (6)

A convective term appears when ξ is not written as a function of the material coordinates. The spatial time derivative vanishes for stationary processes. 2.2. Conservation laws The conservation of momentum given in the current con?guration is σ·? =0 σ·n =t in V on St (7) (8)

2. ALE formulation 2.1. Kinematics The motion and deformation of a continuum can be described in different ways [12]. The state variables are, respectively, written in the material coordinates X or the current coordinates x for the Lagrangian or the Eulerian description. The referential description is used for the ALE method in which the state variables are written in referential coordinates χ. One-to-one mappings between the different domains exist and are functions of time. Therefore, a refer-

where n is the outward normal on the boundary, t the prescribed loads. Body and inertia forces are neglected. In the ?nite element method the equilibrium equations are only weakly enforced. The weak form with the weight function w reads
V

(w? ) : σ dV =

S

w · t d St

(9)

The elastic–plastic material model is given in a rate form with the Jaumann derivative of the stress σ ? , the elasticity and yield tensors E, Y and the rate-of-deformation tensor D: σ = (E ? Y) : D
?

(10)

Therefore the rate of the weak form will be used. The complete rate of the weak form of the ALE equations, which

330

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

{ um }. This expression is added to the stiffness matrix, giving twice as many degrees of freedom and thus larger systems to solve. If the convective displacements are expressed as an explicit function of the material displacements, they can be eliminated from the stiffness matrix. In practice this is limited to relatively simple functions [15]. Assembling and solving these equations is called the ALE predictor (Fig. 1). In the corrector part the state variables are integrated in time, using the results of the predictor { u1 m} and { u1 g }: ξi = ξ0 + = ξ0 +
t+ t t t t+ t

ξ dt ˙ dt ? ξ
t+ t t

(vm ? vg ) · ? ξ dt

(14)

Fig. 1. ALE methods.

contains as well the material velocity as the grid velocity, can be derived from Eq. (9) as shown in [10,13,14]. ˙ ? (vm ? vg ) · ? σ + σ tr (Lg )] dV w? : [?Lg · σ + σ =
S

V

˙ dS ? w·t

This integral consists of two parts, a material part and a convective part. The calculation of the convective part is treated in more detail in Section 3. From the stresses at the integration points the nodal reaction forces can be calculated. These internal forces should be in equilibrium with the external forces. A number of iterations can be required to reach a certain convergence threshold. A disadvantage of this method is that the convection must be taken into account at each iteration, which is time-consuming. 2.3.2. Decoupled ALE In the decoupled ALE method, also called operator split ALE, each time step is subdivided into two parts which are solved separately [17–19]. First a normal updated Lagrangian step is carried out, discretising and solving Eq. (11) for vg is equal to vm . Now only the material displacement is calculated in the predictor part: [Km ]{ um } = { Fext } (15)

S

w · [(vm ? vg ) · ? σ ]n dS (11)

+

S

w · t J S dS

with the gradient of the grid velocity Lg and the grid derivative of the Jacobian J S . 2.3. ALE methods Some different approaches to solve the ALE equations (11) are shown in Fig. 1. The coupled and the decoupled ALE formulations are discussed here together with a semi-coupled ALE formulation, which has similarities with both of the other approaches. This semi-coupled ALE formulation will be used for the simulations in this article. 2.3.1. Coupled ALE A coupled ALE formulation is used for simulations of forming processes by Liu et al. [13,15,16]. In this method Eq. (11) is discretised in space and time and is therefore the best possible solution of the ALE problem. This discretisation leads to the next set of equations (with stiffness matrices [K]) in which the material and convective terms are solved simultaneously: [Km ]{ um } + [Kg ]{ ug } = { Fext } [A]{ um } + [B]{ ug } = {C} (12) (13)

and in the corrector part only the material increment of the state variables is integrated in time: ξL = ξ0 +
t+ t t

˙ dt ξ

(16)

The predictor and corrector are repeated until equilibrium is satis?ed (Fig. 1). Now the material displacement and the state variables after the Lagrangian ξ L step are known, the grid displacement is determined using the material displacement and the new values of the state variables. Finally, the state variables at the integration points of the new grid ξ n+1 at t n+1 have to be calculated from the values ξ L in the integration points of the old grid, also at t n+1 . This is an interpolation problem as no time change is involved, but it can also be written as a convection problem:
n+1 n+1 ξ n+1 (xg ) = ξ L (xm ?

uc )

(17)

To solve this expression, extra equations (Eq. (13)) are needed which describe the relation between { ug } and

The advantage of the decoupled ALE method compared to the coupled ALE method is the ?exibility in de?ning the new grid, as the latter is determined after the Lagrangian step

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

331

is completed. The decoupled method is also more ef?cient since the grid displacement and the convective increment are calculated only once per time step. A drawback of the decoupled method is that equilibrium is not checked at the end of the step. Equilibrium can be violated if the convection is not performed accurately. This leads to unbalance forces, which are carried as a load correction onto the next step. Large unbalance forces can lead to possible divergence. 2.3.3. Semi-coupled ALE The semi-coupled ALE method (Fig. 1) can be seen as a coupled formulation, in which the ALE predictor is replaced by the UL predictor [7]. As the convective terms are neglected a slower rate of convergence is expected compared to the coupled method. An Eulerian stiffness matrix (setting { ug } = {0} in the predictor part) may yield better results in the case of ?ow problems. With the UL predictor the material displacement increment is calculated as in the decoupled method. The meshing step is equal to the decoupled method, but has to be performed at every iteration. This is less ef?cient, but the ?exibility is preserved. In the corrector part the convective increment is calculated using the values of the state variables at the beginning of the step,
n+1 n ξ C (x g ) = ξ 0 (xm ?

the only 3D convection scheme currently implemented in the ?nite element code DiekA. The Molenkamp test has been used to investigate the accuracy of some different 2D convection schemes [18,20–22]. From these results, it can be concluded that the WLGS scheme is not the most accurate scheme in 2D. Hence, the development of 3D versions of schemes other than WLGS would be an interesting topic for further research. 3.1. Weighed local and global smoothing The WLGS scheme was ?rst described in [7,23]. The main characteristics of this scheme are summarised in this section. Some changes were made to the original scheme to avoid cross-wind diffusion. It was shown that this method has similarities with the Lax–Wendroff schemes with added arti?cial diffusion used in ?nite difference [21]. The different steps in this scheme are elaborated for trilinear hexahedrals with shape functions N k , with the local coordinates r, z, points are located at √ The integration √h ∈ [?1, 1]. √ (±(1/3) 3, ±(1/3) 3, ±(1/3) 3). The Courant number is de?ned as the difference in local coordinates ( r, z, h) between the old and new integration points, divided by the element size (2 × 2 × 2 in local coordinates). Three Courant numbers Cr , Cz and Ch are calculated for the different directions at every integration point. The integration point values are then averaged to one value for the whole element Cr = 1 8
8 ip=1

uc )

(18)

followed by the integration of the material increment: ξ n+1 = ξ C +
t+ t t

˙ dt ξ

(19)

Equilibrium is checked at the end of the time step and the iteration loop is repeated until the convergence criterion is ful?lled. This leads to a more robust algorithm, compared to the decoupled method, which is especially important in contact problems. Therefore, this formulation has been chosen for the simulations in this article as it combines robustness with ?exibility.

| rip | 2

(21)

2 + C2 + C2 . The total Courant number is then C = Cr z h The WLGS scheme uses a smoothed continuous ?eld as an intermediate step. The integration point values are extrapolated from the integration points to the nodes using the shape functions: 8 k N k (rip βr , zip βz , hip βh )ξip

3. Transfer of state variables In case of an elastic–plastic material model the history dependent state variables ξ (e.g. strain, stress or damage), which are known at the integration points, have to be calculated at the new position of the integration points. This convective increment has to be calculated at each step or iteration, depending on the ALE method used (Eqs. (14), (17) and (18)). It can be written as a convection or as an interpolation problem. The convective velocity vc is assumed to be constant during a time step: ξg =
t+ t t

ξnode =
k=1

(22)

?vc · ? ξ dt ≈ ? uc · ? ξ

(20)

The convective increment is calculated with the weighed local and global smoothing (WLGS) scheme [7], which is

√ The local smoothing parameter β ∈ [0, 3] is introduced to in?uence the extrapolation. β can be a ?xed number or a function of the Courant numbers in the different local direc√ tions β = f(Cr , Cz , Ch ). A value of β = 3 means normal extrapolation. Least-squares smoothing [24] implies for linear elements that the mean value of the integration points is used as a nodal contribution (β = 0). The nodal values are determined by averaging the contributions of the surrounding elements. A continuous ?eld is obtained by interpolation of the average nodal values using the shape functions. This ?eld is much smoother for β = 0 than for β = √ 3. The values of the state variables in the new integration points can be calculated using the constructed continuous

332

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

?eld:
8 new = ξip k=1 ? k N k ((r, z, h)new )ξnode

(23)

α = min(1.64C1.2 , 1). The amount of smoothing increases with an increasing Courant number. Furthermore in case no convection takes place, the integration point values are not changed (C = 0 ? α = 0). 3.2. Test problem The in?uence of the global and local smoothing parameters α and β is illustrated with a test problem. An initial distribution of a state variable ξ is given to a block, with dimensions 1 × 2 × 5. The average nodal values of the initial distribution depend on the value of β, as shown in Fig. 2(a) and (b). Now this initial distribution will be convected with the material. A material displacement of 0.1 per step is prescribed in the direction of the longest edge (arrow in Fig. 2) and the grid displacement is kept zero. Therefore the material ?ows through the block and the initial distribution, which moves with the material, also. The exact solution is a translation of the initial distribution. Due to numerical errors the distribution after some steps will deviate from the exact solution. In Fig. 2 the results after 25 steps are depicted for different values of α and β. In Fig. 2(d) no local or global smoothing is applied, which leads to oscillations. Adding global smoothing without local smoothing will not improve the results (Fig. 2(f)). Also, local smoothing without global smoothing (Fig. 2(c)) gives unstable behaviour. Only a combination of local smoothing

Another possibility is to use the continuous ?eld to calculate the convective increment and add this to the old integration point value:
8 new old ξip = ξip + k=1 8 ? k N k ((r, z, h)new )ξnode

?
k=1

? k N k ((r, z, h)old )ξnode

(24)

The ?rst method (Eq. (23)) is very diffusive and the second method (Eq. (24)) shows spurious oscillations [7]. Therefore, a weighed combination of these two methods was proposed, which was called global smoothing [7]:
8 new = ξip k=1 8 old + (1 ? α) ξip ? k=1 ? k N k ((r, z, h)old )ξnode ? k N k ((r, z, h)new )ξnode

(25)

The weighting factor, α ∈ [0, 1], is a function of the Courant number. According to the results of some numerical experiments a good choice for the weight factor proved to be

Fig. 2. Test problem for the convection scheme. Results for different values of α and β.

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

333

and global smoothing leads to stable, but also diffusive results (Fig. 2(e)). Smaller steps and elements will improve the results [7]. 3.3. In?ow conditions In the simulation of stationary processes a boundary exists, where material ?ows into the domain (vc · n < 0). At these boundaries, in?ow conditions are needed for the state variables and displacements. The calculated average nodal values for the continuous ?eld are overruled by these prede?ned values. Otherwise the convection is a pure extrapolation, which gives small errors at the in?ow boundary. These errors will add up during the simulation and the steady state will never be reached. 3.4. Cross-wind diffusion Another problem during convection is cross-wind diffusion: gradients normal to the streamlines level out in regions where only convection takes place. This is illustrated again with the previous test problem, but now an in?ow condition for the state variable is used. A value of one is assigned to the state variable at the nodes of the element in the left upper corner and zero to the other nodes in the in?ow surface. The initial distribution is given in Fig. 3(a). The distribution after 75 convection steps is illustrated in the other three pictures of Fig. 3. From Fig. 3(b) it can be concluded that √ β = 0 leads to much cross-wind diffusion. Choosing β = 3 solves the cross-wind diffusion problem, but also gives unstable results (Fig. 3(c)). As was seen in the previous test problem local smoothing is needed for a stable convection scheme. Therefore β is made dependent on the ?ow direction, so that the local smoothing is only applied in the ?ow direction, but not perpendicular to it (Eq. (26)). The factor GLF, added to make it possible to keep a little

smoothing perpendicular to the ?ow direction, is generally chosen between 0.9 and 1: βr = (1 ? Cr ) 1 ? Cr C · |GLF| · √ 3, βz = . . . (26) The result of this orthotropic local smoothing is shown in Fig. 3(d). No cross-wind diffusion or instabilities can be seen in this ?gure. Hence it can be concluded that with orthotropic smoothing in grid aligned ?ows, cross-wind diffusion is effectively suppressed and enough smoothing is applied in the ?ow direction to keep the solution stable. 4. Contact modelling Contact elements are used to describe the contact between the sheet and the tools. These elements are based on a penalty formulation and a Coulomb friction law is used [25]. The tools are modelled as rigid bodies. The surfaces of these tools can be described by an analytical function or approximated with a ?nite element mesh. One face of the hexahedral contact element is connected to the bulk material, the opposite face is connected to the rigid tool. The position of the nodes connected to the rigid tools is determined by a projection algorithm. After the projection the penetration gap and penalty can be calculated. However, for some applications the deformation of the tools is an important issue. This means that the deformation of the tools has to be taken into account. In that case the relocation of the nodes on the surfaces of both contacting bodies cannot be done independently for each body, as also the contact elements in between the contacting bodies have to remain well shaped. It was shown that the deformation of the tools could be calculated using these contact elements for stationary processes [26,27].

Fig. 3. Results for three different values of β, α = 1.64C1.2 .

334

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

5. Fracture modelling By itself the ALE method cannot be used to describe a complete fracture process from initiation of a crack till complete separation of the parts, since this requires a change of the mesh topology. However, the ALE method can be used to model crack propagation. The growth of an initially present crack can be followed to a certain extent without changing the mesh topology. The nodes along the crack front should have a critical value for some fracture criterion. If this is not the case, then the nodes on the crack front have to be relocated, such that these nodes ful?l the fracture criterion. No automatic crack adaptation procedure has been implemented yet. Therefore the nodes on the crack front in the slitting simulation of Section 8.2 are spatially ?xed. This means that the crack propagates at a constant speed.

? The mesh management must be ef?cient, as it is carried out at each Newton–Raphson iteration. ? It must be stable and robust. ? The mesh management procedure should be able to maintain or produce mesh re?nements in critical areas. The determination of the new mesh is carried out in the following order: (1) If a crack front is present, its new position has to be determined ?rst. (2) Next, all surface nodes are put back in the x-direction, keeping cross-sections in one yz-plane. The new y- and z-coordinates are calculated with a convection scheme, using the coordinates of neighbouring nodes on a grid line in the x-direction. (3) The surface nodes are redistributed on the surface in the yz-plane, to keep or create re?nements at the locations, where the largest deformations are found. (4) When the new positions of all the surface nodes are known, the new positions of the internal nodes can be calculated. They must follow the movements of the surface nodes in order to preserve a good element shape. (5) Finally, a new contact search is performed. The mesh management procedures will now be explained in more detail. The material coordinates after the predictor part, i.e. the coordinates at the beginning of the step plus the material displacement increments, are denoted by xm . The new grid coordinates, determined using the mesh management procedures, are written as xg . The convective displacement is then given by xm ? xg . 6.2.1. Mesh management in the x-direction The geometry is meshed in such a way that the initial mesh on the surfaces is regular, with grid lines that almost coincide with the x-direction. During the simulation this surface mesh is kept regular. The complete mesh is kept ?xed in the x-direction (the main ?ow direction), but the grid follows the free surfaces in the yz-plane. The determination of the new nodal positions of points on a grid line is an interpolation problem, which has the same characteristics as a convection problem. Therefore convection schemes will be used to calculate the new nodal positions. First-order upwind schemes are known to be diffusive. A more accurate description of the boundary can be obtained with higher order interpolation or convection schemes. Hence, a Lax–Wendroff scheme with van Leer limiters is used here (Eq. (27) and Fig. 4), which shows no spurious oscillations [28]:
i i i i?1 1 yint = ym ? C(ym ? ym )? 2 C(1 ? C)[ψ(ri+1/2 ) i+1 i i i?1 ×(ym ? ym ) ? ψ(ri?1/2 )(ym ? ym )]

6. Meshing and mesh management In this section the meshing and the mesh management part of the ALE formulation is treated. These are closely related in an ALE method as the initial and the ?nal geometries are described with the same mesh topology. The strategy used for the relocation of the nodes in simulations of stationary metal forming processes is explained. 6.1. Initial mesh For the simulation of metal forming problems the use of hexahedral elements is to be favoured above other types of elements. A drawback of the use of hexahedrals is that the generation of a hexahedrals-only mesh can be dif?cult, as a robust and ?exible mesh-generator for an arbitrary domain still does not exist. Already in the de?nition of the initial mesh, the expected ?nal geometry must be anticipated in a ALE simulation. Therefore, it is chosen to create a mesh with a constant topology in every cross-section perpendicular to the x-axis (the ?ow direction). All nodes in a cross-section have the same x-coordinate, but it is allowed to use unstructured meshes in a cross-section. This results in regularly meshed free surfaces and “2D” meshes in every cross-section, which simpli?es the mesh management and allows the use of more accurate convection schemes (orthotropic smoothing, Section 3.4). 6.2. Mesh management The mesh management must ful?l the following requirements: ? The domain must be described accurately. No material should be gained or lost during the simulation, except through in- or out?ow boundaries. ? The element shape must be kept regular.

(27)

i is calculated from y , the An intermediate y-coordinate yint m known y-coordinates after the predictor part. The same procedure is applied to calculate the intermediate z-coordinate.

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

335

using the described scheme. The results show that the maximum decreases a little after 200 steps. Therefore, it can be concluded that the scheme is stable but also shows some diffusion. The Lax–Wendroff scheme is derived for an equidistant grid. Therefore, it has been adapted a little to handle re?nements. However, the deviations are small when neighbouring elements do not differ too much in size. The upwind scheme is used (?rst two terms in Eq. (27)) if not enough points are available for the limited Lax–Wendroff scheme. 6.2.2. Redistribution in the yz-plane A second relocation step in the yz-plane is performed when the relocation of the nodes in the x-direction is completed for all surface nodes. This is needed to keep the initial re?nements or to create new re?nements at locations were large gradients are found. After relocation in the x-direction, the 3D-mesh again consists of 2D cross-sections in the yz-plane. Therefore a method, developed by Ponthot [29,30] for 2D simulations, can be applied here. The mesh in the yz-plane is now subdivided into poles (P1), master lines (ML1) connecting the poles and regions (R1) formed by the master lines (Fig. 6). Straight master lines are used for internal master lines or for master lines that are suppressed in one direction. Preferably, the poles are chosen such that they do not need an extra convective displacement. The nodes on the master lines are repositioned along these master lines using a weight function. The simplest form of the weight function contains the element length only, which results in an uniform element distribution along the master line. A factor for the initial re?nement can be included to maintain the re?nement during the simulation. These simple weight functions can be very effective. More sophisticated weight functions can be composed by incorporating geometrical information, such as the angle between segments,

Fig. 4. Convection of nodal coordinates in the x-direction.

The van Leer limiter ψ(r) stabilises the Lax–Wendroff scheme. In case of large gradients (ψ(r) → 0) this scheme degenerates to the upwind scheme. The in?uence of the limiter is small for small gradients, as the original Lax–Wendroff scheme is retained for ψ(r) → 1: ri+1/2 =
i+1 i ym ? ym r + |r | ψ(r) = 1 + |r | i ? yi?1 ym m

,

ri?1/2 =

i?1 ? yi?2 ym m i ? yi?1 ym m

, (28)

The Courant number is determined directly from the material i is displacement in the x-direction. The new x-coordinate xg equal to the initial x-coordinate: C=
i ? xi xg m i ? xi?1 xm m

(29)

The application of this scheme in the x-direction is illustrated with a test contour. The initial sine shaped contour is shown in Fig. 5. This contour is moving in the x-direction with a velocity of 0.1. The nodes on the contour are adapted

Fig. 5. Contour of test problem.

336

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

For some internal nodes it is chosen to relocate them in proportion to the grid displacement of two other nodes, for example the nodes on the surfaces. Using this method it is easier to maintain the initially applied re?nements.

7. Solving the FEM equations The solutions for subsequent time steps in transient calculations that converge to a steady state are very similar. Therefore, the calculations can be speeded up if the predictor part of the ?rst iteration is skipped. Instead, the converged solution of the last step may be used as an initial guess of the material displacement increment. This proved to be unstable when the solution converged in one iteration per step. In this case only the corrector part is carried out, which is very fast, but the Newton–Raphson iterations will diverge in one of the subsequent steps. Therefore, at least two NR-iterations per step are required, which means that the solution will always be adapted. Furthermore, this method stabilises the solution process, except for abrupt changes in the material displacement increment, e.g. when points come into contact with the tools. Also line-searches are used, in the second and next NR-iterations, to improve the stability of the simulation. For large 3D simulations, as shown in this paper, the amount of memory is insuf?cient for the use of direct solvers. Therefore iterative solvers have to be used. The used solvers, GMRES or BiCGSTAB with a SSOR preconditioner, provide acceptable answers using less time and memory. Examples of the use of iterative solvers in forming processes are given in [31].

Fig. 6. Mesh management in yz-plane.

or error indicators such as the derivatives of the state variables. The new positions of the nodes along the master lines are determined such that the weight function is equally distributed along the curve. For the determination of the new positions of nodes on a master line the Lax–Wendroff convection scheme will be used again with the coordinates yint and zint . A simple interpolation on a line through the new position of the poles is used to calculate the new coordinates of nodes on a straight master line. 6.2.3. Internal nodes The internal nodes can be relocated using a smoothing procedure. The new coordinates of an internal node are calculated as the mean value of all coordinates of the neighbouring nodes. If the new coordinates of the neighbour node are already known (e.g. nodes on a master line), these new values are used. This is a simple form of the Laplacian smoothing technique: ? ? ? N N 1 ? k k? xg = xg + xm (30) N ?
k=1 k=N +1

8. Applications The method described in this paper has been applied to two stationary forming processes: shape rolling [32] and slitting [10]. Simulation results are shown and some details of the mesh management procedure are given. 8.1. Shape rolling The principle of the studied shape rolling process is illustrated with Fig. 7. A strip of sheet metal is clamped to

One loop over all internal nodes has proved to be suf?cient to reach an accurate solution. This method can depend on the sequence in which the nodes are treated, but that effect was not apparent in the simulations. This simple algorithm gives good results for convex domains, but problems can arise near concavities, where the elements can be inverted. A proper subdivision of the domain into convex subdomains solves this problem. A drawback of smoothing is the limited in?uence of the mesh on the master lines on the mesh in the interior of the region. In spite of re?nements on the master line, the smoothing algorithm will produce equally sized elements in the interior of a region. Extensions of the smoothing algorithm (Eq. (30)) are possible by taking into account the element shape (angle, size, etc.) or an error indicator in the weighting of the contribution of each node in the averaging function. These methods are more time-consuming and therefore not always attractive for ALE formulations, since the smoothing is performed at each iteration.

Fig. 7. A shape rolling process.

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341 Table 1 Material and process data Slitting Material σ0.2 Flow curve Friction coef?cient Width × thickness sheet Clearances Radii slitters/roll Stainless steel 272 MPa σy = 1228(0.023 + ε)0.4 0.1 10 mm × 1 mm hor = 0.05 mm; vert = 0 mm 150 mm Rolling

337

Aluminium 83 MPa σy = 1 + 170(0.026111 + ε)0.2 0.11 40 mm × 3 mm – 63.5 mm (maximum)

the ?xed lower die. This strip is deformed by the upper die, which rolls over the lower die. The shapes of both tools determine the shape of the product. This process can be modelled as a stationary process, when the material in the dotted box is followed. In this way the lower die translates, the upper die only rotates and in- and out?ow boundaries are present. Dimensions and properties are given in Table 1. The initial and steady state meshes are shown in Figs. 8 and 9. During the simulation, as the material ?ows through the mesh, the mesh is adapted in order to follow the free surfaces. The spread of the strip can be seen clearly. The mesh topology of a cross-section can be seen in Fig. 9. In each cross-section some poles are de?ned at the four edges and in the middle (Fig. 10). The poles are connected by master lines. The master lines ML1a and ML1b (ML2a and ML2b, etc.) are chained into a new master line ML1, as long master lines are preferred. The element size on each master

Fig. 9. Cross-section of the mesh of the rolling simulation.

Fig. 10. Mesh management in yz-plane.

line is controlled by a re?nement factor to keep the smallest elements at the edges, where the reduction is largest. The internal nodes in the dark area of Fig. 10 are repositioned proportional to the surface nodes, while the nodes in the light area are repositioned by the smoothing procedure. The unstructured mesh at both edges is created to be able to describe the barrelling deformation at the edge. Information about the deformation of the material cannot be obtained from the mesh in an ALE simulation. Therefore material points with the same initial z-coordinate before and after rolling are shown in Fig. 11. This deformation pattern

Fig. 8. Initial and steady state mesh of the rolling simulation.

Fig. 11. Material ?ow during rolling.

338

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

Fig. 12. Hydrostatic pressure (×100 MPa) during rolling.

Fig. 15. Initial and steady state mesh of the slitting simulation.

Fig. 13. Equivalent plastic strain during rolling.

8.2. Slitting Slitting is a sheet metal cutting process with circular blades. It is used in industry to split wide coiled sheet into narrower width or for edge trimming of rolled sheet. This process is schematically drawn in Fig. 14. The sheet is ?rst plastically deformed (A–B). Continued deformation leads to a ductile fracture (B–C) and complete separation of the parts (D). Some parameters, e.g. the horizontal and vertical clear-

has been con?rmed by experiments [32] and is explained by the maxima of the hydrostatic pressure (Fig. 12). It can be seen in Fig. 13 that the cross-wind diffusion is minimal, as the equivalent plastic strain in the material does not change after passing the rolls.

Fig. 14. The slitting process.

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341

339

Fig. 16. Mesh in a cross-section of the slitting simulation.

Fig. 18. Deformation of the material during slitting.

ance (Table 1), can be used to in?uence the slitting process in order to produce nicely sheared edges. Due to symmetry only one of the multiple cuts has to be modelled. Slitting is a stationary process and can thus be modelled in the same manner as the rolling process. However, the stationary crack growth makes it more complex. The initial mesh (Fig. 15) contains already a crack front. The cracked surfaces are treated as any other free surface. For the moment the crack front is taken ?xed in space, which means that the crack is growing at a constant pace. The mesh management procedures have to be extended to make it possible to predict the correct location of the crack front. From the initial and steady state meshes in Figs. 15 and 16 can be seen that the mesh management procedures smoothen some sharp edges in the free surface of the initial mesh and maintain the re?nement in the shear zone. The subdivision of the nodes in a cross-section into poles, curves and regions is given in Fig. 17. The poles at the edge of the tools (P2, P3, P6 and P7) are relocated such that they remain at these edges and on the surface of the material. The surface nodes in the dark regions do not get an extra convective displacement in the xy-plane, as the deformations are small in this area. The internal nodes in the dark region move proportional to the lower and upper surfaces to make it possible to follow the rotations of the material (bending). The light area is subdivided into two convex regions R2 and R3 by ML5. The internal nodes in these regions are relocated with the smoothing algorithm.

Fig. 19. Hydrostatic pressure (×1000 MPa) during slitting.

The shearing deformation of the material is shown in Fig. 18, which gives the y-coordinate of a material particle in the undeformed state. The hydrostatic pressure and the equivalent plastic strain are presented in Figs. 19 and 20. It can be seen that the hydrostatic stresses and strains, which are important variables for the onset of fracture, are very high around the crack front. These results agree with the results found in the literature for 2D cutting simulations [10,22].

Fig. 17. Mesh management in yz-plane before fracture.

Fig. 20. Equivalent plastic strain during slitting.

340

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341 [7] J. Huétink, On the simulation of thermomechanical forming processes, Ph.D. Thesis, University of Twente, 1986. [8] P. Schreurs, F. Veldpaus, W. Brekelmans, Simulation of forming processes, using the arbitrary Eulerian–Lagrangian formulation, Comput. Meth. Appl. Mech. Eng. 58 (1986) 19–36. [9] H. Huisman, J. Huétink, A combined Eulerian–Lagrangian three-dimensional FE analysis of edge rolling, J. Mech. Work. Technol. 11 (1985) 333–353. [10] H. Wisselink, Analysis of guillotining and slitting, ?nite element simulations, Ph.D. Thesis, University of Twente, 2000. [11] J.L.F. Aymone, E. Bittencourt, G.J. Creus, Simulation of 3D metal-forming using an arbitrary Lagrangian–Eulerian ?nite element method, J. Mater. Process. Technol. 110 (2001) 218–232. [12] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice-Hall, Englewood Cliffs, NJ, 1969. [13] W. Liu, H. Chang, J.-S. Chen, T. Belytschko, Arbitrary Lagrangian– Eulerian Petrov–Galerkin ?nite elements for non-linear continua, Comput. Meth. Appl. Mech. Eng. 68 (1988) 259–310. [14] J. Wang, M. Gadala, Formulation and survey of ALE method in nonlinear solid mechanics, Finite Elements Anal. Des. 24 (1997) 253–269. [15] H. Bayoumi, M. Gadala, J. Wang, Numerical simulation of metal forming processes, in: J. Huétink, F.P.T. Baaijens (Eds.), Simulation of Materials Processing: Theory, Methods and Applications, Proceedings of the NUMIFORM’98, Balkema, Rotterdam, 1998, pp. 103–108. [16] M. Gadala, J. Wang, Simulation of metal forming processes with the ?nite element method, Int. J. Numer. Meth. Eng. 44 (1999) 1397– 1428. [17] D. Benson, An ef?cient, accurate, simple ALE method for nonlinear ?nite element programs, Comput. Meth. Appl. Mech. Eng. 72 (1989) 305–350. [18] H.C. Stoker, Developments of the arbitrary Lagrangian–Eulerian method in non-linear solid mechanics, applications to forming processes, Ph.D. Thesis, University of Twente, 1999. [19] F.P.T. Baaijens, An U-ALE formulation of 3D unsteady viscoelastic ?ow, Int. J. Numer. Meth. Eng. 36 (1993) 1115–1143. [20] G. Rekers, Numerical simulation of unsteady viscoelastic ?ow of polymers, Ph.D. Thesis, University of Twente, 1995. [21] P.N. van der Helm, J. Huétink, R. Akkerman, Comparison of arti?cial dissipation and limited ?ux schemes in arbitrary Lagrangian Eulerian ?nite element formulations, Int. J. Numer. Meth. Eng. 41 (6) (1998) 1057–1076. [22] D. Brokken, Numerical modelling of ductile fracture in blanking, Ph.D. Thesis, Eindhoven University of Technology, 1999. [23] J. Huétink, P.T. Vreede, J. van der Lugt, Progress in mixed Eulerian–Lagrangian ?nite element simulation of forming processes, Int. J. Numer. Meth. Eng. 30 (1990) 1441–1457. [24] E. Hinton, J. Campbell, Local and global smoothing of discontinuous functions using a least squares method, Int. J. Numer. Meth. Eng. 8 (1974) 461–480. [25] J. Huétink, P.T. Vreede, J. van der Lugt, The simulation of contact problems in forming processes using a mixed Eulerian–Lagrangian ?nite element method, in: A.A. Balkema (Ed.), Numerical Methods in Industrial Forming Processes, Proceedings of the NUMIFORM’89, Rotterdam, 1989, pp. 549–554. [26] J. van der Lugt, A ?nite element method for the simulation of thermo-mechanical contact problems in forming processes, Ph.D. Thesis, University of Twente, 1988. [27] H.H. Wisselink, J. Huétink, Tool deformation during the shape rolling of stator vanes, in: M. Pietzryk, Z. Mitura, J. Kaczmar (Eds.), Proceedings of the Fifth ESAFORM Conference, 2002, pp. 371–374. [28] R. Akkerman, Euler–Lagrange simulation of nonisothermal viscoelastic ?ows, Ph.D. Thesis, University of Twente, 1993. [29] J.P. Ponthot, Ef?cient mesh management in Eulerian–Lagrangian method for large deformation analysis, in: E.G. Thompson, R.D. Wood, O.C. Zienkiewicz, A. Samuelsson (Eds.), Numerical Methods

9. Conclusions The proposed ALE formulation for stationary processes has been applied successfully to rolling and slitting. It is a suitable method to calculate the steady state of such processes, as a complete 3D remeshing can be avoided. Also steady crack growth and the deformation of the tools can be incorporated [27]. The 3D initial mesh is chosen such that all cross-sections perpendicular to the x-axis have an equal topology, which fact has been used for the mesh management procedure. The relocation of the nodes is split into two successive 2D problems, ?rst in x-direction and next perpendicular to this in the yz-plane. The results show that good element shapes and initial re?nements can be preserved during the entire simulation. A suf?ciently ?ne mesh is needed in order to describe the movement of the free surfaces accurately in regions with a large curvature. However, some compromise had to be made between accuracy and the calculation times, as for such 3D problems the total number of elements is large (rolling = 12 144 elements; slitting = 9300 elements). The used algorithm for the convection of the state variables is stable and shows very little cross-wind diffusion, but is diffusive in the streamline direction.

Acknowledgements This research was carried out under number IOPC.94.704.UT.WB of the innovative research programme of the Dutch Ministry of Economical Affairs and under Project number MC1.98056 in the framework of the Strategic Research Programme of the Netherlands Institute for Metals Research in The Netherlands.

References
[1] A.E. Tekkaya, Fully automatic simulation of bulk metal forming processes, in: J. Huétink, F.P.T. Baaijens (Eds.), Simulation of Materials Processing: Theory, Methods and Applications, Proceedings of the NUMIFORM’98, Balkema, Rotterdam, 1998, pp. 529–534. [2] C.J.M. Gelten, Application of automatic remeshing and contact to 3D forging simulations, in: J. Huétink, F.P.T. Baaijens (Eds.), Simulation of Materials Processing: Theory, Methods and Applications, Proceedings of the NUMIFORM’98, Balkema, Rotterdam, 1998, pp. 535–540. [3] E.N. Dvorkin, E.G. Pet?cz, An effective technique for modelling 2D metal forming processes using a Eulerian formulation, Eng. Comput. 10 (4) (1993) 323–336. [4] D. Benson, A mixture theory for contact in multi-materials in Eulerian formulations, Comput. Meth. Appl. Mech. Eng. 140 (1997) 59–86. [5] T. Hughes, W. Liu, T. Zimmermann, Lagrangian–Eulerian ?nite element formulation for incompressible viscous ?ows, Comput. Meth. Appl. Mech. Eng. 29 (1981) 329–349. [6] J. Donea, S. Guiliani, J. Halleux, An arbitrary Lagrangian Eulerian ?nite element method for transient dynamic ?uid structure interaction, Comput. Meth. Appl. Mech. Eng. 33 (1982) 689–723.

H.H. Wisselink, J. Hu? etink / Journal of Materials Processing Technology 148 (2004) 328–341 in Industrial Forming Processes, Proceedings of the NUMIFORM’89, Balkema, Rotterdam, 1989, pp. 203–210. [30] J.P. Ponthot, The use of the Eulerian–Lagrangian FEM with adaptive mesh, applications to metal forming simulation, in: Computational Plasticity, Fundamentals and Applications: Proceedings of the Third International Conference, 1992, pp. 2269–2280. [31] A.H. van der Boogaard, J. Huétink, A.D. Rietman, Iterative solvers in forming processes, in: J. Huétink, F.P.T. Baaijens (Eds.), Simu-

341

lation of Materials Processing: Theory, Methods and Applications, Proceedings of the NUMIFORM’98, Balkema, Rotterdam, 1998, pp. 219–224. [32] H.H. Wisselink, J. Huétink, M.H.H. van Dijk, A.J. van Leeuwen, 3D FEM simulations of a shape rolling process, in: A. Habraken (Ed.), Proceedings of the Fourth International ESAFORM Conference on Material Forming, 2001, pp. 843–846.


相关文章:
更多相关标签: