Shock-induced phase transition of iron studied with phase field method

To develop phase field method in high strain rate loading, it derived the governing equations of phase field model in the foundation of entropy functional, and the model adopted a diffusion interface, which made automatic tracking on interface possible. Added the order parameter to characterize the different phases of materials and simulated the soften and phase transition in shock-induced loading. According to the conservation theorems and principle of entropy increasing, we established the governing equations of shock-induced phase transition in metal dielectric. Then the work discretized phase field model by using high order difference scheme, introducing the artificial viscosity, and obtained the numerical solution of the phase field model. Compared the numerical value with the theoretical solution, the accuracy of the developed phase field model was verified. Based on the model described above, the work study the α→ε polymorphic phase transition of iron-base alloy induced by one-dimensional shock wave. Considering the spatial and temporal distribution of each parameter field and the Hugoniot data for iron, it further analyzed the coupling effect of wave and phase transition by tracking the changes of pressure and speed in each unit. The simulation results fitted well with the experimental results, which indicated that the phase field model developed was successful, and it would help people to understand the microstructural evolution of metallic materials under the shock wave.


INTRODUCTION
In the study of solid-solid phase transition in the metal on the impact loading, it is representative and focused attention wildly that the α ↔ ε polymorphic phase transition of iron base material.Tang et al. (Ref.[1]) studied the phase transition and the spalling characteristics induced during shock loading and unloading in an α phase Fe-based alloy, and they found that abnormal spalling phenomena are believed to be related to the shocked α → ε phase transition, and a possible reason for the abnormal spalling is discussed also from the interaction process of stress wave.Zhang et al. (Ref.[2]) conducted light gas gun plate impact experiments on a FeMnNi alloy and phase transition and spallation occurred.Experimental results indicated that the spallation strength of this alloy didn't increase distinctly after phase transition.Liu et al. (Ref.[3]) simulated α-iron phase transition process by using discrete element method (DEM), and they also conducted numerical simulations on α-iron and ___________________________________ *Corresponding Author: 499146173@qq.compolycrystalline iron by DEM, which combined with undiffused two-phase transition model, thermodynamic consistent free energy function, and phase transition kinetics of relaxation equation.
The phase field method is a powerful tool for studying the microstructure evolution of materials, and it has been fully developed in recent decades.In the 1990s, Penrose and Fife (Ref.[4]) introduced a phase field model based on entropy functional.It was a phase field model, which was consistent with classical thermodynamics.Wang et al. (Ref.[5]) deduced the general phase field model of pure substance based on entropy function.Bi et al. proposed the general phase field model of binary alloy.Their work laid the foundation for the unified phase field model.In 2000, the latest polycrystalline phase model proposed by Kobayashi et al. (Ref.[6]) Two multidimensional order parameters were adopted, one of which represented the direction of crystallography, and the other represented the important orientation of grain, thus dealing with grain rotation and migration of anisotropic interfaces.
The martensite phase transition of iron-base alloy under shock wave loading is the main objective of this paper, so we need to establish a shock wave propagation model, which is suitable for the phase field method.In this study, Murnaghan state equation and internal energy function are adopted, which has higher accuracy than Gruneisen equation, and the introduced interface gradient term is used to derive the control equation.So that the phase field model that can reflect the process of shock wave loading in metal materials well would be produced.In Sec. 2 and Sec. 3 we determine the entropy functional form and introduce the gradient phase gradient of density gradient to establish the phase field model of iron shockinduced phase transition.Then the finite difference method is adopted to carry out numerical calculation, and the evolution equation is solved by dispersing the time and space of the partial differential equations.Finally, we use Fortran to write the program, and calculate the case that pressure is 20GPa.In Sec.IV the influence of the viscosity coefficient and interface parameters on the wave propagation is studied.We also study the variation rules of phase change pressure threshold, and relaxation time of phase transition under different pressures.

PHASE FIELD MODEL
To establish a phase field model based on entropy functional thermodynamic consistency, we construct the entropy functional of the system firstly, and make sure that it obeys the dynamics equation and satisfies the entropy increase principle.For different problems, the gradient items introduced in entropy functional are different, resulting in different forms of control equations.When phase transition happening, phase stability is controlled by external factors and internal factors, and the influence of these factors (or field variables) on the phase transition of materials is mainly reflected in the places where the field quantities change greatly.Therefore, we should introduce the gradient of internal energy and density in order to study its contribution to system entropy.Considering the purpose of this article, which is researching the impact of the metallic iron phase change process, in addition to the introduction of containing density gradient item, we also introduce the order parameter of gradient ∇.When the values of  is 0 or 1, it denotes respectively in the two-different homogeneous phase, the initial phase or in the transition phase.When its value is between 0 and 1, it denotes in a two-phase interface.We assume that the quality of the metal M, momentum P, energy E and system entropy S satisfy the following form in the volume Ω().

𝑀𝑀 = � 𝜌𝜌 𝑑𝑑𝑑𝑑 𝛺𝛺(𝑡𝑡)
(1) In order to establish the phase field model of shock wave loading in iron, we need to derive some relational preparation according to the laws of fluid mechanics and thermodynamics, including a relationship between the introduced interface gradient items.Then according to the law of conservation of mass, the momentum theorem, the law of conservation of energy, and system entropy increase principle, we got the final control equation.
where the parameters were taken as follows.
In the process of shock-induced phase transition in metal, the state equation of solid phase is usually applied to Gruneisen state equation, but in this paper, a more accurate Murnaghan state equation was used to describe the shock-induced phase transition of iron.The EOS is Where  0 and  are fitting constants, which are related with the isothermal bulk modulus,  0 is the initial volume.The entropy function and the internal energy expression could be obtained respectively.
The change rule of mixing phase in impact phase transition was also studied.So further considering two-phase interface transition zone, we put forward the hypothesis that specific volume and internal energy could satisfy the following relations in the whole area, Where ℎ() is a function of , and its values should be satisfied that ℎ(0) = 0, ℎ(1) = 1, ℎ ′ (0) = ℎ ′ (1) = 0, and when 0 <  < 1, the values of function meet requirement that 0 < ℎ() < 1.
Then the kinetic control equations were discretized and solved numerically, so that we obtained the distribution of each field variable in the system.In view of the simplicity and flexibility of finite-difference method and its applicability, and the program is concise and easy to use computer for calculation, so we decided that the finite difference method was used to discretize the control equation of the phase field model of shock-induced phase transition in iron.To simplify the calculation during the discrete process, the numerical solution would be solved only from one dimensional discrete format point of view.We used the McCormack method for dispersion.This method is an explicit difference algorithm and uses two-step estimation correction process.The forward difference is used in the estimation step, and backward difference is used in correction step, which make it have second order precision.The McCormack method avoids the calculation of the second derivative, and that simplifies the calculation.Otherwise we noticed that the parameters varied dramatically near the shock wave front, to get a smooth and stable solution, we need to add a specific artificial viscosity to the McCormack method.
After defining the initial condition and boundary condition, the local threshold condition, which controls the phase change evolution, was also considered during the calculation.Local threshold condition means that when Gibbs free energy difference of a lattice exceeds that of activating phase transition, phase transition will be triggered immediately.If the free energy difference does not reach the phase change activation energy, then the phase transition condition is not met, and meanwhile the value of  remain the same.

RESULT FOR 𝜶𝜶 → 𝜺𝜺 TRANSITION
The phase field model obtained in Sec. 2 was used to study the impact phase transition in metal medium.We assumed a metal iron of 10mm in length, and there was a constant loading which is 20GPa on the left side.Under the stress condition, the iron occurred  →  phase change.We used Intel Visual Fortran 2013 to write the calculation program.Through the visualizations, it could be observed and analyzed the distribution of the density, temperature, and internal energy of the shock wave in the iron.The parametric field distribution in Fig. 1 showed that the density, temperature and internal energy had two square wave structures.The value of  where the first wave located correspondly was zero, which meant that the wavefront was the wavefront of shock wave, and the physical quantities on both sides of the interface were abrupt.The phase transition did not happen yet.When the value of  was between 0 and 1, the density, temperature and the internal energy jumped again, and a second square wave occurred.That indicated the impact phase transition would cause the parameter to change again.Because shock wave traveled faster than phase transition, the shock wave would pass to the right edge soon.If the shock wave reached the right boundary, the parameters would oscillate violently, and the shock wave would reflect.
In general, a typical double wave structure will appear in the system under the condition of 20GPa loading.The shock wave front is separated from the phase change interface, and the velocity of the shock wave front is faster than the phase transition interface.This conclusion is consistent with result obtained by experiment and theoretical analysis.As the phase transition interface goes further, we can infer that the value of  will turn from 0 to 1 in the right area, which means all iron have transformed from  phase to  phase.We can see the propagation process of impact and phase transition interface in metal iron more directly in Fig. 2.

ANALYSIS 4.1. Principal Hugoniot
The Impact adiabatic curve, also called Hugoniot, is the focus of the analysis of thermodynamic equilibrium state (Ref.[7]).In order to verify the reliability of the phase field model we established, the results obtained in the previous chapter are compared with the experimental results.It can be seen from the Fig. 3 that the calculated results are in error with the experimental data and theoretical values.This is because the mixing phase is in a stable state during the phase transition, so that the mixed phase platform is less obvious.Overall, the anastomosis was good, and it is proved that the phase field model we established can reflect the impact compression and phase transition of iron.
Then considering the relation of shock wave velocity  and particle velocity after wave   under different pressures, it can be seen from the Fig. 4 that the results obtained through the simulation of the phase field model, which was established by us, are in good agreement with the experimental results.The calculation results indicate that when the impact pressure is less than the phase change pressure threshold, the particle velocity is also small.At this time, there is only one wave in the medium, and the velocity of shock wave is approximately linear with the velocity of wave.When the impact pressure increases to exceed the phase change pressure threshold, the double wave structure occurs in the medium (shock wave and phase change wave).At this time, the velocity of the shock wave increases slowly, and it is approximately maintaining a steady wave velocity propagating within the medium.Meanwhile, the phase change wave velocity increases rapidly with the increase of particle velocity.When the pressure reaches 36GPa, the phase change wave velocity is equal to the shock wave velocity.Thus, it can be inferred that if the impact pressure continues to increase, the phase change wave velocity will catch up with the shock wave velocity.The phase change wave and the shock wave merge into a stable shock wave, and they propagate forward together.Next, we further explore the phase change pressure threshold under different pressure.According to the difference of Gibbs free energy needed for the transition, we can calculate a critical pressure to satisfy local threshold condition of free energy at certain temperature.We call this critical pressure phase transition stress threshold.In table 1, the temperature and critical pressure value of impact phase transition is given, which is the relation between temperature and phase transition pressure threshold.As seen from the table, as the temperature increases, the threshold of phase transition pressure decreases, which is consistent with the experiment, and it proves the correctness of the numerical calculation.
Fig. 5 shows the part of phase diagram of the iron under quasi-static loading conditions.It can be seen that there is some error in the result, but the trend and the data under higher temperature are consilient better.It proves the phase field model we established can more accurately describe the phase transition process in iron medium.

Phase Transition Relaxation Time
Relaxation time refers to the time required for the system to reach a stable state from an unstable state, namely the time required to achieve the thermal equilibrium.For metallic martensitic transformation, the physical substance of relaxation time is nucleation time.In this paper, the relaxation time means the time for  phase converting to  phase completely.We took two space lattices to study relaxation time of phase transition by observing the value of .
For the relaxation time of phase transition in iron under 20GPa, the first observation made by Forbes (Ref.[12]) was that it was less than or equal to 50ns.But Barker revised the data to about 180ns by observing with VISAR.You can see the calculation results of this article from Fig. 6.After 60 to 80 time steps (2ns per step), the value of  changed from 0 to 1, which means that phase transition was complete.This results are consistent with the experimental results.Then the influence of different loading pressure on the relaxation time was further considered.It can be seen from Fig. 7 that the pressure has an obvious effect on the phase transition relaxation time.As the impact pressure increases, the phase change relaxation time decreases, indicating that the rate of phase transition increases.This conclusion is consistent with the results obtained by Boettger (Ref.[13]).

Interface Parameters
According to the specific form of the control equation and the thermodynamic parameters determined, we have selected certain gradient coefficient and interface parameters for numerical calculation.These coefficients are closely related to the propagation of shock wave and phase transition interface.Therefore, the value of each gradient coefficient would be adjusted, and we would study the influence of its value on the phase transition interface.And the rationality of the value we have used in the calculation would also be verified.
First of all, keeping other parameters constant, we research effects of parameter , which is related to the phase change rate, on shock wave propagation characteristics.
Different  almost had no effect on shock wave propagation in metallic iron.However, it could be seen that the phase transition diffusion interface have changed significantly.With the decrease of , the starting point of the phase transition was almost the same.But the end point of phase transition moved backward, namely the diffusion interface thickness of phase transition increased.This indicated that the two-phase mixing region increased, and the rate of change for parameters gradient on both sides of the interface slowed down.It also showed that decrease of  made the phase transition rate go down, phase transition relaxation time become longer.If we continue to reduce , there will appear unreasonable results.So the values should be controlled in a reasonable range.
Different values of   had little impact on shock wave propagation in metallic iron.However, the parameters near the phase change diffusion interface were affected.It could be seen from Fig. 10.(a) that there was no obvious change in the thickness of the phase change diffusion interface, and the starting and ending point of phase transition didn't move either.But as   value increasing, the value of  increased at the same cell in the mixed phase.It indicated there were more iron in  phase, and completion of phase transition increased.It could be seen from Fig. 10.(b) that as   values increasing oscillation of phase transformation zone was reduced, but there was a relatively stable extremal.The influence of   on shock wave and phase interface was almost the same.
The influence of different values of   on the propagation process and phase transition interface was similar with that of   .The difference was the oscillation degree of shock wave front and phase transition zone did not change with the variation of   , and platform zone after shock wave had little influence.

CONCLUSION
Through the study on the influence of viscosity coefficient and interface parameters on the shock wave propagation and phase transition, the rule was summarized, and the rationality of the value was verified.Secondly, the calculation results obtained by numerical simulation in the previous chapter were analyzed, combined with the existing theoretical data and experimental results.Then we studied the change rule of phase transition pressure threshold with temperature.Finally, through value of  at space lattice to study phase transition relaxation time.From various angles we confirmed the validity of the calculation and the reliability of phase field model we established.The results show that the phase field model can accurately describe the real physical process of metallic iron impact transition.
The following conclusions can be obtained through analysis.a) All parameters, such as ,   ,   et al, are closely related to shock wave propagation, and the value of each parameter can affect the phase transition interface and the shock wave front.Therefore, it is necessary to rationally value to describe the phase transition process by establishing a good phase field model; b) The increase in temperature can result in a decrease in the phase transition pressure threshold; c) When the shock pressure is less than the phase transition pressure threshold, the velocity of the shock wave is approximately linear with the velocity of the particle after wave.When the impact pressure exceeds the phase transition pressure threshold and no more than 36GPa, the shockwave velocity increases slowly, and approximately maintains a steady wave velocity propagating within the medium.The phase change wave velocity increases rapidly with the increase of particle velocity (that is, the increase of impact pressure); d) Increasing of pressure can make phase transition rate increase, and phase transition relaxation time decrease.

Fig. 6 .
Fig.6.The curve of value of  with time in two lattices under 20GPa.

Fig. 7 .
Fig.7.The values of  with time in certain lattice under different pressures.

Table . 1
. The relationship between temperature and pressure threshold.