Computational study on microstructure evolution and magnetic property of laser additively manufactured magnetic materials
Min Yi Bai-Xiang Xu Oliver Gutfleisch





Computational study on additive manufacturing of magnetic materials        



M. Yi et al. 

M. Yi , B.-X. Xu, O. Gutfleisch Institute of Materials Science, Technische Universität Darmstadt, Darmstadt 64287, Germany  

              Tel.: +49 6151 16-22922

              Fax: +49 6151 16-21034

              yi@mfm.tu-darmstadt.de; xu@mfm.tu-darmstadt.de     

 
 


Computational study on microstructure evolution and magnetic property of laser additively manufactured magnetic materials
    Min Yi Bai-Xiang Xu Oliver Gutfleisch
    Received: July 10, 2018 / Accepted: 
=========================================================================================================================




IAdditive manufacturing (AM) offers an unprecedented opportunity for the quick production of complex shaped parts directly from a powder precursor. But its application to functional materials in general and magnetic materials in particular is still at the very beginning. Here we present the first attempt to computationally study the microstructure evolution and magnetic properties of magnetic materials (e.g. Fe-Ni alloys) processed by selective laser melting (SLM). SLM process induced thermal history and thus the residual stress distribution in Fe-Ni alloys are calculated by finite element analysis (FEA). The evolution and distribution of the γ-Fe-Ni and FeNi_3 phase fractions were predicted by using the temperature information from FEA and the output from CALculation of PHAse Diagrams (CALPHAD). Based on the relation between residual stress and magnetoelastic energy, magnetic properties of SLM processed Fe-Ni alloys (magnetic coercivity, remanent magnetization, and magnetic domain structure) are examined by micromagnetic simulations. The calculated coercivity is found to be in line with the experimentally measured values of SLM-processed Fe-Ni alloys. This computation study demonstrates a feasible approach for the simulation of additively manufactured magnetic materials by integrating FEA, CALPHAD, and micromagnetics.







§ INTRODUCTION

Fe-Ni permalloys are typical soft magnetic materials with extraordinary magnetic, mechanical, and electrical properties <cit.>. Due to their low coercivity, high magnetoconductivity, high permeability, and moderate saturation magnetization, they are of great interests for applications in electromagnetic devices, including transformers, sensors, and electric motors <cit.>. In order to realize these applications, suitable manufacturing techniques have to be identified since they significantly affect the magnetic properties. In the past, numerous conventional manufacturing methods such as sintering, thermal spraying, ball milling, and magnetron sputtering have been used to obtain the desirable performance of Fe-Ni alloys. Nevertheless, within the scope of these methods, the direct consolidation of different types of powders into bulk magnetic components with magnetism preserved is always challenging.
Moreover, these conventional methods may lead to the decrease of magnetic properties due to the excessive grain growth under a low-speed heating and cooling. They are also weak in producing precise magnetic components with complex shape and geometry.

Selective laser melting (SLM), as a typical additive manufacturing (AM) technique, enables the quick production of complex shaped three-dimensional (3D) parts directly from metal powders. Up to now, a large number of studies about SLM-AM or electron-beam-AM have been focused on structural materials with mechanical properties as the focus, such as aluminium alloys <cit.>, Ti-Al-V alloys <cit.>, Ni-based superalloys <cit.>, stainless steel <cit.>, etc. In contrast, the application of SLM-AM to functional materials is still in its infancy. Nonetheless, SLM-AM undeniably provides a promising route for breaking the bottlenecks of traditional techniques to fabricate complex shaped functional and miniaturized magnetic devices or systems directly from metal powders. Ongoing efforts have been devoted to the production of magnetic materials by SLM-AM.
The initial work was carried out on the SLM processing of magnetic Fe-Ni alloy by Zhang et al. <cit.>. Depending on the composition and processing, Fe-Ni alloy can be either a soft magnetic material or in the L1_0 phase as a rare-earth-free alternative for permanent magnets <cit.>, thus making Fe-Ni alloy a very interesting material. Later, Moore et al. fabricated magnetocaloric La(Fe,Co,Si)_13 geometries by SLM <cit.>.
However, after the early work <cit.> in 2012 and 2013, it is found from the literature survey that few studies followed. Only lately in 2016-2018, studies continue with focus on SLM processed magnetic materials such as Fe-Si alloy <cit.>, Fe-Si-Cr alloy <cit.>, Fe-80%Ni permalloy <cit.>, Fe-30%Ni alloy <cit.>, Ni-Fe-V and Ni-Fe-Mo permalloys <cit.>, Fe-Co-1.5V soft magnetic alloy <cit.>, permanent magnets including NdFeB <cit.> and AlNiCo <cit.>, ect. Electron beam melting (EBM) is also tried to produce MnAl(C) magnets <cit.>.
Apart from the SLM and EBM based AM technique, other 3D printing technologies without high energy input and high temperature, such as binder jetting and material extrusion, are recently applied to the production of polymer-bonded magnets <cit.>.
These experimental studies reveal the notable effect of AM process on the microstructure and magnetic properties of magnetic materials, and provide insight into the challenges for the design and control of magnetic properties by AM.

Despite of these recent experimental efforts, no literature is found about the modeling and simulation of the fabrication of magnetic alloys by SLM-AM. Almost all computational studies are dedicated to the structural materials by SLM-AM  with a focus on the temperature, microstructure, residual stress, strength and ductility <cit.>, possibly driven by the related experimental contributions which are continuously flourishing. As for fabricating magnetic alloys by SLM-AM, numerical simulations are also essential for the optimization of SLM-AM processes without intensive and expensive trial-and-error experimental iterations, as well as for the understanding of underlying physical phenomena which are difficult to observe experimentally.

In this work, taking magnetic Fe-Ni alloy as a model material, we attempt to computationally predict the microstructure evolution and coercivity of SLM processed magnetic materials through the integration of finite element analysis (FEA), CALculation of PHAse Diagrams (CALPHAD), and micromagnetic simulations. Temperature history and distribution were calculated by FEA within the framework of heat transfer. By using the temperature information as the input, thermomechanical simulation by FEA were performed to get the residual stress distribution. Furthermore, integrating temperature results with CALPHAD output resulted in the temporal evolution of liquid, γ-Fe-Ni phase, and FeNi_3 phase. Finally, by incorporating the residual stress into the magnetoelastic energy of micromagnetics, the magnetic hysteresis and coercivity were calculated.
It is anticipated the computational study could provide a possible general routine or procedure for enlarging the process understanding of the underlying physical mechanisms in the SLM processed magnetic materials.





§ THERMAL ANALYSIS

The fabrication of magnetic Fe-Ni alloy by direct SLM processing of powders is illustrated in Fig. <ref>. An Fe-Ni alloy substrate with a dimension of 360 μm × 300 μm × 250 μm is chosen for the additional layer-by-layer growth of new Fe-Ni alloy layers. One of the most important features of SLM-AM is the complex temperature history generated by the laser irradiation. Predicting the temperature history forms the foundation for the subsequent simulation of residual stress, microstructure, and magnetic property. Using the commercial FEA code ABAQUS <cit.>, here we design and implement a non-linear transient thermal 3D model to obtain the laser induced global temperature history.

The governing equation for the energy balance of heat transfer in the SLM process is given as

    k(ϕ,T)T_,ii = ρ(ϕ,T)C_p(ϕ,T)dT/dt,

in which T is the temperature, ρ is the material density, C_p is the specific heat capacity, t is the time, and k is the thermal conductivity. The initial condition for Eq. (<ref>) is

    T(x_i,t_0)=T_0,

in which T_0 is the ambient temperature 300 K. The temperature of the substrate bottom surface is set as constant T_0=300 K. On other surfaces, the thermal flux includes convection part q_con and radiation part q_rad which can be given as

    q_con=h_c(T)(T-T_0)

and

    q_rad=σ_sbϵ_sb(T^4-T_0^4),

respectively. In Eqs. (<ref>) and (<ref>), h_c is the temperature dependent convective heat transfer coefficient, T is the temperature of the corresponding surface, σ_sb is the Stefan-Boltzmann constant, and ϵ_sb is the surface emissivity.

In Eq. <ref>, ϕ is a field variable to indicate the material state, i.e., whether the material has ever gone beyond the liquidus temperature T_l.
Each element stores its temperature T and ϕ. We set ϕ=0 for a powder state and ϕ=1 for a bulk state. ϕ is designed to realize the irreversible melting process from powder to bulk state by the subroutine USDFLD of ABAQUS. The powder elements are initialized with ϕ=0. ϕ is changed from 0 to 1 upon melting and will retain 1 afterwards, i.e. the fused material can never go back to powder. ϕ of the substrate elements is initialized and always remains as 1.

The material parameters ρ and C_p are determined by the CALPHAD approach which is capable of predicting thermodynamically consistent properties. In the CALPHAD model, the Gibbs free energy per gram of one phase in a multicomponent system can be expressed as

    G_1g(P,T)=∑_i c_i G_i^0 + RT∑_i c_i  lnc_i + G^excess,

in which c_i is the composition of element i in the multicomponent system, G_i^0 the Gibbs free energy of pure element i, R the gas constant, and G^excess the excess Gibbs energy of mixing. Once the Gibbs energy is obtained from CALPHAD data, all other thermodynamic properties can be derived. For example, ρ and C_p can be calculated as

    ρ(T)=1  g/V_1g(P,T)=1  g/(∂ G_1g/∂ P)_T

and

    C_p(T)=-T (∂ G^2_1g/∂ T^2)_P,

as shown in Fig. <ref>(b) and (c). The enthalpy per gram can be derived as

    H_1g(T)=G_1g - T (∂ G_1g/∂ T)_P.

The latent heat L due to the change in enthalpy Δ H_1g of the system during the entire solid-liquid phase change can be calculated from Eq. (<ref>) as 291.14 J/g. It should be mentioned that according to the phase diagram of Fe-Ni alloy <cit.>, the solidus temperature T_s of permalloy with the composition around Fe_20Ni_80 is about 1,709 K, which is only 0.2 K lower than T_l. Therefore, a constant L is taken here. For the system with a wide temperature region of solid-liquid mixture, Δ H_1g is a function temperature and the latent heat effects are often included in the temperature dependent effective specific heat <cit.>.

The thermal conductivity of powder (k_p) is usually very small and here is assumed to be 1% of that of bull material (k_b) and slightly increase with temperature before melting. The powder density is assumed to be half of the bulk density. The powder specific heat is set the same as the bulk one. The temperature dependent material parameters used for thermal analysis are shown in Fig. <ref>(a)–(c).

The interaction between the top surface and laser beam is simulated by a moving surface heat flux with a Gaussian distribution, i.e.

    q_a=2η P_a/π R_a^2exp[-2∥𝐫-𝐫_0(v_a,t) ∥^2/R_a^2],

where η is the powder bed absorption coefficient with an assumed value of 0.5, P_a is the laser power, R_a is the laser beam radius, 𝐫 is the coordinate of the point in the material, and 𝐫_0 is a function of laser moving speed, and v_a is the coordinate of laser beam center. The moving laser heat flux is dependent on the scanning strategy and can be realized by the subroutine DFLUX of ABAQUS. If not specified in the following, the laser beam parameters are chosen as P_a=100 W, R_a=50 μm, and v_a=0.4 m/s, according to the experimental work <cit.>. This laser parameter may be difficult to be realized in industrial applications. Here we limit ourselves to the feasibility of the proposed computational scheme and will not focus on its application to the real industrial AM at the current stage.





Fig. <ref> shows the thermal history results for the single-track scan along the middle line perpendicular to y axis. The powder layer is 50 μm thick and the substrate/powder model is discretized by hexahedral FE meshes, as displayed in Fig. <ref>(a).
Fig. <ref>(b) presents the temperature profile around the laser center at t=0.6 ms, as well as the temporal evolution of temperature at three FE nodes, e.g. one on the powder surface, one in the powder interior, and one in the substrate/deposit interface. It can be found that the powder at all these three nodes is quickly heated up to temperature above the melting point and then gradually cools down to the room temperature. The melting of node III in Fig. <ref>(b) ensures the good connection between the substrate and the deposited layer. By measuring the slope of the line connecting the maximum temperature of the peak to the temperature at 5 ms, an average cooling rate in the order of 10^5 K/s can be obtained. Analogously, the average heating rate can be estimated to be in the order of 10^6 K/s. These estimated rates indicate the feature of fast heating and cooling during SLM.

By examining the temperature profile at a certain time (e.g. t=0.6 ms), the molten pool geometry can be obtained, as shown in Fig. <ref>(c). The length and depth of the molten pool is estimated as 226 and 54 μm, respectively. The molten pool appears like a comet tail, whose asymmetry could be attributed to the laser movement, as well as the temperature and material state dependent thermal conductivity. In front of the laser, the material is in powder state with low thermal conductivity, thus leading to slow heat transfer and high temperature gradient. On the contrary, in rear of laser the bulk material state possesses higher thermal conductivity and wide temperature distribution. A similar melt pool geometry is also reported in literature on nonmagnetic materials <cit.>.

In contrast to the single-track scan for a strip-like material, the in-plane and out-of-plane multi-track scans are simulated to build a one-layer and multi-layer bulk material, respectively. The associated thermal results are presented in Fig. <ref>. 
During the multi-track scanning process, an idle time of 5 ms between the completion of one track and the beginning of the subsequent track is assumed. The idle time is demonstrated to be important <cit.>, but its optimization is out of the scope here. 
Cyclic heating and cooling is remarkable during the multi-track scanning process, as shown in Fig. <ref>(b) and (d). The temporal evolution of temperature at the interfacial nodes between the first and second scanning track (nodes 3198 and 7844 marked in Fig.<ref>(a)) indicates four heating-cooling cycles. Especially, these two material nodes experience notable melting in the first scan and remelting in the second scan. The remelting means that during the second scan the molten pool can extend to the previously deposited track, resulting in good inter-track bonding. During the third and fourth scans, although these two nodes do not melt again, heating and cooling with a temperature change around 500 K still occurs and may raise debonding and thermal fatigue issues.
For the out-of-plane layer-by-layer multi-track scan, the continuous addition of powder is considered by using the element deactivation and activation, i.e. successive discrete addition of new elements into the new scanning track at the beginning of each time step. By using the heat accumulation effect in SLM process <cit.> in which the heat stored in the previous layer affects the next processing layer and induces overheating, the laser power can be varied layer by layer, i.e. large powder for the initial layers and small power for the subsequent layers. Fig. <ref>(c) and (d) presents the typical thermal results in the case of 100 W for the first layer, 60 W for the second layer, and 50 W for the other five layers. The temperature evolution of nodes (marked in Fig. <ref>(c)) at the surface of each layer in Fig. <ref>(d) shows melting and remelting process, indicating the possibility of inter-layer bonding and the integration of deposited layers into a bulk material.


§ MECHANICAL ANALYSIS

The mechanical analysis is subsequently performed independently, since it is reasonable that the mechanical response has a neglectable effect on the thermal history, and the thermal and mechanical analyses are weakly coupled. The analysis is based on the thermal history dependent quasi-static mechanical model, which takes the above thermal results as thermal loads. The governing equation for the stress equilibrium is

    σ_ij,j = 0,

in which σ_ij is the stress. For the mechanical boundary condition, the rigid body motion is restricted and the substrate bottom surface is free to deform. The mechanical constitutive law can be given as

    σ_ij = C_ijkl(ϕ,T)ϵ_kl^e,

where the elastic tensor C_ijkl can be expressed by temperature dependent elastic modulus (Fig. <ref>(e)) and a Poisson ratio of 0.33 for the isotropic material behavior considered here. The total strain is decomposed into elastic strain ϵ_ij^e, plastic strain ϵ_ij^p, and thermal strain ϵ_ij^T, i.e.

    ϵ_ij = ϵ_kl^e + ϵ_kl^p+ ϵ_kl^T.

Thermal strain is given by ϵ_ij^T=α(ϕ,T)(T-T_0)δ_ij in which α is the thermal expansion coefficient, T_0 is the initial temperature, and δ_ij is the Kronecker delta.
α is calculated by the CALPHAD approach through the Gibbs free energy in Eq. (<ref>), i.e.

    α(T)=1/31/V∂ V/∂ T = 1/31/(∂ G_1g/∂ P)_T∂^2G_1g/∂ P ∂ T,
 
as shown in Fig. <ref>(d). For the computation of plastic strain, the linear isotropic hardening model and von Mises yield criterion are used. The yield function is computed as

    f(σ_ij,σ_Y^0,σ_Y^h) = σ_mises - (σ_Y^0 + σ_Y^h).

In Eq. <ref>, σ_mises is the von Mises stress calculated from the stress tensor σ_ij. σ_Y^0(ϕ,T) is the initial yield stress without equivalent plastic strain, as shown by the left curve in Fig. <ref>(f). σ_Y^h represents the hardening and linearly correlates with the equivalent plastic strain ϵ_e^p through the hardening coefficient E^h(ϕ,T) (right curve in  Fig. <ref>(f)), i.e. σ_Y^h=E^hϵ_e^p. The plastic strain is computed by combing the yield criterion in Eq. <ref> and the Prandtl–Reuss flow rule.





Fig. <ref> gives the results of residual stress in the case of single-track scan. The distribution of residual stress σ_11 along the scanning direction is shown in Fig. <ref>(a). It is apparent from the contour plot that compressive σ_11 appears in the deposited layer and tensile σ_11 in the substrate region close to the deposited layer. Fig. <ref>(b) shows the distribution of through-thickness residual stresses and plastic strain ϵ^p_11 (along the line 0 marked in Fig. <ref>(a)). It can be seen that both the stress components σ_11 and σ_22 are compressive in the deposited layer, but change from tensile to compressive in the substrate. The stress component σ_33 through the thickness is relatively small.
Furthermore, the residual stress distribution along x direction in the midplane is examined in terms of the 4 lines defined in Fig. <ref>(a).
It can be seen from Fig. <ref>(c) that σ_11 gradually changes from compressive along line 1 to tensile along line 4.
The compressive σ_11 on the free surface of the deposit are caused by the steep temperature gradient, i.e. the expansion of the hotter top-layer material is prohibited by the underlying material with much lower temperature. In addition, the thermally induced plastic strain should be responsible for the residual stress; because pure elasticity with homogeneous material parameters under no external constraint will not generate residual stress after cooling down to a uniform temperature. The distribution of plastic strain ϵ^p_11 in the deposited layer, as shown in Fig. <ref>(b), also favors the compressive σ_11 after cooling down to the room temperature. The tensile stresses in the substrate/deposit interface can be attributed to the cooling down of the molten material <cit.> and the self-balance of the whole structure. Generally, compressive residual stresses in the top part of the deposit are favorable for increasing the load resistance and preventing crack growth. But tensile residual stresses in the bottom part of the deposit are disadvantageous since they could reduce the load resistance and accelerate crack growth.

The residual stress distribution in the multi-track scan is in Figs. <ref> and <ref>. For the in-plane four-track scan in Fig. <ref>, the stress distribution in Fig. <ref>(b) and (c) is similar to that in the single-track case in Fig. <ref>. But the tensile stress σ_11 in the substrate is lower.
For the out-of-plane layer-by-layer multi-track scan in Fig. <ref>, the residual stress is even much lower. The through-thickness stress distribution in Fig. <ref>(c) shows a average residual stress around 50 MPa, much smaller than that in the single-track and in-plane multi-track scan. The reason could be related to the partial relief of stress under reheating and cooling during by the subsequent laser scanning for depositing the adjacent layers.
The cyclic heating and cooling during the multi-track scan also result in cyclic stress history, as depicted by Fig. <ref>(d) and  Fig. <ref>(c). It can be seen that all stress components are almost zero when the liquidus temperature is reached. Most importantly, σ_11 and σ_22 evolution in Fig. <ref>(d) and  Fig. <ref>(c) manifest that the marked nodes experience somewhat cyclic tension and compression along x and y directions. Since the marked nodes are in the interface of adjacent layers, the cyclic tension and compression could weaken the interface bonding, or even lead to interface failure.



It should be mentioned that in the above thermal and mechanical analysis, the calculation methodology for the magnetic FeNi material is similar to that for the conventional alloys. No special treatment is proposed to deal with the magnetic contribution to the temperature and stress/strain. This is an approximation which is reasonable due to the following two aspects. Firstly, the influence of magnetic properties of FeNi on the heat-transfer thermal analysis is negligible. Secondly, the magnetostrictive coefficient of FeNi is in the order of 10^-7 to 10^-6, which is so small that the effect of magnetization on the stress/strain can be neglected when compared to the effect of thermal expansion <cit.>. So the mechanical analysis can be performed by using the similar method for conventional alloys. However, if one deals with giant magnetostrictive materials (e.g. Terfenol-D with 10^-3 order of magnitude of magnetostriction), the magnetization contribution in the stress/strain calculation cannot be ignored <cit.>.



§ MICROSTRUCTURE EVOLUTION

Microstructure plays a critical role in the property of the products processed by SLM-AM and is required to be predicted if possible. Here we attempt to predict the microstructure evolution in Fe-Ni permalloy during SLM process by using the thermal history and the temperature-dependent phase fraction estimated from CALPHAD. The CALPHAD method is capable of predicting not only thermodynamical properties for material design, but also microstructural evolution through comprehensive physical models of materials processing <cit.>. For example, recently CALPHAD has been combined with phase-field simulation <cit.> and heat-transfer simulation <cit.> to predict the process-phase relationships.





Based on the thermodynamic data of Fe-Ni alloy in the CALPHAD software Thermo-Calc <cit.>, we can predict the phase at any temperature for a given Fe-Ni composition. Since in present work we are interested in the magnetic Fe-Ni permalloy whose composition is around Fe_20Ni_80, only the Thermo-Calc calculated results around Fe_20Ni_80 are shown in Fig. <ref>. The phase distribution of permalloy region in Fig. <ref> is partitioned into four parts by four curves f_1 (T), f_2 (T), f_3 (T) and f_4 (T), which can be either fitted by piecewise-smooth functions or directly used as scattered data from Thermo-Calc output. For numerical study here, we extract scatter data from these four curves and interpolate f_i(T) values at any temperature. The peak temperature is around 787 K and the corresponding mole fraction of Ni (x_Ni) is around 0.72.

As a first attempt and for simplicity, here we only consider FeNi_3 and γ-Fe-Ni phases. The real experimental case should be more complicated, whose comprehensive modeling cannot be achieved within one step and will be continuously explored in the near future. In this way, we can find from Fig. <ref> that γ-Fe-Ni phase exists at temperature between 761 K and 1709 K for x_Ni=0.8. Above 1709 K, only liquid exists. In the region bounded by curves f_2(T) and f_3(T), only FeNi_3 exists. When x_Ni is 0.72, single phase appears at any temperature and the stoichiometric is the same as the powder. In the shadow region bounded by curves f_1(T) and f_2(T), and curves f_3(T) and f_4(T), the coexistence of FeNi_3 and γ-Fe-Ni phases occurs. For the phase-coexistence region with 0.72<x_Ni<0.90, by using the lever rule, we can calculate the FeNi_3 phase mole fraction Ψ_FeNi_3 and the γ-Fe-Ni phase mole fraction Ψ_γ-Fe-Ni by applying functions f_1(T) and f_2(T). Similarly, for the case 0.52<x_Ni<0.72, Ψ_FeNi_3 and Ψ_γ-Fe-Ni can be calculated from functions f_3(T) and f_4(T). Then if x_Ni of the initial powder is given, at any time t the level rule reads


    Ψ_FeNi_3[T(t)] =  {[ f_1[T(t)]-x_Ni/f_1[T(t)]-f_2[T(t)],   0.72<x_Ni<0.90;;    x_Ni-f_4[T(t)]/f_3(T)-f_4[T(t)],   0.52<x_Ni<0.72; ].

and

    Ψ_γ-Fe-Ni[T(t)] =  {[ x_Ni-f_2[T(t)]/f_1[T(t)]-f_2[T(t)],    0.72<x_Ni<0.90;; f_3[T(t)]-x_Ni/f_3[T(t)]-f_4[T(t)],    0.52<x_Ni<0.72. ].

 
By using Fig. <ref>, Eqs. <ref> and <ref>, and the temperature evolution available at each FE node from the thermal analysis, the concurrent and spatially varying microstructure evolution could be predicted. Taking the initial powder with composition x_Ni=0.8 as an example, Fig. <ref> shows the temporal evolution of powder, liquid phase, and the phase fraction of γ-Fe-Ni and FeNi_3 in the case of single-track scan. At t=0.6 ms (Fig. <ref>(a)), the laser beam is on and the temperature gradient is very large. In front of liquid phase (black color), powder is still there and no phase forms. Due to the high temperature around the liquid phase, the already deposited layer behind the liquid phase and the substrate region close to the liquid phase possess only γ-Fe-Ni phase.
Because of the neglectable difference between the liquidus and solidus temperatures, it can be seen from Fig. <ref>(a) that the region with liquid-solid coexistence is almost unobservable and the interface between liquid and solid is sharp. This result is different from the previous work on stainless steel 316L with a solidus-liquidus temperature difference around 120 K calculated from CALPHAD <cit.>. Meanwhile, the region with coexistence of γ-Fe-Ni and FeNi_3 phases is also very narrow and only exists in substrate. When the time goes to 0.9 ms, the laser beam is off. As presented in Fig. <ref>(b) at 3.85 ms, the temperature is lower than 793 K, no liquid phase remains, and the deposited layer completely covers x direction. In the deposited layer, there are wide regions where γ-Fe-Ni and FeNi_3 phases coexist. More specifically, Fig. <ref>(c) and (d) depicts the phase fraction of γ-Fe-Ni (Ψ_γ-Fe-Ni) along the surface node path marked in Fig. <ref>(b) at various times. It can be found from Fig. <ref>(c) that within 1.3 ms, γ-Fe-Ni (Ψ_γ-Fe-Ni=1), liquid, and powder exist along the node path, but the interface between them is extremely narrow. During the cooling down process in Fig. <ref>(d), Ψ_γ-Fe-Ni is found to be between 0 and 1 in a wide region along the node path, indicating the obvious coexistence of γ-Fe-Ni and FeNi_3 phases in the deposited layer.
In order to clearly show the mixture of FeNi_3 and γ-Fe-Ni phases, the time in Fig. <ref>(d) is selected so that the corresponding temperature falls into the region of phase coexistence.
Fig. <ref> also indicate that after cooling down to 300 K, the final state is the ordered intermetallic compound FeNi_3 phase. For the composition Fe_20Ni_80, the additional Ni element is present with the formation of solid solution in FeNi_3 phase. The conclusion with final FeNi_3 phase also agrees with the experimental X-ray diffraction results from the literature <cit.>, in which it is demonstrated that FeNi_3 phase is readily identified and the impurity like other Fe-Ni intermetallic compounds or Fe-Ni simple substance cannot be found in samples.





Fig. <ref> shows the phase evolution during the in-plane multi-track scan. At t=6.5 ms in Fig. <ref>(a), the laser beam is around the center of the second track and the powder state remains in front of the melt pool. A narrow interface between γ-Fe-Ni and FeNi_3 is found in both the substrate and the first-track deposit.  During the cooling down of the whole sample at t=23.3 ms in Fig. <ref>(b), γ-Fe-Ni/FeNi_3 interface is significantly broadened and a large region in the deposited layer contains two phases. It is obvious that in the contour plots in Fig. <ref>(a) and (b), the subsequent laser scan induces phase changes in the previously deposited layer. More specifically, Fig. <ref>(c) depicts the temporal evolution of phase fraction at the edge node marked in Fig. <ref>(a). It can be seen that this material node changes from powder to liquid during the first-track scan, and then experiences cyclic phase changes between γ-Fe-Ni and FeNi_3 during the subsequent scans. The cyclic phase change is originated from the cyclic temperature history as discussed in Fig. <ref>(b).

Similarly, Fig. <ref> presents the phase evolution during the out-of-plane multi-track scan. As shown in Fig. <ref>(a), during the deposition of the first layer at t=0.6 ms, liquid, powder, γ-Fe-Ni and FeNi_3 all exist and the interfaces among them are very sharp. At the beginning of depositing the third layer (t=17.8 ms), the previously deposited two layers are γ-Fe-Ni and a wide region with the coexistence of γ-Fe-Ni and FeNi_3 appears in the substrate. During the cooling process (t=55 ms) after the completion of seven-layer deposition, a layered region possessing both γ-Fe-Ni and FeNi_3 emerges in the built vertical structure. The temporal evolution of phase fraction at the edge node (as marked in Fig. <ref>(a)) in Fig. <ref>(b) shows that the material node experiences cyclic phase changes.

It should be noted that in the above analysis, the liquid-solid coexistence is ignored and the solidification be assumed to occur instantly due to the neglectable difference between liquidus and solidus temperatures from the equilibrium phase diagram. However, nonequilibrium solidification process could happen in the real case. In order to deal with the nonequilibrium solidification, we utilize the Scheil-Gulliver model <cit.> which has been implemented in the Thermo-Calc package. The model assumes perfect mixing in the liquid and no diffusion in the solid phase. As introduced by Scheil in 1942 <cit.>, the partition coefficient (κ) during solidification can be defined as the ratio of the local composition of the solid phase C_S to that of the liquid phase C_L, C_S/C_L as determined from the phase diagram. If C_0 is the starting composition and f_S is the fraction of solid, C_S can be obtained through the following equation

    (C_L-C_S)df_S=(1-f_S)dC_L.

The solution of Eq. (<ref>) is C_S=κ C_0 (1-f_S )^κ-1. The calculated fraction of solid (f_S) in both equilibrium and nonequilibrium solidification for FeNi powder with composition x_Ni=0.8 is shown in Fig. <ref>. It can be found from Fig. <ref>(a) that within a temperature range of 300–2000 K, the difference between the equilibrium result and the nonequilibrium result is hardly observable. Minor difference only appears around the liquidus and solidus temperatures. As shown in the inset of Fig. <ref>(a), the liquid-solid coexistence region is about 0.2 K in the equilibrium result, and is about 0.9 K in the nonequilibrium result. The nonequilibrium result gives more accurate prediction on the fraction of solid within the temperature range of 0.9 K (from 1709.36 to 1708.46 K in Fig. <ref>(a)). However, the deviation of the above equilibrium results from the nonequilibrium result only exists in this small temperature range of 0.9 K. In addition, Fig. <ref>(b) indicates that during the nonequilibrium solidification, microsegregation including Ni depletion and Fe enrichment in γ-Fe-Ni occurs, but it is extremely weak.



Solid phase transformation during the SLM processing induced heating and cooling is also important, i.e. the transformation between γ-Fe-Ni phase and FeNi_3 phase for the FeNi alloy with a composition x_Ni around 0.8. It should be mentioned that the results on the transformation between γ-Fe-Ni and FeNi_3 at low temperatures in Fig. <ref>(d), Fig. <ref>(c), and Fig. <ref>(b) are obtained by using the phase diagram and lever principle without the consideration of nonequilibrium effect. This is an extremely rough approximation. More accurate or quantitative modeling and computation of heating and cooling induced phase transformations between γ-Fe-Ni and FeNi_3 phases with nonequilibrium effect are highly recommended, but are still challenging. On the one hand, in contrast to the disordered γ-Fe-Ni solid solution, FeNi_3 is a ordered phase whose modeling is difficult due to the complexity of involved sublattices and the lack of mobility parameters. On the other hand, phase-field modeling provides the feasibility of simulating microstructure evolution during solid phase transformation <cit.>, but it requires the comprehensive thermodynamic and kinetic data of both disordered γ-Fe-Ni and ordered FeNi_3 phases which are unfortunately not readily available at the current stage. Future efforts have to be made towards modeling microstructure evolution and nonequilibrium effect in heating and cooling induced phase transformations between γ-Fe-Ni and FeNi_3.


§ MAGNETIC PROPERTIES

As a functional magnetic material, Fe-Ni permalloy is featured by its magnetic properties, in contrast to the structure materials  with the merit of strength and ductility. Therefore, the focus of property prediction here is different from that in SLM processed structural materials.
Generally speaking, the phase transformation during laser scanning process will influence the intrinsic magnetic properties of the formed phases, such as Curie temperature T_c, saturation magnetization M_s, magnetocrystalline anisotropy K_a, and magnetostriction λ_s. Since magnetic properties are concerned only when the temperature is below T_c (above which a paramagnetic state exits), the Curie temperature curves are important and thus are plotted in the phase diagram, as shown in Fig. <ref>. It can be found that FeNi_3 possesses higher T_c than γ-Fe-Ni. It means that if γ-Fe-Ni is transformed into FeNi_3 in the phase coexistence region, T_c of the system will be increased. But the quantitative calculation of T_c in this phase coexistence region is still challenging. Taking the composition x_Ni=0.85 as an example, the formed phases will be magnetic only below 815 K. In addition, M_s and λ_s of both γ-Fe-Ni and FeNi_3 with the same composition x_Ni are similar <cit.>. As typical soft magnetic materials, both γ-Fe-Ni and FeNi_3 have extremely small magnetocrystalline anisotropy and thus negligible K_a values. Therefore, under the fixed composition x_Ni around 0.8, the phase transformation below the Curie temperature curves in Fig. <ref> will not significantly influence the intrinsic magnetic properties M_s, K_a, and λ_s.
In contrast, magnetic coercivity is an extrinsic magnetic property, which is a very important indicator for the assessment of magnetic materials and is the manifestation of these intrinsic magnetic properties combined with extrinsic factors such as shape, stress/strain, etc. To this end, we attempt to calculate the magnetic properties in the SLM processed Fe-Ni alloy by using micromagnetic simulations, with the consideration of SLM induced residual stress.

In the micromagnetic framework <cit.>, the spatial distribution of the magnetization is described as 𝐦M_s in which 𝐦 is the unit vector of the magnetization direction and M_s is the saturation magnetization. Micromagnetic model is in the framework of continuum theory that handles magnetization processes on a length scale that is small enough to resolve the transition of the magnetization within domain walls but large enough to replace the atomic magnetic moments by a continuous function of position. The free energy density of a magnetic body is defined as

    E^tot = E^exch + E^ani + E^demag + E^Zeeman + E^m-ela.

In Eq. (<ref>), the exchange energy

    E^exch=A_e∇𝐦^2

is related to the gradient of 𝐦 and contributes to the domain wall energy, with A_e denotes the exchange parameter which is taken as 13 pJ/m for Fe-Ni Permalloy <cit.>. The magnetocrystalline anisotropy energy

    E^ani = {[ K_a(m_1^2m_2^2+m_2^2m_3^2+m_3^2m_1^2),                       cubic anisotropy;                      K_a [1-(𝐦·𝐮)^2 ],                    uniaxial anisotropy ].

physically means that a magnetic material is said to have magnetocrystalline anisotropy if it takes more energy to magnetize it in certain directions than in others. K_a is the magnetocrystalline anisotropy constant and 𝐮 is the unit vector along the anisotropy direction. For Fe-Ni Permalloy, K_a is usually taken as zero. The demagnetization energy

    E^demag=-μ_0/2M_s𝐦·𝐇_d

in which μ_0 is the vacuum permeability constant and 𝐇_d is the demagnetization field generated by the magnetic body itself. The Zeeman energy

    E^Zeeman=-μ_0M_s𝐦·𝐇_ex

is the energy of magnetization in an external magnetic field 𝐇_ex.

All these above energy terms are only related to the magnetization. The most important term considered here is E^m-ela, i.e. the magnetoelastic energy originated from the coupling between magnetization and stress/strain, which for a polycrystalline with isotropic magnetostriction can be given as <cit.>

    E^m-ela=-3/2λ_s(∑_i=1^3σ_iim_i^2 + ∑_i≠ j^3σ_ijm_i m_j ),

where λ_s is the magnetostriction coefficient. It can be seen from Eq. (<ref>) that the residual stress σ_ij induced by the SLM process contributes to the total free energy and thus will affect the magnetic hysteresis and coercivity.

!t
Conjugate gradient method for energy minimization and hysteresis calculation in micromagnetics

 




The minimization of total energy E^tot with respect to 𝐦 can be realized by the conjugate gradient method. The conjugate gradient method for the hysteresis calculation is given in Algorithm <ref>, which is implemented in OOMMF code <cit.> in the framework of finite difference method (FDM). The details of the numerical implementation are concisely described here. 
Through the discretization by FDM, the energy in Eq. (<ref>) at an external magnetic field 𝐇_ex can be rewritten as

    E^tot(𝐦̂,𝐇̂_ex)=1/2𝐦̂^T𝒞𝐦̂-1/2𝐇̂_d^Tℳ𝐦̂-𝐇̂_ex^Tℳ𝐦̂.

In Eq. (<ref>), the magnetization direction vectors 𝐦^I at each cell I of a FDM mesh are gathered into a vector 𝐦̂∈ℝ^3K (K is the total number of FDM cells), i.e. 

    𝐦̂=[ m_x^1,m_y^1,m_z^1,…,m_x^K,m_y^K,m_z^K ]^T.

In the similar way, the external magnetic field and the demagnetizing field at the nodes or cells are gathered into 𝐇̂_ex and 𝐇̂_d, respectively. 
Since the magnetization magnitude does not change and only the magnetization direction varies (i.e. 𝐦^I=1), a renormalized vector 𝐦̂ has to be used during the numerical calculations and is correspondingly defined as

    𝐦̂=[m_x^1/𝐦^1,m_y^1/𝐦^1,m_z^1/𝐦^1,…,m_x^K/𝐦^K,m_y^K/𝐦^K,m_z^K/𝐦^K]^T.

The sparse matrix 𝒞 contains grid information associated with the exchange, anisotropy, and magnetoelastic energies. The matrix ℳ accounts for the local variation of the saturation magnetization M_s within the magnet. The effective magnetic field at each cell I is calculated as

    𝐇^I=-1/μ_0M_s^I∇_𝐦^IE^tot(𝐦̂,𝐇̂_ex).

It should be mentioned that in contrast to the normal conjugate gradient method which directly includes the energy gradient (e.g. 𝐇^I in Eq. (<ref>)) in the calculation of search direction, here the cross product ℍ^I=𝐦^I×𝐇^I×𝐦^I is used, as shown in Algorithm <ref>. By using the constraint 𝐦^I=1, ℍ^I is simplified as 𝐇^I-(𝐇^I·𝐦^I)𝐦^I, i.e. it only represents the effective magnetic field component perpendicular to 𝐦^I. The application of this cross product is physically reasonable, since the field parallel to 𝐦^I cannot induce the magnetization change according to the Landau–Lifshitz–Gilbert equation <cit.>. Following the similar notation for 𝐦̂, this cross product vectors ℍ^I at each cell I of a FDM mesh are gathered into a vector ℍ∈ℝ^3K, i.e.

    ℍ = [(𝐦^1×𝐇^1×𝐦^1)^T,…,(𝐦^K×𝐇^K×𝐦^K)^T]^T.

In order to compute the value of β in Algorithm <ref> for the update of the search direction in the conjugate gradient method, both the Fletcher–Reeves method <cit.> and Polak–Ribiere method <cit.> are numerically implemented. It should be mentioned that for the Polak–Ribiere method, the vector ℍ_j is needed. While for the Fletcher–Reeves method, only the scalar values ℍ_j^Tℍ_j needs to be saved between iterations, which reduces the memory requirement and thus is chosen for the micromagnetic energy minimization in this work. Repeating the minimization under different external magnetic fields 𝐇̂_ex^(i) to give the corresponding 𝐦̂^(i) (Algorithm <ref>) will result in the magnetic hysteresis from which the magnetic properties can be calculated.

The magnetic parameters of Fe-Ni alloy with three different compositions used for micromagnetic simulations are listed in Table <ref> <cit.>. The nominal magnetostriction of Fe_20Ni_80 is zero. In order to investigate the possible effect of tiny deviation from the nominally zero magnetostriction, small λ_s (± 1.0× 10^-7) is also considered for the case of Fe_20Ni_80. Fe_15Ni_85 and Fe_25Ni_75 possess negative and positive magnetostriction, respectively. The micromagnetic model is meshed by cubic cells with a size of  10 μm × 10 μm × 5 μm.  Since the mesh for stress calculation and micromagnetic simulation is different, the stress distribution in the micromagnetic model is obtained by performing interpolation of the nodal stress from the previous mechanical analysis through a Delaunay triangulation of the scattered data <cit.>.






The magnetic hysteresis is calculated by applying external magnetic field along the z direction. Fig. <ref>(a)-(d) shows typical hysteresis curves for the case of Fe_20Ni_80 and Fe_15Ni_85 with and without SLM processing. Zooming in around the original point leads to the coercivity observable. When there is no SLM processing, both Fe_20Ni_80 and Fe_15Ni_85 possess very small coercivity around 0.5 mT, as depicted in Fig. <ref>(a) and (c). This is expected for the magnetically soft Fe-Ni alloys. For Fe_20Ni_80, we check the influence of small magnetostriction coefficient on its coercivity, as shown in Fig. <ref>(b) and Fig. <ref>(a). It can be seen that if λ_s is set as small values ± 1.0× 10^-7, SLM processing with different laser power can only lead to a coercivity change within 0.05 mT. On the contrary, for Fe_15Ni_85 which has a negative λ_s of -5.8× 10^-6, SLM processing with 100 W laser can increase the coercivity to ∼ 1.5 mT, as presented in Fig. <ref>(d).
Fig. <ref>(b) shows the calculated coercivity for different alloy compositions and laser power. It can be found that the coercivity of SLM processed Fe_15Ni_85 and Fe_25Ni_75 is increased to ∼1.7 and ∼0.8 mT, respectively. Within the laser power P=100-130 W, the coercivity is found to only slightly increase with P. The reason may be that the stress magnitude and distribution are not significantly changed within P=100-120 W. It should be mentioned that the experimentally measured coercivity of SLM processed Fe-Ni alloys is around 0.06-4 mT <cit.>. Our simulation results on the coercivity are in accordance with these experimental measurement.
Apart from the coercivity values, the remanent magnetization (μ_0 M_r) and magnetic domain structure are also affected by the SLM process.
Fig. <ref>(e) and (f) shows the magnetic configuration at the remanent state of the sample in (c) and (d), respectively.
For Fe_20Ni_80 (Fig. <ref>(a) and (b)) and Fe_15Ni_85 without SLM processing (Fig. <ref>(c)), μ_0 M_r is as low as ∼1 mT. Accordingly, the remanent magnetic configuration is composed of large-area magnetic domains perpendicular to z direction or along the negative z direction, as shown in Fig. <ref>(e).
However, after Fe_15Ni_85 is processed by SLM, μ_0 M_r is enhanced to ∼5 mT (Fig. <ref>(d)), and magnetic domains along the positive z direction occupy larger areas (Fig. Fig. <ref>(f)). These results computationally confirm that SLM processing could affect the coercivity, remanent magnetization, and magnetic domain structure in Fe-Ni alloys.




§ SUMMARY AND OUTLOOK

In conclusion, we have integrated FEA, CALPHAD output, and micromagnetics to demonstrate the first attempt for the computational evaluation of the microstructure evolution and magnetic coercivity of SLM processed magnetic Fe-Ni alloys.
The flowchart and coupling schemes of the integrated framework are summarized in Fig. <ref>, including heat-transfer model, mechanical model, CALPHAD, and micromagnetic model. Both the heat-transfer and mechanical models are numerically solved by finite element method. The micromagnetic model and thus the magnetic hysteresis are calculated through the energy minimization by the conjugate gradient method within the finite difference framework. These models are coupled through the information transfer among them. For example, the temperature history T(x,t) from the heat transfer model can be used as input for the CALPHAD and mechanical model. In detail, finite element (FE) nodal values of field variable ϕ(x,t) and temperature T(x,t) from the finite element simulation of heat transfer model are mapped to the CALPHAD model by a Python script to obtain the temporal and spatial phase fraction. Meanwhile, FE nodal values of ϕ(x,t) and T(x,t) are also imported to the mechanical model for the calculation of transient stress σ(x,t) and strain ϵ(x,t) by finite element method. In turn, CALPHAD can provide thermodynamically consistent parameters for the Gauss quadrature points in finite element simulations, such as specific heat C_p(T), density ρ(T), and latent heat L for the heat transfer model and thermal expansion coefficient α(T) for the mechanical model. Furthermore, the stress fields from the mechanical model can be input to the micromagnetic model for the calculation of magnetic properties, including coercivity, remanent magnetization, and magnetic domain structure. For transferring the stress fields from the finite element method to the micromagnetic model which is solved by finite difference method, we use a Delaunay triangulation of the scattered FE nodal values to perform the linear interpolation. Overall, by using the temperature-dependent material states/parameters and the step-by-step element activation of powder mesh, finite element simulation of the heat-transfer model is performed to calculate the temperature distribution and evolution during the SLM processing of Fe-Ni alloys.
With the thermal history from finite element simulations as input, thermomechanical analysis is performed by using the elastic-plastic material constitutive model.  By integrating the thermal history and CALPHAD output, the evolution and distribution of liquid, powder, FeNi_3 and γ-Fe-Ni during the SLM process are calculated. 
Finally, micromagnetic simulations which treats the residual stress as the magnetoelastic energy are carried out to calculate the magnetic property of SLM processed Fe-Ni alloy.
By using this computational framework, the melting pool geometry and the cyclic thermal history are identified. The cyclic tension and compression are confirmed in the interface of two neighboring tracks, which could degrade the interface bonding. It is found that the material firstly changes from powder to liquid and then experiences cyclic phase changes between γ-Fe-Ni and FeNi_3. SLM process is found to obviously enhance the remanent magnetization and increase the coercivity of Fe_15Ni_85 and Fe_25Ni_75 to 0.8-1.7 mT. The calculated coercivity is shown to agree with the experimental values.

While we have shown a promising method for the simulation of additively manufactured magnetic materials by integrating FEA, CALPHAD, and micromagnetics, the work here is a first attempt and lots of issues have to be thoroughly considered in the near further. As an initial work, here we limit ourselves to the working principle and potential feasibility of the proposed computational scheme, and do not focus on the accurate prediction of industrial or real AM process at the current stage. Several issues in this work related to the real AM process have to be deliberated and resolved in the next step, as listed in the following:

(1) The AM processability of FeNi alloys by SLM is an open question. Here we hypothesize defect-free processing in our simulation. The low processability related effects from fluid dynamics, surface tension and gas flow have to be considered, with more experimental information on the processability.

(2) The temperature-dependent material parameters in Fig. <ref>, which are important for the thermal and stress analysis, are not accurate, due to the lack of sufficient experimental data. It will be good to make more efforts on experimental measurements and obtain thermodynamic property information of magnetic alloys form CALPHAD database. Including the possible vaporization may be also necessary for the accurate prediction of the peak temperature. The phase change and composition variation induced stress should also play a role.

(3) Apart from the weak microsegregation of elements by using the Scheil-Gulliver model, the phase fraction is taken as the main indicator for the microstructure evolution in this work. However, in the real experimental case, microstructure would be more complicated. Other phases except for γ-Fe-Ni and FeNi_3 may also exist and spatial element microsegregation and composition inhomogeneity can occur. There are also other general issues in the prediction of  microstructure evolution during SLM process, such as the non-equilibrium state, fast cooling induced solute trapping, effect of free energy from magnetic contribution, etc. Phase-field simulation using CALPHAD information of magnetic materials <cit.> could be a viable methodology for predicting the microstructure evolution during the SLM processing, and deserves our future efforts.

(4) The dependence of magnetic properties on the phase fraction and microstructure is very complicated and not involved in this work.
As for the prediction of magnetic properties of SLM-processed magnetic materials in the real case, more microstructure information (e.g. phase distribution, grain boundaries, grain orientation, porosity, surface roughness, etc.) except for residual stress should be considered.

(5) Even though the calculated coercivity is shown to be in line with the experimental one for the laser additively manufactured magnetic FeNi, the simulation work here lacks the experimental validation on several points such melt pool size, temperature, stress, microstructure, etc. Collaborative experimental work has to be carried out in order to make the computational framework fully convinced.



Acknowledgements The support from the German Science Foundation (DFG YI 165/1-1 and DFG XU 121/7-1), the Profile Area From Material to Product Innovation – PMP (TU Darmstadt), the European Research Council (ERC) under the European Union´s Horizon 2020 research and innovation programme (grant agreement No 743116), and the LOEWE research cluster RESPONSE (Hessen, Germany) is acknowledged. The authors also greatly appreciate their access to the Lichtenberg High Performance Computer of Technische Universität Darmstadt.



cm-num



