Modeling of the continuous casting process of steel via phase-field transition system. Fractional steps method

1.   Introduction

We consider a material in a region $\Omega \subset {{\mathbb{R}}^{n}}$, $n\le 3$, which may be in either of two phases, e.g., solid and liquid (see Figure 1). Let us denote by

the reduced temperature distribution, where $\theta(t, x)$ represent the temperature of the material and $\theta_M$ is the melting temperature (the temperature at which solid and liquid may coexist in equilibrium, separated by an interface).

In the following we will describe the framework of our problem. So, let us consider the interface as a continuous region, more vast (in which the liquid can coexist with the solid) and of finite thickness, in which the change of phase occurring continuously.

The following nonlinear parabolic system

with the non-homogeneous Cauchy-Neumann boundary conditions

and with the initial conditions

represents the mathematical model called the phase field transition system, introduced by G. Caginalp (see ^[3] and references therein) to model the transition between the solid and liquid phase in melting/solidification process to a matter occupying a region $\Omega$, while:

● $u(t, x)$ - represents the reduced temperature distribution in $Q$;

● $\varphi(t, x)$ - is the phase function (the order parameter) used to distinguish between the states (phases) of material which occupies the region $\Omega$ at every time $t\in[0, T]$;

● $C_\rho = \rho V$; $\rho$ - the density, $V$ - the casting speed;

● $\ell, k, \alpha, \xi, h$ are physical parameters representing, respectively: the latent heat, the thermal conductivity, the relaxation time, the measure of the interface thickness, the heat transfer coefficient;

● $s_{_\xi} = {{m[S]_{\!E}}\over{2\sigma}}T_{\!E}$, a bounded and positive quantity, expressed by positive and bounded physical parameters: $m = \int\limits_{-1}^1\bigl(2F(s)\bigr)^{1\over 2}ds$, $F(s) = {1\over 4}(s^2-1)^2$, $[S]_{\!E}$ - the entropy difference between phases per volume, $\sigma$ - the interfacial tension, $T_{\!E}$ - the equilibrium melting temperature (see Caginalp & Chen ^[4]);

● $w_1(t, x)\in W^{1-{1\over 2p}, 2-{1\over p}}_p(\Sigma)$, $p\ge 2$, - is a given function: the temperature of the surrounding at $\partial\Omega$ for each time $t\in[0, T]$ (can also be interpreted as boundary control);

● $w_2(t, x)\in W^{1-{1\over 2p}, 2-{1\over p}}_p(\Sigma)$ - is a given function;

● $u_0, \varphi_0\in W^{2-{2\over p}}_p(\Omega)$, with $k{\partial\over\partial\nu}u_0 + hu_0 = w_1(0, x)$ and $\xi{\partial\over\partial\nu}\varphi_0 = w_2(0, x)$;

The model (1.1)–(1.3) represents a refinement of the classical Stefan problem (see ^[21,23,24]) in two phases by adding a new nonlinear parabolic equation, derived from the Euler-Lagrange equations for the free energy in Landau-Ginzburg field theory. This new mathematical model reflects more accurately the physical phenomenon of solidification, like: superheating, supercooling, etc.

Different other nonlinearities capable to come out the complexity of the physical phenomena (the effect of surface tension, separating zone of solid and liquid states, etc) have been proposed by several authors (see Cârjă, Miranville & Moroşanu ^[5], Kenmochi & Niezgódka ^[8], Miranville & Moroşanu ^[9], Moroşanu ^[16], Moroşanu & Motreanu ^[20], Penrose & Fife ^[22] and Temam ^[24]). The general nonlinear term in Moroşanu & Motreanu ^[20], is (possibly) non-convex and non-monotone and cover a large class of nonlinearities, including the known cases as well as other new relevant situations. Moreover, different types of boundary conditions on $\Sigma$ can be associated to (1.1) (see Moroşanu ^[16] for more details).

2.   Well-posedness of solutions to the nonlinear system (1.1)–(1.3)

In the present section we will investigate the solvability of the first boundary value problems of the form (1.1)–(1.3) in the class $W^{1, 2}_p(Q)$, $p\ge 2$.

The main result of this section establishes the dependence of the solution $u(t, x)$, $\varphi(t, x)$ in the nonlinear parabolic system (1.1)–(1.3) on the terms $w_1(t, x)$, $w_2(t, x)$ in the right-hand side of (1.2).

Theorem 2.1 Problem $\rm(1.1)–(1.3)$ has a unique solution $(u, \varphi)$ with $u\in W^{1, 2}_p(Q)$ and $\varphi\in W^{1, 2}_\nu(Q)$, where $\nu = min\{q, \mu\}$, $q\ge p\ge 2$. In addition $(u, \varphi)$ satisfies

where the constant $C$ depends on $|\Omega|$ $\rm($the measure of $\Omega$$\rm)$, $T$, $n$, $p$, $q$ and physical parameters, but is independent of $u$, $\varphi$, $w_1$ and $w_2$.

Moreover, given any number $M_d > 0$, if $(u_1, \varphi_1)$ and $(u_2, \varphi_2)$ are solutions to $\rm(1.1)–(1.3)$ for the same initial conditions, corresponding to the data $w^1_1, w_1^2$, $w^1_2, w^2_2\in$ $W^{1-{1\over 2p}, 2-{1\over p}}_p(\Sigma)$, respectively, such that $\|\varphi_1\|_{L^\nu(Q)}$, $\|\varphi_2\|_{L^\nu(Q)}\leq M_d$, then the estimate below holds

where the constant $C$ depends on $|\Omega|$, $T$, $M_d$, $n$, $p$, $q$ and physical parameters, but is independent of $u_1, u_2$, $\varphi_1, \varphi_2$, $w^1_1$, $w^2_1$, $w_2^1$ and $w^2_2$.

Proof. The basic tools in the analysis of the problem (1.1) (see ^[16] and references there in) are the Leray-Schauder degree theory, the $L^p$-theory of linear and quasi-linear parabolic equations, as well as the Lions and Peetre embedding Theorem, which ensures the existence of a continuous embedding $W^{1, 2}_p(Q)\subset L^\mu(Q)$, where the number $\mu$ is defined as follows

The proof of Theorem 2.1 was given in Moroşanu ^[16] noting that there formulation differs from this by certain physical parameters, which implies different values for the constant $C$ in (2.1) and (2.2). Moreover, corresponding to different boundary conditions (including nonlinear and nonhomogeneous boundary conditions), similar results were proved in Cârjă, Miranville & Moroşanu ^[5] and Miranville & Moroşanu ^[9].

Corollary 2.2 Under hypotheses $\rm H_0$ and $\rm H _2$ in ^[20] the problem $\rm(1.1)$ possesses a unique solution $(u, \varphi)\in W^{1, 2}_p(Q)\times W^{1, 2}_p(Q)$.

Proof. Let $w^1_1 = w^2_1 = w^1_2 = w^2_2 = w$ in the Theorem 2.1. Then (2.2) shows that the conclusion of the corollary is true.

3.   Approximating scheme

The aim of this section is to use the fractional steps method in order to approximate the solution of the system (1.1)–(1.3), whose uniqueness is guaranteed by Corollary 2.2. To do that, let us associate to the time-interval $[0, T]$ the equidistant grid of length $\varepsilon = \frac{T}{M}$, for any integer $M\ge 1$. Then, the following approximating scheme can be written in order to approximate the solution of the nonlinear boundary value problem (1.1)–(1.3):

where $Q^\varepsilon_i = [i\varepsilon, (i+1)\varepsilon]\times\Omega$, $\Sigma^\varepsilon_i = [i\varepsilon, (i+1)\varepsilon]\times\partial\Omega$, $i = 0, \cdots, M-1$, and $z(t, \varphi^\varepsilon_-(i\varepsilon))$ is the solution of the Cauchy problem

computed at $s = \varepsilon$, for $i = 0, \cdots, M\!-\!1$. Here $\varphi^\varepsilon_+(i\varepsilon) = \lim\limits_{t\downarrow i\varepsilon}\varphi^\varepsilon(t)$ and $\varphi^\varepsilon_-(i\varepsilon) = \lim\limits_{t\uparrow i\varepsilon}\varphi^\varepsilon(t)$.

For later use, we set: $W = L^2(0, T;H^1(\Omega))\cap W^{1, 2}([0, T];(H^1(\Omega))')$.

Definition 3.1. By weak solution to the problem (1.1)–(1.2) we mean a pair of functions $(u, \varphi)\in W\times W$ which satisfy (1.1)–(1.2) in the following sense

for all $(\psi, \zeta)\in L^2(0, T;H^1(\Omega))\times L^2(0, T;H^1(\Omega))$, together with the initial conditions (1.3). In (3.4) and (3.5) we have denoted by the symbol $\int\limits_Q$ the duality between $L^2(0, T;H^1(\Omega))$ and $L^2(0, T;(H^1(\Omega))')$.

The main result of this section is the following

Theorem 3.2. Assume that $u_0, \varphi_0\in W^{2-{2\over p}}_p(\Omega)$ with $k{\partial\over\partial\nu}u_0 + hu_0 = w_1(0, x)$ and $\xi{\partial\over\partial\nu}\varphi_0 = w_2(0, x)$. Let $(u^\varepsilon, \varphi^\varepsilon)$ be the solution of the approximating scheme $\rm(3.1)–(3.3)$. Then, for $\varepsilon\rightarrow 0$, we have

for any $t\in [0, T]$, where $u^*, \varphi^*\in W^{1, 2}([0, T];L^2(\Omega))\cap L(0, T;H^2(\Omega))$ is the weak solution to the problem $\rm(1.1)–(1.3)$.

Proof.(see ^[2]) The proof is based on compactness methods. As a matter of fact it turns out from Theorem 3.2 that if $u_0, \varphi_0\in L^2(\Omega)$, then the weak solution $(u^*(t), \varphi^*(t))$ of the system (1.1)–(1.3) is a strong solution, i.e., it is absolutely continuous in $t$ on $[0, T]$ and satisfies a.e. the system (1.1)–(1.3). So Theorem 3.2 can be also viewed as a constructive way to prove the existence in (1.1)–(1.3).

The result in Theorem 3.2 remains true by replacing the boundary condition (1.2) with $k{\partial\over\partial\nu}u+hu = w(t)g(x)$ and $\xi{\partial\over\partial\nu}\varphi = 0$ (see ^{[6,7,10,11,12,13,14,15,16,17,18,19]}).

The Cauchy problem (3.3) has the solution

and then the general algorithm to compute the approximate solution by means of fractional steps method consist in the following sequence ($i$ denotes the time level)

Begin algfrac

  $i: = 0\ \rightarrow u^{\varepsilon, 0} = u_0, \ \varphi^{\varepsilon, 0} = \varphi_0$ from the initial conditions (1.3);

  ${\tt For}$ $i: = 0$ ${\tt to}$ $M-1$ ${\tt do}$

    Compute $z(\varepsilon, \varphi^\varepsilon_-(i\varepsilon, x))$ from (3.7);

    $\varphi^\varepsilon_+: = z(\varepsilon, \varphi^\varepsilon_-(i\varepsilon, x))$;

     Compute $\varphi^{\varepsilon, i+1}, u^{\varepsilon, i+1}$ solving the linear system (3.1)-(3.2);

  ${\tt End-for}$;

${\tt End}$.

A comparison between the fractional steps method and the standard iterative Newton method can be found in Moroşanu ^[16].

4.   Numerical model. The 2D case

The finite element method (fem in short) is a general method for approximating the solution of boundary value problems for partial differential equations. This method derives from the Ritz (or Gelerkin) method (see Axelson & Barker ^[1]), characteristic for the finite element method being the chose of the finite dimensional space, namely, the span of a set of finite element basis functions.

The steps in solving a boundary value problem using fem are:

P0. (D) The direct formulation of the problem;

P1. (V) A variational (weak) formulation for problem (D);

P2.     The construction of a finite element mesh (triangulation);

P3.     The construction of the finite dimensional space of test

    function (called finite element basis functions);

P4. (V$_{nn}$) A discrete analogous of (V);

P5.     Assembly of the system of linear equations;

P6.     Solve the system in P5.

The finite element method is used in the sequel in order to deduce the discrete state equations. A conceptual numerical algorithm of fractional step type is then formulated to approximate the weak solution corresponding to (3.1)–(3.2), that is:

together with the initial conditions

By $(\cdot, \cdot)$ we have denoted the scalar product in $L^2(\Omega)$.

Let $\varepsilon = T/M$ be the time step size. We assume that $\Omega\subset {\mathbb{R}}^2$ is a polygonal domain. Let $T_r$ be the triangulation (mesh) over $\Omega$ and $\bar\Omega = \cup_{K\in T_r}K$, and let $N_j = (x_k, y_l), j = \overline{1, nn}$, be the nodes associated to $T_r$. Denoting by $V_{nn}$ the corresponding finite element space to $T_r$, then the basic functions $\{b_j\}_{j = 1}^{nn}$ of $V_{nn}$ are defined by

and so

For $i = \overline{1, M}$, we denote by $u^i$ and $\varphi^i$ the $V_{nn}$ interpolant of $u^\varepsilon$ and $\varphi^\varepsilon$, respectively. Then $u^i, \varphi^i\in V_{nn}$ and

where $u^i_l = u^\varepsilon(t_i, N_l)$, $\varphi^i_l = \varphi^\varepsilon(t_i, N_l)$, $i = \overline{1, M}$, $l = \overline{1, nn}$ are the unknowns to be computed.

Using in addition an implicit (backward) finite difference scheme in time, we introduce now the discrete equations corresponding to (4.1)–(4.2) as follows (see ^{[7,12,13,14,15,16]}) for more explanations)

where $u^i_l$ and $\varphi^i_l$, $l = \overline{1, nn}$, are the vectors of unknowns for time level $i$.

From the initial conditions (1.3) we have

and then from (4.6) we get (see (4.3)–(4.4))

The numerical algorithm to compute the approximate solution by fractional steps method can be obtained from the following sequence (again, $i$ denotes the time level)

${\tt Begin}$ algfrac_fem

  $i: = 0\ \rightarrow$ Compute $u^0_l, \varphi^0_l$, $l = \overline{1, nn}$ from (4.7);

  Choose $w^{i, l}_1 = w_1(t_i, N_l), \ w^{i, l}_2 = w_2(t_i, N_l), \ N_l\in\partial\Omega, \ i = \overline{0, M\!-\!1}, l = \overline{1, nn}$;

  ${\tt For}$ $i: = 1$ ${\tt to}$ $M$ ${\tt do}$

    Compute $z_l = z(\cdot, N_l),$ $l = \overline{1, nn}$ from (3.7);

    $\varphi^{i-1} : = z_l,$ $l = \overline{1, nn}$;

   Compute $u^i_l,$ $\varphi^i_l$, $l = \overline{1, nn}$, solving the linear system (4.5);

  ${\tt End-for}$;

${\tt End}$.

The convergence result established by Theorem 2.1 guaranty that the approximate solution computed by the conceptual algorithm algfrac_fem is in fact the approximate solution of the nonlinear parabolic system (1.1)–(1.3).

5.   Industrial implementation

The aim of this section is to present an industrial implementation of conceptual algorithm algfrac_fem established in the preview section (in fact an implementation of the numerical model stated by the linear system (4.5) (see P5)).

From the thermo-kinetics point of view, the solidification and cooling, as well as the simultaneous heating in a continuous casting process of steel, represents a very complex problem of non-stationary heat and mass transfer (see ^[25]). To solve at present such a problem it is impossible without numerical models of the temperature field and computer technology. In the sequel we will present in short a continuous casting machine, including the key phenomena of interest and a new numerical model that we have used in the settlement of the problem mentioned above.

The continuous casting process in the metallurgy. In a modern steel casting machine (its essential features are illustrated in Figure 2), the molten metal is tapped from a ladle into a copper mold (crystallizer). Here, the water-cooled walls of the mold (the primary cooling zone) extract heat what leads to solidify a shell that contains the liquid pool. Below the mold (the secondary cooling zone), the product is supported by rollers and is cooled down by water sprays that extract heat from the surface, and, eventually, the core becomes fully solid when the metallurgical length increase at 12÷13m. After the end of secondary cooling zone the product is cooled only by radiation (traiber). Finally, the continuous-casting product must be cut into the optimum lengths (cutting) to achive a maximum yield of metal.

The application of the numerical model (4.5) to the continuous casting process, requires experimental research and measurements of operational parameters at MTC2 from Mittal Steel S.A. Galaţi, as well as laboratory research. So, the most important input data in order to do this are (in round bracket we have written the value used by our Matlab program to do the numerical simulations):

● the casting speed ($V = 12.5$ mm/s);

● physical parameters:

   $\circ$ the density ($\rho = 7850$ kg/m$^3$),

   $\circ$ the latent heat ($\ell = 65.28kcal/kg$),

   $\circ$ the relaxation time ($\alpha = 1.0e+2$),

   $\circ$ the length of separating zone ($\xi = .5$),

   $\circ$ the coefficients of heat transfer ($h = 32.012$),

   $\circ$ $T = 44s$;

● the boundary conditions: ($w_1(t_i, N_l), \ N_l\in\partial\Omega, \ \ i = \overline{0, M\!-\!1}, l = \overline{1, nn}$, in the primary cooling zone (see Figure 4)), ($w_1(t_i, N_l), \ N_l\in\partial\Omega, \ \ i = \overline{0, M\!-\!1}, l = \overline{1, nn}$, in the secondary cooling zone (see Figure 8)) and $w_2(t_i, N_l) = 0, \ N_l\in\partial\Omega, \ \ i = \overline{0, M\!-\!1}, l = \overline{1, nn}$;

● dimensions of cristallizer (650 x 1900 x 220), in mm;

● the casting temperature ($u_0 = 1530 ^0$C);

● the termal conductivity $k(u)$:

k(u) = [20 100 200 300 400 500 600 700 800 850 900 1000 1100 1200 1600;

    1.43e-5 1.42e-5 1.42e-5 1.42e-5 1.42e-5 9.5e-6 9.5e-6 9.5e-6 8.3e-6...

    8.3e-6 8.3e-6 7.8e-6 7.8e-6 7.4e-6 7.4e-6].

Numerical experiments. In Figure 3 it can be seen the number of nodes associated to the mesh $T_r$ in the $x_1$ and $x_2$ – axis directions of one half of a rectangular profile. Considering the symmetrical heat removal from the continuous casting (CC) according to the vertical symmetry axis of the rectangular profile, only a half of the cross-section is used in the computation program.

The numerical model (4.5) uses the temperatures $w(t, x, y), t\in [0, T]$, $(x, y)\in \partial\Omega$ measured by the thermocouples; the values are illustrated in the Figure 4.

Figures 5–7 represents the approximate solution $u^i$, $\varphi^i$ (see (4.3), (4.4)), corresponding to different moments of time ($i = 1$, $i = 5$ and $i = M$).

A close examination of the Figures 5–7 tell us the dimension of the solid and liquid zone resulting by runing the Matlab computation program developed on the basis of the conceptual algorithm algfrac_fem.

Moreover (see ^[17]), the shape of the graphs shows the stability and accuracy of the numerical results obtained by implementing the fractional steps method.

The most interesting aspect that we can observe analysing the Figures 6–7 are the supercooling and superheating phenomenon.

6.   Conclusions

● The solidification model that we have considered in this work consist in a nonlinear system of two parabolic differential equations ^[4]. This new mathematical description of the real phenomenon reflects more accurately the physical aspects, like: superheating, supercooling (see Figures 7 and 11, for example), the effects of surface tension, separating zone of solid and liquid states, etc.

● From numerical point of view, the main difficulty in treating the phase field transition system (1.1) is due to the presence of the nonlinear equation corresponding to phase function $\varphi$. Thus it is intensely motivated the work in finding more efficient algorithms in order to compute numerically the solution of such system. A scheme of fractional steps type is considered in this sense. This numerical method avoids the iterative process required by the classical methods (e.g., Newton's type approaches) in passing from a time level to another. Numerical tests show that the fractional steps method is faster (CPU-time spent is very small) and the stability and accuracy are higher (see ^[16,17]) than the Newton's methods.

● The distribution of the temperature and the thickness of the solidifying shell, calculated with the numerical model (4.5) obtained following this technique, show that it really is (see Figures 5–7). New fundamental material properties can also be extracted by analysing the implementation of the numerical model 4.5 (see Figure 7).

This model is able to simulate the temperature field of a CCM (Continuous Casting Machine) as a whole or any of its parts. In addition, the program elaborated may be used for different slab profiles. The industrial implementation of the numerical model enable us the analysis of the temperature field of the slab when it passes through primary, secondary and traiber zone. The Figures 9–11 display the calculated isotherms in the secondary cooling zone, while the Figure 12 display the temperature curves on the mobile part, corresponding to the final time level $t^M$.

● In order to refer the continuous casting process, the following boundary conditions was considered: ${\partial\over\partial\nu}u+hu = w_1(t, x, y)$, where $h$ is the heat transfer coefficient and the given function $w_1(t, x, y)$ represents the temperature of the surrounding at $t\in [0, T]$ and $(x, y)\in\partial\Omega$, $\Omega\subset {\mathbb{R}}^2$.

Generally, the numerical method considered here can be used to approximate the solution of a nonlinear parabolic equation (system) containing a general nonlinear part.

● The numerical solution of the phase field transition system of solidification, approximated by this numerical scheme, can be considered as an admissible one for the corresponding boundary optimal control problem (see ^[16]), formulated in order to study the optimization of the continuous casting process. The numerical results presented in Figures 9-11 illustrate the accuracy of the numerical model (the influence of the density of the net: M = 20, nkx1 = 10, nky1 = 10; M = 20, nkx1 = 20, nky1 = 10; M = 40, nkx1 = 20, nky1 = 10). So, it is strongly motivated to investigate in further the numerical stability of this new numerical model taking into account all parameters (steel casting parameters, physical parameters, net parameters, etc.)

A detailed discussions on the errors produced by the fractional steps method, illustrate the influence of time and space parameters as well as of all physical parameters (see ^[17]).

Conflict of interest

The author declares that there is no conflicts of interest in this paper.