Local radial basis function collocation method for solving thermo-mechanics of hot shape rolling of steel

Local radial basis function collocation method for solving thermo-mechanics of hot shape rolling of steel V International Conference on Computational Methods for Coupled Problems in Science and Engineering COUPLED PROBLEMS 2013 U. Hanoglu and B. Šarler LOCAL RADIAL BASIS FUNCTION COLLOCATION METHOD FOR SOLVING THERMO-MECHANICS OF HOT SHAPE ROLLING OF STEEL U. HANOGLU*AND B. ŠARLER†* *University of Nova Gorica Laboratory for Multiphase Processes Vipavska 13, SI-5000 Nova Gorica, Slovenia email: umut.hanog[email protected], www.ung.si † Institute of Metals and Technology Laboratory for Simulation of Materials and Processes Lepi pot 11, SI-1000 Ljubljana, Slovenia email: [email protected], www.imt.si Key words: Hot Shape Rolling, Thermo-Mechanical Modeling, Meshless Methods, Local Radial Basis Function Collocation Method, Steel. Abstract. The aim of this paper is to demonstrate the suitability of the novel Local Radial Basis Function Collocation Method (LRBFCM) [1] in a coupled thermo-mechanical problem of hot shape rolling of steel. The physical concept of such a large deformation problem is based on a two dimensional traveling slice model [2], which assumes deformation and heat flow only in the perpendicular direction to rolling. The solid mechanics is, respectively, based on the steady Navier’s equation and the thermal field on the transient heat conduction equation. The displacement and traction boundary conditions are assumed in the mechanical model and Dirichlet and Neumann boundary conditions in the thermal model, both specific for hot shape rolling. The solution procedure is based on local collocation on a five noded influence domain with multiquadrics radial basis functions, augmented with the first order polynomials. The steel used in the calculations is assumed to have an ideal plastic behavior which obeys von Misses flow rule, defined by effective stress (,,)T    in terms of effective strain  , effective strain rate  and temperature T . The LRBFCM results of hot shape rolling of steel for a continuous 5 stand rolling mill in Štore Steel company are presented for the case of rolling of a rectangular billet with initial dimension 80 x 95 mm to a circular bar with diameter of 60 mm. The advantage of the meshless method is in accuracy and straightforward node generation, that does not require any polygonisation. The paper presents one of the increasingly emerging examples of the use of the LRBFCM in industrial applications. 1 INTRODUCTION Hot shape rolling of steel, that usually follows continuous casting of billets and blooms, provides different shapes of long products for various applications such as automotive 116 U. Hanoglu and B. Šarler. 2 industry, construction, shipbuilding, railway, etc. The demand is constantly changing due to newly introduced designs and the production should be quickly adaptable. Respectively, computational modeling of the continuous casting and hot rolling gains its importance for better understanding, control, and better insight into these processes. It helps to improve the quality, productivity, safety and environmental impact of the production. The principal goal of rolling simulation is to connect the process variables, such as the rolling speed, temperature, rolling stand geometry, to the temperature, strain, strain rate, stress field in the billet, its microstructure with static and dynamic recrystallization, and calculation of the rolling torque and power. The modeling of rolling started with Hitchcock [3] where he solved the problem of the roll deformation. The current state of the rolling technology can be perceived from [4] and rolling modeling from [5]. The majority of the simulations in solid mechanics are done by using Finite Element Method (FEM) [6]. This method requires meshing as pre-processing which might be problematic and time consuming, especially in case of rolling, where the shape is constantly and drastically changing. In the last decade, meshless numerical methods [7] started to represent an appealing alternative to the classical numerical methods, such as FEM. Meshless method is a numerical technique that uses a set of arbitrary distributed nodes, both on the boundary and within the computation domain, to represent the solution of physical phenomena. The main feature of meshless methods is omission of the polygonalisation between the nodes which can be remarkably demanding, particularly in realistic 3D geometrical situations. One of the simplest meshless methods, able to solve the fluid flow problems and solid mechanics problems is Local Radial Basis Function Collocation Method (LRBFCM). This method was first developed in [8] for elasticity problems and in [1] for diffusion problems. The idea behind this method is to approximate the function and its derivatives locally over a set of neighboring nodes using RBFs [9] and to use collocation for determining the expansion coefficients. The method has been recently applied to numerous scientific and engineering problems, connected with fluid mechanics [10-14] as well as solid mechanics [15-17]. In this paper the simulation of hot shape rolling is performed by using LRBFCM. 2 PHYSICAL MODEL Figure 1: Scheme of slice model of hot shape rolling 117 U. Hanoglu and B. Šarler. 3 2.1 Slice model assumptions A scheme of the hot rolled billet is shown in Figure 1 with sketch of traveling slices and computational nodes. The temperature and the deformation field of the slice can be computed from the known time dependent boundary conditions. The slice time is on the other hand associated with the position in the rolling mill. The Cartesian coordinate z p measures the length from the beginning of rolling. The billet geometry in z direction is assumed straight and the thermal and mechanical interactions in rolling direction are neglected. The z p coordinate can thus be considered parabolic, while the x p and y p coordinates are elliptic. In this way all fields at a given z p coordinate depend only on the slice history, including its cooling intensity and deformation through the rolls, as a function of time. The slices form at the z start p longitudinal coordinate of rolling and travel in the direction of the z i base vector with the rolling speed roll V , (see Figure 1). For calculating the cooling intensity of the slice as a function of time, a connection between the z p coordinate of the rolling mill and the slice history t is needed () (') ' ; () () start t roll zstart roll roll z t zt Vtdtp Vt t     Vi ,   ( ') ( ') / roll start V t St St , (1) where start t represents the initial time of a slice and () start St and   'St represent the initial and the cross-section of the slice at time 't , respectively. In case when the rolling speed is constant and there is no deformation, we obtain the following simple connection between the z p coordinate of the rolling mill and the slice history     / start z z start roll tz t p p V . A quadrant symmetry of the geometry and fields of the slice with domain  and boundary  are assumed and the north-east quadrant (I-st quadrant) is computationally coped with, described by Cartesian coordinates x p and y p and base vectors x i and y i . 2.2 Thermal model The governing equation of the thermal model is   p T c kT Q t      , (2) with , , ,, , p c TtkQ  standing for density, specific heat, temperature, time, thermal conductivity and heat generation rate due to deformation, respectively. The boundary temperatures of the travelling slice were obtained from the Robin type energy balance equations in point p       ,, ,; R Tt k hT t T t          pp pp n (3) 118 U. Hanoglu and B. Šarler. 4 at slice boundaries  , with h standing for the heat transfer coefficient, R T  standing for the reference temperature. On the two quadrant symmetry axis, h is set to 0. The heat source due to the deformation is calculated as Qd   , (4) where Q is the heat generated due to plastic work which is defined in terms of effective stress  and effective strain  . 2.3 Mechanical model The main aim of the mechanical model is to calculate the displacement field of the slice due to deformation by the roll, in order to get the new shape of the slice. The boundary is divided into natural u  and essential t  part ut   . The governing equation of the mechanical model is TLσ+b 0 , (5) where L is the derivative operator, σ is the vector of stresses, and b is the body force vector, considered b0 . Two dimensional plane strain model is assumed. The essential boundary conditions are considered between the slice and the roll and natural boundary condition are assumed at the boundary parts with no contact as well as along the symmetry lines. The natural boundary conditions are defined as ;, , i ii j ij i n n ij xy     , (6) where in i n is component of unit normals and i  is the prescribed shear stress component, equal to 0 . The essential boundary condition is described as ; , ii u u i xy  , (7) where i u is the displacement and i u is the prescribed displacement due to the action of the roll. The material is assumed to be ideally plastic which means that all the energy put into material turns into plastic deformation at yield stress. The yielding of a solid material is defined by von Mises flow rule. A non-linear plastic modulus P H is used instead of elastic parameters to relate the stresses with strains 2 3 ij ij p ij sH   , (8) where ij s is the stress deviator tensor, p H   is the plastic modulus [18] which is the derivative on the effective stress (   32 ij ij    ) – effective strain (   23 ij ij    ) curve. In metal deforming process, the stress deviator tensor can be assumed to be equal to the 119 U. Hanoglu and B. Šarler. 5 components of the stress tensor. The strain vector ε can be written in terms of displacement vector as ε = Lu . (9) Therefore, the strong formulation of deformation problem gives two individual balance equations in each principle direction on a 2D slice in terms of displacement ,, i ui xy which are 2 22 22 21 0 33 y xx pp x y xy u uu HH p p pp            (10) 22 2 22 21 0 33 yy x pp y xy x uu u HH p pp p            (11) 3 SOLUTION PROCEDURE 3.1 Solution strategy The coupled thermo-mechanical simulation is structured in the following way. First, nodes are generated on the un-deformed slice and temperature and deformation fields are calculated. Afterwards, the nodes are newly generated on the new, deformed shape of the slice, based on elliptic node generation [19]. The calculated values are interpolated to the new position of the nodes and the procedure is repeated, by taking into account that the calculated values of the old slice serve as initial values for the next slice. If there is a roll contact, the new shape is calculated with the mechanical model. 0 ; , iii p p u i xy (12) Afterwards, the temperature is calculated whether there is a contact or not. We assume that the time discretization is made in explicit Euler way.   00 T T t kT Q        (13) A detailed description of the LRBFCM solution of the thermal slice model is given in [11]. The slice is discretized into N  domain and N  boundary nodes, in total NN N    nodes. A five noded sub-domain is associated with each of the nodes. The initial fields are interpolated on each of the subdomains by using collocation with radial basis functions, augmented with first order polynomials. In order to get the results for the field  , it needs to be approximated with interpolation function   n  p and coefficients n      1 P NN nn n       pp , (14) 120 U. Hanoglu and B. Šarler. 6 where the interpolation function is defined by scaled multiquadrics radial basis functions inside five noded influence domains ( 5N   ), and first order polynomials ( 3 P N ). A flowchart of the thermomechanical simulation is given in Figure 2. Figure 2 : Flowchart of the thermo-mechanical simulation             2 22 22 max max 1 2 mea 2 mea / /, 1, , , n x xn y yn n nx ny pp x pp y c px py           p pp p (15) where max x , max y , mea x , mea y represent maximum distance between the five nodes in x i and y i directions, and mean position of the five nodes in x i and y i directions, respectively. c is set to 32. The calculation of the coefficients n  is elaborated in [20]. The derivatives of the functions, needed in the solution, are calculated from the derivatives of the radial basis functions 121 U. Hanoglu and B. Šarler. 7     1 p NN nn n pp           pp ,     22 22 1 ;, p NN nn n xy pp             pp , (16) The solution of the thermal model requires inversion of a matrix of the size P NN   for each of the subdomains and no global matrix is formed. 3.2 Solution of mechanical model The mechanical model is solved by expressing the displacements by radial basis functions     1 ,, P NN nn n u xy          pp (17) The governing equation for node m becomes 2 22 22 1 11 21 0 33 pPP NN NN NN mn mn mn pm xn pm xn yn n nn x y xy HH p p pp                       (18) 2 22 22 1 11 21 0 33 pPP NN NN NN mn mn mn pm yn pm xn yn n nn y xy x HH p pp p                       (19) The coefficients ;, n xy    are expressed with the unknown values of the displacements ;, m u xy    and a global 22NN sparse matrix is formed for calculation of the unknown displacements. The assembling of the respective global matrix is described in [21]. The heat generation rate is calculated as    11 1 2( ) l l ll m m mm l mll Qtt        (20) where  is Taylor-Quinney parameter, defining the ratio of mechanical work turning into heat, m is node number and l is slice number. 3.3 Elliptic node generation After each deformation step, the nodes on the boundary are redistributed with equal distance. After that, the internal nodes need to be distributed as orthogonal as possible with its neighboring nodes. The new arrangement of the nodes is done in two steps. First, Trans Finite Interpolation (TFI) [19] is used for aligning of the nodes in accordance with the boundary nodes. This is followed by Elliptic Node Generation (ENG) [19] which repositions them in an iterative way to make them more orthogonal. 100 iterations are allowed for this purpose. The performance of described LRBFCM has been assessed for large deformations by comparison with the commercial FEM code DEFORM [16]. 122 U. Hanoglu and B. Šarler. 8 4 NUMERICAL EXAMPLE In the numerical example a continuous rolling schedule with 5 rolling stands is considered. The parameters, used in the simulation are given in Table 1 and the characteristics of the rolling mill are given in Table 2. The result are shown in terms of displacement field of an 95 x 80 mm initial billet shape, rolled into circular bar with radius 30r mm. Figure 3 represents the increase of the rolling speed due to the deformation, Figure 4 represents the node positions at the start of rolling and after each of the rolling stands. Figure 5 represents the displacement field after each of the rolling stands. Table 1: Thermal and mechanical parameters used in calculations where 4.16A , 0.23b , 0.214c R is the gas constant and K T is the absolute temperature Heat transfer coefficient to air air h 20 2 W mK Heat transfer coefficient to roll roll h 16000 2 W mK Thermal conductivity of steel k 29 W mK Specific heat of steel p c 630 J kgK Initial rolling temperature fur T 1250 °C Initial rolling speed roll V 0.76 m/s Ambient temperature am T 25 °C Roll temperature roll T 500 °C Taylor-Quinney parameter  0.9 - Time step dt 10-4 s Effective stress    38000 , , exp bc K TA RT      GPa Table 2: Rolling schedule used in the simulation. Stand Number Position (m) Type 1 0.45 OVAL 85 (H) 2 3.45 OVAL 85 (H) 3 6.05 ROUND 40 (V) 4 9.65 OVAL 55 (H) 5 12.65 ROUND 30 (V) 123 U. Hanoglu and B. Šarler. 9 Figure 3: Velocity of slices in the rolling direction. Figure 4: Initial shape of the slice and displacement vectors of the slice at the exit of each of the five rolling stands. 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Position [m] Velocity [m/s] 124