UDK 621.74.047:519.61/.64 ISSN 1580-2949 Original scientific article/Izvirni znanstveni članek MTAEC9, 48(2)281(2014) SIMULATION OF CONTINUOUS CASTING OF STEEL UNDER THE INFLUENCE OF MAGNETIC FIELD USING THE LOCAL-RADIAL BASIS-FUNCTION COLLOCATION METHOD SIMULACIJA KONTINUIRNEGA ULIVANJA JEKLA POD VPLIVOM MAGNETNEGA POLJA NA PODLAGI METODE KOLOKACIJE Z RADIALNIMI BAZNIMI FUNKCIJAMI Katarina Mramor1, Robert Vertnik2 3, Božidar Sarler1'3 4 1CO BIK, Tovarniška c. 26, 5270 Ajdovščina, Slovenia 2Štore Steel, Železarska c. 3, 3220 Štore, Slovenia 3University of Nova Gorica, Vipavska 13, 5000 Nova Gorica, Slovenia 4Institute of Metals and Technology, Lepi pot 11, 1000 Ljubljana, Slovenia katarina.mramor@cobik.si Prejem rokopisa - received: 2013-12-30; sprejem za objavo - accepted for publication: 2014-01-13 The initial results obtained with the local-radial basis-function collocation method (LRBFCM) for a two-dimensional (2D), simplified model of continuous casting of steel with an externally applied magnetic field are presented. The multiphysics model is composed of turbulent-solidification equations (mass, momentum, energy, turbulent kinetic energy, turbulent dissipation rate) and Maxwell's equations that are numerically solved for the non-uniform node arrangement. The numerical procedure is structured using the explicit time stepping and local collocation with multiquadric radial basis functions (MQ RBF) on the overlapping five-node subdomains. The pressure-velocity coupling follows the fractional-step method (FSM) and the convection is treated with adaptive upwinding. The novel LRBFCM has been already verified in several benchmark test cases, such as the natural convection in a cavity with a magnetic field, the lid-driven cavity, and the flow over the backward-facing step with a transverse magnetic field. Keywords: continuous casting of steel, turbulent flow, solidification, magnetohydrodynamics Namen članka je predstavitev prvih rezultatov s poenostavljenim 2D-modelom za kontinuirano ulivanje jekel pri zunanjem magnetnem polju, izračunanih z metodo kolokacije z radialnimi baznimi funkcijami. Numerični model združuje Navier-Stokesove in Maxwellove enačbe, ki jih rešujemo na neenakomerni porazdelitvi točk. V numeričnem postopku je uporabljena eksplicitna časovna shema na pettočkovnih poddomenah, na katerih so uporabljene multikvadrične radialne bazne funkcije. Sklopitev tlaka in hitrosti se rešuje z metodo delnih korakov. Omenjeno metodo smo poprej verificirali za izračun različnih preizkusov, kot so naravna konvekcija v kotanji z magnetnim poljem, kotanja z vsiljenim tokom ter tok preko stopnice v kanalu z magnetnim poljem. Dobljeni rezultati kažejo dobro ujemanje z drugimi numeričnimi postopki, kot je npr. metoda končnih volumnov. Ključne besede: kontinuirano ulivanje jekla, turbulenten tok, strjevanje, magnetohidrodinamika 1 INTRODUCTION stand and further improve the process. A number of different numerical models3 have so far been used in the The production of continuously cast steel1 has greatly simulations of the problem. They include the finite- expanded in recent years. Continuous casting of billets, volume method (FVM),4-8 the finite-element method blooms and slabs is the most common process of steel (FEM)9 and some more advanced meshless methods like production.2 The continuously growing demand for cast the local-radial basis-function collocation method steel has fueled the need to produce the steel of even (LRBFCM).10 better quality. Although the process of continuous casting of steel is cost efficient, exhibiting a high yield and good quality of the products, it can be further improved 2 GOVERNING EQUATIONS by introducing the electromagnetic (EM) field into it. The EM force, which is a result of an applied magnetic Jhe sy.!tem of governin8 equation^ that describes the turbulent heat transfer and the fluid flow as well as the magnetic field in the continuous casting of steel, is based field, affects the velocity and temperature fields. By adjusting the magnetic field, the amount of defects in the material can be significantly reduced. on.the Rey"oldIs time-averaging approach to m«delinng a The magnitudes of the velocity, the temperature and turbulent flow^ and a mixture-continuum model, first magnetic fields are all crucial to the final quality of a introduced by Bennon and Incropera.12 The system con- product. As all these quantities are difficult or impossible sists of five time-averaged mixture equations: to measure, numerical models help us to better under- V v = 0 (1) p ^dV+pV( vv) = -Vp + v[( n L + n t)[Vv+(Vv) ^ ]]- 2 ß L(1-/L)2 - 3 P^k - Cf3 (v - v s)+Pß T 8(T - T/)+Fm p ^dhh+pV(vh ) = V(XVT)+pV{vh - fS v S h S - (2) -/l v L h l) + V /l ^^^ Vh L (3) dk P +P^( vk ) = V ß L + p,. Vk +^k + G k - (4) n (1-/l)2 , -pe+pD - ß L ——^— k C/l3 ß L +: ß t ^ 1 Ve +pE - de p dt + ^V( ve) = V (1-fL)^ r T e -ß L C/' ^ + e fl(Pk + CseGk ) - C2e f2 Pe] ^ (5) where v is the velocity of the mixture, p = pS= pL is the density of the mixture, assumed to be constant; t stands for the time and p for the pressure; ^t is the turbulent viscosity and ^l is the dynamic viscosity; k represents the turbulent kinetic energy and K is the permeability of the porous matrix. /3t, g, T, and Tref represent the thermal-expansion coefficient, the gravitational acceleration, the temperature and the reference temperature, respectively. Fm stands for the Lorentz force, h for the enthalpy andX for the thermal conductivity. /s, /l, vs, vl, hS, and hL represent the solid-volume fraction, the liquid-volume fraction, the velocity of the solid phase, the velocity of the liquid phase, the enthalpy of the solid phase and the enthalpy of the liquid phase, while Vt is the turbulent kinematic viscosity. Symbols e and C stand for the dissipation rate and the morphology constant of the porous media. Symbols Pt, Pk, Oe, cie, fi, cie and /2 are the closure coefficients of the turbulence model. Pk, Gk, D and E are the shear production of the turbulent kinetic energy, the generation of turbulence due to the buoyancy force, the source term in the k equation and the source term in the e equation, respectively. The closure relations defined by Abe, Kondoh and Nagano13 are used in the present work. A detailed description of the closure coefficients, source terms and damping functions are given in the paper by Sarler et al14. The Lorentz force is defined as: Fm = j X B (6) where j and B are the current density and the magnetic-flux density. Maxwell's equations are used to calculate the current density: where o is the fluid electric conductivity, ^ is the fluid electric potential and B is the externally applied magnetic field. The induced magnetic field is assumed negligible in comparison with the applied magnetic field. 2.1 Initial and boundary conditions The governing equations are strongly coupled. For the steady solution of a problem it is, therefore, very important how the initial conditions are chosen in order to minimize the required iterations to reach the steady state. Five different boundaries are chosen: the inlet, free surface, wall, outlet and symmetry of the present model. A detailed description of the initial and boundary conditions are given in the article by Sarler at al14. The model of the domain is presented in Figure 1 and the computational domain is depicted in Figure 2. 3 SOLUTION PROCEDURE The solution of the governing equations is obtained with the LRBFCM employing the explicit time stepping. The fractional-step method (FSM)15 is used to solve the pressure-velocity coupling and an upwinding scheme is used to stabilize the highly convective situation.16 The solution procedure begins by calculating the initial Lorentz force (equation (6)). Afterwards, the intermediate velocity is calculated, without a pressure gradient. The pressure is then calculated from the Poisson equation14 by assembling and solving a sparse matrix.17 The calculated pressure gradient is afterwards used to correct the intermediate velocities. After the solution of the velocity field, the equations for the turbulent kinetic energy and dissipation rate are solved. This is followed Figure 1: Simplified 2D model of continuous casting of steel Slika 1: Poenostavljeni 2D-model kontinuirnega ulivanja jekla j = o(-V0 +v X B) (7) Figure 2: Scheme of the computational domain Slika 2: Shema ra~unske domene by the solution of the enthalpy equation. The temperature is calculated from the enthalpy, using the temperature-enthalpy constitutive relation.1718 Finally, the turbulent viscosity, velocity, temperature, turbulent kinetic energy and dissipation rate are updated and the solution is ready for the next time step. The LRBFCM is structured in the following way: Approximation function 6 is represented on each of the subdomains as a linear combination of the radial basis functions (RBFs) as: MI q( i p „) =X i ^ i(i p „) i y i (8) where i^i, ryi, and M represent the RBF shape functions, centred in points /p„, the expansion coefficients, and the number of shape functions, respectively. The most commonly used RBF is the multiquadric RBF:19,20 i(P) = 4 irf (P) + c' (9) where c stands for the dimensionless shape parameter, set to 32 in all the calculations, and the distance between the nodes: irf(p) = ^ x - x, \ l imax / + y - Ti V l y imax J (10) A subdomain is formed around each of the calculation points consisting of the iM - 1 nodes nearest to node ip„. For the purpose of this work, five-node overlapping subdomains are used. By considering the collocation condition of: i 6(i p„) = 6i(i,„) (11) a linear system of equations is obtained. To solve the partial differential equations (PDEs), the first and second derivatives of function i6(p) have to be calculated: dM dj 7i 6(p) i ^ i (p) i y i dz i=1 dx' (12) where index ' = 1,2 is used to denote the order of the derivative and x = x,y. A detailed explanation of the solution procedure is given in the paper by Vertnik et al.17,18 The schematics of the discretization scheme is shown in Figure 3. 4 NUMERICAL IMPLEMENTATION The sparse-pressure matrices, used for the solution of the pressure and the flux-density Poisson equations, are solved by applying the Pardiso routine and Intel Math Kernel Library 11. The OpenMP Library is used for the parallelization. The post processing is performed in PGPlot, Gnuplot 4.4 and Octave 3.6.1. The results were calculated on an HP Proliant DL380 G/7 server running on 64 bit MS Windows. 5 NUMERICAL EXAMPLES The numerical procedure has so far been verified on the following benchmark test cases: the lid-driven cavity,21 the natural convection in a cavity,22 and the backward-facing step.23 The lid-driven test case was used is scaled with iXimax and iyimax, the scaling parameters of subdomain i in the x and y directions, respectively (Figure 3). Figure 3: Discretization scheme. F, Q, ixiMAX and iyiMAx represent the boundary, domain and scaling parameters in the x and y directions, respectively. The circles represent the boundary nodes, whereas the black dots represent the domain nodes. Slika 3: Diskretizacijska shema. F, Q, ixiMAX in iyiMAX ozna~ujejo rob, obmo~je ter skalirna parametra v smeri x in y. Robne to~ke so ozna~ene s krogci, obmo~ne to~ke pa s pikami. =1 to verify the coupling between the mass and momentum equations.21 In the natural-convection case, the energy-conservation equation was added.22 For the backward-facing step problem, the boundary conditions for the inflow and outflow problems were tested. The natural-convection and backward-facing-step problems were first tested without a magnetic field and then with an external magnetic field. The results of all of the test cases are in good agreement with the reference results, calculated with a commercial code or obtained from several published sources. Since the results of our calculations, mentioned in the above history of test cases, are in good agreement with the reference results, the method is subsequently applied to the simplified problem of continuous casting of steel with a magnetic field, as defined in the paper by [arler et al.14 5.1 Continuous-casting geometry and material properties The simplified 2D continuous-casting model geometry is shown in Figure 1 and the elements of the discretization are presented in Figure 2. The computational domain coincides with a half of the longitudinal section of the billet, taken to be 1.8 m long and 14 cm wide. The SEN diameter is 3.5 cm, the mold height is 0.8 m and the EM-coil height is 10 cm. The material properties of steel are temperature and steel grade dependent. However, for the purpose of the present simplified model, constant values are used. The values are given in Table 1. Table 1: Simplified material properties of steel Tabela 1: Poenostavljene snovne lastnosti jekla Property Value P 7200 kg/m3 1 30 W/(m K) Cp 700 J/(kg K) Ts 1680 K Tl 1760 K hm 250000 J/kg 0.006 Pa s ßT 1 • 10-4 K-1 C 1.6 • 108m-2 The results of the numerical simulation are represented in the following sections of the paper. 5.2 Convergence of the method and a comparison with the reference results for the continuous-casting process First the convergence of the method was tested for the node arrangements with 20 220, 50 951, 73 940, 100 089, 131 452, 165 426 non-uniformly arranged nodes. The convergence was tested for the velocity and the temperature fields. The vertical and horizontal components of the velocity were compared for three different vertical cross-sections: just before the application of the mag- netic field, just after the application of the magnetic field and at the end of the computational domain, as can be seen in Figure 4 (left). The temperature was verified at two different vertical directions, the one at the center of the domain and the one at the outside wall, as shown in Figure 4 (right). As can be seen in Figures 5 to 12, the smallest appropriate node arrangement is 100 089. The results were also compared with the results obtained with the commercial code based on the FVM (Fluent23). The smallest appropriate amount of the finite volumes in the commercial code for reaching a reasonable convergence was 169 169. The agreement of the velocity profiles with the commercial code (the solid line denoted with F) is similarly good for the horizontal as well as vertical velocities. The largest differences occur at the positions closer to the inlet. At the positions closer to the end of the computational domain, where a fully developed flow and a partial solidification take place, the agreement between the results obtained with the LRBFCM and FVM is excellent. In the early stages of the flow, slightly larger differences can be observed. The differences between the commercial-code results and the results obtained with the in-house built LRBFCM are reasonably small. The exact cause for the differences is Figure 4: Left: the positions of the velocity-field comparison lines; hI = -0.8 m, hII = -0.9 m, hIII = -1.8 m. Right: the positions of the temperature-field comparison lines; vI = 0.0 m (centreline), vII = 0.07 m (surface). Slika 4: Levo: Položaj linij, kjer je primerjano hitrostno polje; hI = -0,8 m, hII = -0,9 m, hIII = -1,8 m. Desno: Položaj linij, kjer je primerjano temperaturno polje; vI = 0,0 m (sredina), vII = 0,07 m (površina). not known; however, several reasons why the results are not identical might be identified: the solution procedure, e.g., the energy equation in the commercial code is solved iteratively and the temperature is obtained directly from the equation, whereas in our method, the enthalpy equation is solved first and the temperature is obtained by solving the enthalpy-temperature relationship. Another reason might lie in the differences between the methods, e.g., in the commercial code a second-order upwind technique is used, whereas in our case an adap- Figure 5: Horizontal-velocity profiles at vertical position hI = -0.8 m, at the top boundary of the magnetic field Slika 5: Profili horizontalne hitrosti na mestu hI = -0,8 m na zgornji meji magnetnega polja Figure 8: Horizontal-velocity profiles at vertical position hII = -0.9 m, at the bottom boundary of the magnetic field Slika 8: Profili horizontalne hitrosti na mestu hII = -0,9 m na spodnji meji magnetnega polja Figure 6: Horizontal-velocity profiles at vertical position hI = -0.8 m, at the top boundary of the magnetic field Slika 6: Profili horizontalne hitrosti na mestu hI = -0,8 m na zgornji meji magnetnega polja 2.010-^ 0.0-10° -2.010^ I -4.010'^ 7 -6.010"^ -1.0 10"^ -1.210 ,-3 m ...... - ■io-m • 100089 " F — ■ S09S1 1314M • 73S40 16M26 • 1 1 x[m) Figure 9: Horizontal-velocity profiles at vertical position hIII = -1.8 m, at the end of the computational domain Slika 9: Profili horizontalne hitrosti na mestu hIII = -1,8 m na koncu računske domene Figure 7: Horizontal-velocity profiles at vertical position hII = -0.9 m, at the bottom boundary of the magnetic field Slika 7: Profil horizontalne hitrosti na mestu hII = -0,9 m na spodnji meji magnetnega polja -0.01 -0.015 •0.02 ^ -0.025 1 -0.03 ^ -0.035 -0.04 -0.045 -0.05 - T-1-1-1 - - 20220 100089 > F — 50951 131452 • - 7a94o 165426 • I 1 1 1 1 oödööödo x(m] Figure 10: Horizontal-velocity profiles at vertical position hIII = -1.8 m, at the end of the computational domain Slika 10: Profili horizontalne hitrosti na mestu hIII = -1,8 m na koncu računske domene 1800 1795 1790 1785 g 1780 H 1775 1770 1765 1760 1755 20220 ' 'S 50951 " M ' 73M0 s - 100069 fl M ' 131<52 * jB ■ 1ftM26 J8 F /i - ........ CO -0.01 -0.015 -0.02 -0.025 -0.03 -0.035 -0.04 -0.045 -0.05 - 1 1 1111 - - B=0 6=0.26 T - 8=0.0026 T » 8=2.6 T > - B=0.D2B T 1 1 □ 1 1 1 1 CO o O) o