Stress concentrator at the plate-like inclusion tip as an enhancement factor of diffusion flux

The work aims to evaluate stress-assisted effects in the diffusion growth of plate-like inclusions. We estimate the correction for diffusion flux onto plate-like inclusion tip accounting for the mechanical stress produced by phase transition. Also, the external stress effect on the stress-concentrator and upper-bound flux is evaluated. The stress factor results in the shift of solubility measured as the mean concentration in solution at precipitation. The impact of the stress concentrator was estimated for different thermodynamic systems: θ′-phase (Al2Cu) growth in Al-Cu alloy, zirconium hydrides growth, and perlitic transformation. The stress factor could play a sufficient role in nonstationary conditions and is negligible in steady-state approximation. Stress factor can also have substantial effects on the inclusion orientation if an external stress is present: rigid inclusion lengthening parallel to the external stress, while stacks of inclusions can behave differently depending on nucleation barriers competition for shear and hydrostatic stress components.


Introduction
Diffusional growth of second phase particles in metals is an important and well-studied problem in material science. Important applied issues include modeling nucleation and growth kinetics, effects of external stresses and texture on the orientation and shape of inclusions, and others. For this purpose, numerical methods such as phase-field and cluster dynamic are widely applied. In the first case, the diffusional problem is solved in the computational domain, and physical parameters of the material are included in the empirical functional, ascribing free energy to the material locally as a function of the diffusant concentration and its gradient. In the cluster dynamic approach, analytical equations define the material's parameters for nucleation, diffusion flux, and interaction of inclusions.
The plate-like shape of some inclusions in solid bodies is defined mainly by orientation relationships of two crystal lattices and volume dilatation anisotropy in phase transition [1]. A model for the diffusional growth of ellipsoidal inclusion was developed in [2]. Authors of [2] assumed the constant ratio of the ellipsoid axes and diffusant concentration at the interface, ignoring processes related to forming the actual inclusion's shape. The consistent approach to modeling plate-like inclusion should consider the different growth rates of the tip and flat sides. Nucleation of the new surface step (i.e., thickening) is less favorable than the motion of existing ones (resulting in lengthening). Apart from the non-coherent interface formation, nucleation of the new step demands excess energy due to interaction with transformational dislocations with opposite signs. Both mechanical interaction and interaction of transformational dislocation with the concentration field of the diffusant [3] defines the inclusion's tip shape (equivalent curvature radius). The thickness-to-length ratio depends on the energy threshold of nucleation a new interface dislocation. However, the described picture is not universal. Some plate-like inclusions microscopically contain ordered clusters or stacks of smaller precipitates, for example, hydrides in zirconium [4]. The material's ductility under different types of loading depends on the stacks morphology: their orientation, length, and continuity.

Diffusion flux to the inclusion
We consider plate-like precipitates (lamella) long enough to view them in a two-dimensional approximation. The lamella grows due to diffusional flow to its tip, which we determine as the linear sink and a concentrator of mechanical stresses, see figure 1. The diffusion flux onto side faces of the inclusion is neglected for the following reasons. Firstly, solubility near the side faces is higher than near the tip due to interface effects prominent for plate-like inclusions. Secondly, the authors of [5] demonstrated that diffusant concentration near the side is lower than far from the inclusion due to the compressive stress generated by the transformational strains.

Stationary approach
Since the problem is considered cylindrical symmetry, we must choose the outer cutoff radius R max (see figure 1). The diffusion flux depends in a logarithmic way on this parameter; therefore, the distance between precipitates tips can be determined with satisfactory accuracy (see figure 2).
The diffusional task in radial symmetry accounting for the stress field near the concentrator was solved in [7,10] for diffusion flux to the crack tip with tension applied normal to crack plane. The diffusional flux density j is [7,10]: where С-an atomic fraction of the diffusant in the metal matrix (hydrogen in zirconium alloy in [7,10]), D-diffusion coefficient, Ω-the atomic volume of a metal atom, R-gas constant, T-absolute temperature, μ-chemical potential: V-the excess volume per diffusant atom. Diffusion flux per unit length of the tip can be found from the equation (1) in cylindric geometry applying the substitution ¢ = -D C C p V kT exp ( ) as (according to [10]): where R min -inner cutoff radius (see figure 1), C(R min ) is understood here as solubility limit for precipitation, C(R max )-diffusant concentration in solid solution far from the sink. As can be seen from equation (3), the stress field's role is associated with the diffusant concentration change by the factor -D p V kT exp ( )which shifts the solubility equilibrium. If we neglect the mechanical stress field produced by transformation strain (p(r)=0), the equation (3) will transform to the classical sink strength of parallel lines: The mechanical stress field near the plate-like rigid inclusion tip with the thickness h and strain e^due to volume dilatation in the perpendicular direction to the precipitate plate (see figure 1) is the same as produced by the edge dislocation with the Burgers vector e =b h perpendicular to the precipitate plate. The pressure in the planestrain approximation (σ zz =υ·(σ xx +σ yy )): G-shear modulus, υ-Poisson ratio, Equations (3) and (4) depend on both outer R max and inner R min cutoff radiuses. There are two possible lengths that can be chosen as the R min : the inclusion thickness, if it does not emit dislocations, and the plastic zone    (7) and (8): where A=π(1+υ)/4≈1. We can now find the value of integral in equation (3), taking into account equations (5), (6) and (9): We expand the function under integral equation (10) in a series assuming R max ? R min .
Then the flux from equation (3) transforms to: In the current approach, the microscopic stress yield should be applied, which is due to the emission threshold of dislocation misfit loops or other dislocations by the inclusion.

External stress influence
External stresses can influence the orientation of plate-like precipitates in different thermodynamic systems. The work done by external stress during precipitation depends on orientation if we consider uniaxial tension and plate-like inclusion with a single component of the self strain e^(anisotropic dilatation in a more general way). For rigid plate-like inclusion in a matrix loaded in tension perpendicular to its plane, the work done per unit precipitate volume is e ŝ , ext while the work for the same inclusion tensioned parallel to its plane is zero. This assumption is the base, for example, for some models of hydrides reorientation in zirconium alloys and ascribed to nucleation barrier reduction [11][12][13].
However, the hypothesis about hydride reorientation due to different orientation relationships is not the only explanation known in literature [14]. Some other researchers consider the hydride orientation to be the result of micro-hydrides with identical orientation relationships but stacked at different angles to each other depending on external stress [15,16]. The reasons why precipitates can be self-organized in different orders are currently not well understood. In this section, we will apply our approach to estimate the external stress factor on the growth rate of inclusions with different orientations.
The stress field around the rigid ellipsoidal cylinder in the elastic matrix under external load was found in [17] and described in [18]. This result can be applied to find hydrostatic and shear components of the stress field and, thus, diffusion fluxes to different sides of the inclusion with the given aspect ratio. For thin elongated inclusion, equations simplify and become in some way similar to the stress field around the crack [19,20] in coordinates x′, y′, where x′ coincides with the main axis and zero point at the center of the inclusion (see figure 3): Here a is a half-length of inclusion.
The prime symbol marks external stress in the coordinate system x′, y′; parameter κ is equal to k u = -3 4 for plane-strain and k u u = - (14) is valid even for the plastic area around the inclusion in plane strain when υ→0.5 and κ→1, which was confirmed by finite element calculations in visco-plastic approximation [18]. The presence of the inclusion does not disturb the shear stress equal to s ¢ ¢ ¢ x y in coordinates x′, y′. If we consider the finite thickness and finite elastic modulus of the inclusion, the solution will expand in a series for the ratio of matrix and inclusion elastic modulus [19].
Let us consider two rigid plate-like inclusions without transformation strains under uniaxial external tension s . ext The first is oriented along with the tensile stress (s s s . According to equation (13), hydrostatic stresses ratio near the tips is It means that the inclusion aligned with the external tension grows at lower concentrations (faster) than the inclusion oriented normally.
To evaluate maximum shear stress, one must choose the direction of its maximum (ahead the inclusion in . Spatial distribution for stress in x′, y′ is [21]: x a sin 2 1 2 cos 2 cos 3 2 16 where f is the angle between the line connecting the specific point and the inclusion tip and the x′ axis. Differentiating an expression by variable u=cos(f/2) and equating it to 0 resulted in a quadratic equation for u 2 with two roots: with '+' sign stands for maximum value, which is denoted f¢. The value ofΞin equation (16) and the value of equation (17) at different κ are shown in figure 4. Shear stress has its maximum at the angle∼50-60 degrees with the inclusion's plane. The hydrostatic component at this angle is f¢ 4 cos 2 · ( ) / or 10%-15% lower than its maximum value at f = 0. Therefore, it creates conditions for the precipitation of new precipitates in this direction if the nucleation barrier is more sensitive to the shear stress and shear transformation component of mechanical stress than to hydrostatic stress.
It is fruitful to compare it with the previously analyzed self-stress field of rigid inclusion with only self-strain e^, resembling the stress of equivalent edge dislocation. Shear and hydrostatic components of the edge dislocation stresses are relatively narrow and in antiphase: shear component has a maximum at 90 degrees, while the hydrostatic one at 0 degrees. Suppose the nucleation of new hydride precipitates in zirconium is the autocatalytic process where both shear and hydrostatic stresses reduce the nucleation barrier via mechanical interaction with dilatation and shear transformation strains [22] for rigid critical inclusion with volume V crit :  and at f = 90 degrees otherwise. In the axisymmetric approximation of the diffusion problem, the lower sink radius is limited by the cutoff radius where the yield condition is reached: Parameter κ is around 2 for Poisson ratio 1/3 in plane-stress, 1.66 for plane-strain, and 1 for the start of material plastic deformation. According to equation (20) and monotonic behavior ofΞon κ ( figure 4), hydrostatic stresses at the inclusion tip decrease in the following order: plastic deformation, plane-strained elastic, and plane-stressed elastic material.

Nonstationary approach
We can estimate the limiting R max, which does not depend on the precipitates interspacing. The growing lamella eating out the diffusant profile at a nonstationary stage as: is the maximum concentration of the diffusant in the metal. The higher supersaturation leads to the lower R max . Suppose the inclusion grows with the velocity v. In that case, the domain near the tip gets an additional diffusant amount equal to the equation (4).
There are three parallel processes at every time step dt: the inclusion tip movement for dl, diffusion flux to the tip, and the diffusant profile eating out: m a x m a x m i n С · · · ( ( ) ( )) ( ) J stac -the stationary flux according to equation (4). Then, from equation (23): The growth velocity has an upper limit which is determined by the diffusion flux and the task geometry, R max ≈D/v max , and the maximum flux: The same flux is equal to: From equation (26): where r max is the distance between lamellas when the growth velocity reaches its maximum value v max : Note that the equation (28) is valid only for R max <r max , i.e., if neighbor lamellas do not influence each other's growth. The function in equation (28) is plotted in figure 5, which shows that growth velocity reaches its maximum value at R max ≈0.75·r max ; R min was assumed to equal to 0.01·r max for the estimation.

The model application
3.1. θ′-phase growth in Al-Cu We applied the diffusion task solution (equation (12)) to simulate θ′-phase (Al 2 Cu) growth kinetic in Al-4 wt.% Cu alloy studied at different temperatures after quenching in experiments [23]. The θ′-phase inclusion is thin discs growing due to the copper diffusion flux. The phase transition from α-Al (solid solution of Cu in Al) to θ′phase is accompanied by the volume dilatation approximately equal to 5%. The temperature scenario in experiments [23] is quenching from the homogenization temperature to the annealing temperature, which varied. We simulate the precipitates growth during annealing, assuming the diffusion flux to the inclusion equal to the equation (12). The average diameters as functions on the annealing time are shown in figure 6 for temperatures 225°C, 250°C, and 275°C. The model parameters values used in this simulation are collected in table 1. The volume concentration of the θ′-phase nucleus was set in modeling equal to 1.3·10 20 m −3 . This value corresponds to the value (2±0.7)·10 20 m −3 , measured in experiments [24] with a similar temperature scenario (at the end of quenching. The dependence of the average precipitates diameter on experiment time is shown in figure 6; parameter values used in the modeling are in table 1. Simulation for every temperature scenario was performed in three assumptions: accounting for the own stress field (equation (12)), without accounting, and maximal theoretical possible growth velocity (equation (25)). The stress factor accelerates the growth rate up to two times and shifts the growth rate close to the theoretical limit. The ratio R max /r max changes in a wide range during cooling and reaches the value 1 at the highest supersaturations. The θ′-phase thickness was assumed to be constant during the precipitate growth and equal to 10 nm, which is the average value measured in [23] 1 . Figure 6. The size of θ′-phase (Al 2 Cu) inclusions in Al-4 wt.%Cu alloy during annealing at different temperatures after quenching. 1 The thickness kinetic in principle could be simulated if we determine two different solubilities near the inclusion tip and its plane surface, as it was done in simulation of anisotropic growth of the hydride near the crack tip in [33]. The solubility difference at different inclusion surfaces arises due to anisotropy of the surface energy and volume dilatation as described in Introduction. 8

Zirconium hydrides growth
Hydrogen behavior in zirconium alloys significantly impacts the industrial use of this material in water-cooled nuclear power plants. As the secondary product of the water corrosion, hydrogen can penetrate the metal and be accumulated during operation. If hydrogen concentration in solid solution exceeds the solubility limit, the hydrides precipitate in the metal matrix and produce embrittlement risks [34]. The stable hydride phase, which usually arises in practical applications, is FCC δ-ZrH 1.6 [34]. Hydrides precipitate in clusters or stacks as flat dendrites. The hydride morphology is named 'corn flakes' in [35], where it was studied with the synchrotron X-ray tomography technique.
An experimental study [36] of hydrides growth in samples made from Zr-1.0Nb-1.0Sn-0.1Fe with 250 ppm of hydrogen is simulated in this section. The temperature scenario of the experiments is cooling from 400°C with a constant cooling rate of 2°C min −1 down to T quench . After reaching T quench, the temperature was dropped fast, leading to 'freezing' the hydride structure. As a result, new hydrides formed in quenching were much smaller and visually separable from previously formed hydrides.
The simulation of experiments [36] was performed based on the application of equation (12). The used numerical algorithm was the same for modeling Al-Cu thermodynamic system in section 3.1, but with appropriate parameters values shown in table 2. Model-to-experiment comparison can be seen in figure 7; we consider it as satisfactorily. The cooling rate in experiments [36] was slow, and the thermodynamic system is close to equilibrium at every moment. Therefore, the account of own stress has no significant effect (both approximations give effective diffusion sink). Figure 8 shows the result of another experiment simulation [37] with hydrogen-charged zirconium alloy but in substantially different conditions. Specimen made of zirconium-based alloy Zircaloy-4 with 600 ppm of hydrogen was subjected to fast cooling from 570°C to 400°C with the cooling rate of 100°C min −1 followed by annealing at 400°C. The hydrogen concentration in solid solution was measured online during cooling and annealing using the synchrotron X-ray diffraction method [37]. Before cooling started, all hydrogen was in a solid solution, while at the final temperature, around half of hydrogen precipitate as hydrides in the metal matrix.
The concentration of the hydride nucleus was set 10 11 m −3 in the simulation of experiments [37]. That is the relatively low value corresponding to the beginning of the nucleation process (the nucleation itself was simulated nowhere in this paper). Results of the modeling with and without accounting for the effect of self-stress are shown in figure 8. Relative supersaturations are estimated in cases (a) with self-stress, (b) without it, and (c) assuming the maximal diffusion flux equals 1.27, 1.35, and 1.07. The difference may be small; however, it can substantially differ in new inclusions' nucleation conditions. The difference between lines may not be significant. It is better illustrated with the ratio of fluxes J stress to J-with and without account for self-stress shown in figure 9. In the beginning, diffusion acceleration reaches 1.5-2 times. As a more substantial sink results in faster relaxation of supersaturation, after cooling the ratio J stress /J becomes less than 1. The maximum value of R max /r max is around 0.65 during the calculation, but much less for most of the calculation time, which means that the interface velocity is far from the theoretical limit, unlike the conditions of experiments [23] with the Al-Cu alloy described above.

Austenite decomposition
Microstructure and, therefore, mechanical properties of carbon steels depend on austenite (γ-phase of iron) decomposition kinetic, which itself is governed by carbon diffusion. Perlitic transformation of austenite Figure 7. Zirconium hydrides length when cooling from 400°C to T quench room temperature with constant cooling rate 2°C min −1 following quenching from T quench to room temperature. The process direction is marked by narrows. proceeds by nucleation and growth within the γ-phase two new phases (structure called 'perlite'): ferrite (αphase) and cementite (iron carbide Fe 3 C). This process is usually modeled as the growth of parallel lamellas with thickness h separated by distance λ as it is done analytically by Zener [41] or in phase-field modelings (for example, in [6]). The schematic view on the perlite front in austenite is shown in figure 10. Zener and Hillert developed lamella growth velocity in [42] as:    figure 11 together with experimental data [43] and [44]. Parameters used for the assessment are collected in table 3, the ratio h=λ/20 was applied. This value was chosen so that the stress in front of the lamella tip agrees with the calculations performed in [6]. Figure 11 shows that the stress factor is substantial only in the narrow temperature range near the eutectoid temperature, increasing the perlite front velocity up to three times due to a shift of the diffusant concentration and, consequently, precipitation temperature. For lower temperatures, the self-stress does not influence significantly on the process rate.

Conclusion
This paper presented an estimation of the effects of self-stress at the tip of a rigid plate inclusion and external stresses on diffusion sink efficiency. The approach was based on an existed analytical solution for the diffusion equation and can be implemented to the existed calculation tools or applied independently. Based on the simulations for precipitation in systems Al-Cu, Zr-H, and Fe-C, it was demonstrated that the stress factor could play a sufficient role in nonstationary conditions. The perlite growth velocity can differ up to three times if selfstresses are taken into account. For θ′phase (Al 2 Cu) in Al-4wt.%Cu alloy, after quenching, the self-stress accelerates growth up to two times. The effect, in this case, is close to the maximal possible estimation in nonstationary approximation, where inter-precipitate distance limits the available flux. Modeling of hydrogencharged zirconium alloys during fast cooling (100°C min −1 ) revealed that the calculated maximal relative supersaturation of hydrogen was equal to 1.27 and 1.35 in conditions of the simulated experiment with and without accounting for the self-stress factor. It can lead to a substantial difference in nucleation conditions and could change the nucleation frequencies for many orders. Stress factor can also have substantial effects on the inclusion orientation if an external stress is present. For rigid inclusion, it promotes lengthening parallel to tension axis. For a stack of inclusions growing due to nucleation of new micro-precipitates, the external stress field can provide conditions for growth with an angle of 50-60 degrees to the tension axis. The stack behavior depends on nucleation barriers competition for shear and hydrostatic stress components.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Funding
The reported study was funded by RFBR, project number 19-32-60031.