3D Model of Carbon Diffusion during Diffusional Phase Transformations

The microstructure plays a crucial role in determining the properties of metallic materials, in terms of both their strength and functionality in various conditions. In the context of the formation of microstructure, phase transformations that occur in materials are highly significant. These are processes during which the structure of a material undergoes changes, most commonly as a result of variations in temperature, pressure, or chemical composition. The study of phase transformations is a broad and rapidly evolving research area that encompasses both experimental investigations and modeling studies. A foundational understanding of carbon diffusion and phase transformations in materials science is essential for comprehending the behavior of materials under different conditions. This understanding forms the basis for the development and optimization of materials with desired properties. The aim of this paper is to create a three-dimensional model for carbon diffusion in the context of modeling diffusional phase transformations occurring in carbon steels. The proposed model relies on the utilization of the LBM (Lattice Boltzmann Method) and CUDA architecture. The resultant carbon diffusion model is intricately linked with a microstructure evolution model grounded in FCA (Frontal Cellular Automata). This manuscript provides a concise overview of the LBM and the FCA method. It outlines the structure of the developed three-dimensional model for carbon diffusion, details its correlation with the microstructure evolution model, and presents the developed algorithm for simulating carbon diffusion. Demonstrative examples of simulation results, illustrating the growth of the emerging phase and affected by various model parameters within particular planes of the 3D calculation domain, are also presented.


Introduction
The microstructure plays a crucial role [1] in forming the mechanical properties of materials.Materials engineers strive to manipulate the microstructure [2] through processes such as heat treatment, mechanical processing, and other methods to achieve the desired material characteristics.In cases of metals, ceramics, polymers, and other materials, the microstructure significantly influences their strength [3], hardness, plasticity, brittleness, and many other mechanical properties.One of the fundamental phenomena that affects the development of the microstructure in materials is phase transformation [4].
Phase transformations in materials depend on a material's chemical composition and refer to changes in the state of matter or internal structure of a material in response [5] to changes in conditions, such as temperature and pressure.
Phase transformations in steels can be divided into two main groups [6]: diffusional transformations [7] and diffusionless transformations [8].Diffusionless transformations are usually much faster [9] than diffusional transformations, and they occur under extreme cooling (and heating) or deformation conditions [10].On the other hand, diffusional transformations require longer times (especially at low temperatures [7]), allowing for the diffusion of atoms and a more gradual change in the crystalline structure.
for the solid-state phase transformations taking place in multicomponent steels that contain substitutional alloying elements.Lahiri presented a review on phase-field modeling of phase transformations in multicomponent alloys [25], where he noticed that, in view of the versatility and quantitative accuracy of the discussed PF models, it can be anticipated that simulations of microstructure evolution using such quantitative PF models will continue to accelerate the design and development of novel multicomponent alloys in the future.Lv et al. [26] found that numerical uncertainty and parameter sensitivity have a major impact on the simulation results in all kinds of phase-field models.
FCA and the LBM are highly effective approaches for simulating various phenomena and processes in the field of materials science.Recently, Svyetlichnyy analyzed the computational time complexity of Frontal Cellular Automata (FCAs) and compared them with Classic Cellular Automata (CCAs) [27], including phase transformation [28].FCA has a very wide range of applications in the context of modeling various processes and phenomena occurring in materials.This method, among others, has been used successfully for the modeling of recrystallization [29], additive manufacturing [29,30], and rolling [31].Recently, parallel computing modeling is also being realized on a larger scale; for example, Zhu et al. [32] realized the three-dimensional multi-phase-field simulation of a eutectoid alloy based on an OpenCL parallel.
This paper presents a 3D model that focuses on the diffusion of carbon in diffusion phase transformations.This model is one of the most important components of a sophisticated system designed to model such transformations.Built upon two modeling methods, FCA and the LBM, the 3D model incorporates the utilization of hybrid computing systems (GPU and CPU) for calculations.Fundamental details regarding the FCA method and the LBM are presented.The structure and primary assumptions of the model are also demonstrated.This paper presents some modeling results of the transformation obtained by this model.Parallel computing is performed on NVIDIA graphics cards with CUDA architecture and with visualization using Open GL.

Lattice Boltzmann Method
The LBM originated from a variant of the cellular automata method known as the lattice gas automaton (LGA) [33].The basis of the method is to solve the Boltzmann transport equation: where f (x,u,t) represents the particle distribution function, x and ξ denote phase space variables (coordinate and velocity), t signifies time, F stands for external force (e.g., gravity), m is the mass, u represents macroscopic velocity, Ω is the collision operator, and e represents lattice velocities.Particles are distributed according to velocities and directions, as follows: where R is the gas constant, D represents the space dimension, ρ signifies gas density, and T denotes temperature.
Along with the discretization of space and time, the LBM uses the discretization of velocities by dividing them into a limited number of directions and setting the velocities according to the lattice space.The velocities ξ are discretized to a set of such discrete velocities e. Equation ( 1) is approximated for characteristic velocities e determined by the grid and by directions using trapezoidal integration for spatial variables and an explicit Euler scheme for time derivatives.The result of such a discrete approximation is the lattice Boltzmann equation (LBE), which presents an evolution of components of the distribution function as follows: Equation ( 3) is the foundation of the LBM.Discretization of two-or three-dimensional space consists of imposing a regular orthogonal grid on a uniform lattice of dimensionless length (the distance between neighboring nodes) equal to one dx = dy = dz = 1.The parameter dt related to the dimensionless time step is also set to unity dt = 1.A fundamental aspect of the system involves selecting a set of velocities e, known as the velocity model, for example, the D2Q9 model, where D indicates dimensionality (D2-two-dimensional model) and Q indicates the number of velocities.In one-dimensional space, the primary models used are D1Q2 or D1Q3, while in two-dimensional space, models such as D2Q4, D2Q5, or D2Q9 are commonly employed.For three-dimensional problems, models like D3Q6, D3Q7, D3Q15, D3Q19, and others find application.
The important element of this equation (Equation ( 3)) is the collision operator Ω i .In general, this operator is a complex nonlinear integral.An approximation of this integral can be fulfilled by linearization near a local equilibrium point.Several approximations are developed: single-time, two-time, and multi-time relaxation schemes and a cascade scheme.The first and simplest single-time relaxation scheme was proposed by Bhatnagar, Gross, and Krook and was called the BGK collision operator: where τ-relaxation time and f eq -equilibrium distribution function.
The macroscopic variables, density ρ and velocity u (through the momentum ρu) (Equation ( 5)), are calculated as follows: The LBE (Equation ( 3)) is solved in two stages: collision and streaming.The postcollision form of the initial distribution function is then as follows: The streaming is a direct transfer of components of distribution functions to appropriate neighboring nodes: Boundary conditions can be considered as a separate operation or in conjunction with a streaming.
Figure 1 shows the general modeling algorithm using the LBM.
priate neighboring nodes: Boundary conditions can be considered as a separate operation or in conjunction with a streaming.
Figure 1 shows the general modeling algorithm using the LBM.A detailed description of the LBM can be found elsewhere [33].

Frontal Cellular Automata
The frontal cellular automaton is a specific type of computational algorithm based on cellular automaton models.The name "frontal" comes from one of its most widespread applications, namely, grain growth during the microstructure evolution, where the boundary of a growing grain is considered as the front of changes occurring in the structure.The main differences in FCA from classical CA are the reversal of information sending (the cell sends information to its neighbors instead of receiving information from the neighborhood), the introduction of an additional state into the automaton, and the organization of the cells into lists.As a result, only small part of all cells is used in calculations in one time step, every cell is maintained only once in the entire simulation, the time complexity is reduced, and simulation time decreases significantly [27].
Three areas can be separated during grain growth in the phase transformation process (Figure 2).Transformation does not begin in the first area, and cells remain in the initial state.The transformation ends in the second area: the cells are in the final state, and this state will not change.The third area is between the first two: it separates them, and the transformation occurs in this area.Only the third area is utilized in FCA calculations.Both the first and second zones should be excluded from the calculations in the current step because changes in the cell states within these zones are not expected and cannot occur.A detailed description of the LBM can be found elsewhere [33].

Frontal Cellular Automata
The frontal cellular automaton is a specific type of computational algorithm based on cellular automaton models.The name "frontal" comes from one of its most widespread applications, namely, grain growth during the microstructure evolution, where the boundary of a growing grain is considered as the front of changes occurring in the structure.The main differences in FCA from classical CA are the reversal of information sending (the cell sends information to its neighbors instead of receiving information from the neighborhood), the introduction of an additional state into the automaton, and the organization of the cells into lists.As a result, only small part of all cells is used in calculations in one time step, every cell is maintained only once in the entire simulation, the time complexity is reduced, and simulation time decreases significantly [27].
Three areas can be separated during grain growth in the phase transformation process (Figure 2).Transformation does not begin in the first area, and cells remain in the initial state.The transformation ends in the second area: the cells are in the final state, and this state will not change.The third area is between the first two: it separates them, and the transformation occurs in this area.Only the third area is utilized in FCA calculations.Both the first and second zones should be excluded from the calculations in the current step because changes in the cell states within these zones are not expected and cannot occur.The use of FCA (Frontal Cellular Automaton) instead of conventional cellular automata reduces simulation time not only in 2D models but also in 3D CA, where it is more noticeable.This is because significant spatial areas are excluded from the calculations at each time step, and only a thin layer participates in the computations.In classical CA, each step requires practically the same computational effort, and the time taken for its deter- The use of FCA (Frontal Cellular Automaton) instead of conventional cellular automata reduces simulation time not only in 2D models but also in 3D CA, where it is more noticeable.This is because significant spatial areas are excluded from the calculations at each time step, and only a thin layer participates in the computations.In classical CA, each step requires practically the same computational effort, and the time taken for its determination remains unchanged throughout the entire simulation.In contrast, the duration of one time step calculations in FCA depends on the number of participating cells, and it exhibits significant variability, but it always remains a small fraction compared to the duration in classical CA.Each cell in FCA participates in calculations only once, in one time step, during the entire simulation, whereas in classical CA the cells participate in calculations in all time steps.
Details about FCA can be found elsewhere [27].

3D Carbon Diffusion Model
A multiscale model, which contains the 3D carbon diffusion model based on the LBM and a microstructure evolution model based on FCA, is shown in Figure 3.This multiscale model was adapted for parallel calculations on GPU with the CUDA architecture.Phase separation, carbon diffusion in different phases, and calculation of the phase fraction in cells on the phase's boundaries are the additional elements implemented into the standard LBM modeling scheme.These elements are shown in Figure 4 as a common block named 'calculation of phase (carbon diffusion)'.Phase separation, carbon diffusion in different phases, and calculation of the phase fraction in cells on the phase's boundaries are the additional elements implemented into the standard LBM modeling scheme.These elements are shown in Figure 4 as a common block named 'calculation of phase (carbon diffusion)'.
The diffusion model is based on the 3D Fourier equation [34]: where c-carbon concentration and D-carbon diffusion coefficient.Phase separation, carbon diffusion in different phases, and calculation of the phase fraction in cells on the phase's boundaries are the additional elements implemented into the standard LBM modeling scheme.These elements are shown in Figure 4 as a common block named 'calculation of phase (carbon diffusion)'.The Bhatnagar-Gross-Krook (BGK) collision operation used in the LBM [33] (Equation ( 4)) is applied for the calculations.
The equilibrium distribution function is calculated as follows: where c is the carbon concentration.
where c-carbon concentration and D-carbon diffusion coefficient.
The equilibrium distribution function is calculated as follows: where c is the carbon concentration.
Determination of the equilibrium distribution functions f eq ; 6.
Collision-calculation of f out ; 7.
Streaming-determination of f in ; 8.
Boundary conditions.

Simulation Results
A comprehensive multidimensional analysis was carried out on the three-dimensional model of carbon diffusion.The established fundamental linkage with the microstructure evolution model, utilizing the FCA method, enabled, among other things, the observation of transformation rates.Selected outcomes of the three-dimensional carbon diffusion modeling are showcased for specified parameter values within the model.
The first two of the modeling variants concern the modeling of carbon diffusion using the developed algorithm in a system composed of one ferrite grain located in the middle of the space.Two variants of the FCA neighborhoods were used for the simulations: von Neumann and Moore.The growth of the grain and changes in carbon concentration were modeled.The simulations were performed utilizing the D3Q19 LBM scheme (Figure 5).The modeling relies on the subsequent fundamental parameters: n x × n y × n z = 128 × 128 × 128, ∆t = 1, τ = 1, ∆x = 1, ∆y = 1, ∆z = 1.Boundary conditions using the bounce-back method were implemented at the edges.The initial concentration for all nodes in the calculation grid was set to 0.6.Concentration in ferrite (α)-0.1,austenite (γ)-0.9, and in interface (c iface )-0.5.The computations were executed utilizing the CUDA parallel programming architecture on graphics cards, namely, NVIDIA GeForce GTX 1060 and NVIDIA GeForce GTX 2080Ti (NVIDIA, Santa Clara, CA, USA).
The first variant with the von Neumann neighborhood is presented in Figure 8.It shows, on the selected planes (x= n x /2, y = 1, 2, . .., n y , z = 1, 2, . .., n z ) of the 3D model space, the growth of a ferrite grain (Figure 8a) and the carbon concentration distributions (Figure 8b).Different blue intensities represent different concentration values, from smaller values in light blue to larger values in dark blue.Figure 9a shows the second variant with the Moore neighborhood.It presents the grain growth for the identical selected plane as the variant in the von Neumann neighborhood.Distributions of carbon concentrations are presented in Figure 9b.The grain growth rate is directly related to the carbon diffusion process, and one of the main driving forces of the growth of the new phase is then the difference in carbon concentrations between the phases.Figure 9a shows the second variant with the Moore neighborhood.It presents the grain growth for the identical selected plane as the variant in the von Neumann neighborhood.Distributions of carbon concentrations are presented in Figure 9b.The grain growth rate is directly related to the carbon diffusion process, and one of the main driving forces of the growth of the new phase is then the difference in carbon concentrations between the phases.Figure 9a shows the second variant with the Moore neighborhood.It presents the grain growth for the identical selected plane as the variant in the von Neumann neighborhood.Distributions of carbon concentrations are presented in Figure 9b.The grain growth rate is directly related to the carbon diffusion process, and one of the main driving forces of the growth of the new phase is then the difference in carbon concentrations between the phases.In the next presented modeling results, the microstructure evolution with five randomly located grains was analyzed.The research was carried out for another selected plane (x = 1, 2, . .., n x , y = n y /2, z = 1, 2, . .., n z ) in the 3D modeled space.The basic parameters used for the modeling were as follows: a number of nodes-n x × n y × n z = 128 × 128 × 128, ∆t = 1, ∆x = 1, ∆y = 1, ∆z = 1, the bounce-back method for boundary conditions, and the Moore neighborhood.The initial concentration for all nodes in the calculation grid was set to 0.6.Concentration in ferrite-0.1,austenite-0.9,and in interface-0.5. Figure 10 shows the grains' growth and the distribution of the carbon concentration in different stages.
In the next presented modeling results, the microstructure evolution with five randomly located grains was analyzed.The research was carried out for another selected plane (x = 1, 2, …, nx, y = ny/2, z = 1, 2, …, nz) in the 3D modeled space.The basic parameters used for the modeling were as follows: a number of nodes-nx × ny × nz = 128 × 128 × 128, t = 1, x = 1, y = 1, z = 1, the bounce-back method for boundary conditions, and the Moore neighborhood.The initial concentration for all nodes in the calculation grid was set to 0.6.Concentration in ferrite-0.1,austenite-0.9,and in interface-0.5. Figure 10 shows the grains' growth and the distribution of the carbon concentration in different stages.The influence of the initial carbon concentration in interface nodes c iface on grain growth and carbon concentration distribution is shown in Figure 11.The size and shape of the grains and the carbon concentration distributions are shown for the same time step tstep = 400.Identical values as in the previous variant were used for the foundational parameters, and the same plane was taken into account.The growth of six randomly distributed grains of the new phase was considered.
The initial carbon concentration in interface nodes directly affects the transformation rate.The difference in concentration between the interface nodes and the neighboring austenite nodes affects the calculation of the boundary velocity.The determination of the volume of the transformed material in the interface node is directly related to the boundary velocity and affects the changes in the new phase fraction.Ferrite grains grow faster under the condition of a lower c iface .
The influence of the initial carbon content c in on grain growth and carbon concentration distribution is shown in Figure 12.The size and shape of the grains and the carbon concentration distributions are shown for the same time step tstep = 700.The base parameters retained the same values as in the previous variant, and the analysis was conducted on the same plane.The calculations were realized with the use of the von Neumann neighborhood.Initial carbon concentration directly affects the transformation rate.Ferrite grains grow faster under the condition of lower initial carbon content in austenite cin.

Summary
Phase transformations, especially diffusion phase transformations in carbon steels, play a crucial role in forming material properties.Numerical modeling remains pivotal in this domain, and newly developed models employ diverse approaches and methods to simulate microstructure formation after transformation.This paper presents the developed 3D carbon diffusion model for the simulation of carbon diffusion during diffusional phase transformations, which is based on the LBM.A detailed presentation of the developed carbon diffusion calculation algorithm is provided.Examples of the simulation results of carbon diffusion and the growth of the new phase in the selected planes of the 3D modeling space for different simulation conditions are also presented.The simulations Initial carbon concentration directly affects the transformation rate.Ferrite grains grow faster under the condition of lower initial carbon content in austenite c in .

Summary
Phase transformations, especially diffusion phase transformations in carbon steels, play a crucial role in forming material properties.Numerical modeling remains pivotal in this domain, and newly developed models employ diverse approaches and methods to simulate microstructure formation after transformation.This paper presents the developed 3D carbon diffusion model for the simulation of carbon diffusion during diffusional phase transformations, which is based on the LBM.A detailed presentation of the developed carbon diffusion calculation algorithm is provided.Examples of the simulation results of carbon diffusion and the growth of the new phase in the selected planes of the 3D modeling space for different simulation conditions are also presented.The simulations used the D3Q19 LBM scheme and the von Neumann and Moore neighborhoods.The developed carbon diffusion model incorporates CUDA parallel computing on GPUs.A fundamental connection was established between the evolved carbon diffusion model and the microstructure evolution model.The developed carbon diffusion model at the next stage will also be connected to the developed 3D heat flow model during diffusional phase transformations.This will be important in the achievement of a comprehensive 3D model of carbon diffusion and heat flow during diffusional phase transformation.

16 Figure 2 .
Figure 2. FCA zones: 1-cells in the initial state, 2-cells in the final state, 3-cells participating in computations.

Figure 2 .
Figure 2. FCA zones: 1-cells in the initial state, 2-cells in the final state, 3-cells participating in computations.

16 Figure 3 .
Figure 3.A 3D model of carbon diffusion (based on LBM) and its integration with the microstructure evolution model (based on FCA) (red arrows symbolize the interaction between these models).

Figure 3 .
Figure 3.A 3D model of carbon diffusion (based on LBM) and its integration with the microstructure evolution model (based on FCA) (red arrows symbolize the interaction between these models).

Figure 3 .
Figure 3.A 3D model of carbon diffusion (based on LBM) and its integration with the microstructure evolution model (based on FCA) (red arrows symbolize the interaction between these models).

Figure 4 .
Figure 4. Carbon diffusion in the calculation of phase in the LBM algorithm.Particular colors represent various stages of the cyclically realized algorithm.The diffusion model is based on the 3D Fourier equation[34]:

Figure 4 .
Figure 4. Carbon diffusion in the calculation of phase in the LBM algorithm.Particular colors represent various stages of the cyclically realized algorithm.

Figure 5 .
Figure 5. D3Q19 LBM lattice.The numbers represent particular directions, and the colors mark the perpendicular planes.

Figure 5 .
Figure 5. D3Q19 LBM lattice.The numbers represent particular directions, and the colors mark the perpendicular planes.

Figure 6 .
Figure 6.Moore (a) and von Neumann (b) neighborhoods.The cell under consideration is shown in burgundy and the cells in a given neighborhood are shown in orange.

Figure 6 .
Figure 6.Moore (a) and von Neumann (b) neighborhoods.The cell under consideration is shown in burgundy and the cells in a given neighborhood are shown in orange.

Figure 7
Figure7shows an algorithm that is used to model carbon diffusion during the transformation.

Figure 6 .
Figure 6.Moore (a) and von Neumann (b) neighborhoods.The cell under consideration is shown in burgundy and the cells in a given neighborhood are shown in orange.

Figure 8 .
Figure 8. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the von Neumann neighborhood.

Figure 9 .
Figure 9. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the Moore neighborhood.

Figure 8 .
Figure 8. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the von Neumann neighborhood.

Materials 2024 , 16 Figure 8 .
Figure 8. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the von Neumann neighborhood.

Figure 9 .
Figure 9. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the Moore neighborhood.

Figure 9 .
Figure 9. FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase utilizing the Moore neighborhood.

Figure 10 .
Figure 10.FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase grains utilizing the Moore neighborhood.The influence of the initial carbon concentration in interface nodes ciface on grain growth and carbon concentration distribution is shown in Figure11.The size and shape of the grains and the carbon concentration distributions are shown for the same time step tstep = 400.Identical values as in the previous variant were used for the foundational parameters, and the same plane was taken into account.The growth of six randomly distributed grains of the new phase was considered.

Figure 10 .
Figure 10.FCA states (a) and carbon concentration distributions (b) at various growth stages of the new phase grains utilizing the Moore neighborhood.

Figure 11 .
Figure 11.FCA states (a) and carbon concentration distributions (b) for different values of ciface.The initial carbon concentration in interface nodes directly affects the transformation rate.The difference in concentration between the interface nodes and the neighboring austenite nodes affects the calculation of the boundary velocity.The determination of the volume of the transformed material in the interface node is directly related to the boundary velocity and affects the changes in the new phase fraction.Ferrite grains grow faster under the condition of a lower ciface.The influence of the initial carbon content cin on grain growth and carbon concentration distribution is shown in Figure12.The size and shape of the grains and the carbon concentration distributions are shown for the same time step tstep = 700.The base parameters retained the same values as in the previous variant, and the analysis was conducted on the same plane.The calculations were realized with the use of the von Neumann neighborhood.

Figure 11 .
Figure 11.FCA states (a) and carbon concentration distributions (b) for different values of c iface .Materials 2024, 17, x FOR PEER REVIEW 14 of 16

Figure 12 .
Figure 12.FCA states (a) and carbon concentration distributions (b) for different values of cin.

Figure 12 .
Figure 12.FCA states (a) and carbon concentration distributions (b) for different values of c in .