Computer Modeling of Grain Structure Formation during Quenching including Algorithms with Pre-and Post-Solidiﬁcation

: Simulation of the grain growth process, as a function of steel heat transfer conditions, is helpful for predicting grain structures of continuous cast steel products. Many authors have developed models based on numerical methods to simulate grain growth during metal solidiﬁcation. Nevertheless, the anisotropic nature of grain structures makes necessary the employment of new mathematical methods such as chaos theory, fractals, and probabilistic and stochastic theories of simulation. The problem is signiﬁcant for steelmakers to avoid defects in products and to control the steel microstructure during the continuous casting process. This work discusses the inﬂuence of nodal solidiﬁcation times and computer algorithms on the dynamic formation of the chill, columnar, and equiaxed zones including physical phenomena such as nucleation and grain growth. Moreover, the model incorporates pre-nucleation and pre-growth routines in the original algorithm. There is a description of the inﬂuence of the mathematical parameter criteria and probabilities over the grain morphology obtained after solidiﬁcation. Finally, an analysis of these algorithms elucidates the differences between these structures and those obtained from models considering only the solidiﬁcation.


Introduction
Flemings [1] was one of the first who studied the grain structures formed in steel products. This author identified three different grain zones: chill, columnar, and equiaxed. These structures are formed during solidification and are known as primary structures. The grain structures' amorphous morphology is a complex problem that the Euclidean geometry cannot describe [2][3][4]. Fortunately, nowadays, the development of Monte Carlo methods [1,[3][4][5][6], fractals, cellular automaton [7][8][9][10][11][12], and stochastic [9,[13][14][15][16] procedures in addition to deterministic methods make it possible to represent it computationally. In this way, many authors [4,10,12,17,18] have been developing models to describe the solidification process. The dendritic evolution is one of the most important topics [3,10,12,17,18], since this structure is present in steel products obtained by the continuous casting process and other metal industrial casting technologies. This structure contributes significantly to the post-manufactured components' properties, performance, and quality [6,7,13]. Dendrites are part of either columnar or equiaxed grains, depending on the local thermal and solute fields. Other authors have used cellular automata models and non-structured meshes to show dendrites' evolution as a function of time [5][6][7][8]13,18].
The understanding of the physical mechanism for the formation of specific structures and the transition zones from chill to columnar transitions (CCTs) and from columnar to equiaxed (CET) is crucial for controlling grain structure formation during solidification [4,6,[9][10][11]. The solidification models developed [1][2][3][4][5]11,12,[14][15][16][17][18][19][20][21] can be classified as deterministic or stochastic. Deterministic models take averaged quantities and equations solved on a macroscopic scale. The CCT and the CET regions can be predicted by tracking the movement of the columnar front and calculating the growth of equiaxed grains in the undercooled liquid in front of it. A vital issue in deterministic models is the selection of suitable criteria for determining the position where the equiaxed grains block the columnar front and cause the CET and CCT transition.
On the other hand, the stochastic models developed aim to track the nucleation and growth of each grain without assumptions regarding the grain morphology [9,[13][14][15][16]21]. The evolution of the computed grain shape is a function of the local thermal behavior and a random command included in the algorithm. Hence, the CCT and CET regions may be determined based on whether the average final grain shape in a casting appears more columnar or equiaxed.
The original algorithms were modified to evaluate the formation times for every transition zone (t CCT ) and (t CET ) values and the criteria to define the borders between the chill, columnar and equiaxed zones [22][23][24][25][26]. In this work, the evolution of each grain is solved using a nested procedure that analyzes each node (n t I,J ) at every step time (t + ∆t); and random routines simulate nucleation and growth processes. The rules related to the solidification times were included in the same basic algorithm to understand the grain morphology.
During the solidification process of any metal or alloy, compositional partitioning occurs, which results in a no homogeneous solid phase, structure, and composition. The dendritic structures that appear during casting are also not homogeneous, and their fineness influences the mechanical properties. Therefore, the need for microscopic modeling without losing the macroscopic viewpoint is advisable [4,10,12,17,18]. The simulation of dendritic evolution requires a microscopic scale, and some authors have developed models to represent it. The model developed here uses a cellular automata representation in combination with stochastic techniques to create a grain structure [22][23][24][25][26].
Many authors have also developed models to simulate the solidification of some alloys [4,[10][11][12]14,19,20]. The solute rejected from the solid-liquid interface develops a concentration boundary layer in the liquid near the interface and causes a lowering of the adjacent liquid temperature and constitutional super-cooling [1,8,12,13,18,19]. This thermodynamically unstable state creates a rough solid-liquid interface, resulting in a mushy zone. The mushy zone lies between the solid and liquid zones and promotes stable solidification by accepting the rejected solute and heat over temperatures [4,6,14,17,18,21]. Therefore, to discuss the microscopic structure and composition of the solidified material, it is important to clarify the formation of the mushy zone in which the microstructures develop. With this objective, the establishment of a microscopic solidification model is essential. The methodology of this work consists of two steps: First building software resources based on the cellular automata approach and second the use of heat transfer model to feed its results as inputs to the former aiming at the micro-macro structures of the solidified metal.

Grain Formation and Transition Zones
In the simulator, the assumption is that the steel volume has homogeneous physical and chemical properties; and the mushy zone is calculated for each node as a function of the corresponding solidification times. Moreover, the mushy zone is a particular period for each node in the present work. When the analyzed node remains between (t sol I,J ) and (t liq I,J ); then, the remaining time in the mushy zone can be calculated directly using Equation (1) [3,8,12,18,19].
The physical frame consists of steel flowing through a nozzle connecting the tundish and a water copper mold in the continuous casting process. Once in the mold, the heat transfer from the liquid steel to the mold walls leads to a solidified shell, adopting its shape, holding the liquid core. The dwell time of the steel inside the mold is crucial to guarantee that shell thickness is strong enough to contain the liquid steel pressure and the friction forces inside the core [5,11,12,14,17,18]. A thin shell can result in dangerous situations such as breakouts or probable product damage. Due to the rapid cooling rate, there is a chill zone over the billet surface while descending in the mold. A significant number of solidified crystals appear quickly at the beginning of the solidification. These points grow randomly until the impingement limits a small size with other crystals. A thin chill zone decreases the heat transfer flow; consequently, a thin columnar zone grows following the opposite direction of the heat flux vector. This solidification pattern is critical because a large volume of liquid steel remains inside the billet core. However, damages, such as reheating or breakouts, and undesirable risk situations can occur when the latent heat inside the billet flows towards the billet surface.
In the simulators developed [12,14,[17][18][19]24], a subroutine evaluates the distance from every analyzed node to the billet surfaces and finds the shortest time difference (dt I,J ) regarding the nearest neighbors used to give a preferential probability of growing. However, the final decision for the growth is possible through a random routine based on Monte Carlo methods. The result is a cellular automaton model that shows the grain growth printing pixels on the computer screen with colors according to a numeric code. The assignment of the numerical code represents every grain, and the numerical values are assigned to the corresponding nodal position; then, the stored information is in a new temporal 2D computer array.
Grains in the equiaxed zone are randomly oriented. These are bigger than those in the chill zone because the solidification speed is also slow, and they have enough time to grow with their neighbors. In the simulator, time differences between the analyzed nodes (dt I,J ) help to define the borders between the chill, columnar and equiaxed zones by direct comparison with the time criteria defined by the user (t CCT ) and (t CET ).

Nucleation and Grain Growth
Nucleation and growth joint kinetic processes are necessary to simulate solidification of metallic materials. Some models create stationary simulations of grain structures [4,5,11,14,17]. These works used programmed sequences for fixing vertexes in one step as a function of random distribution and joined each other according to their availability to create a grain structure. Unfortunately, these models do not relate to the physical phenomena [15,16,19]. These models were the first used by the beginners of computational simulation in material science and were the pioneers in the field. Nevertheless, it is necessary to include nucleation and growth as shown in case A in Table 1, to show the evolution of grain morphology during the solidification using a step-by-step algorithm and establishing a relationship with the steel thermal behavior. Steel remains a long time in the mushy zone near the billet core. Here, pre-solidification is very important as a function of the liquid and solid fractions assigned with the nomenclature (X liq I,J ) and (X sol I,J ). The inclusion of these processes (cases B and C of Table 1) makes more complex the algorithm, but different phenomena can be appreciated during solidification. The developed computing routines in this work allow treating and representing the solidification evolution computationally in a dynamic sense. This computational procedure allows to simulate more sophisticated structures [10,16].
Nucleation and grain growth are the two fundamental processes involved in solidifying metallic materials. The algorithm includes these processes if the computing and the solidification times of any analyzed node are the same ((t + ∆t) = t sol I;J )). Under this unique condition, steel will solidify. Nucleation and growth occur simultaneously at each step time (t + ∆t) during simulation. Nevertheless, the final decision if a pivoting node of position (I, J) nucleates or grows randomly involves the probabilities (P Nuc ) and (P Growth ) for each grain zone [11,18,20]. A definition of a new nucleated crystal fulfills the condition (N ≤ P Nuc ). Where N is a number randomly generated, respecting the condition (P Nuc + P Growth = 100%). The new grain is identified by assigning a particular color in the new mesh and displayed on the computer screen. On the other hand, the growth process involves including a neighboring nodal analysis and the definition of probabilities for all the un-solidified neighbors surrounding the pivoting node; in the same way, a random routine selects the growing direction.
Nucleation and grain growth are the essential processes to simulate the evolution of the grain structure in the chill zone. Nevertheless, the steel remains a long time in mushy, and the solidification speed decreases for the nodes near the core. Consequently, the morphology in the columnar and equiaxed zones is different [11,14,17,18]. Here, a probability of pre-nucleated or pre-growth with the pivoting node complements the algorithm to obtain a more complex and closer similarity to an actual structure.

Pre-Solidification Algorithms
A pre-nucleated node can appear for the condition ((t + ∆t) < t sol I,J ) if the solid fraction X sol I,J in a node is large enough to be considered as a potential node for nucleation or growing unit. The solid fraction increases at each step time (t + ∆t) as the simulation continues while the liquid fraction (X liq I,J ) decreases. Thus, the probability for a node to be pre-nucleated also increases; but the final decision if the node is pre-nucleated or not is also randomly selected. As the computational time loop advances, many points can be pre-nucleated and grow in columnar and equiaxed zones. Nevertheless, the pre-growth process is considerably different from the growth explained in Table 2. The pre-growth procedure uses inverse criteria in comparison with growth. Therefore, the programming conditions for searching an available neighbor must be logically modified. A pre-nucleated node can grow with an un-solidified neighbor but not neighboring a solidified one because this sentence does not have a logical meaning. Moreover, a solidified neighbor cannot return to be un-solidified because it is already a part of any other grain. Table 2. Description of the processes involved in solidification (the sub-indexes i and j indicate the neighbors' positions).

Including the Only Solidification
Including Pre-Solidification To generate a simulated grain structure using the computational algorithms is necessary to define the values of the parameters shown in Table 3. The probabilities for nucleation, growth, pre-nucleation, and pre-growth must be defined for the three zones (i.e., chill, columnar, and equiaxed) [19,21]. The sum of the probabilities for all the nucleation and growth processes must be equal to (100%). Nevertheless, this condition is not necessary to define pre-solidification routines because pre-nucleation and pregrowth are not developed sufficiently during a simulation [9][10][11][13][14][15][16]21]. In Table 3, values with zero define the chill zone, meaning that there is neither a pre-nucleation nor pre-growth in the chill zone. In the columnar and equiaxed zones, there is 35% and a 5% probability that a node was nucleated or growth, respectively; the remaining 60% assumes the probability of an event of no-occurrence, and it is calculated as the difference of the previous values. Then, a minimum probability for being pre-nucleated or pre-growth (5%) is assigned, if possible; the same probabilities are defined for the central equiaxed zone. A node solidifying when (t sol I,J = t + ∆t) has only two options: nucleating separately (a new grain borns) or growing with a previously solidified neighbor. However, this condition is not necessary for pre-solidification. Here an un-solidified node prenucleates or pre-grows or remain in the mushy zone without change during the simulation before fulfilling the condition (t sol I,J = t + ∆t). Figure 1 shows the structures of a pair of grains obtained using the same input information and the parameters for defining the transition zones, assuming only solidification routines for the simulation [22][23][24]. Here, the three different zones can be easily identified. Nevertheless, the grains in the columnar zone only grow from the chill zone towards the billet core without interruption because only one nucleation sequence performs. In consequence, these grains are very large. Some grow until they meet other grains near the corners or impinge with the grains in the equiaxed zone. Finally, the grains in the central zone are bigger and randomly oriented in the billet core, according to an equiaxed grain morphology. Figure 2a,b show a pair of close-ups near the billet corner. The solidification routine simulates the structure using the input data for cases (1) and (2) in Table 3. The chill zone in Figure 2a is shorter than that in Figure 2b because the criterion used to define the transition zone (dt CCT ) fixes it to be shorter. In these pictures, the grain sizes in the chill zone for case (1) were smaller than those in case (2) because the probability of growth was higher in this case. Near the corner, the columnar grains initially nucleated grow until impinging with the grains growing perpendicularly. Figure 3a,b show the dynamic process, including the probabilities for pre-solidification. The simulations correspond to cases (2) and (3) in Table 3, respectively. In Figure 3a, some grains begin to appear in the columnar zone when only pre-nucleation works. In Figure 3b the pre-nucleated nodes were pre-growth towards the billet surfaces (against the original growth direction of the nucleated grains). In consequence, the growth of the nucleated grains can be interrupted when their leading fronts impinge with the pre-nucleated grains. data for cases (1) and (2) in Table 3. The chill zone in Figure 2a is shorter than that in 2b because the criterion used to define the transition zone ( ) fixes it to be sho these pictures, the grain sizes in the chill zone for case (1) were smaller than those (2) because the probability of growth was higher in this case. Near the corner, the nar grains initially nucleated grow until impinging with the grains growing perpe larly. Figure 3a,b show the dynamic process, including the probabilities for pre-sol tion. The simulations correspond to cases (2) and (3) in Table 3, respectively. In Fig some grains begin to appear in the columnar zone when only pre-nucleation wo Figure 3b the pre-nucleated nodes were pre-growth towards the billet surfaces ( the original growth direction of the nucleated grains). In consequence, the growth nucleated grains can be interrupted when their leading fronts impinge with the p cleated grains.   ups near the billet corner. The solidification routine simulates the structure using the input data for cases (1) and (2) in Table 3. The chill zone in Figure 2a is shorter than that in Figure  2b because the criterion used to define the transition zone ( ) fixes it to be shorter. In these pictures, the grain sizes in the chill zone for case (1) were smaller than those in case (2) because the probability of growth was higher in this case. Near the corner, the columnar grains initially nucleated grow until impinging with the grains growing perpendicularly. Figure 3a,b show the dynamic process, including the probabilities for pre-solidification. The simulations correspond to cases (2) and (3) in Table 3, respectively. In Figure 3a, some grains begin to appear in the columnar zone when only pre-nucleation works. In Figure 3b the pre-nucleated nodes were pre-growth towards the billet surfaces (against the original growth direction of the nucleated grains). In consequence, the growth of the nucleated grains can be interrupted when their leading fronts impinge with the pre-nucleated grains.

Assumptions for the Computational Simulation
The developments of these simulations were under the assumptions explained be low. It is important to point out that the simulations had as the initial conditions the dat file provided by a heat transfer model previously published [22][23][24]. The specific castin conditions of application are in the Appendix A.

•
The nodes represent a squared area forming a squared mesh that discretizes one vol ume of the steel; • There were no grains at the beginning of the simulation ( = 0); this means compu tationally that the cellular automata initialized or returned to their original conditio (all the numerical values were cleaned ( , = 0)). This operation was necessary t avoid errors and numerical conflicts with returned or false values from previous sim ulations; • It was possible to identify when a node belonged to a part of the grain using a logica comparison between the original and corresponding values with the last iteration When the sentence ( , ≠ 0) was true, the pivoting node was already a part of an grain (due to the color assigned to the node).
The following assumptions were applied when only solidification routines were ex ecuted: • An un-solidified node would become a solidified node if the sentence ( , = + ∆ ) was verified as accurate; • When the sentence ( , = + ∆ ) took place, an instantly solidified node could grow with a new un-solidified neighbor; • The leading front in chill, columnar, or equiaxed zones was assumed flat for the nu cleated grains as shown in Figure 3a; and there was only one single solidification front; • The leading solidification front was assumed as unique, flat, and continuous becaus it was a function of the solidification times ( , , , or , ) with the actua simulation time ( + ∆ ); The assumptions used in programming the pre-solidification routines were: • A neighbor of a pre-nucleated node in columnar or equiaxed zones had three option during each step time. Firstly, the node was a pre-growth one if it was an un-solidi fied neighbor. Second, it could be pre-nucleated separately, or third, it could remai

Assumptions for the Computational Simulation
The developments of these simulations were under the assumptions explained below. It is important to point out that the simulations had as the initial conditions the data file provided by a heat transfer model previously published [22][23][24]. The specific casting conditions of application are in the Appendix A.

•
The nodes represent a squared area forming a squared mesh that discretizes one volume of the steel; • There were no grains at the beginning of the simulation (t = 0); this means computationally that the cellular automata initialized or returned to their original condition (all the numerical values were cleaned (C I,J = 0)). This operation was necessary to avoid errors and numerical conflicts with returned or false values from previous simulations; • It was possible to identify when a node belonged to a part of the grain using a logical comparison between the original and corresponding values with the last iteration. When the sentence (C I,J = 0) was true, the pivoting node was already a part of any grain (due to the color assigned to the node).
The following assumptions were applied when only solidification routines were executed: • An un-solidified node would become a solidified node if the sentence t sol I,J = t + ∆t was verified as accurate; • When the sentence t sol I,J = t + ∆t took place, an instantly solidified node could grow with a new un-solidified neighbor; • The leading front in chill, columnar, or equiaxed zones was assumed flat for the nucleated grains as shown in Figure 3a; and there was only one single solidification front; • The leading solidification front was assumed as unique, flat, and continuous because it was a function of the solidification times (t sol I,J , t liq I,J or t mushy I,J ) with the actual simulation time (t + ∆t); The assumptions used in programming the pre-solidification routines were: • A neighbor of a pre-nucleated node in columnar or equiaxed zones had three options during each step time. Firstly, the node was a pre-growth one if it was an un-solidified neighbor. Second, it could be pre-nucleated separately, or third, it could remain unsolidified depending on its actual solid fraction and the probabilities previously defined; • When, the leading front of a pre-solidified grain met the front of the growing grain, the original pre-growth direction was randomly changed if there was at least one available neighbor for growing.
The results in Figure 3a include the pre-nucleation routine. Consequently, a new growing front developed for each new grain; this was the result originating from the performance of more than one nucleation sequence, but the nodes began to grow only when the leading solidification front reached the pre-nucleated node. Multiple grains began to grow against the original solidification front when the simulations included the pre-growth routine as shown in Figure 3b. Figure 4 shows some pre-nucleated nodes in the equiaxed central zone. These nodes began a randomly growth, creating multiple growing fronts. These grains begin to impinge with the grains on the columnar zone forming a non-flat and a non-homogeneous transition CET. However, the grains in the billet core impinged with other neighbors growing isotropically. Nevertheless, the grain leading fronts were multi-directional (without a preferential growth direction). The new pre-nucleated grains appeared as a function of the increment on the solid fraction (X sol I,J ) during the simulation process. Here, the equiaxed zone included the pre-growth phenomena. However, this subroutine worked when the value of the random number generated agreed with the defined probability. When a node was under a pre-solidified process, the latent heat inside was conducted isotropically towards the nearest available un-solidified neighbors, decreasing their solid fractions (X sol I,J ) and increasing the final solidification times neighbors. The solidification front was not flat. There were multiple advancing fronts as new solidified grains appeared (each new grain generated a new advanced front). A neighboring analysis was applied to each pre-nucleated node in the columnar and equiaxed zones. Here, only the nodes in the grain border could grow. The nodes inside the grain were ignored, assuming no participation. un-solidified depending on its actual solid fraction and the probabilities previously defined; • When, the leading front of a pre-solidified grain met the front of the growing grain, the original pre-growth direction was randomly changed if there was at least one available neighbor for growing.
The results in Figure 3a include the pre-nucleation routine. Consequently, a new growing front developed for each new grain; this was the result originating from the performance of more than one nucleation sequence, but the nodes began to grow only when the leading solidification front reached the pre-nucleated node. Multiple grains began to grow against the original solidification front when the simulations included the pregrowth routine as shown in Figure 3b. Figure 4 shows some pre-nucleated nodes in the equiaxed central zone. These nodes began a randomly growth, creating multiple growing fronts. These grains begin to impinge with the grains on the columnar zone forming a non-flat and a non-homogeneous transition CET. However, the grains in the billet core impinged with other neighbors growing isotropically. Nevertheless, the grain leading fronts were multi-directional (without a preferential growth direction). The new pre-nucleated grains appeared as a function of the increment on the solid fraction ( , ) during the simulation process. Here, the equiaxed zone included the pre-growth phenomena. However, this subroutine worked when the value of the random number generated agreed with the defined probability. When a node was under a pre-solidified process, the latent heat inside was conducted isotropically towards the nearest available un-solidified neighbors, decreasing their solid fractions (Xsol I,J,) and increasing the final solidification times neighbors. The solidification front was not flat. There were multiple advancing fronts as new solidified grains appeared (each new grain generated a new advanced front). A neighboring analysis was applied to each pre-nucleated node in the columnar and equiaxed zones. Here, only the nodes in the grain border could grow. The nodes inside the grain were ignored, assuming no participation.

Mathematical and Computational Aspects
The times in the mushy region (t mushy I,J ) were compared and classified to establish the transition zones CCT and CET according to the time criteria shown in Table 3. The thickness of the zones (d chill ), (d col ), and (d equi ) were measured from the billet surface, and their dimensions were directly related to the criteria (t CCT ) and (t CET ). If the value of the criterion (t CCT ) is reduced, the width of the chill (d chill ) zone will also be shorter. The limits of CET and CCT are rounded near the corners and straight over the billet surface, spite that the billet shape is square, as shown in Figure 5a-c, due to the cooling rate near the corners. The computational algorithm used sentences "if" to compare the simulation time with the values stored in the arrays (t liq I,J ), (t sol I,J ), and (t mushy I,J ). Figure 5a-c show the resulting profiles for the values classification using the following criteria: (t CCT = 2 s) and (t CET = 12 s). Here, it is possible to appreciate that the best fit is possible when using (t mushy I,J ). The use of (t liq I,J ) gives a good approach but provides a very thin chill zone and long straight transition zones near the billet corner. These features do not agree with the real grain zones. On the other hand, the use of (t sol I,J ) gives a rounded region between the CET; but the chill zone is nearly nonexistent.
Metals 2022, 12, x FOR PEER REVIEW 9 of 18 Figure 4. Simulation of pre-nucleated nodes in the central equiaxed zone; here, the CET began to form, while multiple pre-solidified grains began to grow in the liquid steel pool inside the billet core.

Mathematical and Computational Aspects
The times in the mushy region ( , ) were compared and classified to establish the transition zones CCT and CET according to the time criteria shown in Table 3. The thickness of the zones ( ), ( ), and ( ) were measured from the billet surface, and their dimensions were directly related to the criteria ( ) and ( ). If the value of the criterion ( ) is reduced, the width of the chill ( ) zone will also be shorter. The limits of CET and CCT are rounded near the corners and straight over the billet surface, spite that the billet shape is square, as shown in Figure 5a and ( = 12 ). Here, it is possible to appreciate that the best fit is possible when using ( , ). The use of ( , ) gives a good approach but provides a very thin chill zone and long straight transition zones near the billet corner. These features do not agree with the real grain zones. On the other hand, the use of (tsol I,J) gives a rounded region between the CET; but the chill zone is nearly nonexistent.   Figure 5a-c, evidencing the size of the chill, columnar, and equiaxed areas after considering these parameters establishing the limits between them. These figures demonstrate that using ( , ) as the classification criterion yielded the best approach to estimate the CCT and CET region sizes.   Figure 5a-c, evidencing the size of the chill, columnar, and equiaxed areas after considering these parameters establishing the limits between them. These figures demonstrate that using (t mushy I,J ) as the classification criterion yielded the best approach to estimate the CCT and CET region sizes.

Procedure and Validation for Multi-Scale Models
The criteria for the estimations of ( ) and ( ) were useful for analyzing an orig inal squared mesh (100 × 100 nodes) obtained initially after the heat removal calculation Nevertheless, the same classification profiles are possible for an enhanced-finer mesh. Th model yielded the same geometrical results for the length of every zone ( , ) by a neighboring analysis nested in the original algorithm and applied t the original or an enhanced-finer mesh. Equations (2)-(6) help obtain this purpose. Here the sub-indexes (I) and (J) belong to the pivoting node; sub-indexes (i) and (j) belong t the nearest surrounding neighbors; (nn) is the number of available neighbors of ever analyzed node; the subindex (ave) is used to indicate the averaged value again of ever analyzed node. These data are particular for each node in the array, and the procedure i nested to update it at every step time; the maximum for a 2D model is eight for nodes i central positions, but it is minor for nodes over the billet surface. Moreover, it is necessar

Procedure and Validation for Multi-Scale Models
The criteria for the estimations of (t CCT ) and (t CET ) were useful for analyzing an original squared mesh (100 × 100 nodes) obtained initially after the heat removal calculation. Nevertheless, the same classification profiles are possible for an enhanced-finer mesh. The model yielded the same geometrical results for the length of every zone (d chill , d col , and d equi ) by a neighboring analysis nested in the original algorithm and applied to the original or an enhanced-finer mesh. Equations (2)-(6) help obtain this purpose. Here, the sub-indexes (I) and (J) belong to the pivoting node; sub-indexes (i) and (j) belong to the nearest surrounding neighbors; (nn) is the number of available neighbors of every analyzed node; the subindex (ave) is used to indicate the averaged value again of every analyzed node. These data are particular for each node in the array, and the procedure is nested to update it at every step time; the maximum for a 2D model is eight for nodes in central positions, but it is minor for nodes over the billet surface. Moreover, it is necessary to include the corresponding sums and comparisons statements in the nested algorithm. In the same way, the solution of Equations (5) and (6) Another option for validation is to find the shortest time difference (dt I,J ) using Equations (7)-(9) separately in nested loops. Here the sub-index (n) is identifying every corresponding neighbor, and a comparison process yields the shortest times, though the same geometrical results for the distances (d chill , d col , and d equi ) were obtained again.
dt liq I,J,n = t liq I,J − t liq i,j dt mushy I,J,n = t mushy I,J − t mushy i,j The shortest times differences (dt sol I,J min , dt liq I,J min , and dt mushy I,J min ) are stored and compared to find the nearest neighbors' status in the 2D model. During neighboring analysis, if two or more leading fronts impinge; one of the following cases occurs:

•
A border between the grains arises if the leading fronts have different numerical codes; • If there are three or four leading fronts, there is a vertex formation; • If the leading fronts are identical in a code number, the grains coarsen, forming a more complex grain. Figure A1a,b show the flow charts for the simulation considering only solidification routines and including pre-solidification routines, respectively. The neighboring process consists of a pair of nested subroutines in the main algorithm; neighboring 1 helps analyze the surrounding neighbors to the pivoting node (I, J) to find the most probable node to grow with or without previous nucleation. Whether neighboring procedure is executed again but to re-distribute the latent heat in time units. Figure A2 shows the procedure how solidification times are saved while thermal behavior is calculated. Figure 7a,b prove that different grain structures are possible because of the inclusion of random selections. In these figures, the probabilities, according to Table 3, influence the changes in the grain size. Only two computational arrays in 2D are required to re-build the entire solidification process again. The use of these arrays let save the information for the changes from liquid to mushy and from mushy to solid (t liq I,J ) and (t sol I,J ) [22][23][24]. Metals 2022, 12, x FOR PEER REVIEW 12 of 18  (2) and (3) in Table 3.
The casting speed and the heat removal strongly influence the steel solidification, as shown in Figure 8a-d. Figure 8a shows that high slopes (with a high tendency towards the billet core) are near the lateral billet surfaces due to the high cooling rate. Nevertheless, the slope describes a slow solidification speed near the billet core. This behavior is typical of the zones where steel remains longer times in the mushy region. Thus, steel thermal behavior can be described as a function of the nodal temperature and simulation time using Equation (10). Figure 8a corresponds to a ¾ cut view of in-built steel thermal behavior, while Figure 8b shows the structure of Figure 7a over-placed on Figure 8a, resulting in the identification of the exact moment when every grain was formed. Figures 8c and  8d, show a ½ cut view and a color contrast view to analyze the sizes and the formation of the CCT and CET. Figure 9 shows the grain growth evolution during the solidification process. Here, it is possible to appreciate the formation of every zone and how the algorithms modify the grain morphology.   (2) and (3) in Table 3.
The casting speed and the heat removal strongly influence the steel solidification, as shown in Figure 8a-d. Figure 8a shows that high slopes (with a high tendency towards the billet core) are near the lateral billet surfaces due to the high cooling rate. Nevertheless, the slope describes a slow solidification speed near the billet core. This behavior is typical of the zones where steel remains longer times in the mushy region. Thus, steel thermal behavior can be described as a function of the nodal temperature and simulation time using Equation (10). Figure 8a corresponds to a 3 4 cut view of in-built steel thermal behavior, while Figure 8b shows the structure of Figure 7a over-placed on Figure 8a, resulting in the identification of the exact moment when every grain was formed. Figure 8c,d, show a 1 ⁄2 cut view and a color contrast view to analyze the sizes and the formation of the CCT and CET. Figure 9 shows the grain growth evolution during the solidification process. Here, it is possible to appreciate the formation of every zone and how the algorithms modify the grain morphology.  (2) and (3) in Table 3.
The casting speed and the heat removal strongly influence the steel solidification, as shown in Figure 8a-d. Figure 8a shows that high slopes (with a high tendency towards the billet core) are near the lateral billet surfaces due to the high cooling rate. Nevertheless, the slope describes a slow solidification speed near the billet core. This behavior is typical of the zones where steel remains longer times in the mushy region. Thus, steel thermal behavior can be described as a function of the nodal temperature and simulation time using Equation (10). Figure 8a corresponds to a ¾ cut view of in-built steel thermal behavior, while Figure 8b shows the structure of Figure 7a over-placed on Figure 8a, resulting in the identification of the exact moment when every grain was formed. Figures 8c and  8d, show a ½ cut view and a color contrast view to analyze the sizes and the formation of the CCT and CET. Figure 9 shows the grain growth evolution during the solidification process. Here, it is possible to appreciate the formation of every zone and how the algorithms modify the grain morphology.

Conclusions and Comments about Results
The influence of the programmed algorithms and solidification speed over the dynamic evolution of the grain structure on the chill, columnar, and equiaxed zones was presented. The conclusions derived from this study are as follows:

•
The most appropriate parameter to be compared is the (t mushy I,J ) to obtain the appropriate dimensions for the grain zones (d chill , d col , and d equi ); • The transition zones, CCT and CET, of the grain structures obtained were easily identified. These zones were continuous when only solidification routines were used, but they did not become regular when pre-solidification routines were included. • Near the corner, the shape of the grains tended to be perpendicular to the billet surface, because the algorithm developed takes into count the perpendicular distance to surfaces. This assumption provides a more accurate shape for the grains in this region according to heat removal; • The model can simulate grain growth on improved meshes to resolve the cellular automaton better and provide a multi-scale validation; • The pre-solidification routines must not be included in the chill zone due to the rapid solidification; • The inclusion of pre-solidification routines does not generate flat transition zones, and more complex grain morphologies can be obtained. However, grains are very fragmented. Thus, its inclusion for simulation must be considered unnecessary because a similar morphology can be obtained, reducing the solidification criteria; • Although pre-solidification was not physically proved, the inclusion of the subroutines in an algorithm with only solidification routines resulted in a better approach to the actual grain morphology for steel billets. Funding: This research study was funded by the Consejo Nacional de Ciencia y Tecnología (CoNa-CyT), SNI, the Instituto Politécnico Nacional (IPN), and the Universidad Autónoma de Coahuila (UAC).
Data Availability Statement: Not applicable.

Acknowledgments:
The authors wish to thank the Autonomous and Technological Institute of Mexico (ITAM), the Consejo Nacional de Ciencia y Tecnología (CoNaCyT), the Instituto Politécnico Nacional (IPN), and the Universidad Autónoma de Coahuila (AUC) for their technical and financial support.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The heat transfer model employs the following operations of the caster. The temperature field and the solidification times across the square section of in the billet work as inputs in the software of the automata cellular shown in Figure A1. Figure A2 shows the previous saving processes to obtain solidification times during heat removal calculation. Tables A1-A4 show the casting conditions, temperatures and chemical composition of the steel, and some parameters of the secondary cooling system, respectively.     Figure A2. Scheme of the previous savings processes to obtain solidification times during heat removal calculation.