Lobachevskii Journal of Mathematics Vol. 13, 2003, 3 – 13

©R. Dautov, R. Kadyrov, E. Laitinen, A. Lapin, J. Pieska¨, V. Toivonen

R. Dautov, R. Kadyrov, E. Laitinen, A. Lapin, J. Pieska¨, V. Toivonen
ON 3D DYNAMIC CONTROL OF SECONDARY COOLING IN CONTINUOUS CASTING PROCESS

ABSTRACT. In this paper a 3D-model for simulation and dynamic control of the continuous casting process is presented. The diffusion convection equation with multiphase transition is used as a simulation model. The developed model is discretized by finite element method and the algebraic equations are solved using pointwise relaxation method. Two different type of methods are used to control the secondary cooling, namely PID and optimal control method. The numerical results are presented and analyzed.

________________

Key words and phrases. Stefan problem, continuous casting process, optimization problem, FEM, implicit schemes.

1. Introduction

In the continuous casting process the molten steel is poured into a bottomless mold which is cooled with internal water flow. The cooling in the mold extracts heat from the molten steel and initiates the formation of the solid shell. The shell formation is essential for the support of the slab after mold exit. After the mold the slab enters into the secondary cooling area in which it is cooled by water mist sprays. The secondary cooling region is usually divided into cooling zones in which the amount of the cooling water can be controlled separately.

The control of cast cooling is of central importance in continuous casting process because it has a considerable influence on formation of cracks and other defects which can be formed in cast material. The cast should be cooled down according to a certain temperature field which depend on e.g. steel quality, cast shape, casting speed. Accurate knowledge of temperature field and liquid pool length is also important especially when using soft reduction or ”near final shape” casting. Many numerical models for simulation of the casting process have been developed in recent years [6, 9, 10]. Some of them have been applied to control and optimize casting process [4, 5, 11]. To our knowledge, all real-time industrial control applications are based on two dimensional models. In many cases two dimensional models are sufficient for control purpose. However, nowadays the highly automated and instrumented casting machines in steel factories allows the use of more sophisticated simulation and control models.

Our aim is to use the developed 3D-model in dynamic process control. We have considered PID and optimal control methods to control the secondary cooling. The two methods are very different from each other. The PID algorithm is quite simple and computationally inexpensive. However, the use of PID method is limited to the control of the surface temperature. Moreover, the PID algorithm contains experimental tuning parameters. On the other hand the optimal control method can be very complicated. The optimal control method minimizes a cost function which is constructed by means of metallurgical cooling criteria. Several different cooling criteria can be used, e.g. the maximum length for the liquid pool, the maximum reheating and cooling rates on the slab surface, the minimum and maximum temperature at the unbending point and the maximum and minimum temperature on the surface along the casting machine [6]. The development work of our optimal control model is still ongoing. Therefore we use quite simple cost function in our numerical example.

2. Mathematical formulation of state problem

Let Ω ⊂ R2 be a rectangular domain [0,LX] × [0,LY ] and V = Ω × [0,LZ] be a 3D domain. The boundary Γ = ∂V consists of the parts:

Γ0= Ω ×{0}, ΓN= {(x,y) ∈ ∂Ω : x = 0 ∨ y = 0}× [LM,LZ], ΓS= {(x,y) ∈ ∂Ω : x ⁄= 0 ∧ y ⁄= 0}× [0,LZ) ∪ Ω ×{LZ}, ΓM= {(x,y) ∈ ∂Ω : x = 0 ∨ y = 0}× [0,LM],

where LM is the length of the mould, LZ is the length of the strand and LX,LY are the width and thickness of the calculation domain.

PIC
Figure 1: The calculation domain.

We define the temperature T = T(x,y,z,t), dependent enthalpy function H(T) and the Kirchoff’s temperature K(T) by

H(T) = ρ∫ 0T c(ξ)dξ + L(1 − f s(T)),K(T) = ∫ 0T k(ξ)dξ,

where ρ, c(T), k(T), and L are density, specific heat, thermal conductivity and latent heat, fs (T) is the solid fraction.

The mathematical model of the continuous casting problem can be written as:

∂H(T) ∂t + v∂H(T) ∂z − ΔK(T) = 0 in V× (0,tf], T = T0 on Γ0 × (0,tf], ∂K(T) ∂n + h(T − Tw) + σε(T4 − T ext4) = 0on Γ N × (0,tf], ∂K(T) ∂n = 0 on ΓS × (0,tf], ∂K(T) ∂n = Q on ΓM × (0,tf], T(x,y,z, 0) = T0  in V. (1)

Here n is the unit vector of outward normal to ∂V, h is the heat transfer coefficient, v is the casting speed, and T w ,Text are known temperatures. The σ is the Stefan-Boltzmann constant and ε is the emissivity. The cooling capacity Q in the mould is known constant and tf is the simulation time.

3. Mesh approximation

We partition Ω into a set of quadrilateral finite elements as shown on figure 2.

PIC

Figure 2: Partition of Ω on quadrilateral finite elements.

Step sizes are smaller near left and down boundaries of Ω, and all quadrilaterals are rectangular except those elements near the round corner. The 1D-domain [0, LZ] is divided by nz mesh points in z-direction. We decompose our 3D domain V into a set of prismatic 3D finite elements with quadrilaterals in its cross section. Since the constructed mesh is Cartesian product of 2D mesh in Ω and 1D mesh on [0,LX], basis functions of finite element space will also be a product of 2D and 1D basis functions.

We discretize our problem in time by using semi-implicit approximation, more precisely, we approximate the term ∂ ∂t + v(t) ∂ ∂zH by using the characteristics of the first order differential operator [2, 7, 8]. Namely, if (x,y,z,t) is the mesh point on the time level t we choose z˜ = z −∫ t−τtv(ξ)dξ and approximate this term by:

∂ ∂t + v(t) ∂ ∂zH ≈ 1 τ H(x,y,z,t) − H(x,y,z˜,t − τ) . (2)

We denote by H˜(x,t − τ) = H(x,y,z˜,t − τ). If z˜ < 0 then we set H˜(x,t − τ) = H(x,y, 0,t − τ).

The weak solution of the semi-discrete problem is defined by integral identity

∫ V(H(T) −H˜)∕τηdx∫ V∇K(T) ⋅∇ηdx + ∫ ΓN{h(T − Tw) + σε(T4 − T ext4)}ηdΓ N + ∫ ΓMQηdΓM = 0 (3)

for all trial functions η ∈ V 0 = {η ∈ H1(V) : η = 0 on Γ0}. Here H1 (V) is the Sobolev space of first order and T = T(x) = T(x,y,z,jτ), j ≥ 1.

Let V h be a finite element approximation of the space H1 (V), V h T0 = {v h ∈ V h(V) : vh = T0h on Γ0} with T 0h being the V h -interpolant of T 0 and obviously defined subspace V h0. Let further ϕ1(x),…,ϕn(x), be the basis functions in V h. Thus, the function

Th(x) = ∑ i=1nT iϕi(x)

is the finite element approximation of T(x), T = (T1,…,Tn)t is the vector of the unknown nodal values of the function T h (x). Below the notations are: H(T) = (H(T1),…,H(Tn))t, K(T) = (K(T1),…,K(Tn))t and

fh(Th) = ∑ i=1nf(T i)ϕi(x)

for every function f(T), which depends on T(x).

Approximation of the semi-discrete problem (3) by a finite element scheme is defined as

∫ VHh(Th) −H˜h τ ηdx + ∫ V∇Kh(Th) ⋅∇ηdx + ∫ ΓN{h(Th − Twh) + σε((T4) h − (Text4) h)}ηdΓN + ∫ ΓMQhηdΓM = 0 (4)

for all trial functions η ∈ V h0. Discrete problem (4) is equivalent to the system of nonlinear algebraic equations:

M H(T) −H˜ τ + AK(T) + B(h)(T − Tw) + B(T4 − T ext4) + DQ = 0, (5)

where M,A,B(h),B and D are the n × n matrices with entries

mij= ∫ Vϕjϕidx, aij = ∫ V∇ϕj ⋅∇ϕidx, bij(h)= ∫ ΓNhϕjϕidΓN,bij = σε∫ ΓNϕjϕidΓN, dij= ∫ ΓMϕjϕidΓM.

More precisely, let L(l,k) and K(s,t) be the global indices of any two mesh points, where k and t indicate the node number in z-direction and l,s in xy-plane, respectively. By definition of mass matrix M we have

mLK = ∫ VϕLϕKdx = ∫ 0LZ ϕk(z)ϕt(z)dz∫ Ωhϕl(x,y)ϕs(x,y)dxdy = mktzm lsΩ.

Therefore matrix M is equal to Kronecker product of standard 1D mass matrix Mz and 2D mass matrix MΩ. Similarly we can define for the stiffness matrix A the entries

aLK = mktza lsΩ + a ktzm lsΩ

and for matrices B,B(h) and D

bLK= σεmktzm lsb, bLK(h)= hkmktzm lsb, dLK= mktzm lsb.

Temperature dependent enthalpy H and Kirchoff transform K are defined as piecewise linear and continuous functions. For the solution of system (5) the modified SOR-method [3] for points in V is used and for nonlinear part on ΓN Newton-Raphson method is used with one inner iteration.

We note that to implement our method it is sufficient to construct only 2D- and 1D-matrices. Therefore the computational efficiency of our model is very good. Also the memory allocation requirements are not so strong than in the case of ordinary 3D-brick elements.

4. Control of secondary cooling process

The secondary cooling region is divided into cooling zones (see figure 3) in which the values of heat transfer coefficients, hi , can be controlled separately in each cooling zone.

PIC

Figure 3: Schematic presentation of 3D-cooling zones for quarter of the slab

Two different control methods are used for optimizing the secondary cooling process on the boundary of the steel slab, namely PID and optimal control method.

4.1. PID control. In PID control the water cooling in each secondary cooling zone is done separately and indepentendly of each other. Heat transfer coefficient in each cooling zone i are calculated using the equation

hi(t) = P ΔTi(t) + I ∫ 0tΔT i(s)ds + DdΔTi(t) dt , (6)

where ΔTi(t) = Ti(t) − Titar(t) is the difference between calculated temperature T and target temperature Ttar at time t and P,I,D are experimentally known tuning parameters. After discretization of integral and derivative we obtain the recursive form [1]

hik = h ik−1 + P1 + D τ ΔTik −1 − τI + 2D τ ΔTik−1 + D τ ΔTik−2, (7)

with constant time step τ.

4.2. Optimal control. In optimal control method our aim is to minimize a cost function which is constructed by means of metallurgical cooling criteria. We can formulate our optimal control problem in the following way: Find h∗ ∈ U ad such that

J(T(h∗)) = min h∈UadJ(T(h)), (8)

where Uad = {h(t)∣hmin ≤ h(t) ≤ hmax} and T(h) is the solution of the state equation (5).

Our cost function, J, has the form

J(T) = 1 2 ∫ LMLZ (T − Ttar)2dz.

To solve the optimal control problem (8) we use the following gradient method. For given initial guess h0 ∈ U ad we find hi+1 ∈ U ad for i = 0, 1,… by the formula

hi+1 = hi + ρ ipi, (9)

where pi is the antigradient of the cost function I(h) = J(T(h)) at the point hi and ρi is an iterative parameter. To find I′(h) we use the method of Lagrange multipliers. Namely, we define Lagrange function

L(T,λ,h) = (Ψ(T,h),λ) + J˜(T), (10)

where

Ψ(T,h) := M H(T) −H˜ τ + AK(T) + B(h)(T − Tw) + B(T4 − T a4) + DQ

and

J˜(T) = 1 2 ∑ i=1nzl˜ i(Ti − Titar)2

is the approximation of the cost function J and by (⋅, ⋅) we mean inner product in Rn. We denote the length of the i-th element by l ˜ i. Differentiation with respect to variable T gives

1 τM(H′(T)⋅∂T)+A(K′(T)⋅∂T)+B(h)⋅∂T +B(4T3 ⋅∂T),λ+J˜′(T)⋅∂T = 0

and we derive a system of linear algebraic equations

1 τH′Mλ + K′Aλ + B(h)λ + 4T3Bλ = −J˜′. (11)

After solving the adjoint state problem (11) the descent direction

pj = 0,(hj = 0 ∧ Ij′ > 0) ∨ (h j = hmax ∧ Ij′ < 0) −Ij′,otherwise,

I′(h) = ∂Ψ ∂n (T,h),λ + ∂J˜ ∂h (T) = (B(1)(T − Tw),λ).

is found.

Choose a h0 ∈ U ad and natural number n. Given hi ∈ U ad we calculate the descent vector pi and the upper bound for the iterative parameter ρi :

ρmax = max{ρ ≥ 0 : hi + ρpi ∈ U ad}

The optimal step size in the equation (9) is still unknown. We use the following algorithm to determine ρopt:

Step 1:
ρ0 = ρmax 2n+1

If J(T(h(ρ0))) < J(T(h(0))) GOTO Step 2;

Else ρopt = 0 END;

Step 2:

For i = 1, 2,...,n

ρi = (2i+1 − 1)ρ 0

ρopt = ρi

If J(T(h(ρi))) < J(T(h(ρi−1))) CONTINUE;

Else ρopt = ρi−1 END;

Our aim is to use the optimal control method in dynamic casting process. At each iteration the calculation of the solution of the discretized state problem (5) is needed, which is very time consuming. Therefore, we use n = 7 in our numerical calculations.

5. Numerical example

Here we present some numerical results for solving continuous casting problem. The secondary cooling region is divided into eight cooling zones. The essential input parameters are presented in the table 1.

  
  
  
Half width of the slab Lx = 0.675m
Half thickness of the slab Ly = 0.105m
Length of the slab Lz = 32m
Length of the mould LM = 0.9m
Average cooling in mould Q̂ = 1229kWm−2
Latent heat L = 261kJkg−1
Density ρ = 7312kgm−3
Solidus Ts = 1763K
Liquidus Tl = 1793K
Incoming steel temperatureT0 = 1808K
Water temperature Tw = 293K
External temperature Text = 293K
Emissivity ε = 0.8 ⋅ 10−3
Stefan-Boltzmann constantδ = 5.67 × 10−8Wm−2K−4
Radius R = 0.005m
Casting speed v = 1.4m∕min


Table 1: Input parameters for continuous casting problem.

In our first numerical test we compare our 3D-model with existing 2D continuous casting model used at Rautaruukki steel factory in Finland. The cooling is assumed to be constant along each cooling zone during the casting process. Thus, the values of heat transfer coefficients hi are also constant. The values of hi are presented in the table 2 and the calculated core and surface temperatures are shown in the figure 4.

zonelength [m]h[kWm−2K−1]
1 0.47 0.4950
2 1.46 0.3708
3 1.62 0.3257
4 1.90 0.2783
5 3.84 0.1918
6 5.76 0.1175
7 5.18 0.1233
8 8.20 0.0995


Table 2: Zone lengths and values of h in each cooling zones.

PIC PIC

Figure 4: Calculated core and surface temperatures for 2D-model (on the left) and 3D-model (on the right) for piecewise constant h.

We are also interested in calculated lengths of liquid pool (longitudinal length of the liquidus isotherm in the centre of the slab), ’stable’ liquid pool ( T l − 2K-isotherm) and mushy pool (Ts-isotherm). Calculated pool lengths for 2D- and 3D-methods are shown in the table 3.

2D-method3D-method
Liquid pool 11.20m 10.55m
Stable liquid pool 16.88m 16.76m
Mushy pool 22.42m 22.63m


Table 3: Calculated pool lengths.

Our results show that calculated temperature fields of the two models are very similar but a difference in calculated lengths of liquid pools appear. From the control point of view the information about temperature field and the pool lengths is important because internal cracks of final product arise during solidification process. Moreover, soft reduction technology requires exact information about liquid and mushy pool lengths.

PIC PIC

Figure 5: Calculated temperatures using 3D-model with PID-control. Left graph shows target and calculated surface temperatures on the centreline of wide face. Right graph shows calculated corner temperatures.

In our second numerical test the PID and optimal control algorithms were also tested to optimize the secondary cooling i.e. heat transfer coefficients. The results are shown in figures 5 and 6. We can see from the figure 5 that PID control algorithm adjusts water cooling in the way that the difference between target and calculated temperatures at the end of each cooling zone is minimized. The optimal control method on the other hand minimizes the cost function J. From the figure 6 we see that the difference between target and calculated surface temperature is minimized on the whole length of the slab.

PIC PIC

Figure 6: Calculated temperatures using 3D-model with optimal control. Left graph shows target and calculated surface temperatures on the centreline of wide face. Right graph shows calculated corner temperatures.

References

[1]   S. Bouhouche, Contribution to quality and process optimisation in continuous casting using mathematical modelling, Doctoral Dissertation, der Technischen Universita¨t Bergakademie Freiberg, Germany, (2002).

[2]   J. Jr. Douglas and T.F. Russel, Numerical methods for convection-dominated diffusion problem based on combining the method of characteristic with finite element or finite difference procedures, SIAM J. Numer. Anal., V. 19., (1982), pp. 871–885.

[3]   C. M. Elliot and J. R. Ockendon, Weak and variational methods for moving boundary problems. Pitman Advanced Publishing Program, Boston, (1982).

[4]   B. Filipic, B. Sarler, Evolving parameter settings for continuous casting of steel, in Proceedings of the 6th European Congress on Intelligent Techniques & Soft Computing, Aachen, Germany, (1998), pp. 444–449.

[5]   M. Jauhola, E. Kivela¨, J. Konttinen, E. Laitinen, S. Louhenkilpi, Dynamic secondary cooling model for continuous casting machine, in Proceedings of ATS 93 Steelmaking, Paris, France, (1994), pp. 110–111.

[6]   E. Laitinen, On the simulation and control of the continuous casting of steel, Doctoral Dissertation, University of Jyva¨skyla¨, Finland, (1989).

[7]   E. Laitinen, A. Lapin, J. Pieska¨, Mesh approximation and iterative solution of the continuous casting problem, ENUMATH 99 - Proceeding of the 3rd European conference on Numerical Mathematics and Advanced Applications, ed. by P. Neittaanmki, T. Tiihonen and P. Tarvainen, World Scientific, Singapore, (2000), pp. 601–617.

[8]   E. Laitinen, J. Pieska¨, Comparison of upwind and characteristic schemes for solving multiphase diffusion-convection equation, Computer assisted Mechanics and Engineering Sciences, 7, Institute of Fundamental Technological Research, Polish Academy of Sciences, Warsaw, Poland, (2000), pp. 421–426.

[9]   S. Louhenkilpi, E. Laitinen, R. Nieminen, Real-time simulation of heat transfer in continuous casting, Metall. Trans., Vol. B24, (1993), pp. 685–693.

[10]   J.F Rodrigues, F. Yi, On a two-phase continuous casting Stefan problem with nonlinear flux, Euro. J. App. Math., Vol. 1, (1990), pp. 259–278.

[11]   C.A Santos, J.A. Spim Jr., M.C.F. Ierardi, A. Garcia, The use of artificial intelligence technique for the optimisation of process parameters used in the continuous casting of steel, Applied Math. Modelling., 26, (2002), pp. 1077–1092.

KAZAN STATE UNIVERSITY, DEPARTMENT OF APPLIED MATHEMATICS, 18, KREMLEVSKAYA ST., KAZAN, RUSSIA, 420008

E-mail address: rdautov@ksu.ru

KAZAN STATE UNIVERSITY, DEPARTMENT OF APPLIED MATHEMATICS, 18, KREMLEVSKAYA ST., KAZAN, RUSSIA, 420008

E-mail address: rkadyrov@ksu.ru

UNIVERSITY OF OULU, FINLAND.

E-mail address: ejl@sun3oulu.fi

KAZAN STATE UNIVERSITY, DEPARTMENT OF APPLIED MATHEMATICS, 18, KREMLEVSKAYA ST., KAZAN, RUSSIA, 420008

E-mail address: alapin@ksu.ru

UNIVERSITY OF OULU, FINLAND.

E-mail address: jpieska@sun3.oulu.fi

UNIVERSITY OF OULU, FINLAND.

E-mail address: valtteri@mail.student.oulu.fi

Received October 15, 2003