Norton-Hoff model for deformation of growing solid shell of thin slab casting in funnel-shape mold

A funnel-type mold is commonly used to provide necessary clearance for the submerged entry nozzle in the thin slab casting (TSC). The partially solidified shell is subjected to the mechanical deformations, which can lead to the defects formation and, as a results, to a breakout. Traditionally, the results of the flow simulation, performed by the finite volume method (FVM), are fed to the external package for the finite element analysis of stress and strain. A coupled model was assembled using “creeping solid” approach by blending the Norton-Hoff viscoplastic stress for the solidifying shell with the Newtonian viscous stress of the liquid melt. The FVM was used to combine both liquid and solid stress models within a single solver. The iterative procedure based on the improved both side diffusion method was introduced to treat the nonlinear relation between the viscoplastic stress and the strain rate. The modeled shell thickness was verified by previously published breakout measurements and the simulation results. Temperature distribution, obtained during the TSC simulation, dominantly corresponds to the viscoplastic range. Developed numerical approach is robust and has direct industrial application.


List of symbols
Operators r Gradient operator (m -1 ) r T Transpose of gradient operator (m -1 ) rÁ Divergence operator (m -1 ) r 2 Laplace operator (m -2 ) a b Outer product of two vectors A : B Double inner product of two tensors dev A ð Þ Deviatoric part of a tensor tr A ð Þ Trace of a tensor The original online version of this article was revised: assignment of affiliations was not correct. Greek symbols _ e Strain rate tensor (s -1 ) _ e vp Viscoplastic strain rate (s -1 ) _ e vp Viscoplastic strain rate tensor (s -1 ) _ e vp eq ; _ e eq Von Mises equivalent strain rate (s -1 ) k Thermal conductivity (W m -1 K -1 ) k ' Liquid thermal conductivity (W m -1 K -1 ) k s Solid thermal conductivity (W m -1 K -1 ) K Perzyna flow rule coefficient (Pa -1 s -1 ) l ' Dynamic viscosity of liquid (Pa s) l vp Apparent viscoplastic viscosity (Pa s)

Introduction
During the solidification, the properties of the alloys, according to the pioneering work of Flemings [1], drastically vary at the two-phase region, also called the mushy zone. In the continuous casting (CC), the defects are formed and accumulated under the mechanical strain during the casting process and can result in a breakout due to the variation of the alloy properties. The course of the initial solidification in the CC mold, the interaction of the growing shell with the turbulent melt flow and its behavior during a withdrawal are one of the key features for the final product quality in the CC process [2][3][4][5][6].
The extended reviews on the methodology and practical applications of the numerical modeling of the continuous casting are presented throughout two decades starting from early two-thousands by Thomas [7] and Atkinson [8] and supplemented by recent works by many researchers [2,4,5,[9][10][11].
In early two-thousands, Thomas [7] stressed in his review that the enhanced heat flux at the exit of the mold could be an indication of the mold taper pressing to strongly against the shell. Additionally, the following factors, such as shell thickness at mold exit, the metallurgical length and the corresponding cooling practice are important, since the slab becomes susceptible to cracks due to the off-corner depression, reheating or non-uniform shell growth.
According to Ramirez Lopez et al. [5], the numerical modeling became a key tool to investigate the lubrication concepts in the continuous casting, giving insights into the heat transfer and solidification guided by the flow with an aim to predict defect formation, mentioned by Flemings [1] and Thomas [7]. It is promoted that by the oversimplification of the numerical models, the key issues can be missed. The numerical modeling achievements supplemented with the industrial observations should improve the knowledge and the conceptual point of view on modeling techniques.
Novel modeling techniques are applied nowadays to resolve the transient multiphase phenomena during CC, e.g., meshless approach employed by Vertnik et al. [12] to tackle the turbulent flow and solidification at the presence of the electromagnetic steering. Reduction in the problem dimensions allows to perform the investigation of the solidification during casting in shorter terms [13][14][15][16] or combine them with the industrial measurements to define proper casting strategy to avoid defects formation [17,18]. Current authors suggested a novel single mesh approach to simulate the heat transfer during CC including convection in the liquid slag to predict the parasitic solidification in the submerged entry nozzle (SEN) region [19]. The scale adaptive modeling of the multiphase mold flow with application to the analysis of the opened eye formation was recently published elsewhere [20,21].
As a funnel-type mold is used to provide necessary clearance in the narrow cavity region for the SEN, the (partially) solid shell of the thin slab casting is much stronger subjected to the mechanical deformations than that of conventional CC. Traditional approach for the funneltype mold simulation, as presented in Vakhrushev et al. [22], provides steady-state distribution of the solid shell velocities. The full elasticity of the solid structure and neglecting its rheology variations due to the cooling are restrictive assumption. However, the steady-state approximation cannot tackle the transient phenomena that occur during continuous casting. The shell compression and tension are the provenance of the cracks and the hot tearing formation [23]. As recently revealed, the interaction of the highly conductive shell with the induced electric current during electromagnetic brake leads to the complex topology of the Lorentz force and to the complete alternation of the melt flow pattern [24][25][26].
The aim of the current work is to track the transient mechanical behavior of the solidifying shell during continuous casting in a funnel-type thin slab mold by employing a transient viscoplastic stress model. Traditionally, a decoupled one-way approach is employed in the computational fluid dynamics (CFD) community for complex full-size geometries. The flow simulation results are typically obtained with the finite volume method (FVM). Next, they are fed to an external package for the finite element analysis (FEA) [27][28][29]. Substantial contribution in this field was done in previous studies [30][31][32][33][34] by developing the arbitrary Lagrangian-Eulerian (ALE) formulation in the context of the alloy casting. The method includes a solution of the governing conservation laws for mass, momentum, and energy on the Eulerian mesh. The coupled solution is reflected via elasto-viscoplastic stress model to the corresponding mesh motion algorithm. To model the continuous casting, such an approach requires several steps: firstly, the thermomechanical analysis is performed by treating the mush as a non-Newtonian fluid; next, the initial and boundary conditions for the two-phase simulation are extracted from the first calculation by defining the point of the coherency of the two-phase core ([ 65% of solid) and continuously crowing mesh in the frame of the ALE approach with the casting speed using special buffer zone as well.
Due to the complex cross-methodological (FVM/FEA or Lagrangian/Eulerian) coupling, such approaches are limited to 2D or axisymmetric geometries, reducing the studies mostly to the steady-state cases. The currently presented model is an attempt to overcome these limitations.
Additional motivation of the present study is to verify a single-phase volume-averaging approach for the application to model a viscoplastic stress in the two-phase region during alloy solidification. On the contrary to the previously applied multiphase Eulerian-Eulerian volume-averaging model for the equiaxed crystal sedimentation [35][36][37] and the twin-roll casting process [38][39][40], the socalled mixture model is very robust, numerically stable and fast in sense of the calculation time required. As an advantage of the presented work, a full 3D model is considered including an SEN, a funnel-type thin slab mold (879 mm in height) and strand part (4000 mm along casting direction including mold). The boundary layers were introduced to correctly predict near-wall flow and heat flux. A highly turbulent flow in the SEN (with maximum velocities of 3 m s -1 ) drastically limited the efficiency of the time marching algorithm to integration step of 10 -4 s. Taking all arguments into account, the multiphase approach was found to be inappropriate for transient fullscale industrial simulation.
In the present study, a coupled model is assembled using viscoplastic approach in the frames of the ''creeping solid'' hypothesis. In this method, as the solid fraction increases specifically above a certain (coherency) limit, a considerable strain resistance of the solidifying alloy develops. The mush represents a columnar skeleton with liquid in the inter-dendritic mushy zone [41], which can sustain high stresses, while the solid and the liquid become strongly coupled. That is usually reflected in the viscoplasticity theory by an abrupt jump of the effective viscosity.
In current work, the viscoplastic Norton-Hoff type stress is applied to model a complex process of the thin slab continuous casting, and the first results are discussed. The blending of the stress tensor for the solidifying shell is done with the viscous liquid stress based on the solid fraction distribution. The advantage of the proposed approach is using finite-volume formulation within a single solver for both solid shell mechanical behavior and the strong melt flow. The latter is typically excluded in most studies on the mechanical stress calculations by simplifying to a plugtype flow. In the present study, the flow phenomenon is fully considered. The iterative procedure is introduced for the treatment of the nonlinear relation between the viscoplastic stress and the strain rate tensor.
The classical thermomechanical models were developed for the semi-solid state, providing the constitutive equations for porous metallic materials saturated with liquid. The difference between the cases lies in the solid-phase permeability and the material globular [41] as well as the dendritic morphology [42][43][44]. Some of the models are also restricted to the isothermal conditions only. The great endeavors were done by Fachinotti et al. [31,45] by revising constitutive temperature and composition-dependent models and combining with the Lagrangian-Eulerian approach of Bellet and Fachinotti [30] to apply it to the actual casting processes.
According to the book of Rappaz et al. [46], the metal alloys perform as an elastic material till 1/3 of their absolute melting temperature. Then, the creeping occurs corresponding to the elasto-viscoplasticity of the material, and approximately, above 2/3 of the melting temperature, a fully viscoplastic behavior is observed. The viscoelasticity of an alloy represents rate-dependent unrecovered deformations of the material. Rate dependence in this context corresponds to the creeping of the solid according to the rate of load has been applied. Temperature distribution, obtained during the thin slab casting simulation, lies dominantly within the viscoplastic range slightly catching transition to the elastoplastic behavior.
The numerical results comparison showed good agreement with the previously obtained and experimentally verified ones. The detailed description of the numerical model and the simulation results are presented in the next sections; finally, the discussion and the conclusions are made.

Numerical model
Based on the previous experience of developing a numerical Eulerian-Eulerian model to simulate behavior of the packed bed of equiaxed crystals [35][36][37] and the macrosegregation during the twin roll casting [38][39][40], the viscoplastic stress model was applied to the thin slab casting by reducing it to the ''mixture'' formulation, which is more robust for the industrial scale simulations.
The following assumptions were made to model the viscoplastic behavior of the solidifying shell during its withdrawal: (i) An averaged ''single''-phase approach is employed considering the simulated media to be a mixture of the liquid and solid. (ii) Main material properties such as density, specific heat and thermal conductivity are constant across both phases, and moreover, they are temperatureindependent; thermo-solutal buoyancy is ignored. (iii) Temperature and pressure field equilibrium is assumed between liquid and solid phases. (iv) Total effective stress incorporates both the Newtonian and the viscoplastic one according to the phase fraction in the control element of the computational grid; thermomechanical stress is ignored.

General equations of motion
The continuity and momentum equations for the incompressible liquid become: To perform the coupling between the solid and liquid phase, the total stress tensor is divided into deviatoric and hydrostatic parts, representing shear and volumetric forces correspondingly. Next, a volume averaging of all terms is considered: Omitting the volume averaging notation hereinafter for all further derivations, the deviatoric part of the total stress tensor can be expressed as a blending of the viscoplastic stress for the solid and the Newtonian one for the liquid melt: where the condition f s þ f ' 1 is valid. The deviatoric part of the liquid stress tensor is where _ e represents a symmetric part of ru: Norton-Hoff model for deformation of growing solid shell of thin slab casting in funnel-shape mold 91 To clarify the reformulation in right side of Eq. (5), one should revive the property of the _ e from Eq. (6) by recalling the continuity condition of incompressible fluid Eq. (1) and by simple manipulation. Consistently, it can be derived that In other words, the strain rate is originally a deviatoric (or traceless) tensor for incompressible media. That will be later used in this work for the corresponding substitutions in further derivations, e.g., in the equivalent strain rate formulation. However, for the rheological models including the mush compressibility, such as reported by Fachinotti et al. [31], the non-hydrostatic terms will reflect the so-called bulk viscosity and the corresponding two-phase region compression. Here, the authors left this topic for the future detailed investigations. As a concluding assumption of the presented approach, the pressure balance at the solid-liquid interface is taken, thereby leading to a single ''pressure'' field for both phases. To avoid confusion, the ''pressure'' is regarded as a spherical tensor of the volumetric forces defining compression-depression in the material regardless of its ''liquid'' or ''solid'' nature.

General elasto-viscoplasticity at high temperatures
General solids show the elastic-viscoplastic behavior. The solid becomes purely viscoplastic close to the solidus temperature when the yield and hardening are negligible. Viscoplastic behavior is referred as the Ostwald-de-Waele law or more commonly the Norton-Hoff law [47].
With the rising temperatures, the strength properties of steel undergo significant weakening and become strainrate-dependent, showing creeping behavior. At low strains, the solid stress can be split in two components by decoupling elastic and viscoplastic parts as formulated by Norton for the uniaxial stress [48]: The difference r elvp À r stat is a so-called extra-stress.
Below that threshold, the material acts as the elastic one. Above the limit, an irreversible portion of deformation happens, characterized by a _ e vp . The classical power law in Eq. (8) was introduced by Norton [48] including two parameters and m. The von Mises equivalent stress is introduced as: where dev r elvp À Á ¼ r elvp À 1 3 tr r elvp À Á I. The von Mises equivalent strain rate for the incompressible media by recalling Eq. (7) is The multiaxial analogue of Eq. (8) is obtained by generalizing it to the form of Applying the overstress model and formalization of Perzyna [49,50], the elasto-viscoplastic flow rule is obtained [46], where its viscoplastic part corresponds to with n ¼ dev r elvp À Á (see Eq. (9)) and the viscoplastic consistency K vp . By manipulating left-and right-hand side of Eq. (12) to extract the von Mises invariants according to Eqs. (10) and (9), one obtains: After the nonlinear function K defined in Eq. (12) is plugged into Eq. (13), the extra-stress becomes which allows to derive an expression for the consistency parameter in Eq. (11): Combining the extra-stress from Eq. (14) with the Perzyna formulation (see Eq. (12)) gives a compact form of the elasto-viscoplastic flow rule by reducing it to Summarizing, in this section, the rheological model reflecting the solidified alloys behavior at a high temperature range (above 1/3 of the liquidus) is described. Two model parameters, K vp and m, are obtained experimentally and will be described later.

Pure viscoplasticity without static threshold
With continuously growing temperature at some critical point, the plasticity limit, r stat , cannot be defined. In this regime, the unrecoverable deformations occur when the material undergoes even the weak loads. Such a rheological phenomenon is regarded as the pure viscoplasticity, which is valid, according to the book of Rappaz et al. [46], at temperature higher than 2/3 of their melting point.
The pure viscoplasticity law can be directly derived from the elasto-viscoplasticity formulation as presented in Eq. (8). The latter is reduced by excluding the static plastic threshold, assuming r stat 0, to the pure power law. The flow rule Eq. (12) together with its compact form Eq. (16) can be inverted to the pure viscoplastic constitutive law, referred in the literature as the Norton-Hoff law: where R vp is accepted to be a deviatoric part of the r vp .

Constitutive equation for viscoplastic stress and strain
As it was mentioned here previously, the viscoplastic model parameters are estimated experimentally. In the present study, a constructive equation model is used corresponding to the Model AI reported in Kozlowski et al. [51]. As reported by the authors, it is applicable for a wide range of low strain rates between 10 -3 and 10 -6 s -1 at high temperature range of 1223-1673 K when the carbon content is varied between significantly low (0.005 wt.%) and the high one (1.54 wt.% It is worth reminding that a peritectic steel casting is considered with the 0.06 wt.% C alloy (the same as in the previous study [22]) which fits in the range of Kozlowski experimental work [51].
The original model includes constitutive equation, which expresses the _ e vp as a function of the r vp , temperature and the carbon content of the alloy [51]. The onedimensional flow rule is expressed as an Arrhenius-type equation [45]: where the relation between the Q and R is fixed to a value Q R ¼ 49; 890 K. According to the accepted notation, the Q is defined as a variable parameter (which is generally a case); however, in the present study, it is included in the constant relation and does not vary. The material parameters are described based on T and c as follows [45]: However, Eq. (18) is derived for the uniaxial stress and strain case. To extend it for the multidimensional deformations, its scalar analogue for the von Mises invariants are employed. Omitting the direct flow rule formulation, we switch to a reverted version derived in respect to the equivalent stress: The elasto-viscoplastic extra-stress (see Eq. (14)) is rewritten for pure viscoplasticity without static plasticity threshold as: Simple comparison of Eqs. (21) and (22) reveals that the strain rate sensitivity is expressed via the n as m ¼ 1=n, which results in Finally, the viscoplastic stress law Eq. (22) and the inverted Kozlowski flow rule Eq. (23) are combined to derive the viscoplastic consistency parameter K vp from the constitutive model in a form of Thus, the parameter K vp is now expressed as a function of temperature and the carbon content.
It is crucial to notice that since in Model IA in Kozlowski et al. [51] the stress is measured in MPa and the material constant C is taken in MPa -n s -1 [45,51], thereby the K vp (after some trivial manipulations) gets MPa s 1/n or MPa s m as its units.

Viscoplastic model update algorithm
An important remark should be done for the one-phase formulation used in this study. The equality of the viscoplastic and liquid strain rate ( _ e vp _ e) is assumed since a single velocity field is used for both phases. Thereby, hereinafter, a notation of _ e is used in the completed definition of the blended stress tensor deviator in Eq. (4). Moreover, in the von Mises invariant _ e vp eq ; the viscoplastic superscript is omitted as well as leading to its redefinition as _ e vp eq _ e eq . Generalizing the viscoplastic stress formulation in Eq. (17), the ''effective'' or ''apparent'' l vp is introduced as a nonlinear function of the equivalent strain rate _ e eq : The strain rate characterizes the velocity gradient between the adjoin ''layers'' of the creeping material: the higher the strain rate is, the more velocity change is observed at the neighboring locations. The viscoplastic viscosity depends on this local characteristic of the motion which is referred as the rate-of-load. Since the strain rate is integrated into the ''effective viscosity'' calculation, a coupling procedure is applied to converge the nonlinear terms: (i) the strain rate is estimated using Eq. (6); (ii) according to Eqs. (10) and (25), the _ e eq and the l vp are updated; (iii) a new solution of the flow/heat transfer equation system is integrated based on the latest estimation of the blended stress distribution advancing to the next time step; (iv) repeat procedure by applying under-relaxation on pressure, velocity fields, released latent heat term as well as on the effective viscosity; (v) repeat above steps to reach a convergence criterion for the coupled solution; (vi) proceed to the next time step.
Additional coupling step was to overcome the solver stability issues. Since a strong nonlinear term in the Norton power law can cause the divergence of the solution, the enhanced numerical treatment was applied according to the improved both side diffusion method by Fernandes et al. [56] and more recently extended by Araújo et al. [57]. It employs the artificial diffusion term added to both sides of the momentum equation. The terms are treated separately in the implicit and explicit manner for the left-and the right-hand sides correspondingly. The method originates from the FEA and aims to tighten the velocity/stress coupling. It simultaneously helps avoiding the checkerboard pattern in solution. The improved algorithm is based on the larger computational stencil. Such a discretization is done due to the extended flexibility of the OpenFOAM CFD package.

Heat transfer and solidification model
As it was shown previously by the authors, the latent heat advection plays a dramatic role in build-up of the solid shell thickness and depends on the columnar phase velocities predicted by the model [22,58]. However, in the coupled model for the flow and solidification, newly suggested here, which considers the viscoplasticity of the solidified shell, no information regarding the solid-phase velocities can be obtained due to the single-phase formulation. Assuming the constant q and k, the energy equation is taken in its general advection-diffusion form regarding the total enthalpy: The total enthalpy for the mixture is estimated in the current approach by ignoring the temperature dependency of the C p and the L: By inserting Eq. (27) into Eq. (26), the latent heat advection term is derived as follows: The simulation results, using previously developed formulation [22,58], are considered as a reference case. The corresponding latent heat advection term is: The difference between two models, described in Eqs. (28) and (29), becomes obvious by simple comparison: in the current approach, the u is used, whereas the u s is responsible for the latent heat advection within the frames of the previously used method. The details and outcomes are discussed later.

Modeling thin slab casting
According to the approach, presented in Sect. 2, the coupled numerical model was implemented by the authors as a standalone solver in the open-source CFD package Open-FOAM. The description of the numerical setup, the simulation results and their discussion come in the very next sections.

Viscoplastic model parameters
The basic parameters of the viscoplastic deformation model were selected here based on the Kozlowski constitutive Model IA (see Eqs. (19) and (20)), and the corresponding expression for the K vp is derived here in Eq. (24). Figure 1 shows that the temperature dependency of m and K vp is presented and analyzed for the selected constitutive model used in this study. As it is observed, the m drops during cooling, which corresponds to faster response of the apparent viscoplastic viscosity to small strain rates, as shown in Fig. 2. Corresponding viscoplastic consistency grows during cooling, defining stronger viscoplastic stress magnitude in general. Additionally, the temperature range of the experimental work of Kozlowski et al. [51] is marked in Fig. 1 along with the liquidus and solidus temperatures of the alloy used for the thin slab casting simulation (see Table 1 for the details).
As shown in Fig. 1, the sensitivity parameter varies between 0.136 and 0.170 for the investigated alloy and in the pure viscoplastic temperatures, whereas the K vp fits into the 5-115 MPa s m range. However, to simplify the analysis of the results, the presented studies were limited to the constant viscoplastic model parameters. The fixed strain rate sensitivity m * = 0.138 and the viscoplastic consistency K Ã vp = 44.62 MPa s m were selected corresponding to the reference temperature T Ã % 1170:33 K, which is referred to the transition between elasto-and pure viscoplastic rheology hypothesis according to the book of Rappaz et al. [46].
As shown in Fig. 2, with a decreasing parameter m, the apparent l vp dramatically grows in correspondence to Eq. (25) for the range of the _ e eq below 10 -2 s -1 , and thus the material strongly resists to the small velocity gradients namely the relative deformations in the material. On the contrary, its motion is allowed for significantly higher deformations since its resistance drops and the material starts to creep easier. For the sensitivity parameter m close to identity for the whole rate-of-strain range, the viscoplastic viscosity does not change dramatically, and more relative motion is allowed within the solid phase. However, based on previous studies and literature sources for the metal alloys, the casting range for the parameter m lies between values of 0.1 and 0.2, thereby in the region of rapid response to strain rates.   The strain rates range (expressed for the von Mises invariant _ e eq ), which was used in the experimental work of Kozlowski et al. [51], is marked in Fig. 2 with a cyan stripe. Typical casting range for the strain rate sensitivity is marked in gray according to the literature and to the previously plotted temperature dependency in Fig. 1. The rheological curve of the apparent viscosity l vp ¼ l vp ðm; K Ã vp; _ e eq Þ according to the m * is marked in red. Current model uses, as shown in Fig. 3, the solid fraction/temperature correlation for the industrial alloy (see Table 1) determined by an engineering software IDS [59,60], which is more complex than, for example, corresponding linear relation in Svensson et al. [27].

Numerical grid and model setup
The modeling domain sketch is shown in Fig. 4a; it includes the SEN, water-cooled mold cavity of the funnel type (879 mm in height) and the strand part (4000 mm along casting direction including mold) marked as the secondary cooling zone. The details of the dominantly hexahedral mesh are shown in Fig. 4b. The refined boundary layers are introduced to accurately resolve the extracted wall heat flux, which defines the growth of the solid shell. The melt flow/ two-phase region interaction should also be carefully tracked along with the corresponding viscoplastic stress in the mush during the shell withdrawal. Thereby, the mesh refinement was extended toward the bulk region. Geometry dimensions for the 1726-mm-wide and 72-mm-thick funnel-type mold are shown in Fig. 4c. The corresponding applied heat flux profiles, adopted from the measured values in Camporredondo et al. [61], are gathered in Fig. 4d.
To simplify the analysis and the comparison of the results, the sub-grid unresolved turbulence flow effects are excluded in the present study. Instead, a laminar liquid stress model is employed by the authors on a very fine numerical grid (so-called ''coarse DNS'' approach).
The total grid size used in the study included around 10 million volume elements. The calculations are performed on the HPC cluster using two nodes each with 28 cores (Intel Xeon E5-2690v4 2.60 GHz) and 128 GB RAM. Time integration step in the simulation is 10 -4 s to keep maximum Courant number below 1. For the time integration, the PIMPLE algorithm (a merge of PISO and SIM-PLE) is utilized to stabilize the numerical procedure using under-relaxation to get the required convergence level.

Simulation results
The developed model considering liquid flow, solidification and viscoplastic stress in the two-phase region was used to perform a simulation of the thin slab continuous Fig. 3 Temperature dependency of solid-phase fraction. T sol ¼ 1755:5 K; T liq ¼ 1797:9 K; T cast ¼ 1825 K Fig. 4 Simulation domain and setup details. a 3D geometry of thin slab mold including submerged entry nozzle and slab part; b details of numerical grid; c geometry dimensions (mm); d heat flux along vertical direction of mold adopted from Refs. [22,61] casting process based on the process conditions and geometry shown in Fig. 4a and in Table 1, as well as it was reported previously by Vakhrushev et al. [22].
The flow and solidification simulation results are shown in Fig. 5, and the instantaneous velocity magnitude distribution in Fig. 5a is displayed along with temperature field and shown in Fig. 5b. As it is mentioned previously, the ''coarse DNS'' approach was applied on a numerical mesh with high cell resolution; as shown in Fig. 5a, the highly turbulent flow structures are obtained during simulation, which are further employed for solving coupled heat transfer problem. Figure 5 shows the corresponding solution of the energy equation: the temperatures are scaled between the liquidus and the pouring temperature from Table 1 to provide better contrast and to distinguish the fresh and cold melt regions. The details of the temperature distribution in the vicinity of the narrow walls of the mold are analyzed in the zoomed Fig. 5c where the dark blue region corresponds to the growing solidified shell. The hot fresh melt jet, impinging and going along the mush in both vertical directions, is detected as well, damping the shell growth, as it was previously reported by Wu et al. [58].

Validity of viscoplastic regime in simulated domain
The viscoplastic deformation law was assigned in the presented study for the solidified shell description.
However, a question arises, to which extending it reflects the real mechanical behavior of the cast steel without violating it. Thereby, the validity of the applied Norton-Hoff viscoplastic law is investigated next according to the criterion reported in Rappaz et al. [46]. As mentioned previously, the viscoplastic behavior is postulated to be valid for the alloys above two thirds of their absolute melting temperature. Applying for the presented study, this range corresponds to the temperature values above the T sol Â 2 3 % 1170:33 K. The surface temperature distribution is analyzed and shown in Fig. 6. The color scaling is performed in such a way that the critical viscoplastic regime temperature values (1170-1200 K) are reflected by the blue colors that allows to clearly detect them. As one observes, only limited part of the domain along the corners lies below the viscoplastic temperature range; it starts at 1710 mm below the mold exit (see Fig. 6).
The temperature evolution along centerline of the wide and narrow face of the mold and strand part of the slab is shown in Fig. 7. One can see that at the mold's exit, the cooling curve reaches 1510 K at the wide and 1610 K at the narrow face. At the location of 4 m below meniscus, the surface temperature drops to 1460 K (wide side) and 1450 K (narrow side) correspondingly. Thus, the estimated centerline surface temperatures are significantly higher than a critical value for the pure viscoplastic law. More significant slab surface temperature decrease is detected during performed numerical studies at the strand's corner regions where it goes beyond the viscoplastic limit of 1170.33 K. As referred to the literature [46], the elastic behavior of the material is observed below one third of the melting temperature, which corresponds to 585.17 K in the present study. However, the minimum temperature observed for the simulation domain including 879-mm mold and the 3-m strand in the casting direction is higher than 1000 K. While insignificantly catching the transition from elasto-viscoplastic regime, the simulation results remain dominant in the frames of the pure viscoplasticity hypothesis, thus confirming the utilized stress model to be valid in the present study.

Solid shell thickness
The solid shell thickness was compared (see Fig. 8) between two cases: case I uses the precalculated solid velocities based on a fully elastic traditional approach disclosed elsewhere [22]; case II employs a viscoplastic model from this paper to treat the solid shell. It should be mentioned that the constant material properties are used in both cases. The validity of this assumption is discussed in the next section.
Firstly, the comparison along the wide face of the slab is shown (see Fig. 8a, b). Corresponding numbers on the isolines show the shell thickness in millimeters. It is observed that only slight difference occurs in the casting (vertical) direction. However, it is more pronounced in the transversal direction: in the referenced Case I, the shell is thicker in the central part of the wide face, whereas for the results using new approach, it is thicker near the narrow side.  The deviation between the predicted shell thickness can be described by recalling Eq. (29), based on the previously derived formulations [22,58], and the latent heat advection using the mixture velocity (as in Eq. (28)). Comparing these two formulations, the corresponding dference DS e is obtained in the following form: According to Eq. (30), DS e is proportional to u r . As a result, in the updated model formulation, used here, slightly more of the latent heat (released within the mushy zone during solidification) is advected along casting direction by the dendritic structure as compared to the case of the separately precalculated solid velocities. Thus, the growth of the solid shell is damped, and some deviation in shell thickness is observed, as shown in Fig. 8c, d for the slab's horizontal cross section at 2000 mm below meniscus level. The shell thickness at the corresponding sections A-A in Fig. 8 is 20.5 mm for the model with predefined solid velocities (Fig. 8c) and 18.1 mm for the viscoplastic deformation model (Fig. 8d) with the 13.3% difference. For the narrow face distribution (section B-B), the thicknesses correspondingly are 14.8 and 13.6 mm been measured at the center of the slab's narrow side, showing lower deviation of 8.8% between the models than at the wide side.

Model validation
The presented modeling results for the solidified shell thickness were verified against other simulation works [22,62,63] and the breakout measurements [61].
The solid shell thickness for the cases with the precalculated solid velocities using temperature-dependent properties (Fig. 9a, previous study [22]) and with the constant material data (Fig. 9b) is displayed against the industrial measurements [61] together with the presented viscoplastic model results (Fig. 9c). The mushy zone profiles are plotted along the vertical line located at 750 mm from the wide face center, where the shell thickness is least disturbed by the hot jets. It must be clarified that the proper consideration of the operation delay after the breakout and the partial solidification of the alloy, entrapped in the mush, during the drainage should be explained when comparing to the breakout shell [64].
Considering these uncertainties, the simulated shell thickness shows a good agreement with the measured one for the presented variety of the numerical setups (see Fig. 9). The experimental data lie between the iso-lines of .01 and f s = 0.30. As seen from the plots, the evolution of the shell can be well-tracked. The modeling results presented elsewhere [61][62][63] show the similar tendencies.
It is worth mentioning that a slight difference is observed between the case with temperature-dependent properties (Fig. 9a) and simulations with the constant ones (Fig. 9b, c). Moreover, there is a small deviation in the mushy zone distribution between the Case I in Fig. 9b and Case II in Fig. 9c due to the different withdrawal models, as discussed previously. However, all mismatches are minor. Thereby, the model verification is satisfying; the assumption of the constant material properties can be accepted. The material properties variation based on the temperature field, incorporation of the solid fraction in the constitutive model parameters, and appearance of the bulk viscosity due to the mush compressibility [31] are essential topics for the future studies.
As already shown, the verification data for the CC process mostly come after costly breakouts. However, in the future, we are looking for the experiments such as very recently published by Balogun, O'Malley and co-authors [65] combined with the withdrawal of the semi-solid shell.

Viscoplastic stress in solidified shell
Two viscoplastic stress components acting correspondingly along transversal and casting directions are investigated as shown in Fig. 10. The compression and tension zones in the solidified shell are indicated along the iso-surface of the 95% solid fraction by blue and red colors correspondingly. The iso-surface is enlarged 4 times in thickness direction for better resolution of the results. Strong tension/compression areas are detected especially along the curved funnel part of the mold cavity. It was found that the shell is under tension horizontally at the central and outer parts of the funnel (check for R vp xx \0 in Fig. 10) as well as vertically at the low part of the funnel, after it has been straightened in the parallel part of the mold (blue zone for R vp yy in Fig. 10 Fig. 9 Verification of modeled solid shell thickness (750 mm from wide face centerline) against breakout measurements [61], and cases with precalculated solid velocities using temperature-dependent properties [22] (a) and constant material properties (b), results for viscoplastic solid shell model with constant material properties (c) Fig. 10 Transversal R vp xx (left) and vertical R vp yy (right) components of viscoplastic stress tensor R vp horizontal transition zones from the curved funnel and the parallel mold plates.
As shown in Fig. 10, the critical tensile values in the simulation lie in the range of 250 kPa. The ultimate tensile stress (UTS) for the steel alloys is strongly temperatureand morphology (read composition)-dependent. As reported by Outinen and Mäkeläinen [66] and recently reviewed by Wang and Lui [67], the strength properties of the steel drop from the room temperature level of 500-1000 MPa to 5%-10% at the high temperatures above 1273 K. By extrapolating the experimental data [66,67] up to 1723 K, the lowest range for the UTS becomes less than 0.5-1.0 MPa. That estimation makes the simulation values close to the critical limit. Thereby, at the upper part of the mold, the hot jet impingement and the local remelting of the solidifying skeleton can lead to the strength weakening till the rapture occurs.
However, it is not fully valid to compare the presented results with the real steel properties, since the constant viscoplastic model parameters were used for the apparent viscosity l vp .

A Norton-Hoff-type viscoplastic model was applied
for the two-phase solidifying region to study the key phenomena in the curved mold during thin slab casting. A single FVM framework was used with the enhanced numerical treatment of the nonlinear viscoplastic stress tensor according to the improved both side diffusion method. 2. A full 3D simulation of the industrial thin slab caster was performed including SEN region, the funnel-type mold, and the secondary cooling zone. The melt flow, heat transfer and viscoplastic stresses in the solidifying shell were calculated for the alloy with compositions of C 0.06, Ni 0.1, Mn 0.13, Si 0.1, Cu 0.08, Al 0.035, P 0.015 and S 0.012 within the typical casting range. 3. The viscoplasticity was found dominantly valid in the calculation domain except small regions catching the elasto-viscoplastic transition. The tension and compression of the solid shell was detected along the mold's funnel within the corresponding straightening and bending zones, where the local reheating or nonuniform solidification can cause the crack formation leading later to a breakout. 4. The modeling results were verified against the breakout shell measurements and the numerical studies published elsewhere, showing a good agreement. 5 The described approach can be integrated with the free-surface tracking models, and include the magnetohydrodynamic effects, etc. The aim of the future development is to estimate the breakout risks using modern quality indices and defect detection techniques to be directly applied for the industrial process enhancement.