Evolutionary Mechanism of Solidification Behavior in the Melt Pool During Disk Laser Cladding with 316L Alloy

Abstract

Laser cladding is an emerging environmentally friendly surface-strengthening technology. During the cladding process, the changes in molten pool temperature and velocity directly affect the solidification process and element distribution. The quantitative revelation of the directional solidification mechanism in the molten pool during the cladding process is crucial for enhancing the quality of the cladding layer. In this study, a multi-field coupling numerical model was developed to simulate the coating process of 316L powder on 45 steel matrices using a disk laser. The instantaneous evolution law of the temperature and flow fields was derived, providing input conditions for simulating microstructure evolution in the molten pool’s paste zone. The behavior characteristics of the molten pool were predicted through numerical simulation, and the microstructure evolution was simulated using the phase field method. The phase field model reveals that dendrite formation in the molten pool follows a sequence of plane crystal growth, cell crystal growth, and columnar crystal growth. The dendrites can undergo splitting to form algal structures under conditions of higher cooling rates and lower temperature gradients. The scanning speed of laser cladding (6 mm/s) has minimal impact on dendrite growth; instead, convection within the molten pool primarily influences dendrite growth and tilt and solute distribution.

1. Introduction

The laser cladding technique is a comprehensive surface-strengthening technology that integrates various disciplines such as materials science, physics, chemistry, and optics. By introducing external strengthening materials into the matrix and achieving a sound metallurgical bond through laser irradiation, its performance can be significantly enhanced [1]. Laser cladding technology presents significant advantages in terms of precision, efficiency, and minimal heat-affected zone, rendering it highly promising for applications in parts remanufacturing, surface repair, and surface strengthening. Coaxial powder feeding is a widely employed technique in laser cladding processes, as shown in Figure 1. The protective gas is conveyed to the molten pool during the cladding process through the flow of metal powder, wherein the metal powder rapidly melts and combines with the matrix to generate a cladding layer exhibiting excellent performance [2,3].

During the laser cladding process, the formation and development of the molten pool play a crucial role in determining the final quality of the cladding layer. A comprehensive understanding of the evolution mechanism of molten pools is essential for predicting their formation. The inner part of the molten pool consists of a combination of metal powder and matrix material. Within a coaxial powder feeding system, not all metal powders are fully melted into the melting pool; only those that fall within it actively contribute to forming the cladding layer. Moreover, it is important to consider both energy absorption efficiency within the molten pool and energy absorption by the powder flow itself. Simultaneously, variations in powder quantity influence momentum and mass characteristics within the molten pool, thereby further impacting its convection and solidification behavior. Numerical simulation is a highly effective approach for constructing the mathematical model of molten pool evolution [4]. Currently, extensive research has been conducted on numerical simulation of molten pools. Kovalev [5] investigated the impact of powder outlet geometry and substrate distance on powder flux, temperature, and molten pool size. Lei [6] observed that the shape of the melt pool and collective heat-affected zone during NiCrBSi powder cladding on Ti6Al4V substrates was semi-elliptical. C. Li [7] developed a thermal–elastic–plastic flow multi-field coupling cladding model that considered laser light shielding by powder and energy absorption.

The performance of the cladding layer is determined by its microstructure. Therefore, it is crucial to analyze the microstructure of the cladding layer during laser cladding optimization. A comprehensive understanding of coating material performance and stability can be achieved by investigating microscopic characteristics such as grain size, grain boundary distribution, and crystal structure in the cladding layer. This profound analysis and optimization provide theoretical guidance for obtaining exceptional properties in cladding layers that meet specific requirements across various application fields, thereby promoting widespread adoption and advancement of laser cladding technology. The dendritic crystal structure is the most prevalent form observed during the process of supercooling solidification, and these intricate crystal structures exert significant influences on the mechanical and material properties of the cladding layer. The microstructure in the solidification process can be studied through three primary approaches: deterministic, stochastic, and phase field methods. The deterministic method is a simulation approach for the continuous nucleation model proposed by Oldfield [8]. The instantaneous nucleation model proposed by Hunt [9] considers a specific moment, a given nucleation density, and a growth rate. Randomness methods primarily encompass the Monte Carlo (MC) method [10,11] and the Cellular Automaton (CA) method [12]. The interface model employed in this method is a sharp interface model, which comprises the thermal diffusion equations for both solid and liquid phases. It is crucial to clearly define the boundary condition coupling at the solid–liquid interface while considering factors such as interface mobility, curvature, and heat flux. This aspect poses significant challenges. The phase-field method circumvents these issues by incorporating an ordered parameter ψ (representing the physical state property of each location in the simulated system) into the model. In general, ψ = 0 signifies the liquid phase state, while ψ = 1 indicates the solid phase state. The entire microstructure can be continuously represented by the sequence parameter ψ, which is bounded between 0 and 1.

In the phase field method, the solid–liquid state is represented by a constant value, while the interface corresponds to the region of order parameter variation. Consequently, the energy term describing interface motion is expressed as a function of Δψ. By minimizing the system’s free energy and considering the thermodynamic gradient of global thermal effects, both ψ (order parameter) and its gradient function are employed to describe the overall free energy during solidification. The state of the phase in any region is described by the order parameter ψ, which serves as the phase indicator. The parameters and state variables in the phase field method are dependent on this indicator and change solidification. Consequently, the evolution of the phase field can be effectively captured through an energy functional. Notably, there is no longer a need to trace the position of the interface; instead, it can be calculated using the phase field parameters.

The solidification process drives the system towards a low-energy equilibrium state by facilitating the evolution of microstructure, thereby minimizing the free energy of the system. The phase field equation framework couples the thermal and concentration field equations to calculate the microscopic characteristics of the structure, employing a physically driven energy minimization method [13]. This approach considers the contributions of chemical potential, thermal field, and concentration field to the overall free energy system, enabling calculations of dendrite shape, alloy segregation, and solute capture in the paste region. However, this method relies on phase gradient Δψ, enthalpy gradient ΔT, and solute concentration Δc [14,15], necessitating accurate deduction of a multivariable well potential function. By incorporating only Δψ as an interfacial motion determinant while adding heat transfer and diffusion equations for solute fields into consideration, significant simplification can be achieved.

Langer [16] and Calginalp [17] formulated the initial phase field equation based on solidification dynamics, subsequently enhancing the phase field model to incorporate anisotropy effects. Kobayashi [18] conducted simulations of two-dimensional dendritic growth with varying degrees of anisotropy using the Ginzburg–Landau free energy potential. Wheeler [19], employing principles of thermodynamic consistency and entropy functional theory, simulated the isothermal solidification phase field model for binary alloys. Karma [20] developed a solidification phase field model for dilute binary alloys, introducing the concept of absolute trapping phases while eliminating non-equilibrium effects arising from interface thickness. Beckermann [21], by investigating the phase field model from columnar to equiaxial, successfully developed a multi-scale model that integrates finite element and random analysis techniques for simulating the microstructure evolution of niobium-based superalloys during the solidification process in laser additive manufacturing. They conducted a comprehensive investigation into the mechanism of dendrite nucleation and growth, the phenomenon of niobium (Nb) segregation, and the formation process of leaf phase particles. Li [22] proposed a multi-field coupling model of the molten pool, revealing that 95% of the cladding layer in the molten pool exhibited columnar dendrite morphology. Qin [23] conducted simulations on the microstructure evolution of thin iron binary alloy within the molten pool, demonstrating that columnar dendrites grew from the bottom of the molten pool following the temperature gradient.

The laser cladding manufacturing process significantly differs from the directional casting solidification process (DS). In DS, heat transfer is unidirectional, with temperature gradients occurring solely in the vertical direction. However, due to the intricate geometry of the paste region and the presence of mobile heat sources, laser cladding exhibits local temperature gradients and convective boundary conditions that actively promote dendrite growth. The convection of melt during the solidification of pure metal into shear flow was simplified by Tonhardt [24]. A corresponding phase field model was established to simulate the dendrite growth process under shear flow conditions, and the influence of convection on the growth morphology of dendrites was investigated. Based on Karma’s phase field model, Beckermann and Tong [25] proposed a comprehensive phase field model that incorporates the flow field, enabling simulation of the impact of parameters such as flow velocity and direction on both the advancing velocity and morphological selection of dendrite tips under forced convection. Therefore, the morphology of dendrite growth is contingent upon the temperature gradient and convection within the solidification pool. To accurately predict microstructure evolution, it is imperative to consider the influence of mobile heat sources and convection in the molten pool. The governing equation F(ψ, T, c) for binary alloy also incorporates a coupled convection term. Most of the above studies use a single perspective to consider the influence of various factors on microstructures. The combination of macroscopical and mesoscopic is rarely used for systematic analysis. During the cladding process, the dendrite growth morphology is not only dependent on the temperature gradient and convection in the solidification pool but also influenced by the heat source moving speed. To accurately predict the evolution of microscopic structures. It is necessary to analyze the convection effect in the molten pool and consider the influence of moving heat sources on the solidification process of the molten pool. In this paper, the heat source moving velocity term and convection term were coupled, respectively, in the governing equation of the binary alloy. The dynamic evolution of microstructure in the laser cladding process is revealed more comprehensively.

In this study, a multi-field coupling numerical model for the process of coating 316L powder with a disc laser on a 45-steel matrix was established. The analysis and calculation of predicting the microstructure development in the cladding layer are compared with the experimental results. The microstructure evolution under moving heat sources and different convection directions was investigated. By employing a solidification model, predictions were made for the spacing and orientation of primary dendrite arms as well as secondary dendrite arms. The dendrite evolution and solute enrichment, induced by the flow field, can be observed through the coupling of the solute field equation. It is of paramount importance to conduct an in-depth investigation into the microstructure evolution of the cladding layer.

2. Experimental Configurations

The growth trend of microstructure in phase field simulation is confirmed through laser cladding experiments conducted in this study. A disk laser model, TruDisk 4000, operating at a power of 1600 W and with an optical fiber minimum diameter of 200 μm, emitting a wavelength of 1030 nm, was employed for the experimental setup. The cladding material utilized was 316L powder. The present study focuses on the analysis of a monorail cladding specimen. Figure 2a shows a monorail cladding specimen with a cladding layer length of approximately 20 mm. To facilitate microstructure observation, the component is divided into two sections: the longitudinal section (shown in Figure 2b) for scanning direction microstructure examination and the cross-section (shown in Figure 2b) for top, middle, and bottom molten pool microstructure observation.

3. Modeling of Laser Cladding Process

The molten pool model was established in this section to obtain a three-dimensional estimation of the molten pool, paste zone size, and temperature distribution within the molten pool. As shown in Figure 1, laser cladding employs a coaxial powder feeding technique, wherein the powder is evenly distributed around the laser beam to form multiple circumferential beams. Additionally, a protective gas (argon) is supplied to the cladding head. The formation of a monorail laser cladding pool was described using fundamental physical assumptions and model formulas.

The underlying assumptions are as follows:

• The fluid within the molten pool exhibits incompressible Newtonian behavior;
• Laminar flow characterizes the movement of molten metal within the pool;
• Powder addition to the substrate follows a Gaussian distribution pattern;
• The absorption and reflection of laser energy by powder exhibit a linear relationship with the amount delivered.

3.1. Governing Equation

The model data can serve as the initial data and boundary conditions for solidification analysis, thereby establishing the prerequisites for melt pool simulation to enhance model feasibility. In this study, the molten pool model is simulated and predicted using a combination of the finite element method and finite volume method. Additionally, a coupled analysis of the temperature field and flow field is performed to accurately determine the size and temperature gradient of the molten pool. In addition to an expanded range of phase field equations, the numerical calculation of the phase field method becomes challenging due to the necessity of incorporating a coupled solution for fluid flow.

By considering the conservation equations of mass, momentum, and energy in the cladding process, we derive the governing equation for the temperature field (Fourier’s law) as follows.

$$\rho C_p{\displaystyle \frac{\partial T}{\partial t}}+\rho C_p\mathit{u}\nabla T=\nabla \cdot \left(k\nabla T\right)-{\displaystyle \frac{\partial \Delta H_T}{\partial t}}-\rho \mathit{u}\Delta H_T$$

(1)

where ρ represents the density of the molten material, C_p denotes the heat capacity at constant pressure, T signifies the temperature during the cladding process, u indicates the rate of molten flow, and k stands for thermal conductivity.

The latent heat of the melt is taken into consideration, where ΔH_T represents the latent enthalpy and is expressed as follows:

$$\Delta H_T={\beta }_L{f}_L$$

(2)

where β_L is latent heat. f_L is the mass fraction of the fluid, and the expression is as follows:

$$f_L=\left\{\begin{array}{ll}0& T<T_S\\ {\displaystyle \frac{T-T_S}{T_L-T_S}}& T_S\le T\le T_{L}s\\ 1& T_L<T\end{array}\right.$$

(3)

where T_S and T_L are the solid–liquid temperatures of the melt.

The flow of molten material within the liquid pool is governed by fluid dynamics. Under the assumed conditions, the continuity equation and Navier–Stokes momentum equation are employed to describe the motion of the fluid.

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla \cdot \left(\rho \mathit{u}\right)=0$$

(4)

$$\rho \left[{\displaystyle \frac{\partial \mathit{u}}{\partial t}}+(\mathit{u}\cdot \nabla )\mathit{u}\right]=-\nabla \cdot (pI)+\nabla \cdot \left[{\mu }_L\left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)\right]+D_A$$

(5)

where I is the spatial unit matrix, p is the melt pressure, and μ_L is the dynamic viscosity of the melt.

3.2. Boundary Conditions and Source Terms

During the cladding process, it is not always guaranteed that the powder can be fully absorbed by the molten pool, thus necessitating consideration of surface growth in terms of powder absorption rate [26].

$$V_Z={\displaystyle \frac{2V_f{\eta }_f}{{\rho }_f\pi R_f^2}}\exp \left[-{\displaystyle \frac{2\left({\left(x-V_{s}t\right)}^2+y^2\right)}{R_f^2}}\right]z$$

(6)

where V_f is the powder delivery velocity, η_f is the powder absorption rate, ρ_f is the powder flow density, R_f is the powder flow radius, V_S is the scanning velocity, and z is the unit vector in this direction.

During the solidification process, the size of the melt pool is directly influenced by the velocity of the molten material. As solidification progresses and temperature drops below the liquid phase line, a significant decrease in melt speed occurs. Therefore, the inclusion of an additional drag term in the momentum equation elucidates the observed decrease in velocity:

D_A = −μ_LK_Au

(7)

where μ_L represents the temperature-dependent dynamic viscosity, and permeability can be determined using the Kozeny–Carman equation.

$$K_A=K_v{\displaystyle \frac{{\left(1-f_L\right)}^2}{f_L^3+B}}$$

(8)

The formula incorporates the momentum damping coefficient, K_v, as a constant and B as a minimum value. To prevent division by zero errors, it is essential to restrict the liquid phase fraction f_L within the range of 0 < f_L < 1.

The static force balance of the liquid metal is disrupted by the uneven temperature distribution in the melting pool, leading to fluid flow driven by the temperature gradient. The gas–liquid interfaces between the molten pool surface and gas exhibit momentum expression as described [27].

$$F_{l/s}=f_s\mathit{n}{\kappa }_s-\nabla T{\displaystyle \frac{d{f}_s}{dT}}$$

(9)

where f_s is the surface tension, n is the interface normal vector, and k_s is the interface curvature.

The dominant forces governing the flow of molten metal in a molten pool include gravity, buoyancy, and surface tension. By the Boussinesq approximation principle, the expression for thermal buoyancy resulting from non-isothermal flow can be formulated as follows [28]:

$$F_B=-\rho \mathit{g}{\beta }_T\left(T-T_M\right)$$

(10)

where T_M is the melting temperature of the melt, g is the acceleration of gravity, and 9.8 (m∙s⁻²) is taken.

The Marangoni effect is a phenomenon in which the surface tension of a liquid undergoes flow in response to variations in temperature or concentration. The surface tension of the molten pool decreases with increasing temperature, resulting in elevated shear stress near the periphery of the molten pool, thereby inducing lateral expansion and reduced depth of the molten pool [29]. Assuming a linear relationship between temperature and the decrease in surface tension, γ_b can be expressed as γ_b = γ_b0 − ζ(T − T_M). However, this formulation neglects the influence of solute elements present in the molten pool. Therefore, considering thermodynamic principles, an expression accounting for solute elements is proposed [30]:

$${\gamma }_b={\gamma }_{b0}-\zeta \left(T-T_M\right)-R_{G}T{\Gamma }_S\ln \left[1+k_L{a}_i\exp \left({\displaystyle \frac{\Delta H_0}{R_{G}T}}\right)\right]$$

(11)

where γ_b0 is the surface tension of the melt at time T, ζ is the surface tension constant of the temperature coefficient, R_G is the common gas constant, and Γ_S is the saturated excess of the metal surface, the value of which is 1.3 × 10⁻⁵ mol·m⁻². a_i is the solute element function under temperature gradient, k_L is the segregation entropy correlation constant, and the value is 3.18 × 10⁻³. ΔH₀ is the standard heat of adsorption, with the value 1.67 × 10⁴ J∙mol⁻¹.

3.3. Model Parameters and Mesh Settings

The numerical model was solved using a multi-physics finite element computing platform. A rectangular model of 20 mm × 10 mm × 5 mm was established. The model used symmetric boundaries along the x–z plane. The whole mesh was divided by free tetrahedral mesh. The maximum cell is 1.6 mm. The intermediate area was refined. The maximum unit is 0.25 mm. The minimum unit is 0.004. The mesh partition contains 131,583 domain units, 5302 boundary units, and 252 side units. The geometric model and meshing are shown in Figure 3. The laser cladding process parameters are shown in

Table 1. Numerical simulation of process parameters.

Process Parameter	Value
Laser power P0/(W)	1000
Spot radius r/(mm)	2.5
Scanning speed vx/(mm/s)	6
Powder feed velocity Vf/(g/min)	6
Powder flow radius Rf/(mm)	2.5
Ambient temperature Tamb/(K)	298.15

.

3.4. Numerical Analysis and Experimental Verification

The laser cladding process involves non-equilibrium metallurgical solidification, characterized by “hot and cold” conditions, with transient cladding. The physical property parameters of both the powder and matrix exhibit significant fluctuations in response to temperature variations. The CALPHAD method was employed to analyze the thermal properties of both powder and matrix materials as a function of temperature and to simulate the cladding process. Utilizing a multi-physics coupling computing platform, we establish the boundary conditions, source terms, and mesh division while incorporating temperature change physical property parameters to solve the heat transfer equation, continuity equation, and momentum equation.

The longitudinal cross-section of the predicted molten pool is shown in Figure 4a. The transient molten pool is generated by a moving heat source, with the scanning direction aligned with the X-axis and the surface growth direction of the cladding layer aligned with the Z-axis. Laser scanning penetrates through the matrix to obtain a high-temperature transient molten pool. The width of the laser-scanned molten pool is shown in Figure 4b.

Figure 5a shows the cladding temperature distribution cloud map at 0.3 s. The figure shows that the maximum temperature of the molten pool reaches 2730 K, surpassing both the melting temperature of the matrix and cladding powder. Consequently, a minute molten pool initiates formation on the surface of the matrix. However, due to the brief duration of laser irradiation, there is limited powder infiltration into the molten pool, resulting in a relatively diminutive molten pool. Under continuous laser energy input, the surface temperature of the matrix exhibits a rapid increase, leading to the accumulation of heat within the matrix. Consequently, an expanding heat-affected zone is observed, resulting in the complete formation of the melt pool at 1.0 s. The transition zone of cladding preheating is observed between the melting and non-melting zones in Figure 5b, with an elliptical molten pool forming near the spot’s center, where the maximum temperature reaches 2932 K. Figure 5c,d shows that the maximum cladding temperature of the cladding process is about 2950 K. After laser scanning, the molten pool is cooled and solidified instantly. Obviously, a raised cladding layer is formed behind the laser scanning path.

The laser cladding process induces thermal expansion in the irradiated area of the matrix surface, while the unirradiated region restricts this expansion, leading to internal stress within the matrix and subsequently impacting the quality of the cladding layer. The stress distribution of the matrix at different times is shown in Figure 6. Figure 6a presents a cloud map illustrating the stress distribution in the matrix at 0.3 s. Currently, the heat-affected area is limited, and the maximum matrix stress reaches 277 MPa. Figure 6b shows the cloud image of matrix stress distribution at 1.0 s, revealing that stress levels are low in the center of the spot while circularly distributed stresses dominate its front region, with this area exhibiting maximum stress values. Figure 6c shows the cloud image of matrix stress distribution at 2.0 s, demonstrating an expanding range of stress distribution as laser movement progresses; a long strip-shaped zone of high stress emerges behind laser scanning, reaching up to 294 MPa in terms of matrix stress levels. Figure 6d shows the stress distribution cloud map at 3.0 s. At this time, the stress range of the front end of the spot becomes smaller. The stress region after laser scanning is distributed in a “long strip” shape. The maximum matrix stress is 303 MPa.

Figure 7 shows the cloud diagram of flow field distribution in the molten pool during laser cladding. At the beginning (0.3 s), the flow field in the molten pool was not completely formed, and a small circulation could be seen, and the maximum velocity was 0.23 m/s, as shown in Figure 7a. The formation of Marangoni convection in the melt pool was evident at 1 s, during which time the flow field reached its complete development with a maximum velocity of 0.34 m/s, as shown in Figure 7b. Comparatively, the distribution of the flow field within the molten pool remained essentially unchanged at both 2.0 s and 3.0 s, with only slight variations in flow velocity observed, reaching a maximum value of 0.35 m/s, as shown in Figure 7b–d. According to the numerical simulation results of the flow field, a maximum velocity of 0.35 m/s was selected for analysis of its influence on dendrite growth, with the simulated flow direction being from left to right.

To ensure the accuracy of the input conditions for microstructure evolution in the molten pool, a comparison was made between numerical results of macroscopic and experimental observations of cladding layer morphology, as shown in Figure 8. The figure illustrates the cladding layer, with W representing its width and H representing its height. On the left, we present the numerical model calculation result, where the cladding profile has a width of 5.851 mm and a height of 0.682 mm. On the right, we display the scanning electron microscope observation result, which shows a cladding profile width of 5.838 mm and a height of 0.675 mm. Comparative analysis reveals that the experimental results align closely with the numerical simulation results, demonstrating high consistency in terms of heat-affected zone dimensions with an error margin of only 0.22% for width and 1.09% for height, thus validating the effectiveness of our numerical model.

4. Melt Solidification Simulation

The melt solidification simulation problem involves describing the growth of crystallization in a molten state, originally based on the classical solid phase growth model derived from diffusion theory. The most classical approach is based on the sharp interface model, wherein accurate tracking of the solid–liquid interface motion is necessary to determine changes in its position and shape. However, due to the intricate influence from all directions on dendrite growth in the molten pool, this problem becomes particularly complex [31]. Therefore, a multi-scale method is employed to predict microstructure during dendrite solidification from various directions. The phase field analysis of the supercooled paste region is conducted independently without considering the high-temperature molten pool. To better elucidate the microstructure of the solidification process in the molten pool, the paste region is decomposed into two orthogonal cross-sections. The first section represents the cross-sectional view of the molten pool (Figure 9a), which serves to depict the gradient microstructure within. To address the issue of varying cross-sections, we assume that the interface aligned with the scanning direction exhibits consistent solidification behavior. The subsequent section represents a longitudinal view along the cladding direction (as shown in Figure 9b), which effectively characterizes dendrite growth along the scanning direction of the cladding layer.

When simulating the microstructure of the longitudinal section, we set the bottom and left (distal) boundary to the solid-phase line temperature to simulate the boundary conditions of the moving paste region. Meanwhile, due to the continuous release of latent heat, the paste region remains hot.

According to the Ginzburg–Landau theory of free energy [32], within a closed system V, the expression for free energy incorporates the influence of solute:

$$F={\displaystyle \underset{V}{\int }\left[f(\psi ,c,T)+\frac{1}{2}{\epsilon }^2{(\nabla \psi )}^2+\frac{1}{2}{\delta }_c{(\nabla c)}^2\right]}dV$$

(12)

where f(ψ, c, T) is the free energy density function, c is the solute concentration, T is the temperature driving term, and ε and δ are the phase field gradient term and solute gradient term, respectively.

In the phase field method, ψ is employed as a scalar quantity representing different states for the solid phase, liquid phase, and interface. Within the phase field model, each point’s state value ranges from 0 to 1, where ψ = 1 denotes complete solidification and ψ = 0 represents complete melting. Building upon the interface relation of Gibbs–Tomson theory, Boettinger et al. [33] established a geometric derivation of the isotropic phase field equation by relating the interface normal velocity to temperature term, liquidus slope, and interface curvature. The resulting equation can be expressed as follows:

$${\displaystyle \frac{1}{{\xi }_k}}{\displaystyle \frac{\partial \psi }{\partial t}}=\mathsf{\Gamma }\left[{\nabla }^2\psi +{\displaystyle \frac{\psi (2\psi -1)(1-\psi )}{d_1}}\right]+(T_{eq}-T+m_L(c_l-c_0)){\displaystyle \frac{(\psi (1-\psi ))}{d_1}}$$

(13)

where ξ_k is the dynamic coefficient, T_eq is the equilibrium temperature, and c_l is the liquid phase equilibrium concentration of the binary alloy, c₀ is the initial average concentration of the binary alloy, d₁ is the thickness of the solid–liquid interface, Γ is the Gibbs–Tomson coefficient, and m_L is the liquidus slope of the alloy.

The Equation (13) is thus further refined to incorporate the linear correlation between the interface temperature and the binary phase composition as well as melting temperature: m = (d₁/2Γ) × (T_L − T_S) × (1 − T + m_L(c_l − c₀)). The convergence of the numerical model is enhanced by replacing the linear function with a bounded function: m = (α/π)arctan{γ[T_M − T + m_L(c_l − c₀)]}, where α and γ are constants, and mL represents the liquidus slope. In this simulation, molten 316L is considered a binary solution, with Fe acting as the solvent and Ni as the solute. The concentration of the aggregated second-phase precipitation is assumed to deviate no more than 14% from the composition of the initial equilibrium concentration solution.

The original equation describing the phase field variables was derived from the minimum energy functional of ψ. For most metallic materials, the surface energy and kinetic coefficients of the solid–liquid interface depend on crystal orientation. Due to the significant impact of anisotropy on crystal growth and dendrite morphology, it is imperative to modify the anisotropy within the phase field model. Building upon the models proposed by Kobayashi [18], Boettinger [14], Wheeler [19], and Zaeem et al. [13]’s approach, this study refines the anisotropic formula for solidified alloy materials. In a two-dimensional model, the diffusion coefficient can be expressed as a function of the growth direction to encompass all members of the opposite sex: θ is defined as arctan((∂ψ/∂y)/(∂ψ/∂x)), where ε = $\tilde{\epsilon }$σ(θ) and $\tilde{\epsilon }$(1 + δcos(j(θ − θ₀))). In this equation, $\tilde{\epsilon }$ represents the anisotropic gradient energy coefficient at the interface, δ denotes the intensity of anisotropy, and j signifies the modulus of anisotropy (symmetry order). The crystal structure of the material presented in this paper is that of face-centered cubic (FCC), with a coordination number of 4. The initial orientation θ₀ for each dendritic dendrite is determined either randomly or by established principles. Consequently, the anisotropic modified phase field equation can be expressed as follows:

$${\tau }_0{\displaystyle \frac{\partial \psi }{\partial t}}=\nabla \cdot \left(\left(\begin{array}{cc}{\epsilon }^2& -\epsilon {\epsilon }^{\prime }\\ \epsilon {\epsilon }^{\prime }& {\epsilon }^2\end{array}\right)\nabla \psi \right)+\psi (1-\psi )(\psi -{\displaystyle \frac{1}{2}}+m(T,c))$$

(14)

where τ₀ is the phase field fitting parameter, σ′ = dσ(θ)/dθ.

The advancement of the interface in the phase field equation is contingent upon both the temperature field equation and the solute diffusion equation of the binary solution. Considering the volume control of the interface region in the diffusion phase, assuming a time increment of Δt and a phase indicator increment of Δψ in the phase field, ρLΔψ represents the amount of latent heat produced or consumed. Therefore, by discretizing it in time, the rate of heat change can be expressed as ρL∂ψ/∂t. Consequently, the enthalpy variation in the binary system within Δt time is given by Δh = ρC_pΔT − ρLΔψ, where the negative sign represents the transition of ψ from liquid phase (ψ = 0) to solid phase (ψ = 1) during solidification. The equation governing heat transfer is as follows:

$${\displaystyle \frac{\partial T}{\partial t}}={\alpha }_T{\nabla }^{2}T+{\displaystyle \frac{L}{C_p}}{\displaystyle \frac{\partial \psi }{\partial t}}$$

(15)

where α_T is the thermal diffusion coefficient (α_T = k_T/ρC_p) and k_T is the thermal conductivity.

The conservation equation was formulated by Karma and Boettinger et al. [11], incorporating the weighted value of the phase field equation for average concentration, resulting in the derivation of the evolution equation for average concentration.

$${\displaystyle \frac{\partial c}{\partial t}}=\nabla \cdot \left[\psi D_S\nabla c_S+(1-\psi )D_L\nabla c_L\right]$$

(16)

where D_S is the solid phase diffusivity, D_L is the liquid phase diffusivity, and c is the average concentration of the binary alloy solution. The solid phase of the solute concentration for c_S = k_fc/[(1 − ψ) + k_fψ], the concentration of the solute in the liquid phase for c_L = c/[(1 − ψ) + k_fψ].

Considering the impact of rapid solidification through laser cladding, to address the issue of non-equilibrium solute field distribution, the determination of the distribution coefficient incorporates the influence of the solidification rate [34]. Therefore, k_f = (k_eq + v/v_eq)/(1 + v/v_eq), where k_eq and v_eq represent equilibrium quantities, and the actual solidification rate is v = (∂ψ/∂t)/(∇ψ). The final simplified diffusion equation for a binary alloy’s solute field can be expressed as follows:

$${\displaystyle \frac{\partial c}{\partial t}}=\nabla \cdot \widetilde{D_k}\left[\nabla c+{\displaystyle \frac{(1-k_f)c}{1+(k_f-1)\psi }}\nabla \psi \right]$$

(17)

where ${\tilde{D}}_k=D$_S+ [(D_L − D_S)(1 − ψ)]/(1 + (k_f − 1)ψ).

The Equations (14), (16), and (17) are derived under the assumption of negligible convective velocity in the moving heat source. Hence, to incorporate the scanning velocity in the x direction into the longitudinal section (x-z) and modify the coupling equation of (ψ, T, c) [35], we present the following results:

$${\tau }_0{\displaystyle \frac{\partial \psi }{\partial t}}=-{\displaystyle \frac{\partial }{\partial x}}\left(\epsilon {\epsilon }^{\prime }{\displaystyle \frac{\partial \psi }{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\epsilon {\epsilon }^{\prime }{\displaystyle \frac{\partial \psi }{\partial x}}\right)+\psi (1-\psi )(\psi -{\displaystyle \frac{1}{2}}+m(T,c))+{\tau }_0{v}_x{\displaystyle \frac{\partial \psi }{\partial x}}$$

(18)

$${\displaystyle \frac{\partial T}{\partial t}}-{\displaystyle \frac{L}{2C_p}}{\displaystyle \frac{\partial \psi }{\partial t}}={\alpha }_T{\nabla }^{2}T+v_x{\displaystyle \frac{\partial T}{\partial x}}-{\displaystyle \frac{L{v}_x}{2C_p}}{\displaystyle \frac{\partial \psi }{\partial x}}$$

(19)

$${\displaystyle \frac{\partial c}{\partial t}}=\nabla \cdot \widetilde{D_k}\left[\nabla c+{\displaystyle \frac{(1-k_f)c}{1+(k_f-1)\psi }}\nabla \psi \right]+v_x{\displaystyle \frac{\partial \psi }{\partial t}}$$

(20)

where v_x follows the scanning direction of the cladding (consistent with V_s).

The dimensionality normalization in the above Equations (18)–(20) is as follows: T = (T_actual − T_S)/(T_L − T_S), l = l_actual/l₀, t = t_actual/τ, v_x = (v_actualτ/l₀), where l₀ is the characteristic length scale, τ is the characteristic time scale $l_0^2$/α_T, τ₀ is the characteristic constant, and T_L and T_S are the liquid and solid phase temperatures, respectively.

During the laser cladding process, convection is also induced by temperature gradient, concentration gradient, and surface tension in the molten pool. The flow of molten material can impact solute distribution within the pool and influence microstructure growth. The phase field equation of the microstructure [36] with the addition of the convection term is as follows:

$${\tau }_0{\displaystyle \frac{\partial \psi }{\partial t}}=-{\displaystyle \frac{\partial }{\partial x}}\left(\epsilon {\epsilon }^{\prime }{\displaystyle \frac{\partial \psi }{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\epsilon {\epsilon }^{\prime }{\displaystyle \frac{\partial \psi }{\partial x}}\right)+\psi (1-\psi )(\psi -{\displaystyle \frac{1}{2}}+m(T,c))$$

(21)

$${\displaystyle \frac{\partial T}{\partial t}}+\left(1-\psi \right)\mathit{v}\nabla T={\alpha }_T{\nabla }^{2}T+{\displaystyle \frac{L}{2C_p}}{\displaystyle \frac{\partial \psi }{\partial t}}$$

(22)

$${\displaystyle \frac{\partial c}{\partial t}}+\left({\displaystyle \frac{1-\psi }{1-\psi +k_f\psi }}\right)\mathit{v}\nabla c=\nabla \cdot \widetilde{D_k}\left[\nabla c+{\displaystyle \frac{(1-k_f)c}{1+(k_f-1)\psi }}\nabla \psi \right]-\nabla j_{at}$$

(23)

where v is the velocity of forced convection in the simulation domain and the velocity of the dimensionless liquid phase (ψ = 0). μ is the kinematic viscosity of the fluid. The term “jat” proposed by Karma et al. [20] represents the anti-solute retention phenomenon. The introduction of the absolute trap term aims to prevent the discrepancy between the diffusion velocity of the solute and the simulated interface moving velocity caused by the dendrite growth rate in simulations.

The detailed boundary conditions for the variables ψ, T, and c can be found in

Table 2. Cross-sectional boundary conditions.

Boundary Conditions	ψ	T	c
1	Zero flux	1602.26 K	Zero flux
2	Zero flux	1737.26 K	Zero flux

and

Table 3. Boundary conditions of the longitudinal section.

Boundary Conditions	ψ	T	c
3	Zero flux	1602.26 K	Zero flux
4	Zero flux	1737.26 K	Zero flux
5	Zero flux	T linearly varying from 1602.26 K to 1737.26 K	Zero flux
6	Zero flux	1737.26 K	c0
7	Zero flux	1737.26 K	cL

. Figure 9 shows that the dendrite growth nucleated seeds are distributed along borders 4, 5, and 6 on the longitudinal cross-section. The cross-sectional nucleated seed is defined at boundary 1. Additionally, Figure 10 presents the program flow chart of the phase field simulation.

5. Analysis of Numerical Simulation Results

The finite element method was employed for numerical analysis. The feature length was normalized to l₀ = 25 μm, and the mesh was generated using a free triangle mesh approach with a maximum cell size of l_mesh/l₀ = 0.1. In the longitudinal section, there were 52,918 units, while in the transverse interface, there were 43,906 units.

The evolution of the crystal during solidification is shown in Figure 11. Initially, the crystal exhibits a planar morphology that subsequently transforms into a cellular structure due to the influence of temperature and solute gradients. This cellular crystal then undergoes rapid growth and develops directional branches, ultimately forming a dendritic architecture. It is noteworthy that the longitudinal cross-section exhibits variations in the initial temperature of crystal seeds at different boundaries, leading to divergent microstructure growth processes along the three directions. Specifically, the lower section of Figure 11a remains planar, while the left and upper sections tend to cellular crystallization. In Figure 11b, the lower plane has started to transition into cellular crystals, whereas the left and upper regions have displayed indications of lateral branching. This phenomenon arises due to a significant temperature gradient in the left and upper areas, accelerating crystal growth in these regions.

In Figure 11c, it is evident that the crystal tip bifurcates at the left boundary, giving rise to a thick, filamentous structure resembling algae. The formation of seaweed-like structures typically occurs under conditions characterized by a low-temperature gradient and a high cooling rate. This is due to the instability of the solid–liquid propulsion interface, which leads to the continuous division of dendrite tips and the subsequent development of seaweed-like structures [37]. The formation of this algae-like structure is closely linked to the rapid cooling characteristics during laser cladding. It should be noted that the seaweed-like structure represents only one stage in the growth of dendritic structures. As crystal growth continues, the tip of the split dendrite arm gradually diminishes while secondary dendrite arms begin to grow, ultimately forming a complete dendritic structure, as shown in Figure 11d.

In Figure 12a, we investigate the motion of the laser heat source and observe dendrite growth at the 45,000th time step. Compared with the previous Figure 11d, it is evident that there are no significant alterations in dendrite orientation and tip velocity. This can be attributed to the relatively low scanning speed employed during laser cladding, which exerts minimal influence on the solidification process of the molten pool. However, subtle disparities exist in the perpendicular and scanning directions of dendrite arms. In Figure 12a, it is evident that the left secondary dendrites exhibit marginally higher growth rates, while certain dendrites display a slightly greater abundance of left secondary branches. Moreover, as shown in Figure 12b of the experimental results, it is evident that the columnar crystals exhibit a predominantly perpendicular orientation to the molten pool’s bottom and remain minimally influenced by natural convection. The figure also shows the morphology of adjacent crystals, which transition from flat bases to cellular, columnar, and equiaxed structures with increasing height.

The dendrite morphology at the top of the scanning direction is shown in Figure 12c, revealing specific observations. Firstly, the secondary dendrite arm of the columnar crystal fractures at mark 1, resulting in an equiaxed crystal structure formation. Secondly, markers 2 and 3 indicate interference between adjacent dendrites during growth, which restricts the generation of secondary dendrite arms and consequently leads to fewer and smaller dendrites observed within the dendrite gap region. The numerical simulation results also validate this phenomenon, elucidating the influence of competitive growth restriction that enhances the dendrite’s diameter structure and ultimately manifests as equiaxed crystal morphology. The remarkable concurrence between experimental observations and numerical simulations strongly unveils the intricate microstructural evolution during solidification.

The distribution of solutes in the dendrite structure is shown in Figure 13. Solute segregation is prominently observed at the root of the dendrites while gradually diminishing as dendrite growth progresses toward the tip. By comparing the solute field scanned by the laser heat source (Figure 13b), it is evident that the scanning speed exerts an influence, resulting in a slight increase in the content of solute elements between the dendrites and enrichment of solute along the scanning direction. Based on the phase field profiles shown in Figure 11 and Figure 12, it is evident that within the supercooled region, there exists a pronounced inclination of columnar crystals to grow perpendicular to the interface. The presence of similar dendrite structures can also be observed in cross-sections with higher undercooling, as shown in Figure 14. The orientation of dendrites undergoes changes during their growth and approach towards the melt region, or liquidus. The temperature distribution is predominantly influenced by the magnitude of latent heat released at the interface.

The boundary conditions for applying a fixed velocity are shown in Figure 15a. To mitigate the competition between neighboring crystals during columnar crystal growth, it is necessary to increase the spacing between initial crystals. At a dimensionless velocity of v ≈ 2.00 (equivalent to an actual velocity of 350 mm/s), the growth of cylindrical crystals induces a discernible deflection toward the approaching side, as shown in Figure 15b. In the convective temperature field, the flow effectively washes the oncoming side, facilitating rapid dissipation of latent heat released during solidification. Simultaneously, a portion of the heat is transported from the oncoming side to the downstream side, resulting in elevated temperatures at the dendrite arm’s downstream region, which hinders diffusion. Consequently, there is an accelerated advancement of the solid–liquid boundary toward the oncoming side and a corresponding inclination of dendrites toward this direction. Additionally, the growth rate of secondary dendrite arms is observed to be higher on the oncoming side. In the convective solute field, as shown in Figure 15c, a distinct disparity can be observed between the number and length of secondary dendrite arms on the upstream side compared to those on the downstream side. Within the solute field, the concentration of solute elements is evident in the roots of secondary dendrites [38]. Notably, when contrasted with the downstream side, thinner solute films are present on the upstream side, accompanied by a larger concentration gradient at the solid–liquid interface front and faster growth of dendrite arms.

In 2021, Peng et al. [39] used the phase field method to simulate dendrite morphology and solute distribution of Ni-Cu alloy under different unidirectional and mixed convection. The results show that convection has a great influence on dendrite morphology and solute distribution. Convective action causes the dendrites to tilt towards the onstream side, and the secondary dendrite arm tip grows faster on the onstream side. The concentration of solute elements can be seen in the roots of secondary dendrites in the solute field. The growth rate of dendrite arms on the onstream side is faster. In the above study, the numerical simulation of the flow-down dendrite solute field morphology is highly consistent with the results obtained in this paper, which verifies the validity of the model.

6. Conclusions

(1) In this study, a macroscopic multi-field coupling model for laser cladding was established. A phase field model was proposed to describe dendrite growth during solidification and to describe dendrite orientation and solute distribution during contrastive flow. Quantization reveals the transient evolution of the macroscopic physical field. Among them, the peak temperature field is 2950 K. The peak stress field is 303 Mpa. The peak velocity of the flow field is 0.35 m/s. The input conditions are provided to simulate the microstructure evolution of the paste zone in the molten pool.

(2) The microstructural evolution during the solidification process was predicted. The phase field method was successfully used to reveal the evolution process of crystal microstructure under temperature and solute gradients. The process consists of the transformation of plane-cell crystal, algal structure, and columnar crystal. Accurately predicted the microscopic features of the partial loss of dendritic structure and local fracture caused by competitive growth, as well as the changes in solute distribution.

(3) The influence of heat source scanning speed and convection on solidification microstructure was analyzed. The transverse section of the melt pool at 45,000 steps with different scanning speeds of v_x = 0 and v_x = 0.028 was compared. The quantitative analysis revealed the dendrite morphology under the non-dimensional flow velocity v ≈ 2.00 (actual speed 350 mm/s). This study showed that the scanning speed of the moving heat source has a negligible effect on the dendrite growth and orientation in laser cladding. However, the internal temperature gradient and surface tension-induced convection are the main factors affecting the dendrite growth, orientation, and solute distribution.