Lobachevskii Journal of Mathematics Vol. 13, 2003, 3 – 13 \begin{matrix}c \\ ◯ \end{matrix}R. Dautov, R. Kadyrov, E. Laitinen, A. Lapin, J. Pieskä, V. Toivonen

R. Dautov, R. Kadyrov, E. Laitinen, A. Lapin, J. Pieskä, 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.

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 realtime 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.

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

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 Ω ⊂ R^{2} be a rectangular domain [0, L_{X}] × [0, L_{Y}] and V = Ω × [0, L_{Z}] be a 3D domain. The boundary Γ = ∂V consists of the parts:

Γ_{0} = Ω × {0},
Γ_{N} = {(x, y) ∈ ∂Ω : x = 0 ∨ y = 0} × [L_{M}, L_{Z}],
Γ_{S} = {(x, y) ∈ ∂Ω : x ≠ 0 ∧ y ≠ 0} × [0, L_{Z}) ∪ Ω × {L_{Z}},
Γ_{M} = {(x, y) ∈ ∂Ω : x = 0 ∨ y = 0} × [0, L_{M}],

where L_{M} is the length of the mould, L_{Z} is the length of the strand and L_{X}, L_{Y} are the width and thickness of 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) = ρ( ∫_{0}^{T} c(ξ)dξ + L(1 − f_{s}(T))), K(T) = ∫_{0}^{T} k(ξ)dξ,

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

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

\{ \begin{matrix}\frac{∂H(T)}{∂t} + v \frac{∂H(T)}{∂z} − ΔK(T) = 0 & in V × (0, t_{f}], \\ T = T_{0} & on Γ_{0} × (0, t_{f}], \\ \frac{∂K(T)}{∂n} + h(T − T_{w}) + σϵ(T^{4} − T_{ext}^{4}) = 0 & on Γ_{N} × (0, t_{f}], \\ \frac{∂K(T)}{∂n} = 0 & on Γ_{S} × (0, t_{f}], \\ \frac{∂K(T)}{∂n} = Q & on Γ_{M} × (0, t_{f}], \\ T (x, y, z, 0) = T^{0} & in V. \end{matrix}  (1)

Figure 1. The calculation domain.

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}, T_{ext} are known temperatures. The σ is the Stefan-Boltzmann constant and ϵ is the emissivity. The cooling capacity Q in the mould is known constant and t_{f} is the simulation time.

3. Mesh approximation

We partition Ω into a set of quadrilateral finite elements as shown on figure 2. Step sizes are smaller near left and down boundaries of Ω, and all quadrilaterals

Figure 2. Partition of Ω on quadrilateral finite elements.

are rectangular except those elements near the round corner. The 1D-domain [0, L_{Z}] is divided by n_{z} 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, L_{X}], 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 (\frac{∂}{∂t} + v(t)\frac{∂}{∂z} ) H by using the char-

acteristics 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−τ}^{t} v(ξ)dξ and ap-

proximate this term by:

(\frac{∂}{∂t} + v(t)\frac{∂}{∂z} ) H ≈ \frac{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

\begin{matrix}∫_{V} (H(T) − H̃)/τ ηdx ∫_{V} ∇K(T) · ∇ηdx + ∫_{Γ_{N}} \{h(T − T_{w}) \\ +σϵ(T^{4} − T_{ext}^{4})\} ηdΓ_{N} + ∫_{Γ_{M}} Q ηdΓ_{M} = 0 \end{matrix}  (3)

for all trial functions η ∈ V^{0} = {η ∈ H^{1}(V) : η = 0 on Γ_{0}}. Here H^{1}(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 H^{1}(V), V_{h}^{T_{0}} = {v_{h} ∈ V_{h}(V) : v_{h} = T_{0h} on Γ_{0}} with T_{0h} being the V_{h}-interpolant of T_{0} and obviously defined subspace V_{h}^{0}. Let further φ_{1}(x), …, φ_{n}(x), be the basis functions in V_{h}. Thus, the function

T_{h}(x) = ∑_{i=1}^{n} T_{i}φ_{i}(x)

is the finite element approximation of T (x), T = (T_{1}, …, T_{n})^{t} is the vector of the unknown nodal values of the function T_{h}(x). Below the notations are: H(T) = (H(T_{1}), …, H(T_{n}))^{t}, K(T) = (K(T_{1}), …, K(T_{n}))^{t} and

f_{h}(T_{h}) = ∑_{i=1}^{n} f(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

\begin{matrix}∫_{V} \frac{H_{h}(T_{h}) − H̃_{h}}{τ} ηdx + ∫_{V} ∇K_{h}(T_{h}) · ∇ηdx + ∫_{Γ_{N}} \{h(T_{h} − T_{w h}) \\ +σϵ((T^{4})_{h} − (T_{ext}^{4})_{h})\} ηdΓ_{N} + ∫_{Γ_{M}} Q_{h} ηdΓ_{M} = 0 \end{matrix}  (4)

for all trial functions η ∈ V_{h}^{0}. Discrete problem (4) is equivalent to the system of nonlinear algebraic equations:

M \frac{H(T) − H̃}{τ} + AK(T) + B(h)(T − T_{w}) + B(T^{4} − T_{ext}^{4}) + DQ = 0, (5)

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

m_{ij} = ∫_{V} φ_{j}φ_{i}dx,  a_{ij} = ∫_{V} ∇φ_{j} · ∇φ_{i}dx,
b_{ij}(h) = ∫_{Γ_{N}} hφ_{j}φ_{i}dΓ_{N}, b_{ij} = σϵ ∫_{Γ_{N}} φ_{j}φ_{i}dΓ_{N},
d_{ij} = ∫_{Γ_{M}} φ_{j}φ_{i}dΓ_{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 xyplane, respectively. By definition of mass matrix M we have

m_{LK} = ∫_{V} φ_{L}φ_{K}dx = (∫_{0}^{L_{Z}} φ_{k}(z)φ_{t}(z)dz) ∫_{Ω_{h}} φ_{l}(x, y)φ_{s}(x, y)dxdy = m_{kt}^{z}m_{ls}^{Ω}.

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

a_{LK} = m_{kt}^{z}a_{ls}^{Ω} + a_{kt}^{z}m_{ls}^{Ω}

and for matrices B, B(h) and D

\begin{aligned}b_{LK} & = σϵm_{kt}^{z}m_{ls}^{b}, \\ b_{LK}(h) & = h_{k}m_{kt}^{z}m_{ls}^{b}, \\ d_{LK} & = m_{kt}^{z}m_{ls}^{b}. \end{aligned}

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 2Dand 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, h_{i}, can be controlled separately in each cooling zone. 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

h_{i}(t) = P [ΔT_{i}(t) + I ∫_{0}^{t} ΔT_{i}(s)ds + D\frac{dΔT_{i}(t)}{dt}] ,  (6)

where ΔT_{i}(t) = T_{i}(t)−T_{i}^{tar}(t) is the difference between calculated temperature T and target temperature T^{tar} at time t and P, I, D are experimentally known Figure 3. Schematic presentation of 3D-cooling zones for quarter of the slab

tuning parameters. After discretization of integral and derivative we obtain the recursive form [1]

h_{i}^{k} = h_{i}^{k−1} + P [(1 + \frac{D}{τ})ΔT_{i}^{k} − (1 − τI + \frac{2D}{τ})ΔT_{i}^{k−1} + \frac{D}{τ}ΔT_{i}^{k−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∈U_{ad}} J(T (h)),  (8)

where U_{ad} = {h(t)| h_{min} ≤ h(t) ≤ h_{max}} and T (h) is the solution of the state equation (5). Our cost function, J, has the form

J(T) = \frac{1}{2} ∫_{L_{M}}^{L_{Z}}(T − T^{tar})^{2}dz.

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

h^{i+1} = h^{i} + ρ_{i} p^{i},  (9)

where p^{i} is the antigradient of the cost function I(h) = J(T (h)) at the point h^{i} 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 \frac{H(T) − H̃}{τ} + AK(T) + B(h)(T − T_{w}) + B(T^{4} − T_{a}^{4}) + DQ

and

J̃(T) = \frac{1}{2} ∑_{i=1}^{nz} l̃_{i}(T_{i} − T_{i}^{tar})^{2}

is the approximation of the cost function J and by (·, ·) we mean inner product in R^{n}. We denote the length of the i-th element by l̃_{i}. Differentiation with respect to variable T gives

(\frac{1}{τ}M(H′(T)·∂T)+A(K′(T)·∂T)+B(h)·∂T +B(4T^{3}·∂T), λ)+J̃′(T)·∂T = 0

and we derive a system of linear algebraic equations

\frac{1}{τ}H′Mλ + K′Aλ + B(h)λ + 4T^{3}Bλ = −J̃′.  (11)

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

p_{j} = \begin{cases}0, (h_{j} = 0 ∧ I′_{j} > 0) ∨ (h_{j} = h_{max} ∧ I′_{j} < 0) \\ −I′_{j}, otherwise, \end{cases}

I′(h) = (\frac{∂Ψ}{∂n}(T, h), λ) + \frac{∂J̃}{∂h}(T) = (B(1)(T − T_{w}), λ).

is found.

Choose a h^{0} ∈ U_{ad} and natural number n. Given h^{i} ∈ U_{ad} we calculate the descent vector p^{i} and the upper bound for the iterative parameter ρ_{i}:

ρ_{max} = max{ρ ≥ 0 : h^{i} + ρp^{i} ∈ U_{ad}}

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

Step 1: ρ_{0} = \frac{ρ_{max}}{2^{n+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} = (2^{i+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.

Half width of the slab  L_{x} = 0.675m
Half thickness of the slab L_{y} = 0.105m
Length of the slab  L_{z} = 32m
Length of the mould  L_{M} = 0.9m
Average cooling in mould Q̂ = 1229kW m^{−2}
Latent heat  L = 261kJkg^{−1}
Density  ρ = 7312kgm^{−3}
Solidus  T_{s} = 1763K
Liquidus  T_{l} = 1793K
Incoming steel temperature T_{0} = 1808K
Water temperature  T_{w} = 293K
External temperature  T_{ext} = 293K
Emissivity  ϵ = 0.8 · 10^{−3}
Stefan-Boltzmann constant δ = 5.67 × 10^{−8}W m^{−2}K^{−4}
Radius  R = 0.005m
Casting speed  v = 1.4m/min

Table 1. Input parameters for continuous casting problem.

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.

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 h_{i} are also constant. The values of h_{i} are presented in the table 2 and the calculated core and surface temperatures are shown in the figure 4.

zone length [m] h[kW m^{−2}K^{−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.

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 (T_{s}-isotherm). Calculated pool lengths for 2D- and 3D-methods are shown in the table 3.

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

2D-method 3D-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.

Figure 5. Calculated temperatures using 3D-model with PIDcontrol. 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.

12 R. DAUTOV, R. KADYROV, E. LAITINEN ET AL.

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.

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 Universitä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. Kivelä, 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 Jyväskylä, Finland, (1989).

[7] E. Laitinen, A. Lapin, J. Pieskä, 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. Pieskä, 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.

ON 3D DYNAMIC CONTROL 13

[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