A novel mixed Eulerian–Lagrangian finite-element method for steady-state hot rolling processes

Available online at www.sciencedirect.com

International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

A novel mixed Eulerian–Lagrangian #nite-element method for steady-state hot rolling processes Josef Synkaa;∗ , Alexander Kainzb;1 a

Institute for Industrial Mathematics, Johannes Kepler University, Altenberger Street 69, Linz 4040, Austria b Institute for Robot & Manufacturing Systems, Johannes Kepler University, Linz, Austria Received 4 April 2003; accepted 18 December 2003

Abstract For the #nite-element analysis of steady-state rigid-viscoplastic forming processes in hot rolling of steel, a new and e3cient mathematical model was developed. The method is based on a mixed Eulerian–Lagrangian formulation, using an Eulerian co-ordinate in the rolling direction, while employing Lagrangian co-ordinates in the direction of the thickness and width of the strip. This approach leads to an e3cient algorithm, where the time is eliminated as an independent variable and the velocity components in the Lagrangian directions are replaced by displacement components as independent #eld quantities. Due to this concept, the free surface deformations can be accounted for directly and the problems encountered with pure Eulerian or Lagrangian models now appear with reduced complexity and can thus be tackled more easily. The new model is based on a #nite-element formulation, modi#ed to account for a mixed treatment of velocity and displacement components. ? 2004 Elsevier Ltd. All rights reserved. Keywords: Eulerian–Lagrangian modelling; Rigid-viscoplastic material 7ow; Large deformations; Plane-strain rolling simulation; Computational mechanics; Galerkin-based #nite-element method; Contact and friction; Rigid zone and free surface treatment

1. Introduction For metal production and treatment, the immensely complex and interrelated relevant process parameters require the application of highly sophisticated mathematical models, which have to be veri#ed by actual operational experience. This forms the basis and a vital precondition for manufacturing high-quality products satisfying exact tolerance demands. By performing systematic parameter studies and regression methods, sets of characteristic curves can be obtained, which serve as an input ∗ 1

Corresponding author. Formerly VOEST-Alpine Industrieanlagenbau GmbH& Co, Rolling Mill Technologies, RPH 2, Linz, Austria.

0020-7403/$ - see front matter ? 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2003.12.008

2044

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

for online control and process automation systems. Computer simulations, therefore, not only provide an alternative to or enhancement of experimental studies, but also allow for the prediction of stress– strain relations and distributions in the interior of the slab during the metal forming process. Further, they provide a more detailed insight into fundamental rolling aspects, such as pro#le, 7atness, pro#le actuator ranges and strip spread, to mention just a few. Such computer models, when veri#ed against experiments, can in the long term be used to improve the quality of the rolled product and to control the online process of metal forming. For the latter, fast numerical realisations are essential. Though the #nite-element modelling of steady-state forming processes is not new (see, for example, Ref. [1] or Ref. [2]), the modelling of such processes in an intermediary Eulerian–Lagrangian system —brie7y denoted either as (E–L) or as (E–L) system thereafter—is novel. Due to its set up it combines the advantages of both the Eulerian and the Lagrangian models, but simpli#es their intrinsic problems, such as the appearance of free surface and of the contact problem, respectively, by reducing their complexity: updated Lagrangian formulations, on the one hand, are well suited for problems regarding path dependent material properties and free surfaces, but lead to problems with evolving contacts and suHer from numerical problems, when the mesh is distorted heavily. Moreover, explicit time integrations are necessary in steady-state simulations as well. Eulerian formulations (cf. Ref. [3] or Ref. [4]), on the other hand, are able to cope with large material deformations by elimination of the explicit time dependence, but are less suited for the description of free surfaces. Further, our (E–L) method diHers signi#cantly from the well-known and established “Arbitrary Lagrangian Eulerian” (ALE) #nite-element concepts based on a moving mesh which is independent of the material 7ow: in our approach we perform the calculation of all our quantities of interest only with respect to the intermediary system. This means that no Lagrangian updates and no subsequent non-linear mappings based on convection are required. It is these very aspects in which the ALEand our (E–L) formulation are diHerent. The mixed Eulerian–Lagrangian model requires a systematic transformation of the weak formulation of the governing equations and of its quantities involved, such as the stress-, strain- and strain rate tensors and the velocity #eld, into this intermediary reference frame. The concept used for the determination of the velocity #eld, the mean stress in the strip and the contact pressure between the strip and the work roll (penetration constraint) is based on a so-called mixed variational formulation (see, e.g., [5], [1], or [6]). A representation of the Eulerian strain rate tensor and of the Cauchy stress tensor in terms of the Eulerian–Lagrangian co-ordinates turned out to be most e3cient. The transformation from the Eulerian to the Eulerian–Lagrangian level of description can therefore be regarded as a non-linear mapping of the distorted (true, i.e. Eulerian) strip geometry onto the intermediary reference system in (E–L). The main advantage is that this auxiliary con#guration remains unchanged, when updates of the geometry are performed. As a positive consequence, this renders the streamline integration, which is employed for calculating the actual displacement #elds of the strip in width and thickness direction, to become a simple one-dimensional integral along the Eulerian co-ordinate in the longitudinal direction. The rheological behaviour of the strip is treated by using rigid-viscoplastic constitutive laws within the frame of the Levy–Mises theory. Special care is taken of the contact modelling between strip and work roll, in which the Coulomb friction law and the Tresca law (“sticking friction”) are taken into account. The results for various real-world (hot) rolling processes, especially those for which homogeneous deformation is closely satis#ed, are compared with those from elementary theory based on the one-dimensional KKarmKan–Siebel diHerential equations (cf. the model of Orowan [7]).

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2045

2. The mixed Eulerian–Lagrangian (E–L) concept In order to enhance the modelling and solution properties of the rolling processes under investigation, we frame our new model in a mixed Eulerian–Lagrangian co-ordinate system, brie7y denoted as (E–L) henceforth. In our study we de#ne a co-ordinate system such that the longitudinal, vertical and lateral component of a position vector in this system is respectively aligned with the direction of the length, the thickness or height, and the width of the strip or slab. In all our rolling applications the material is brought in and thus “7ows” in longitudinal direction with prede#ned inlet and outlet conditions, while in vertical and lateral directions no comparable 7ow conditions of that kind are present and the material deformations there are simply due to the forming process. This is accomplished by choosing an Eulerian co-ordinate in the longitudinal direction and Lagrangian co-ordinates in the remaining two directions to account for local deformation eHects directly. This choice will simplify the contact problem which arises with the use of Lagrangian co-ordinates. The free surface problem, which results in the 7ow direction, however, can be circumvented by setting the inlet and outlet region su3ciently far away from the roll gap area (typically 1–2 distances between the roll gap entry and exit locations) to render the deformation eHects negligible at the inlet and outlet sides. This set up is used throughout our study. Thus, a position vector in our ˜ is de#ned as Eulerian–Lagrangian based model, denoted as x, x˜ = (x1 ; X2 ; X3 ):

(1)

Throughout this study, unless otherwise indicated, we will denote quantities related to (E) with small and those related to (L) with capital letters, while quantities related to the mixed system in (E–L) will be indicated with a tilde or hat and, where necessary, also with a superscript m. The position vector, as de#ned by Eq. (1), describes a cross-section of the strip at time t = T at the #xed position x1 in space, while the co-ordinates X2 and X3 denote the respective position of the strip in the reference con#guration at time t = 0. By moving the inlet and outlet su3ciently far away from the roll gap—the area of interest—the formation of a rigid zone and thus a stable and accurate numerical result will be ensured. Due to the spatial co-ordinate employed, the contact problem is reduced by one dimension and, as we will show further below, the calculation of the plastic strain can simply be obtained from the strain-rate by a one-dimensional integration along the streamlines. In analogy with the relations between the Eulerian and Lagrangian levels of descriptions, the displacements between the position vectors in (E–L) and (L), denoted as Uˆ and uˆm , and between the position vectors in (E–L) and (E), denoted as u˜ and u˜m , respectively, are introduced, as depicted in Fig. 1. Therefore, the displacement vectors are de#ned as follows: U(X; t) := x(X; t) − X

and

u(x; t) := x − X(x; t);

(2)

ˆ ˜ U(X; t) := x(X; t) − X

and

˜ t) := x˜ − X(x; ˜ t); uˆm (x;

(3)

˜ t) := x(x; ˜ t) − x˜ u˜m (x;

and

˜ t) := x − x(x; ˜ t): u(x;

(4)

The motion of the particle P in the reference con#guration (in (L)) with position vector X to the location p in the current con#guration (in (E)) is given by x = x(X; t):

(5)

2046

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

(E)

2

F U+dU

(L) dX

U, u

+d

Û+dÛ Û, û

X

dx

(E-L)

x

1 3

Fig. 1. Displacement vectors and deformation gradient tensors in (E), (L) and (E–L).

To ensure that it is invertible, the function x(X; t) in Eq. (5) must be single-valued and continuously diHerentiable at a given time t. Thus the Jacobian of x(X; t) must be non-zero in the reference region:


@xk

= |xk; K | = 0: (6) J =

@XK
This condition can be satis#ed by assuming the motion of a particle between (L) and (E) to be su3ciently smooth. Since according to our de#nition of the intermediary system (E–L), the motion of the particle is subdivided into motions in speci#c co-ordinate directions, these “partial motions” will then also be smooth. Thus, the respective functions of motion between (E–L) and (E) and (E– L) and (L), respectively, will then also be invertible. For our (E–L) position vector, as de#ned in Eq. (1), the #rst equation in Eq. (3) takes the simpli#ed form       x1 (X; t) X1 x1 (X; t) − X1           : ˆ 0 ˜ U(X; t) = x(X; t) − X =  (7)  X2  −  X 2  =   X3 X3 0 The velocity of a particle in the Eulerian–Lagrangian con#guration with regard to the Lagrangian system thus reads   V1 (X; t)   d Uˆ ˆ 0  (X; t) =  V(X; t) := (8)  : dt 0

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2047

Consequently, if we draw our attention only to the temporal change of the position of a particle in our (E–L) setting with respect to its position in (E), we obtain   0   d u˜m ˆ V2 (X; t)  ˜ ˜ (x(X; t); t) = V(X; t) − V(X; t) =  t); t) := (9) v˜ m (x(X;  : dt V3 (X; t) According to the transformation rule for line elements in (L) and (E), viz dx = F · dX, analogous ˆ ˜ F·dX relations between (L) and (E–L) and between (E–L) and (E), respectively, being given as d x= ˜ ˜ and dx = F · d x, can be derived. On combining them, the following relation between the deformation gradients in (E), (L) and (E–L) is obtained: ˆ F = F˜ · F:

(10)

The corresponding Jacobians of the deformation gradient tensors (de#ned as being positive; see also Eq. (6)), i.e. Jˆ := det Fˆ

J := det F;

and

˜ J˜ := det F;

(11)

are then related to each other according to the formula J = J˜ · Jˆ:

(12)

For the transformation of the surface and volume elements between the Lagrangian and current con#guration it can easily be shown that    −1 1 · da = J · dA · F = J · F−T · dA (13) @v

and

@V



 v

@V

1 · dv =

V

J · dV:

(14)

Thereby, the integration domain and its surface (or boundary) are given by v and @v in the Eulerian and by V and @V in the Lagrangian co-ordinate space, respectively. The quantities, da = n · da and dA = N · dA denote the respective surface element vectors in (E) and (L) and are written as the product of the outer normal onto the in#nitesimal surface element and the size of the element. The size of an in#nitesimal volume element in (E) and (L) is respectively denoted as dv and dV , and 1 represents the characteristic (or indicator) function which is one if x describes a point of the integration domain and zero otherwise. In analogy with the de#nitions in (E) and (L) we introduce the surface and volume element vectors in the intermediary system (E–L) as d a˜ = n˜ · d a˜ and d v, ˜ respectively. Then, we can write the transformation rules for surface and volume elements between (E–L) and (L) as d a˜ = Jˆ · dA · Fˆ −1 = Jˆ · Fˆ −T · dA

and

d v˜ = Jˆ · dV

(15)

and, respectively, between (E–L) and (E) as da = J˜ · d a˜ · F˜ −1 = J˜ · F˜ −T · d a˜

and

dv = J˜ · d v: ˜

(16)

2048

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

As discussed earlier, the assumption of a su3ciently smooth function x(X; t) between (E) and (L) implies su3ciently smooth functions between (E), (L) and the intermediary system (E–L). Under ˜ this assumption, the material derivative of a functional f = f(x(X; t); t), de#ned in (E–L), may be written by using the summation convention for index k as



df @f @f = + Vˆ k (X; t) : (17) dt X @t x˜ @x˜k Thereby, the subscripts on the right parenthesis identify the variable held constant partial time derivative. If we apply this de#nition of the substantial derivative to the displacement vector in (E–L), as given by Eq. (7), and use Eq. (8), we obtain the following relation between of the displacement with respect to x˜1 and the velocity components, which is valid processes: vk uk; 1˜ = (k = 2; 3): v1

in taking the our system in the derivative for stationary (18)

The main diHerence between the well-known Arbitrary Lagrangian–Eulerian (ALE) and our mixed Eulerian–Lagrangian formulation can best be shown by considering the mesh velocity in either system as follows (cf. [8–10] for the former). From the material time derivative of a function ˜ and f = f(t; x(t; X)) with respect to (E–L) and (E), respectively, the following f = f(t; x(t; x)) relations can be derived:


@f

@f @x

@f

· + = (19) @t
@x @t
@t
x˜

x

and

x˜




@f

@f @x

@f

+ = : · @t
x @x @t
X @t
X

By subtracting Eq. (20) from Eq. (19), we #nally obtain

@f

@f

@f · c; = +

@t X @t x˜ @x

(20)

(21)

where c denotes the convective velocity, de#ned as the deviation of the material from the mesh velocity, i.e.,

@x

@x

c := − (=V − VM ): (22) @t
@t
X

x˜

It is obvious that the mesh and convective velocities in (E) and (L) are given by VM = 0, c = V and VM = V, c = 0, respectively. In ALE-formulations, the mesh velocity speci#es the corresponding reference frame, i.e. VM = @x=@t|x˜ is prescribed, and the convective velocity is therefore given as c = V − VM , which is non-zero in general. These relations hold true irrespective of whether the underlying process is steady-state or transient. For our mixed Eulerian–Lagrangian formulation, however, when applied to steady-state processes, the nonlinear mapping between (E) and (E–L) is time-independent and thus the mesh velocity remains identically zero, whereas the convective velocity equals the material velocity.

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2049

3. The governing equations In #nite deformation elasto-plasticity the choice and decomposition of the strain and strain rate measures are not as straightforward as in the classical in#nitesimal theory. It can be shown that the rate of deformation tensor D, which is the #nite strain rate measure in the current con#guration, is in general not equal to the sum of the elastic and plastic rate of deformation tensors, De and Dp , D = De + Dp :

(23)

However, in hot metal forming processes the elastic strains are in most cases negligibly small compared to the plastic strains. It is thus reasonable to assume that the elastic deformation gradient, written as Fe = I + e with e usually being of the order 10−3 in metals, may be replaced by the unit tensor (cf. [11], Chapter 7). By ignoring the higher-order in#nitesimal quantities of e , Eq. (23) can then be approximated by (cf. Ref. [12]), D ≈ De + Dp ;

(24)

which is a commonly used relation in the plasticity theory of #nite deformations in metal forming processes. A more detailed discussion on these aspects and a comprehensive investigation on the material objectivity in connection with elasto-plastic models framed in the mixed Eulerian–Lagrangian reference system was performed by Synka et al. [13,21]. Our present model for the simulation of steady-state hot rolling processes is based on the following main assumptions and simpli#cations: on volume constancy, the acceleration forces being negligible compared to the roll separating forces present in the roll gap, on an isotropic rigid-viscoplastic material behaviour and on symmetric conditions at the upper and lower surface of the strip. Further, we assume the friction coe3cient  of being a prescribed phenomenological constant, with the friction being modelled either by a Coulomb or Tresca law (“sticking friction”) or by a truncated law, where the tangential frictional stress is limited by the yield shear stress. Moreover, in the present study we restrict ourselves to isothermal plane-strain conditions without heat transfer between the strip and the surrounding including the rigid work rolls. The strong formulation of the governing equations on domain V for the rigid-viscoplastic model (D ≡ Dp ) may thus be summarised as follows: (div )i = 0;

i = 1; 2;

tr D = 0; f(; z) := 12  :  − @f D = ˙ : @

(25) (26)

y2 (z) = 0; 3

(27) (28)

Thereby, the Cauchy stress tensor  is decomposed into the sum of its deviatoric part,  , and its isotropic mean stress contribution (negative hydrostatic pressure), m , i.e.,  =  + m I:

(29)

2050

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

In the von Mises 7ow rule, as given by Eq. (27), J2 := 12  :  represents the second invariant of the stress deviator, while the colon denotes the dyadic product of two tensors A and B, which is de#ned as A : B = Aij Bij . In our model the hardening parameter z in the yield function f is of the form z = ( S; S˙; !), where ! denotes the temperature at a given location in the material, which is assumed to be known in our applications. The uniaxial equivalent (plastic) strain rate for #nite deformations is de#ned as   S˙ = DS := 23 D : D : (30) S denotes the uniaxial equivalent (plastic) strain and can be obtained by a simple integration along the streamline in the longitudinal (x˜1 -) direction, i.e.,  x˜1 ˙ S(z; x˜2 ) d z: (31) S(x˜1 ; x˜2 ) = v1 (z; x˜2 ) 0 For stationary processes this holds true since in our intermediary system (E–L), the relation between the eHective strain path and the eHective strain rate is, according to Eqs. (17) and (8), of the simple form @S S˙ = v1 : (32) @x˜1 Due to the use of the von Mises 7ow rule, the proportionality factor in the normality rule, given by Eq. (28), can be identi#ed as 

3 DS ; ˙ = 2 S where the eHective Mises stress, denoted by , S is de#ned as  S := 32  :  :

(33)

(34)

The (equivalent) yield stress, y ≡ , S in the von Mises 7ow rule (27) is obtained from the underlying constitutive law for the rolled material. In our applications we utilised either a power-law-type model or the more general Hensel–Spittel formula. The power-law model is of the form y = y ( S˙) := " S˙# ;

(35)

whereas the commonly used Hensel–Spittel formula can be represented as y = y ( S; S˙; !) := $e−m1 ! S m2 S˙ m3 e−m4 S:

(36)

The parameters involved in the constitutive laws depend on the material used in rolling. In our #rst test case, where we utilised a constitutive law of power-law type, we chose " = 150 (N=mm2 ) and #=0:2, a setting which is typically used in the analysis of plane strain rolling of metal slabs (see, for example, the study of Iguchi and Yarita [4]). In all other cases studied, we applied a constitutive law based on the Hensel–Spittel formula for a low carbon steel grade ST37 with constant temperature. Both formulae given above have to be modi#ed to account for the rigid zone eHects as discussed further below. Moreover, since the Hensel–Spittel formula does not provide a good approximation for the yield stress for small values of eHective strain and strain rate, the formula can, for example, be linearised with respect to these quantities. The limiting values for the equivalent strain and strain

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2051

rate used in all our applications were chosen as Smin = 0:01 and S˙min = 0:01, respectively. These settings turned out to yield a qualitatively correct approximation of the deformation properties of the rolled material even in or close to the rigid zones. 4. The weak formulation and two models in (E–L) By applying the method of weighted residuals to the equations of motion and the plastic incompressibility condition, given by Eqs. (25) and (28), with the weighting factors, respectively, denoted as %wi (i = 1; 2) and %m , and integrating the equations over the given domain V , the weak formulation of the governing equations may be written in compact form as

 @ji dvi %wi + tr D%m dv = 0: (37) −' @xj dt V For convenience, let us introduce an auxiliary tensor D(%w) in analogy with the rate of deformation tensor D(v), i.e. we de#ne

1 @%wi @%wk Dik (%w) := : (38) + 2 @xk @xi In analogy with the tensor D, we will now decompose this auxiliary tensor into a so-called deviatoric and volumetric part as D(%w) = D (%w) + tr D(%w)I;

(39)

where tr D(%w) = Dii (%w) denotes the trace of D(%w). Using Green’s theorem (integration by parts), the symmetry of the stress tensor and the auxiliary tensor D(%w), Eq. (27) may be rewritten in the form     ij Dij (%w) dv + 'v˙i %wi dv + tr D%m dv = ti(n) %wi ds; (40) V

V

V

S

where S := @V denotes the surface or boundary of the domain V , ds := da is the size of an in#nitesimal surface element vector, and the surface traction vector t (n) is de#ned as ti(n) = nj ji

or; in tensor form; as

t(n) = nT :

(41)

In all our applications, the surface integral will be divided into an integral over the inlet and outlet region, the symmetry plane, and the upper surface of our domain, where the latter is subdivided into a free surface and a contact area. In 7at rolling processes the acceleration forces are by at least three orders of magnitude smaller than the roll forces prevailing in the contact region. Thus, without signi#cant loss of accuracy the acceleration term in Eq. (40) will henceforth be omitted in our model. A detailed discussion on appropriate boundary conditions for the hot rolling process will be given at the end of this section. Under these assumptions and with the use of the decomposition of the stress tensor, according to Eq. (29), Eq. (40) may now be written as       : D(%w) dv + m I : D(%w) dv + tr D%m dv = t(n) %w ds: (42) V

V

V

S

2052

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

After substitution of the normality and the von Mises 7ow rule, respectively being given by Eq. (28) with Eqs. (33) and (27), and by utilising Eq. (39), the #nal form of the weak formulation of the governing equations in the Eulerian system (E) reads     2 y (z)   m tr D(%w) dv + tr D%m dv = t(n) %w ds: (43)  D : D (%w) dv + S 3 D V V V S Applying the transformation rules between the systems in (E) and (E–L), as presented in Section 2, to Eq. (43), yields the weak formulation of the governing equations in the Eulerian–Lagrangian system (E–L), which may be written in the form V˜ I + V˜ II + V˜ III − S˜ I = 0

(F :=) with



V˜ I :=



V˜



V˜ III := S˜ I :=

2 y (z)  D : D (%w) J˜ d v; ˜ 3 DS 

V˜

V˜ II :=

 S˜

(44)

V˜

m tr D(%w)J˜ d v; ˜ ˜ tr D%m J˜ d v;

and

˜ d s: t(n) %wJ˜ F˜ −T n ˜

(45)

˜ denotes the size of an in#nitesimal surface element vector in (E–L) and Thereby, d s˜ := d a is obtained by using Eq. (16). It should be stated that the rate of deformation tensor D and the auxiliary tensor D(%w) are yet to be transformed to our system in (E–L). Both tensors may due to their de#nition be written in combined form as D(-) with D(-) = 12 (L(-) + LT (-))

(46)

˜ F˜ −1 , and - set equal either to v or to %w. The quantity L(-) := @-=@x transforms as @-=@x=(@-=@x) − 1 ˜ where F˜ := @x=@x denotes the inverse deformation gradient tensor between (E–L) and (E). Depending on the choice of the weighting functions, we consider two diHerent Galerkin-type models in (E–L). In the weak formulation, Eq. (44) with Eq. (45), the dependent quantities can be expressed either in terms of the three unknowns v1 ; v2 and m or, by using relation (18), in terms of the unknowns v1 ; u2 and m . Consequently, the corresponding Galerkin-type weighting functions will be chosen either as (%v1 ; %v2 ; %m ) ∈ Vv1 × Vv2 × Vm ;

(47)

(%v1 ; %u2 ; %m ) ∈ Vv1 × Vu2 × Vm ;

(48)

or as

where Vx denotes the vector space for the unknown x.

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2053

According to their derivation we will refer to these models as natural Galerkin-type model either in (E) or (E–L), or brie7y either as vv-model or as vu-model according to the characteristic unknowns involved. Appropriate vector spaces for the unknowns over the domain . with proper boundary conditions would thus be either [H 1 (.)]2 × L2 (.) or H 1 (.) × H 2 (.) × L2 (.), where H 1 (.) denotes the well-known Sobolev space of functions whose functions and generalized #rst derivatives thereof  belong to L2 (.), the Lebesgue space of measurable functions v de#ned on . for which . |v(x)|2 d x is #nite. H 2 (.) denotes the Sobolev space of functions whose generalized derivatives up to order 2 belong to L2 (.), i.e., H 2 (.) := {v ∈ L2 (.) | Di v ∈ L2 (.); |i| 6 2} and Di v is called the ith generalized derivative of the function v. Before we formulate appropriate boundary conditions for our problem and for both models, we will brie7y discuss their main diHerences: • The vv-model is equivalent to the standard modelling approach used in (E) based on the unknowns v1 , v2 and m , followed by a transformation to the Eulerian–Lagrangian co-ordinate system. The displacement component u2 , arising in the transformation rules and required for a mesh update at the free surface, is not actually available but has to be calculated from relation (18) once the velocities are known. Hence, in contrast to the vu-model, all quantities related to u2 will lag behind in an explicit-type vv-model. • Due to relation (18) between the velocity components and u2 , second-order derivatives of u2 arise in the weak formulation of the vu-model. This requires not only the continuity of the trial function, but also the continuity of its #rst derivative with respect to x˜1 , i.e. a C1 -type continuity for u2 . It necessitates a combination of Hermite and Lagrange interpolation for the set up of shape functions for the unknowns. However, the main advantage of the vu-model is that it directly accounts for the local deformations at the free surface and in the contact region and simpli#es the imposition of displacement boundary conditions, which naturally arise in the formulation of the contact problem and the impenetrability constraint. For constitutive laws, where the yield stress is strain-independent (cf. the power-law type model (35)), the vu-model does not require an outer loop for streamline integration compared to the vv-model, where it is still necessary for the geometry update.

5. Boundary conditions The boundaryof V , denoted as S := @V , is divided throughout this study into #ve distinct regions S1 ; : : : ; S5 with i=1(1)5 SS i = S and Si ∩ Sj = ∅ for i = j ∈ {1; : : : ; 5}. The #ve regions are subsequently de#ned as the inlet and outlet regions, the symmetry plane, the free surface and the contact area, i.e., S1 := Sin , S2 := Sout , S3 := Ssymm , S4 := SF and S5 := SC , respectively. For hot rolling processes, as sketched in Fig. 2 and investigated in our study, the following boundary conditions are implemented: • At the inlet and outlet, boundary conditions according to the metallurgical and technological demands are posed for the back and front tensile forces, TB and TF , respectively, which are present to stabilise the rolling process, and/or for the velocity #eld. In order to ensure numerical

2054

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060 Backup roll

Work roll

TB/2

h0 2

S0 S1

Slab

h1 2

TF /2

Symmetry plane

Fig. 2. Sketch of a rolling process.

• • •

•

convergence, these boundaries have to be located su3ciently far outside the roll gap area. For the vu-model we may further assume that the mesh is yet in an undeformed state. This means that the meshes in (E) and (E–L) coincide at the inlet, which is equivalent to the displacement being identical to zero there. At the symmetry plane in addition to Newton’s third law the condition that all vectors are being mirrored there must hold. The free surface, located at the upper boundary in Fig. 2 and outside the contact region, can be assumed to be free of traction with the 7ow being tangential to it. At the contact region, located between the strip and work roll and bounded by the roll gap entry (S0 ) and exit (S1 ) locations, the eHect of friction has to be taken into account. Since the position of the neutral point or zone is a priori unknown, the frictional stress is approximated by a velocity-dependent frictional stress to eliminate its sudden change of direction at this location. To avoid erroneous formulations for rigid-viscoplastic materials, the friction has to be limited by the maximum shear stress, denoted by ksh (cf. Ref. [2]). Since both the Coulomb and the Tresca sticking friction law are discontinuous at the location of the neutral point (or region), where the slip velocity changes its sign, the velocity-dependent friction law was regularised as proposed by Kobayashi et al. [1]. In the contact zone the material must at no point penetrate the surface of the work roll, which is stated by the impenetrability constraint. In the Eulerian as well as the Eulerian–Lagrangian description the strip has to be #xed to at least one point in the contact region. The formulation is based on the well-known Signorini problem (of elasticity; cf. Ref. [14]). The iterative procedure for its implementation is based on the assumption that there exists a unique shape of the contact zone, which yields the prescribed inlet and outlet thicknesses of the slab, denoted as h0 =2 and h1 =2, respectively, in Fig. 2. This is accomplished by a so-called contact algorithm in which nodes are either removed from the contact surface if traction prevails or added if they penetrate the work roll in this region. In contrast to the vv-model, the implementation of such a displacement-type boundary condition in the vu-model is an easy task.

In the summary of the boundary conditions for both our models in (E–L) the following notation is used: the prescribed back and front tensile forces are denoted by TS B and TS F . The values of the strip entry and exit velocities are given by C0 and C1 with e1 and n being the unit vector in longitudinal direction and the normal vector onto the slab surface, respectively. The surface traction vector t(n)

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2055

is decomposed into its tangential (4t(n) ) and normal (4(n) ) component. The tangential vector onto the roll and the friction coe3cient are, respectively, speci#ed by tr and . The slip velocity vs is de#ned as vs := v − vroll , where v denotes the strip velocity and vroll the roll speed. The prescribed vertical displacement of the work roll is denoted as u2roll . And #nally, the value of the smoothing parameter −4 −3 v in the regularised friction law is chosen in the range of [10 ; 10 ] times the averaged slab velocity, as proposed by Kobayashi et al. [1]. This setting proved to work su3ciently well in all our hot rolling applications. The boundary conditions for both models, where the direction of rolling is aligned with x˜1 , may thus be written as follows: vv-model: TB = TS B

and

v = C 0 e1

on Sin ;

(49)

TF = TS F

and

v = C 1 e1

on Sout ;

(50)

v·n=0

and

12 = 0

on Ssymm ;

(51)

nT = : t(n) = 0 v · n = 0;

and

v·n=0

contact algorithm;

4t(n) = −min{mksh ; −4n(n) }

on SF ;

(52)

and

(53)



vs · tr 2 arctan 7 v

on SC :

(54)

vu-model: TB = TS B ;

v1 = C0 ;

u2; 1˜ = 0

(and u2 = 0)

TF = TS F ;

v1 = C1 ;

u2; 1˜ = 0

and

(u2 = 0; )

u2; 1˜ = 0

t(n) = 0 u2 = u2roll ; 4t(n)

=

and

12 = 0

on Sin ;

12 = 0

on Sout ;

on Ssymm ;

on SF ;

(56) (57) (58)

contact algorithm;

−min{mksh ; −4n(n) }

(55)

and



vs · tr 2 arctan 7 v

(59) on SC :

(60)

It should be noted that the operator D(v), as de#ned in Eq. (46) with - := v, which is linear in v, transforms in the two-dimensional case to D(v1 ; u2 ), which according to Eq. (18) forms a nonlinear (quasi-linear) operator in (E–L) with respect to the unknowns v1 and u2 . Hence, the vu-model introduces a nonlinearity in the governing equations, compared to the vv-model, even for linear constitutive laws. The advantage of using the vu-model rather than the vv-model, however, is that no condition on the velocity #eld has to be posed at the free surface. This lies in the fact that for the

2056

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

vu-model the normal vector is obtained directly from the displacement changes at the free surface, i.e. (u2; 1˜; −1)T nroll := −nslab =  ; 1 + u2;2 1˜

(61)

and, due to relation (18), the velocity #eld naturally vanishes in normal direction. In contrast, this condition on the velocity #eld of being tangential to the free surface, stated by v · n = 0, has to be posed explicitly in the vv-model and can only be achieved iteratively.

6. Solution process In order to evaluate the performance of our novel vu-model in the mixed Eulerian–Lagrangian co-ordinate system compared to the vv-model and to simpli#ed models based on the KKarmKan– Siebel diHerential equations in the case of homogeneous deformation processes, we employed the well-known #nite-element method to our problem. Since the displacements and deformations of the mesh occur in vertical direction only and due to the nonlinear mapping between (E) and (E–L), as described in Section 2, a rectangular domain is obtained in (E–L) in our applications. Hence, a simple rectangular mesh can be used for discretisation in the (E–L) frame of reference (see Fig. 3). Consequently, rectangular #nite elements can be utilised for the discretisation of the domain of interest. In the vv-model the standard #nite element of C0 -continuity with bi-linear or bi-quadratic approximation of the velocity #eld in two dimensions and a constant representation of the mean stress su3ces to obtain a stable and converged solution (cf., e.g., Refs. [1,6]). Due to the appearance of second-order derivatives of u2 in the vu-model, a #nite element of C1 -continuity and Hermite interpolation has to be used in the #nite-element modelling: An element, which uses bilinear and cubic-linear interpolation for the longitudinal velocity and the vertical displacement component, respectively, and a constant approximation of the mean stress turned out to yield a stable and converged solution in this case. In addition, this #nite-element type when applied to diHerent hot rolling situations, ranging from strip to slab rolling, appears to be even more robust than its vv-model counterparts (see Ref. [15]). This #nite-element approach results in a system of non-linear equations, which is then linearised using either the Newton–Raphson method or a double dog-leg strategy (see Ref. [16]). The latter approach was implemented in a separate VAI-internal Eulerian–Lagrangian #nite-element code based on a vv-model approach, which was simultaneously developed to a university-based code for x2

~ x 2 = X2

x1

~ x =x 1

1

Fig. 3. Distorted and rectangular mesh in (E) and (E–L), respectively.

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

2057

both the vu- and vv-model and incorporates a non-circular arc-model, where the elastic work roll deformation is treated by utilising the Jortner in7uence function technique (see, e.g., Ref. [17]). The VAI-developed model is documented in a paper of Kainz and Finstermann [18], and in several unpublished internal technical VAI-reports. The appearance of a rigid zone in the vicinity of the neutral point (or zone) and outside the roll gap region require a special treatment since the vanishing strain rates would yield to zero divisions in the leading term of Eq. (44). Further, the integration scheme must be adopted such as to avoid volume locking eHects. Only then the system of linear equations, solved via special sparse-matrix techniques, including scaling of the diHerent quantities, will yield a stable and su3ciently accurate solution to our problem. For a detailed discussion on the various techniques applied in the numerical solution process, the reader is referred to [15,19,20], respectively.

7. Veri2cation of the vu-model In order to validate our novel vu-model, as outlined above, we also implemented the more standard-like vv-model in our numerical code and tested both against simpli#ed rolling models based on the KKarmKan–Siebel diHerential equations for homogeneous deformation processes (KKarmKan –Siebel, Sims, Orowan; [7]). An industrial application with nearly homogeneous deformation properties is given, for example, by the following setting, as provided by our industrial partner, VOEST Alpine Industrieanlagenbau GmbH & Co. (VAI), Linz: A strip of steel grade ST37 (low carbon) with a thickness reduction from 2 to 1:2 mm using a work roll of radius 350 mm and a prescribed angular velocity of 21:4 rad=s. Concerning the tribologic aspects, a Coulomb friction law limited by the shear yield stress and assuming a friction coe3cient of  = 0:1 was utilised. No tensile forces were applied at the strip inlet and outlet and the temperature #eld inside the strip was assumed to be homogeneous and #xed at a value of 900◦ C. Table 1 summarises the results for some essential rolling parameters, as obtained by utilising the vu-model framed in the Eulerian–Lagrangian co-ordinate space and compared to data attained by solving the KKarmKan–Siebel diHerential equation. For small roll gap aspect ratios, de#ned by the ratio of the mean strip thickness to the contact length, hm =ld , both methods are expected to yield similar results. This is con#rmed by our test results. Line plots for the #eld distributions of interest, the velocity and vertical displacement components, the mean stress m , the Cauchy shear stress components and the yield stress y are depicted in Fig. 4 for the vu-model. Since the results obtained with the vv-model deviate by less than one per cent from those of the vu-model, they are not presented. From the distribution of the longitudinal velocity component and of the mean stress one can easily see that the deformation of the strip is indeed close to homogeneity. The good agreement of the vu-model results with the KKarmKan–Siebel data—brie7y denoted as KS in the #gures—veri#es the validity of our novel vu-model. In Fig. 5 the distribution of the normal stress multiplied by the friction coe3cient, which equals 0.1 in our case, of the velocity-dependent tangential stress and of the ratio between these two quantities in longitudinal direction is presented. From the latter sub#gure it can be seen that the neutral point is signi#cantly shifted toward the roll gap exit, as expected from practical experience, and that the numerical approximation of the sign-function is su3ciently accurate.

2058

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

Table 1 Results for the vu- and KKarmKan–Siebel model for hm =ld = 0:088 Quantity Number of elements in x1 -direction Number of elements in x2 -direction Number of equations Speci#c rolling torque (kN) Speci#c rolling force (kN/mm) Horizontal force residuum (kN/mm) Strip entry velocity (mm/s) Strip velocity at neutral point (mm/s) Strip exit velocity (mm/s) NS1 =S0 S1 Backward slip (%) Forward slip (%) Gap angle (deg) Neutral angle (deg) Number of Newton steps Contact length, projected onto x1 (mm)

vu-model

K–S model

102 13 5652 60.13 7.37 0.007 4936 7490 8215 0.38 34.10 9.68 2.69 1.01 54 16.45

— — — 55.64 7.22 0.002 4895 7489 8158 0.37 34.64 8.93 2.74 1.00 — 16.73

At center nodes: S11

v1 vs. z1 - from top boundary to symmetry plane (cf. u2)

200

10000

0 -200

σ

v1

11

8000 6000 4000

KS

-400 -600 0

5

10

15

20

25

30

35

40

0

5

10

15

20

25

30

35

20

25

30

35

20

25

30

35

z1 S22

v2 vs. z1 - from top boundary to symmetry plane (cf. u2)

500

100

0

σ22

v2

0 -100

-500

KS

-200 -300

0

5

10

15

20

25

30

35

40

-1000

0

5

10

15

z1 S12

u2 vs. z1 - from top boundary to symmetry plane

100

0.4

50

12

0

0

σ

u2

0.2

-50

-0.2 -0.4

0

5

10

15

20

25

30

35

-100

40

0

5

10

15

z1 Effective stress (yield stress)

Mean stress vs. z1 - from top boundary to symmetry plane (cf. u2)

300

0

200

KS

σ

y

σm

500

-500 -1000

100

0

5

10

15

20

z1

25

30

35

0

0

5

10

15

20

Slab length

Fig. 4. Results for v1 ; v2 ; u2 ; m ; y and the Cauchy stress tensor components.

25

30

35

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060 −µτ n vs. z1

−sgn(v s) |τT| vs. z1

100

100 T

−sgn(v ) |τ |

−µτ n KS

50

n

s

50

−µ

2059

0 -50

S0

0

10

S1

20

30

0 −50 −100

40

S0

0

10

S1

20

30

40

−sgn(v ) |τ | / −µτ vs. z s

T

n

1

2 1 0 −1 −2

S

S

0

0

10

1

20 z

30

40

1

Fig. 5. Results at the top boundary.

8. Conclusions In the present paper a new modelling approach for the simulation of steady-state hot rolling processes in a mixed intermediary Eulerian–Lagrangian frame of reference was presented. The transformation of the governing equations and the corresponding weak formulation for rigid-viscoplastic materials with large deformations, as originally given in the Eulerian co-ordinate system, to this intermediary frame of reference were derived and discussed in detail. Comparisons with elementary rolling models based on the KKarmKan–Siebel diHerential equation for a homogeneous deformation case, for which the latter models yield reliable agreement with experimental investigations, were performed in order to validate our new formalism. The results obtained for this strip rolling test case with our novel vu-model were in excellent agreement with the KKarmKan–Siebel data and nearly identical to the vv-model predictions. For other rolling situations, ranging from thick slabs to ultra-thin hot strips, our new vu-model proved to be a valuable extension to existing rolling theories. Further aspects of the numerical methodology of our vu-model and its application to various test cases of practical interest will be subject of a subsequent publication (cf. [15]). Acknowledgements The #nancial support from the Christian-Doppler Society, and from the Federal Ministry for Economics and Labour (BMWA) and the Province Upper Austria in the framework of the Industrial Mathematics Competence Center, as well as from the Industrial Plant Building Company VAI (VOEST-Alpine Industrieanlagenbau GmbH & Co) is gratefully acknowledged. Special tribute should be paid to Dr. Gerhard Finstermann, Head of the VAI-Linz Cold Rolling Department, to O.Univ.-Prof. Dr. Heinz W. Engl and to O.Univ.-Prof. Dr. Klaus Zeman for their valuable contributions.

2060

J. Synka, A. Kainz / International Journal of Mechanical Sciences 45 (2003) 2043 – 2060

References [1] Kobayashi S, Oh SI, Altan T. Metal forming and the #nite-element methods. New York, Oxford: Oxford University Press; 1989. [2] Montmitonnet P, Buessler P. A review on theoretical analyses of rolling in Europe. ISIJ International 1991;31(6): 525–38. [3] Karabin M, Smelser R. A quasi three-dimensional analysis of the deformation processing of sheets with applications. International Journal of Mechanical Sciences 1990;32(5):375–89. [4] Iguchi T, Yarita I. 3-Dimensional analysis of 7at rolling by rigid-plastic FEM considering sticking and slipping frictional boundary. ISIJ International 1991;31(6):559–65. [5] Bathe KJ. Finite element procedures in engineering analysis. Englewood CliHs, NJ, USA: Prentice-Hall; 1982. [6] Zienkiewicz OC, Taylor RL. The #nite element method, 4th ed., Vol. 1. London, New York: McGraw-Hill Book Company; 1994. [7] Orowan E. The calculation of roll pressure in hot and cold 7at rolling. Proceedings of the Institution of Mechanical Engineers 1943;150. [8] Hughes TJR, Liu WK, Zimmerman TK. Lagrangian–Eulerian #nite element formulations for incompressible viscous 7ows. Computer Methods in Applied Mechanics and Engineering 1981;29:329–49. [9] Liu WK, Belytschko T, Chang H. An arbitrary Lagrangian–Eulerian #nite element method for path dependent materials. Computer Methods in Applied Mechanics and Engineering 1986;58:227–46. [10] Liu WK, Chang H, Belytschko T. Arbitrary Lagrangian and Eulerian Petrov-Galerkin #nite elements for nonlinear continua. Computer Methods in Applied Mechanics and Engineering 1988;68:259–310. [11] Khan AS, Huang S. Continuum theory of plasticity. New York, Chichester: Wiley; 1995. [12] Mattiasson K. Continuum mechanics principles for large deformation problems in solid and structural mechanics. Publication 81:6, Department of Structural Mechanics, Chalmers University of Technology, G[oteborg, Sweden, 1981. [13] Synka J. Bericht zur Machbarkeitsstudie “Warmwalzen von Stahlbrammen—Modellierung und Simulation”, Internal Report, Institute for Industrial Mathematics, Joh. Kepler University Linz, Austria, 1997, unpublished. [14] Kikuchi N, Oden JT. Contact problems in elasticity: a study of variational in-equalities and #nite element methods. SIAM Studies in Applied Mathematics. Philadelphia, PA: SIAM; 1988. [15] Synka J, Kainz A, Maringer R, Obereder A. Numerical analysis of a new Eulerian–Lagrangian #nite-element method applied to steady-state hot rolling processes, to be submitted for publication. [16] Dennis J, Schnabel R. Numerical methods for unconstrained optimization and non-linear equations. Englewood CliHs, NJ, USA: Prentice-Hall; 1983. [17] Jortner D, Osterle JF, Zorowski CF. An analysis of cold strip rolling. International Journal of Mechanical Sciences 1960;2:179–94. [18] Kainz A, Finstermann G. A new Eulerian–Lagrangian hybrid #nite element method for the numerical simulation of stationary rolling processes. In: Benyon JH, et al., editors. Modelling of metal rolling processes, Vol. 3. IOM, conference papers, Church House, London, UK, 1999. p. 104–13. [19] Synka J, Kainz A. A novel approach in simulating the hot rolling process. ECMI-Newsletter 2000;28:11–3. [20] Synka J. Steady-state simulation of 7at rolling processes—a novel mixed Eulerian–Lagrangian approach with applications to plane strain problems. PhD thesis, Joh. Kepler University, Linz, Austria, 2002. [21] Synka J, Kainz A. Extension of the concept of material objectivity to mixed Eulerian-Lagrangian reference systems. Acta Mechanica 2003;166(1–4):13–25.