Introduction

The cast-iron foundry industry produces cast parts for automotive, agriculture, transportation, energy, aerospace, manufacturing industry, etc., with 90 pct of manufactured products containing cast-iron parts. All foundry processes generate a certain level of rejected parts, which are irrecoverably lost, being closely related to the type of casting and processes used and the equipment available. As quality demands from end-users of castings increase, it is essential that cast-iron technology move forward, together with green manufacturing as a first step towards sustainability.

Extensive work has been carried out by the foundry community to minimize casting defects such as porosity, slag, and clogging, but few literature reports consider casting metal using more versatile melt control technology.

The flow control system in a tundish has a significant impact on the quality and level of rejection when pouring casts. If the metal does not flow in a consistent stream, casting defects may result from oxidation, air entrapment, and erosion of the pouring mold, among other causes.

In the cast-iron industry, it is usual for molten metal to be transfered from a ladle with a lip-axis pour design.[1] This method is fast and reliable, but inaccurate unless a sophisticated control system is used.[2,3,4] In recent years, our group has been working on technologies that aim to transform the conventional (batch-by-batch) foundry process into a flexible (mold-by-mold) process. This requires that the ladle become a furnace as well, that the tundish be eliminated from the process, and that the molten metal be poured directly into the mold. Moreover, it is convenient for the pouring flow rate to be constant to obtain uniform refilling of the mold. Such a flexible mold-by-mold process is characterized by, firstly, combination of the melting, treatment, and pouring processes into a single cast-iron production cell. Also, this requires integration of an artificial-intelligence-based control system to monitor local structures, phases, and mechanical properties to guarantee high-quality casting in the foundry. Finally, a robot cell will be in charge of the metal finishing process. This also reduces the melt temperature and transport while improving validation of the cast pieces.

In recent decades, extensive effort has been invested in development of CFD simulation methods for application in the casting industry, especially for multiphase flows.[5,6,7] In the particular case of simulations of pour tilt casting, Prakash et al.[8] used the smoothed particle hydrodynamics method to simulate the oxidation process during furnace emptying, and Kuriyama et al.[9] used the software Flow 3D to optimize the pouring velocity in aluminum tilt casting. Pauty et al.[10] carried out two-dimensional (2D) numerical simulations of the liquid transport and thermal convective velocity in a tilting furnace, finding good agreement with water experiments, while admitting that the extension of the method to three-dimensional (3D) flow simulations of molten metal would require a dramatic increase of computational resources and time. Fortunately, advances in computer technologies in recent decades have made this kind of complex 3D fluid flow simulation more affordable; For instance, Davila et al.[11] computed the 3D flow during drainage of a ladle to understand vortex mechanisms in this process.

We present herein a study of the impingement location of liquid poured from a tilting ladle having a lip spout with different curvatures. It has been reported that a curved lip provides better jet control and reduces air entrainment.[12] The purpose of this work is to present a numerical method with low computer power requirements to calculate the liquid trajectory in such a pour tilt casting process. The results could be used to reduce jet dispersion to minimize liquid spilling and reduce sprue size.

Section II develops the theoretical background to calculate the angular velocity required to achieve constant flow rate. The results are the angular velocity and tilting angle vs time. In the next stage, these data are used in CFD simulations. The numerical approach is described in Section III. The simulations were performed using the OpenFOAM® toolbox, which is based on the finite-volume method (FVM). As this is a dynamic process with moving boundaries, the habitual strategy is to use a dynamic mesh. Another option is the FAVOR® technique.[13,14] These two methods are very flexible and allow independent movement of several objects in the domain. However, we propose an alternative method where neither the mesh nor object move, but the gravitational acceleration. This method is simpler and faster than the dynamic mesh or FAVOR® approach, but has the disadvantage that only one moving object is allowed (in this case, the tilting ladle) and the postprocessing is more complicated. The results are presented in Section IV; finally, in Section V, these results are discussed and conclusions drawn.

Theoretical Background

The theoretical calculation of the flow rate poured from the furnace is based on conservation of mass of liquid. The flow rate of liquid is equal to the time derivative of the mass inside the ladle, thus

$$\begin{aligned} \frac{\partial \phantom {}}{\partial t}\int _{\text {VC}}\rho \mathrm {d}V=-\rho q, \end{aligned}$$
(1)

where q is the flow rate through the spout. Since the density \(\rho \) is constant, it can be eliminated from Eq. [1]. Assuming that the free surface of the liquid always remains horizontal, the volume of liquid can be divided into two regions: the volume above a horizontal plane passing through the lower point of the liquid exit surface, \(V_{\rm t}\), and the volume below this plane, \(V_{\rm b}\). The height of \(V_{\rm t}\) is h, so that

$$\begin{aligned} V_{\rm t}=hA_{\text{FS}}, \end{aligned}$$
(2)

where \(A_{\text {FS}}\) is the area of the free surface. Figure 1 depicts a scheme of the pouring liquid with arbitrary tilting angle \(\theta \).

Fig. 1
figure 1

Scheme for pouring liquid for arbitrary tilting angle \(\theta \) (not to scale)

Full size image

Thus, Eq. [1] yields

$$\begin{aligned} \frac{\mathrm {d}V_{\rm t}}{\mathrm {d}t}=-q-\frac{\mathrm {d}V_{\rm b}}{\mathrm {d}t}, \end{aligned}$$
(3)

which can be expressed in terms of the angle \(\theta \) instead of time, thus

$$\begin{aligned} -\omega \frac{\mathrm {d}V_{\rm t}}{\mathrm {d}\theta }=-q+\omega \frac{\mathrm {d}V_{\rm b}}{\mathrm {d}\theta }, \end{aligned}$$
(4)

where \(\omega =-\dfrac{\mathrm {d}\theta }{\mathrm {d}t}\) is the angular tilting velocity of the ladle. Using Eq. [2] in [4] yields

$$\begin{aligned} A_{\text {FS}}\frac{\mathrm {d}h}{\mathrm {d}\theta }+h\frac{\mathrm {d}A_{\mathrm {FS}}}{\mathrm {d}\theta }=\frac{q}{\omega }-\frac{\mathrm {d}V_{\rm b}}{\mathrm {d}\theta }. \end{aligned}$$
(5)

The lower volume, \(V_{\rm b}\), is known to be a cylindrical wedge and can be calculated based on geometric arguments if the spout lip is neglected (i.e., assuming that the ladle is a perfect cylinder). If the free surface does not cut the cylinder (ladle) bottom, the volume can be easily calculated as[15]

$$\begin{aligned} V_{\rm b}=\frac{1}{2}\pi R^{2}\left( H+H_{\rm s}\right) =\pi R^{2}H\left( 1-\frac{R}{H\tan \theta }\right) . \end{aligned}$$
(6)

If, otherwise, the plane cuts the cylinder bottom, the calculation is more complex. In the particular case in which the plane cuts the bottom circle across its diameter, known then as a cylindrical hoof, the volume is one-sixth the volume of the square prism in which the cylinder is inscribed, as shown by Archimedes[16]

$$\begin{aligned} V_{\rm b}=\frac{2}{3}R^{2}H. \end{aligned}$$
(7)

For the general case, the volume has to be calculated by integration, giving

$$\begin{aligned} V_{\rm b}=\frac{R^{3}}{\tan \theta }\left[ \pi \left( \eta _{H}-\eta _{0}\right) -\eta _{H}\arccos \eta _{H}+\frac{1}{3}\sqrt{1-\eta _{H}}\left( 2+\eta _{H}^{2}\right) \right. \nonumber \\ +\eta _{0}\arccos \eta _{0}-\frac{1}{3}\sqrt{1-\eta _{0}}\left( 2+\eta _{0}^{2})\right] , \end{aligned}$$
(8)

where

$$\eta _{0}=\frac{b_{0}}{R}-1,$$
(9)
$$\begin{aligned} \eta _{H}= \eta _{0}+\frac{H \tan \theta }{R}, \end{aligned}$$
(10)

where \(b_{0}\) is the height of the liquid on the opposite side of the cylinder (see Figure 2). For our case of pouring liquid, \(b_{0}=0\) and Eq. [8] reduces to

$$\begin{aligned} V_{\rm b}=\pi HR^{2}-\frac{R^{3}}{\tan \theta }\left[ \pi +\eta \arccos \eta -\frac{1}{3}\sqrt{1-\eta ^{2}} \left( 2+\eta ^{2}\right) \right] , \end{aligned}$$
(11)

where \(\eta =\dfrac{H\tan \theta }{R}-1\). Alternatively, one can use the following expression from Reference 15:

$$\begin{aligned} V_{\rm b}=\frac{1}{3}HR^{2}\left[ \frac{3\sin \phi -3\phi \cos \phi -\sin ^{3}\phi }{1-\cos \phi }\right] , \end{aligned}$$
(12)

where

$$\begin{aligned} \phi =\frac{\pi }{2}+\arctan \left( \dfrac{\eta }{\sqrt{1-\eta ^{2}}}\right) . \end{aligned}.$$

Both equalities, i.e., Eqs. [11] and [12], reduce to the Archimedes theorem, Eq. [7], in the case of the cylindrical hoof (i.e., \(\eta =0\) in Eq. [11]).

Fig. 2
figure 2

Scheme of partially full inclined cylinder

Full size image

The area of the free surface is trivial when the plane does not cut the bottom of the cylinder, being an ellipse with area

$$\begin{aligned} A_{\mathrm {FS}}=\frac{\pi R^{2}}{\sin \theta }, \end{aligned}$$
(13)

but is more complicated when the plane cuts the bottom. In this case, it is a segment of an ellipse, and, according to Hugues et al.,[17] can be calculated as

$$\begin{aligned} A_{\mathrm {FS}}=\frac{R^{2}}{\sin \theta }\left( \psi -\sin \psi \right) , \end{aligned}$$
(14)

where

$$\begin{aligned} \psi =2\arccos \left( -\eta \right) . \end{aligned}$$

The flow rate in Eq. [5] can be estimated using the expression for the flow rate over a circular weir[18]

$$\begin{aligned} q=0.3926C_{\rm d}\sqrt{2g}h^{3/2}D\beta ^{1/2}\left( \sqrt{1-0.2200\beta }+\sqrt{1-0.7730\beta }\right) , \end{aligned}$$
(15)

where \(\beta =h/D\) and \(C_{\rm d}\) is the discharge coefficient, which is also a function of \(\beta \). Nevertheless, a simplified expression[19] is usually applied:

$$\begin{aligned} q=KD^{n_{1}}h^{n_{2}}, \end{aligned}$$
(16)

where \(n_{1}=0.693\), \(n_{2}=1.807\), and \(K=1.598\) when S.I. units are used.

Given the flow rate poured, which is kept constant, Eq. [16] leads to the height of the free surface as

$$\begin{aligned} h=\left( \frac{q}{KD^{n_{1}}}\right) ^{\frac{1}{n_{2}}}. \end{aligned}$$
(17)

Thus, the angular velocity as a function of angle can be calculated from Eq. [5], provided that the flow rate q and the free surface height h are constant,

$$\begin{aligned} \omega \left( \theta \right) =\frac{q}{h\frac{\mathrm {d}A_{\mathrm {FS}}}{\mathrm {d}\theta }+\frac{\mathrm {d}V_{\rm b}}{\mathrm {d}\theta }}, \end{aligned}$$
(18)

where the derivative of the free surface area is calculated from Eqs. [13] and [14], and the derivative of the bottom volume is obtained from Eqs. [6] and [11].

This angular velocity can be expressed in terms of time, instead of angle, using

$$\begin{aligned} t(\theta )=\int _{\frac{\pi }{2}}^{\theta }\frac{\mathrm {d}\alpha }{\omega (\alpha )}. \end{aligned}$$
(19)

Combining Eqs. [18] and [19], one obtains the angular velocity vs time. The dimensionless angular velocity and time can be defined using a characteristic time \(\tau =\dfrac{\pi R^{3}\lambda }{q}\), with a form factor \(\lambda =\frac{H}{R}\), which is the time needed to empty the cylinder at flow rate q,

$$\begin{aligned} \begin{array}{c} \omega ^{*}=\omega \tau \\ t^{*}=\dfrac{t}{\tau }. \end{array} \end{aligned}$$

Figure 3 shows the dimensionless angular velocity, calculated from Eq. [18], against dimensionless time. Figure 4 shows the rotation angle vs dimensionless time.

Fig. 3
figure 3

Theoretical estimation of dimensionless angular velocity vs dimensionless time for constant flow rate

Full size image
Fig. 4
figure 4

Theoretical estimation of rotation angle vs dimensionless time for constant flow rate

Full size image

Numerical Approach

CFD simulations were carried out using the OpenFOAM® toolbox.[20,21,22] This open-source tool provides libraries, solvers, applications, and utilities to solve partial differential equations (PDEs) using the FVM. There are several methods for solving multiphase flows, but the most popular is the volume-of-fluid (VOF) approach,[23] in which a local liquid fraction \(\gamma \) is defined in each cell domain. This liquid fraction defines a cell as completely filled with liquid \(\left( \gamma =1\right) \) or gas \(\left( \gamma =0\right) \), or partially occupied by liquid and gas \(\left( 0<\gamma <1\right) \). The liquid fraction is transported with the flow as

$$\begin{aligned} \frac{\partial \gamma }{\partial t}+\nabla \cdot \left( \gamma \vec {u}\right) +\nabla \cdot \left( \gamma \left( 1-\gamma \right) \vec {u}_{r}\right) =0, \end{aligned}$$
(20)

where the last term on the left-hand side represents interface compression and avoids dispersion of interfaces, keeping a narrow zone of cells with \(\gamma \approx 0.5\). For this purpose, the velocity \(\vec {u}_{r}\) is artificially defined normal to the interface and pointing towards it.[24]

Along with the transport equation for the liquid fraction, the conservation of mass and momentum are solved, i.e.,

$$\begin{aligned} \nabla \cdot \vec {u}= & {} 0, \end{aligned}$$
(21)
$$\begin{aligned} \frac{\partial \rho \vec {u}}{\partial t}+\nabla \cdot \left( \rho \vec {u}\vec {u}\right) =-\nabla p+\nabla \cdot \left( \mu \left( \nabla \vec {u}+\nabla \vec {u}^{T}\right) \right) \nonumber \\ \quad+\rho \vec {g}+\int _{S_{i}}\sigma \kappa \vec {n}\delta \left( \vec {r}-\vec {r}_{i}\right) \mathrm {d}S, \end{aligned}$$
(22)

where \(\vec {g}\) is the gravitational field, \(\sigma \) is the surface tension coefficient, \(\kappa \) is the curvature of the interface, and \(\vec {n}\) is the unitary vector normal to the interface surface \(S_{\rm i}\). The density \(\rho \) and viscosity \(\mu \) in each cell are computed as

$$\begin{aligned} \rho= & {} \gamma \rho _{\rm l}+\left( 1-\gamma \right) \rho _{\rm g},\end{aligned}$$
(23)
$$\begin{aligned} \mu= & {} \gamma \mu _{\rm l}+\left( 1-\gamma \right) \mu _{\rm g}, \end{aligned}$$
(24)

where the subscripts “l” and “g” indicate the liquid and gas phase, respectively.

From the point of view of computational resources and simplicity, it is better to keep the ladle steady and rotate the gravitational field. This requires modification of the solver so that it calculates the tilting angle according to the angular velocity history, as shown in Figure 3, and modifies the direction of the gravitational vector in the momentum Eq. [22].

This also implies that the computational domain is a noninertial frame of reference, hence fictitious forces due to rotation must also be included in the momentum equation. This fictitious acceleration, which has to be included in the right-hand side of the momentum equation, is

$$\begin{aligned} \vec {a}_{\rm f}=-\vec {r}\times \frac{\mathrm {d}\vec {\omega }}{\mathrm {d}t}-2\vec {\omega }\times \vec {u}-\vec {\omega }\times \vec {\omega }\times \vec {r}. \end{aligned}$$
(25)

Although in the typical case in the present study all terms in this equation are much smaller than gravity, they were included in the simulations, except the first one (rotational acceleration), which is typically three orders of magnitude smaller than the other two (Coriolis and centrifugal forces).

With respect to the last term in the momentum Eq. [22], the problem is that, in the VOF method, the interface is not explicitly tracked, so the integral in this term cannot be computed directly. Brackbill[25] overcame this problem by using a continuum surface force model which considers the surface tension force as a volumetric force that acts only on the interface as

$$\begin{aligned} \int _{S_{\rm i}}\sigma \kappa \vec {n}\delta \left( \vec {r}-\vec {r}_{\rm i}\right) \mathrm {d}S\approx \sigma \kappa \nabla \gamma , \end{aligned}$$
(26)

where the curvature \(\kappa \) is computed as

$$\begin{aligned} \kappa =\nabla \cdot \left( \frac{\nabla \gamma }{\left\| \nabla \gamma \right\| }\right) , \end{aligned}$$
(27)

so that it vanishes away from the interface.

The time derivative term is computed using an Euler implicit scheme, which is first order and unconditionally stable. The spatial discretization schemes are second order with Gaussian integration. Since the objective of this study is to simulate the jet trajectory, which is supposed to be not much affected by turbulence structures, no turbulence model was explicitly included in the simulations.

Five ladle geometries, besides the original one, were numerically analyzed. The ladle inner cylinder had height of 0.97 m and base radius of 0.244 m. The original spout geometry is shown in Figure 5. Several lip inclination angles, with respect to the cylinder vertical, were studied for spout curvature of 200 mm. The angles simulated were 0, 10, 20, 30, and 40 deg. The geometries are shown in Figure 6. The ladle rotates around an axis located at y = 0.9 m, z = 1.2 m, with the center of coordinates at the center of the ladle bottom.

Fig. 5
figure 5

Original spout geometry (lengths in mm)

Full size image
Fig. 6
figure 6

Geometries of simulated spouts (lengths in mm) with spout angle of (a) 0 deg, (b) 10 deg, (c) 20 deg, (d) 30 deg, and (e) 40 deg

Full size image

All meshes were created using snappyHexMesh. The meshes have around 95,000 cells, most of them (about 70 pct) being orthogonal hexahedra (Figure 7). The simulations were run on an AMD Opteron 6100 Ubuntu cluster with 64 cores and 64 GB RAM. Each simulation (1000 seconds) required around 25 h with domain decomposition of 16 cores.

Fig. 7
figure 7

Slice of mesh for original geometry

Full size image

Results

As mentioned above, this study focuses on the trajectory of the poured liquid and, especially, on the impingement location of the liquid on the tundish baseline. Figures 8 through 10 show the trajectory of the liquid for the original case and the five geometric modifications for \(t^{*}=0.3\) and \(t^{*}=0.6\). There is a notable difference in the trajectory as a function of spout angle. To quantify this difference for all times against the impingement baseline, further postprocessing calculations must be performed.

Fig. 8
figure 8

Simulated pouring for \(t^{*}=0.3\) and \(t^{*}=0.6\) for original case (a, b) and simulation 1 with spout angle of 0 deg (c, d))

Full size image
Fig. 9
figure 9

Simulated pouring for \(t^{*}=0.3\) and \(t^{*}=0.6\) for simulation 2 with spout angle of 10 deg (a, b) and simulation 3 with spout angle of 20 deg (c, d)

Full size image
Fig. 10
figure 10

Simulated pouring for \(t^{*}=0.3\) and \(t^{*}=0.6\) for simulation 4 with spout angle of 30 deg (a, b) and simulation 5 with spout angle of 40 deg (c, d)

Full size image

The impingement baseline was located 1 m above the tank base, that is, 200 mm below the center of rotation. We define a line with this point and normal to the gravity vector, i.e., dependent on time.

Fig. 11
figure 11

Scheme of rotation of impingement plane as a function of pouring angle: (a) vertical tank, (b) tank tilted by 45 deg, and (c) vertical tank with impingement plane and gravity rotated by 45 deg

Full size image

The sample utility, shipped with the OpenFOAM package, allows computation of the liquid fraction \(\gamma \) along this line. Using a python script, this line was calculated at each time, using the same angle data as the solver, and input to the sample utility. The distribution of \(\gamma \) along this line was thereby obtained for each time. This procedure is shown schematically in Figure 11. From this distribution of \(\gamma \) along the rotated x axis, the position of the maximum value of \(\gamma \) was sought and recorded for each time. The result of this computation for the original (straight) spout is plotted in Figure 12. The size of the points is proportional to the maximum value of \(\gamma \).

Fig. 12
figure 12

Impingement position as function of time for original geometry

Full size image

The impingement position of the poured jet remained more or less stable at \(x\approx \) 0.6 m after the first 300 seconds. However, at the beginning of the process, the dispersion is quite large, with the first impact occurring at around \(x\approx \) 0.3 m. The full extent of the pouring jet on the impingement line is around 170 mm (we consider the extent of the jet to be twice the weighted standard deviation).

Fig. 13
figure 13

Impingement position as function of time for simulation 1 (spout angle 0 deg)

Full size image

When the spout was curved up to 0 deg, the behavior was better, as seen in Figure 13. The dispersion was narrower and closer to the center of rotation, with average position at around \(x\approx \) 350 mm and extension of 190 mm. Figures 14 through 17 show the impingement location of the liquid for simulations 2 to 5.

Fig. 14
figure 14

Impingement position as function of time for simulation 2 (spout angle 10 deg)

Full size image
Fig. 15
figure 15

Impingement position as function of time for simulation 3 (spout angle 20 deg)

Full size image
Fig. 16
figure 16

Impingement position as function of time for simulation 4 (spout angle 30 deg)

Full size image
Fig. 17
figure 17

Impingement position as function of time for simulation 5 (spout angle 40 deg)

Full size image

These results are gathered and plotted in Figure 18 as the mean position of the impingement location and its standard deviation vs spout angle. The average is weighted by the value of \(\gamma \) at each point as well as the standard deviation. The original case (straight spout) is also displayed for comparison.

From these results, it can be concluded that, among the tested cases, the best option is the last one (simulation 5) with spout angle of 40 deg. In this case, the extension of the impingement location is 90 mm.

Fig. 18
figure 18

Impingement location and extension of impact zone for original case and five simulated cases

Full size image

Discussion and Conclusions

Control of the impingement location of a jet poured from a tilting ladle is important in the molten metal manufacturing industry, to avoid spilling liquid outside of the correct location. This work numerically studied its passive control by modifying the angle of the lip spout of the ladle.

Firstly, the angular velocity of the tank required for constant flow rate was calculated, based on geometric considerations and the flow rate over a circular weir. The resulting angular velocity was used as input for a set of CFD simulations to calculate the trajectory of the poured liquid for six cases, viz. the original geometry (straight spout) and five angles (0, 10, 20, 30, and 40 deg) with 200 mm radius of curvature. To simplify the fluid dynamics computations, we proposed to modify the solver to consider variable gravity, rotated at the previously calculated angular velocity. This avoids use of a dynamic mesh, which is computationally expensive. The disadvantage is that the postprocessing is more complicated, since the impingement plane also has to be rotated with gravity.

The results are presented as the \(\gamma \)-weighted average of the impingement location and its standard deviation vs spout angle. These results showed that, among the tested cases, the last one with angle of 40 deg was best, having extension of the impingement location four times smaller than for the original case with a straight spout.