On longitudinal electromagnetic stirring in the continuous casting of steel blooms

Recent work highlighting an anomaly in the modelling of rotary electromagnetic stirring (EMS) in the continuous casting of round steel billets is extended to the case of longitudinal stirring for rectangular blooms. An earlier, still often-cited, model forms the basis of the current analysis, which uses asymptotic methods on the three-dimensional (3D) Maxwell equations and demonstrates how the earlier result for the components of the Lorentz force is but a particular case of a more general form. Time-dependent 3D computations using finite-element methods are also performed to verify the validity of the asymptotic analysis, and the relevance of the results to modulated EMS is noted.

However, recent work by Vynnycky [8] revisited the problem of rotary EMS applied to round-billet continuous casting and found that the method used originally in [2,7] to determine the components of the Lorentz force led to a non-unique solution; this was a consequence of the fact that the normal component of the induced magnetic flux density, rather than the tangential ones, was prescribed as the boundary condition. Moreover, since the normal component was prescribed in models for the case of longitudinal stirring for rectangular blooms also [3][4][5][6], it is natural to expect non-uniqueness in those models too. Furthermore, since the expressions for the components of the Lorentz force given in [3][4][5][6] are still frequently used [9][10][11][12][13][14][15][16][17], it is clear that a resolution of the issue is still timely, especially in view of modern-day interest in modulated EMS [18][19][20]. In this case, magnetic fields of different frequencies are applied and it is the intention that the resulting Lorentz force should have a constant time-averaged part and a time-varying one; it goes without saying that posing the correct boundary conditions for the magnetic field is vital for obtaining meaningful results from modelling. Furthermore, although this paper revisits an earlier EMS model in a similar spirit to [8], the details turn out to be more extensive, but also more subtle. In the case of rotary EMS applied to round-billet continuous casting, the geometry is axisymmetric and considerable progress can be made analytically. In the current paper, on the other hand, the geometry is genuinely three-dimensional (3D), and analytical progress via the use of asymptotic methods can only be made because we identity a suitable asymptotically small parameter in which we can conduct the analysis. In addition, we find that the form of the boundary condition that is responsible for stirring influences the form of the asymptotic expansion that is needed. Last, but certainly by no means least, we verify the accuracy of the analysis by comparison with results from 3D computations.
The structure of the paper is as follows. In Sect. 2, we formulate the problem mathematically, and in Sect. 3 we non-dimensionalize the governing equations. The main body of the analysis, in terms of an asymptotically small parameter, is given in Sect. 4 , whereas Sect. 5 outlines the numerical method used to solve the fully 3D transient problem. Results are given in Sect. 6, and conclusions are drawn in Sect. 7.

Governing equations
We consider, as shown in Fig. 1, a linear travelling stirrer situated at y = 0, −h/2 ≤ z ≤ h/2, 0 ≤ x ≤ l, operating on a molten steel region with cross-section 0 ≤ y ≤ b, −h/2 ≤ z ≤ h/2, as shown in Fig. 2. For simplicity, we neglect unnecessary details such as the existence of a solidified shell around the rim of the cross-section, the presence of a mould wall on the other three sides of the steel, as would occur in mould electromagnetic stirring (M-EMS) or an air gap between the solidified shell and the stirrer, as could occur in M-EMS or strand electromagnetic stirring (S-EMS); instead, the initial focus is on formulating the simplest well-posed problem. Moreover, as suggested by Fig. 2, this constitutes having to consider the air surrounding the other three sides of the principal region of interest.
As already motivated in the introduction, we consider the solution of Maxwell's equations in the magnetohydrodynamic (MHD) approximation, which consist of: -the magnetic field constraint, where B is the magnetic flux density vector; -Ampere's law, Fig. 1 Schematic of the arrangement for the longitudinal stirring of blooms and billets. The figure is not drawn to scale; typically, b and h are more or less equal, and l is around four times larger -Faraday's law, where t is time and E is the electric field; -Ohm's law, where σ is the electrical conductivity of the medium in question, which we take to be constant, and is v the velocity vector, which we will simply set to zero.
Manipulating (2.2)-(2.4), we arrive at In addition, F x , F y and F z are used to calculate the time-averaged components of the Lorentz force-F x ,F y and F z -viā where the integrals are taken with respect to time over one oscillation period, 2π/ω, with ω as the angular frequency of the current, which is related to the frequency of the current, f , by ω = 2π f .

Boundary and interfacial conditions
One of the key points of this paper concerns what should be the boundary condition at y = 0, |z| ≤ h/2 which is in effect responsible for the Lorentz force. Dubke et al. [3] prescribe B y at −h/2 ≤ z ≤ h/2, y = 0, and their results, based on this condition, have been in use ever since; their approach is recapped in Appendix A. In this sense, prescribing the normal component of the magnetic field is consistent with earlier work for rotary EMS acting on a round billet [2,7]. However, Vynnycky [8] recently re-analyzed the same problem and found that prescribing the normal component would lead to a non-unique solution for the components of magnetic flux density vector; on the other hand, prescribing the tangential components would lead to a unique solution. This observation, for a round billet, suggests that applying a boundary condition on the normal component for a rectangular bloom or billet will lead to a non-unique solution also. In turn, this suggests instead that one should prescribe the tangential components of B, and we will therefore set at y = 0, where λ is the wave vector given by λ = π/ p, with p as the pole pitch, i.e. the distance from north to south pole in the inductor. The motivation for choosing the forms given in (2.10) and (2.11) is that they constitute a travelling wave in the x-direction, as was also considered in [3], and also give the most general possible dependence in the z-direction; as a result, we are able later on to relate the Lorentz force directly to what is chosen for α (z) and β (z) , as shown in Sect. 4.1.2.
As regards the interfaces between air and steel at y = b, |z| ≤ h/2 and |z| = h/2, 0 ≤ y ≤ b, we apply the usual electromagnetic interface conditions. These are the continuity of the tangential components of the magnetic field strength and the normal component of the magnetic flux density, which are written as respectively, where t is the unit tangential vector and n is the unit normal vector at these interfaces, and the continuity of the tangential components of the electric field, which is written as In (2.12) and (2.13), [ ] + − denotes the difference in the value of a function in the air (+) and in the steel (−). In the air far from the steel, we should also expect B → 0 as r → ∞, (2.14) where r = x 2 + y 2 + z 2 .

Non-dimensionalization
It is ultimately more instructive to non-dimensionalize the equations, and we do this through and, for k = x, y, z and K = X, Y, Z , where B 0 is the characteristic magnetic flux density, and η s and σ s are the steel magnetic permeability and electrical conductivity, respectively; note also that the corresponding quantities for air, η a and σ a , will also be required later on. Now, for 0 ≤ X ≤ 1, 0 where δ = b/l and Ω = f σ ηb 2 ; however, note that the values for Ω for steel, Ω s , and air, Ω a , are different, but can be calculated from the data in Table 1. Also, for later use, note that and The boundary conditions (2.10)-(2.14) are then: at Y = 0, (3.10)

Analysis
Since δ ∼ 1/4, it is perhaps tempting to consider δ 1 straightaway and to pose regular perturbation expansions for B X , B Y and B Z with δ as the expansion parameter, as would be suggested by Eqs. (3.3) and (3.4). However, the character of the problem turns out to be different to this. First of all, Eq. (3.3) and the boundary conditions in (3.7) imply that the solutions for B X , B Y and B Z will be of the form Thus, in (3.3), the first term will now be of order δΛ, i.e. λb; we see from Table 1 that this non-dimensional quantity is no longer necessarily small (λb ≈ 13.1) and that it does not even contain the length of the stirrer. In addition, the apparently O δ 2 terms in (3.4) are actually O λ 2 b 2 . Thus, to proceed more appropriately, we will first set where Re denotes the real part of a complex number and i = √ −1.
where Δ = λb. In terms of B X , B Y and B Z , the boundary conditions are now: at Y = 0, We can note at this point that it is inconvenient to have an interface condition, (4.5), in terms of E when the rest of the problem is posed in terms of B. We will therefore write out the equations in (4.5) in their full forms in terms of B X , B Y and B Z , which gives: at Y = 1, |Z | = H/2, Now, in the absence of any further simplification for general values of Δ-other than Δ 1, which appears to be not of technological interest anyway-it would be necessary to resort to a numerical solution. Ironically, however, it appears to be simpler to solve (3.3) and (3.4), subject to (3.7)-(3.10), numerically than it is to solve (4.1)-(4.6) numerically, even though the latter is a two-dimensional steady state system of equations, whereas the former is three-dimensional and transient. By "simpler", we mean that there exists commercial software such as the finite element-based Comsol Multiphysics which we will use later anyway, with readily available modules to solve the former, but not the latter, with the main difficulty being the implementation of the interface conditions that arise out of (4.5), i.e. (4.8) and (4.10).
Thus, instead of trying to solve (4.1)-(4.4) and (4.6)-(4.10) numerically, we will first consider these equations analytically for the case Δ 1, which is indeed the regime pertinent to the data in Table 1. Note also that, although the data given are the same as those used in [3], it is quite typical for most casting machines that apply longitudinal EMS; the only significant difference would be that the steel breadth is as small as the width, but this does not affect the value of Δ, which does not depend on the breadth anyway.
In this subsection, we first determine the solutions for B X , B Y and B Z ; thereafter, these are used to constitute respectively. These suggest that, at leading order in ε, which would clearly contradict (4.3) and (4.4) for |Z | < H/2. Thus, since a small parameter multiplies the highest order derivatives in Eqs. (4.11) and (4.12), there must be boundary layers in the problem. The most obvious location for these is at Y = 0 for |Z | < H/2. Anticipating that Y ∼ ε , we set Y = εỸ , upon which (4.11) and (4.12) become respectively. Note that we have replaced Ω in (4.12) by Ω s , since the analysis to follow concerns the steel region.
Observe also that we retain Ω s as a dimensionless O (1) parameter since, on using the values in Table 1, Ω s ≈ 2.8. Setting with the boundary conditions being (4.3) and (4.4), and which must come from matching to the solutions in (4.13), taken at leading order. So, Y , which then clearly satisfies (4.18) for k = Y , indicating self-consistency. We can also observe that (4.20)-( 4.22) will not satisfy (4.9)-(4.10), for which boundary layers where |H/2 ± Z | ∼ ε are required; we omit the details of the analysis for those, as they are of little importance for the main results that follow, although we do return to a discussion of these in Sect. 6.1.
For later use, we need the second and third terms, i.e. at ε 1 and ε 2 , in the expansions also. At ε 1 , we have whence where the prime denotes differentiation with respect to Z . At ε 2 , we have (4.38) Furthermore, although the above analysis has been for the boundary layer near Y = 0 for |Z | < H/2, it is straightforward to construct a composite solution for the entire steel region. This is done in the usual away by summing the inner and outer solutions and subtracting the common part [21]; in view of (4.13), this is effectively the same as replacingỸ by Y/ε in (4.20)-(4.22 ), (4.27), (4.29) and (4.36)-(4.38) above. Note also that Eqs. (4.7)-(4.10) were effectively not used at all in this analysis, although the equations that gave rise to them, (3.9), will be used in the numerical solution that follows later.
Since differentiating with respect to X will bring out a factor Λ, and since δΛ = 1/ε, this suggests the following expansions: Hence, we find, ϕ = 2πτ − λX , An unexpected difficulty now occurs if we haveᾱ = 1,β = 0, which we note is not a pathologically contrived counterexample, but in fact corresponds to the situation originally considered by Dubke et al. [3]: in this case, all three components of J,F andF are zero at leading order, where moreover, the expansion given by Eq. (4.16) turns out not to be of the correct form. We therefore consider this case separately.

4.2ᾱ ≡ 1,β ≡ 0
It is evident from the outset that B Z ≡ 0, as is found by noting that the equations for B Z -(4.12) for k = Z , subject to (4.4) and (4.13)-decouple from the others and have a trivial solution; thus, we only need consider (4.14) and (4.15) for k = X, Y , which are now Setting, instead of (4.16), subject to which leads to and hence and thence

B
(1) With the current choice ofᾱ andβ, we will have J X , J Y ≡ 0, as is evident from (3.5) and the fact that B X and B Y do not depend on Z , and it only remains to determine J Z ; at leading order, this will be of O (ε) , and we will need to set Now, it is evident that F X , F Y ∼ O (ε) , F Z ≡ 0, so that on setting whence, we will find that so, we would need to go to an even higher order in ε 2 to find the leading-order term forF Y . Thus, at ε 4 , we have Similarly, we obtain    (4.95) and thencē Note that, in (4.95), we have written the expansion to two terms, in order to include the first term that provides a non-zero contribution to (4.96); in (4.91)-(4.94), we have written just the first term.

Reconstruction of the earlier solution [3]
It is evident that we can use the results in Sect. 4.1 to obtain a solution for which B y behaves at y = 0 in the same way as the solution of Dubke et al. [3]. If we replace (2.10) by we should now set where B k will now satisfy the same equations as in Sect. 4.2. Omitting all of the intervening algebraic details, we will obtain, instead of (4.91)-(4.95), and thencē Consequently, although the forms in (4.98)-(4.102) are different to those in (4.91)-(4.95), the expressions forF x andF y are identical.

Numerical solution
In order to ascertain the correctness of the asymptotic analysis, Eqs. (2.6) and (2.7), subject to (2.10)-(2.14), were also solved numerically. For this, the Magnetic Fields module of the Comsol Multiphysics software mentioned previously was used. Second-order Lagrangian quadrilateral elements were employed on a mesh having approximately 70,000 elements, corresponding to around 1.4 million degrees of freedom. The computational mesh that was used was decided upon after a grid independence study, although, in the interest of brevity, the results of this are not presented here; the fact that we have obtained good agreement with the asymptotic results, as will be seen in Sect. 6, means that we can be sure that we have computed accurately enough. As regards the size of the computational domain, its outer edge was taken to be a half-cylinder of radius 0.4 m and height 1.2 m. The geometry and the mesh are shown in Fig. 3. For the numerical integration, adaptive time-stepping was used and the convergence criterion at each time step was taken as where (U i ) is the solution vector corresponding to the solution at each time step, E i is the estimated error in the latest approximation to the ith component of the true solution vector, A i is the absolute tolerance for the ith degree of freedom, and R is the relative tolerance; for the computations, R = 10 −3 , A i = 10 −4 for i = 1, . . . , N dof were used, where N dof is the number of degrees of freedom. Note that all computations were carried out on an architecture with an Intel® Core TM i9-7940X CPU @ 3.10 GHz×14 processor with 32 GB RAM. The typical computation time was of the order of 2-3 h; this constitutes running the model for three periods of oscillation, corresponding to 0.06 s, in order to minimize the effect of the initial transient, so as to obtain appropriate data for computing the components of the time-averaged Lorentz force.

Results
In this section, we consider first results for the case when α = B 0 , β = 0, and then for a more general case when both α and β are functions of z.  [3], our asymptotic solution and the full numerical solution. Note that we have plotted the asymptotic solution to two terms, so as to demonstrate how well it agrees with the full numerical solution; since ωt − λx is an even multiple of π (16 in fact) for these values of x and t, we would have had just the zero solution to one term for B x . Moreover, Fig. 4a in particular illustrates that the asymptotic solution is better able to reproduce the numerical solution than the analytical solution from [3]. Also of interest is to see whether the full numerical solution displays the boundary layer indicated by the analysis and to see the behavior near |z| = h/2. This is shown in Fig. 5 which shows the contours of B y for spacings of 0.002 T for 0.002 T≤ B y ≤ 0.02 T, as given by the numerical solution for t = 0.06 s, x = −0.6 m; the boundary layer is evident, as is its breakdown near |z| = h/2. Figure 6a and b compare F x and F y , respectively, as functions of y at z = 0 m, t = 0.06 s, x = −0.6 m, as obtained by the three approaches. This time, the discrepancy between the numerical solution and the solution from [3] is evident in Fig. 6b, whereas the two-term asymptotic solution once again agrees well in both cases with the full numerical solution. Figure 7a and b compareF x andF y , respectively, as functions of y at z = 0 m, x = −0.6 m, as obtained by the three approaches. The results of all three methods agree with each other, and the overall situation is reminiscent of the corresponding situation in rotary EMS, where an earlier model [2,7] provides the correct solution for the timeaveraged Lorentz force components, even though the time-dependent solution is incorrect [8]. Note that although the profile forF y is but the leading term in the asymptotic expansion, it is actually constituted from the third terms in the expansions for B x and B y .
Since it reasonable to expect symmetry about z = 0 m, we take for α and β the simplest possible profiles for which neither is constant, i.e. quadratic powers in z. In particular, we take         -the numerical solution has already achieved periodic behavior, since it agrees well with the analytical solution into which periodicity is already inbuilt; -the frequency of F x is twice that of B x .
We comment also that we have not shown the solution at z = 0, as was done in the case where α = B 0 , β = 0, as it will be identically zero for these profiles of α and β. Moreover, the small discrepancy between the three-term asymptotic solution and numerical solution in Fig. 8 is most likely due to fact that, although the numerical solution is mesh-independent with respect to B X , B Y and B Z , a still finer mesh would be required to achieve mesh-independence for their derivatives; of course, the components of the Lorentz force contain the derivatives. However, as we were already at the operating limit as regards RAM on the computer that was used, we did not attempt a solution on an even finer mesh.
Finally, in Fig. 9, we showF x ,F y ,F z vs. y at z = 0.0875 m, x = −0.6 m. Clearly, these have been computed correctly. Note here that a 3-term expansion is required to obtain a non-zero leading term for all three force components, because of the choice of α and β; in particular,F (1) X = 0 from (4.54) so that the first non-zero contribution toF X is at O ε 2 .

Conclusion
This paper has revisited an early, yet still often-cited, mathematical model for longitudinal electromagnetic stirring in the continuous casting of rectangular steel blooms [3][4][5][6]. As indicated in Appendix A, there are mathematical inconsistencies in the way the model equations were originally posed, leading to the possibility of non-unique solutions; our analysis has resolved this issue. However, it appears that the original model worked well enough because of the parameter regime in which it was used. Furthermore, our analysis identified the key dimensionless parameter, Δ, to be the product of the bloom width, the reciprocal of the inductor pole pitch and π. In typical applications, this parameter is large, which leads, physically, to a boundary layer near the stirrer surface and, mathematically, a singular perturbation analysis in terms of the small parameter ε = 1/Δ. Moreover, it was found that the expressions found originally in [3][4][5][6] were not of the most general form that can occur: in those, the magnitudes of the Lorentz force components depend on the oscillation frequency, whereas in general they do not. In addition, in the most general case, the magnetic flux density components expand in powers of ε, whereas they expand in powers of ε 2 in the original case. The subsequent results of this analysis were compared with those of transient 3D computations of the Maxwell equations using finite element methods and good agreement was found.
The current methodology may now be extended to derive corresponding expressions for the components of the Lorentz force in the continuous casting of slabs, wherein stirring is horizontal, i.e. perpendicular to the casting direction [22]; with reference to Fig. 1, this will mean that a wave of the form cos (ωt − λz) , rather than cos (ωt − λx) , is applied. Morever, and again with reference to Fig. 1, it is the case than b/ h 1, which means the possibility of simplification via a second asymptotically small parameter.
Finally, we point out that determining how the magnetic flux density is perturbed, i.e. at O (ε) or O ε 2 , may be of significance because it will also be perturbed by the melt velocity field. As this perturbation is characterized by the magnetic Reynolds number, R m , which is itself typically small, the size of R m in relation to ε and ε 2 will determine whether the next-order perturbation is due to the velocity field (R m ε or R m ε 2 , depending on the boundary conditions) or the inductor (R m ε or R m ε 2 ).

Appendix A: Solution from Dubke et al. [3]
Before summarizing the solution from [3], we re-cap how it was obtained. All quantities were assumed to be independent of z, and B z was taken to be zero. The governing equations [A-3] and [A-8] in [3] were then taken to be, for y ≥ 0 and −∞ < x < ∞, respectively; for simplicity, since it does not affect the gist of the argument, we have omitted in (A.2) an advection term that was present in [A-8]. A primary objection already at this stage is that there should also be an equation for B x which has the same form as (A.2). Thence, assuming that B y = A (y) e i(ωt−λx) , prescribing B y at y = 0 and assuming that B y → 0 as y → ∞, leads to a unique solution for A (y) from (A.2); the resulting expression for B y is then substituted into (A.1), which is then integrated with respect to x to obtain an expression for B x . However, the correctness of this step is also debatable, since it should introduce a new function of y and t, which was not accounted for in [3]. Furthermore, even if such a function had been (correctly) introduced, it is not clear how it would be determined, since there are no boundary conditions on B x that could help determine it. Now, with φ = ωt − λx − (Imγ ) y, the key expressions are where Im denotes the imaginary part of a complex number and γ = σ s η s ωi + λ 2 , (A.8) so that with F x = σ s ωB 2 0 e −2(Re γ )y 2λ ,F y = σ s ωB 2 0 (Im γ ) e −2(Re γ )y 2λ 2 . (A.10) We may note at this point that the expressions in (A.10) are not identical to those in (4.103), even though the agreement in Fig. 7 was very good. This is explained by noting that, on using the data in Table 1, λ 2 /σ s η s ω ≈ 9.7, which justifies considering λ 2 σ s η s ω; thus, the expressions in (A.9) can be approximated as Re γ ≈ λ, Im γ ≈ σ s η s ω 2λ , rendering (4.103) and (A.10) identical.