# Fileset [Phase-field modeling of solid-state sintering with interfacial anisotropy.pdf](https://mdr.nims.go.jp/filesets/4985cb89-ae69-400b-a8ce-6cde385ca31b/download) ## Creator [Akimitsu Ishii](https://orcid.org/0000-0002-9261-4047), Kyoyu Kondo, Akiyasu Yamamoto, Akinori Yamanaka ## Rights [Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/) ## Other metadata [Phase-field modeling of solid-state sintering with interfacial anisotropy](https://mdr.nims.go.jp/datasets/c9c533e6-6dbe-4fd5-97c3-f6831804055a) ## Fulltext Phase-field modeling of solid-state sintering with interfacial anisotropyMaterials Today Communications 35 (2023) 106061Available online 24 April 20232352-4928/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).Phase-field modeling of solid-state sintering with interfacial anisotropy Akimitsu Ishii a,*, Kyoyu Kondo b, Akiyasu Yamamoto c, Akinori Yamanaka d a International Center for Young Scientists, National Institute for Materials Science, 1–2-1, Tsukuba, Ibaraki 305–0047, Japan b Department of Mechanical Systems Engineering, Graduate School of Engineering, Tokyo University of Agriculture and Technology, 2–24-16, Koganei, Tokyo 184–8588, Japan c Division of Advanced Applied Physics, Institute of Engineering, Tokyo University of Agriculture and Technology, 2–24-16, Koganei, Tokyo 184–8588, Japan d Division of Advanced Mechanical Systems Engineering, Institute of Engineering, Tokyo University of Agriculture and Technology, 2–24-16, Koganei, Tokyo 184–8588, Japan A R T I C L E I N F O Keywords: Phase-field modeling Solid-state sintering Anisotropy Crystal grain orientation A B S T R A C T Sintered structures observed in experiments can consist of faceted crystal grains. To predict the formation of such structures, a new phase-field (PF) model of solid-state sintering that can analyze morphological changes and microstructural evolutions considering the strong interface anisotropies of sintered particles was developed in this study. The developed PF model treats the crystal grain orientation-dependent surface and misorientation- dependent grain boundary anisotropies. Furthermore, this model employs quaternions to calculate the particle rotation, eliminating complicated calculations involved in analyzing the crystal orientation in three-dimensional (3D) simulations. The morphological change in a sintered particle and the grain boundary formulation at triple junctions simulated using the developed model were validated by comparing them with theoretical solutions. The neck growth and densification rates were investigated by performing 3D simulations using two particles with interface anisotropies. The simulation results revealed that neck growth and densification are affected by surface energy and mobility anisotropies. A 3D PF simulation using 200 particles demonstrated that the developed PF model can potentially reproduce sintered structures with faceted particles observed in experiments. The PF model provides a promising simulation for predicting the microstructural evolution of actual materials during solid-state sintering. 1. Introduction Sintering is a manufacturing process that converts powder particles into a solid compact [1]. Various products and materials, including ceramics, metals, hard magnets [2,3], and superconducting materials [4, 5], are produced by sintering. Additionally, sintering has recently attracted attention as a fundamental technology in additive manufacturing [6,7]. The mechanical, electrical, magnetic, and other properties of sintered materials depend on their microstructures. For example, the grain boundary misorientation affects the electromagnetic properties of polycrystalline bulk high-temperature superconducting materials [8], and the crystal orientation alignment contributes to the magnetic properties of anisotropic magnets [9]. Therefore, microstructural evolution during sintering must be controlled to improve the material properties. Numerous experiments have been conducted to understand the morphological and microstructural changes induced by sintering. In recent years, in situ observations of microstructural evolution during sintering have been actively conducted [10–13]. Bah et al. [10] in situ observed a shrunken particle due to grain boundary migration during the sintering of piezoelectric K0.5Na0.5NbO3 particles. Zuo et al. [11] conducted quasi-in situ observations of sintering with Cu nanoparticles and investigated the grain growth and twin formation mechanisms. Takahashi et al. [12] revealed the evolution of density distribution and growth of coarser pores during the sintering of Al2O3 by in situ observation. Phuah et al. [13] evaluated the impact of ultra-high heating rate on densification in sintering through in situ observation of various ceramic nanoparticles. Although these experiments using cutting-edge techniques provide useful information to understand the sintering mechanism, experimental observations alone are insufficient to accurately predict neck growth, pore annihilation, densification, and grain growth, which are essential factors in controlling the microstructural evolution. The microstructural evolution during solid-state sintering was * Corresponding author. E-mail address: ISHII.Akimitsu@nims.go.jp (A. Ishii). Contents lists available at ScienceDirect Materials Today Communications journal homepage: www.elsevier.com/locate/mtcomm https://doi.org/10.1016/j.mtcomm.2023.106061 Received 24 February 2023; Received in revised form 12 April 2023; Accepted 22 April 2023 mailto:ISHII.Akimitsu@nims.go.jpwww.sciencedirect.com/science/journal/23524928https://www.elsevier.com/locate/mtcommhttps://doi.org/10.1016/j.mtcomm.2023.106061https://doi.org/10.1016/j.mtcomm.2023.106061https://doi.org/10.1016/j.mtcomm.2023.106061http://crossmark.crossref.org/dialog/?doi=10.1016/j.mtcomm.2023.106061&domain=pdfhttp://creativecommons.org/licenses/by/4.0/Materials Today Communications 35 (2023) 1060612examined via numerical modeling and simulation using the phase-field (PF) method. The PF method [14] is a powerful numerical methodology for solving free-boundary problems. Because the PF method can analyze the complicated evolution of interfaces without explicitly tracking them, PF simulations of sintering have been used to study microstructural evolution during sintering. For example, Choudhuri et al. [15] used PF simulations to investigate the impact of particle local curvature on neck growth during sintering. To date, various PF sintering models have been developed. Wang [16] was the first to propose a PF model for solid-state sintering that can simulate the rigid-body motions of sintered particles. Rigid-body motion is an important factor to simulate the macroscopic densification of sintered compacts. Several PF models have been developed by incorporating the calculation method for the rigid-body motion proposed by Wang. Biswas et al. [17–19] proposed PF models of solid-state sintering considering elastic energy [17], crystal orientation [18], and thermal conduction [19] as well as rigid-body motion. Zhao et al. [20] developed a sintering PF model that could directly characterize the elastic energy on the neck and contact area of sintered compact to investigate the microstructural evolution in pressure-assisted sintering. Ivannikov et al. [21] reduced the gap between PF simulations and experiments by developing an advection mechanism to calculate the rigid-body motion considered in Wang’s model. Notably, Wang’s original PF model was limited in the number of particles that could be practically analyzed owing to its large computational cost. To overcome this problem, Termuhlen et al. [22] developed a new numerical scheme using a grouping algorithm that enabled three-dimensional (3D) simulations with many particles. Furthermore, Ma et al. [23] helped reduce the computational cost and improved computational efficiency of the Wang model by applying a proper generalized decomposition technique to the 3D PF simulation of solid-state sintering. PF models of sintering have also been developed using a different approach. For example, Hötzer et al. [24] developed a PF model of solid-state sintering based on the concept of grand-potential model [25] and investigated microstructural evolution and densification by performing a 3D simulation of sintering with approximately 25,000 Al2O3 particles. Greenquest et al. [26,27] proposed a sintering PF model based on the grand-potential model and predicted the sintering behavior of UO2 doped with Cr or Mn. Several PF simulations focusing on the sintering process in additive manufacturing have also been reported [28–32]. Abdeljawad et al. [28] used PF simulations to examine the microstructural evolution during solid-state sintering in a direct ink write process. Yang et al. [29,30] constructed a PF model for non-isothermal sintering and performed simulations of selective laser sintering and powder bed fusion. Wang et al. [31] performed a non-isothermal sintering simulation for direct metal laser sintering and clarified the effect of powder size distribution on densification behavior. As is evident, various PF models of sintering have been developed; however, some factors remain to be considered in these models. One such factor is the strong anisotropy of surface energy depending on the crystal grain orientation. Interfacial energy is an important factor for predicting microstructural evolution. In particular, surface energy (that is, surface tension) is related to neck growth, which affects the densification of sintered compacts [1,16]. Most materials have anisotropic surface energies depending on their crystal grain orientation. The anisotropic surface energy changes the shape of isolated crystal grains to an equilibrium shape [33]. Faceted particles are observed in actual sintered material, which can be attributed to the strong anisotropy of the surface energy. As an example, Fig. 1 shows a scanning electron microscopy image of faceted particles of a K-doped BaFe2As2 polycrystalline bulk superconducting material sintered at 900 ◦C for 18 h under ambient pressure. Faceted particles have also been observed in other sintered materials [34–36]. Therefore, PF sintering models must consider the crystal grain orientation of particles and strong anisotropy of the surface energy to accurately predict the microstructural evolution. As mentioned earlier, Biswas et al. [18] developed a PF model for solid-state sintering that can analyze the crystal grain orientation of sintered particles. They performed two-dimensional PF simulations considering the direction-dependent interface diffusion and misorientation-dependent grain-boundary energy anisotropies. However, no PF model of solid-state sintering that can simulate morphological changes in particles owing to the surface energy anisotropy has been reported. In other words, the PF models developed to date cannot predict the formation of the sintered structure shown in Fig. 1. Additionally, to reproduce the structure of sintered compacts, similar to that shown in Figs. 1, 3D simulations must be performed. However, no 3D PF simulation considering crystal grain orientation-dependent anisotropy has been performed. One possible reason for this is that 3D simulations require complicated calculations to analyze changes in the crystal orientation due to rotation when the Euler angles are used for the rotation calculation, as in the conventional model [18]. The PF method enables the simulation of morphological changes in crystal grains to their equilibrium shape that is predicted from the anisotropy of the surface energy. Qin et al. [37] reported that the equilibrium shape of cubic crystals can be reproduced using PF simulations of crystal growth based on the anisotropic surface energy of the cubic crystal. Li et al. [38] experimentally demonstrated that one spherical W carbide particle transforms into a triangular prism shape by performing PF simulations using the strong anisotropic surface energy of the material. The integration of the aforementioned techniques is expected to lead to the development of an evolved PF model of sintering capable of simultaneously considering the rigid-body motion and anisotropies of the surface and grain boundary. In this study, a new PF model for solid-state sintering was developed to predict 3D sintered structures with faceted particles as shown in Fig. 1, which cannot be predicted by conventional PF models of sintering. The developed PF model not only considers the rigid-body motions of sintered particles and the misorientation-dependent grain boundary anisotropy based on the conventional PF models [16,18] but also newly considers the crystal grain orientation-dependent strong surface anisotropy. Moreover, to facilitate 3D PF simulations, a new technique to implement quaternions in the rotation calculation of crystal grain orientation due to the rigid-body motion is proposed. Similar to the previous PF models [16,18], the proposed PF model can simulate the isothermal sintering processes only for single-phase materials. However, our model has significance in that it provides a way to predict 3D sintered structures of actual materials that present interface anisotropy in most cases. In addition, the method of implementing interface anisotropy with quaternions in the PF model of solid-state sintering proposed in this study has the potential to be applied to other PF models of sintering. The remainder of this paper is organized as follows. Section 2 explains the PF model developed in this study. Section 3 describes the conditions for PF simulations conducted in this study. Sections 3.2 and Fig. 1. Scanning electron microscopy image of K-doped BaFe2As2 crystal grains obtained by sintering at 900 ◦C for 18 h under ambient pressure. A. Ishii et al. Materials Today Communications 35 (2023) 10606133.3 explain the definition of the anisotropic properties of the surface and grain boundary. Section 4.1 presents the results of the PF simulation to validate the developed PF model. Section 4.2 investigates the neck growth and densification obtained in the PF simulations, where two particles with anisotropic properties are in contact with each other. Section 4.3 presents the results of the PF simulation using 200 particles. The results demonstrate that the developed PF model is promising for predicting the actual sintered structure. Finally, Section 5 presents the conclusions of this study. 2. Formulation 2.1. Phase-field model The time evolution of microstructural evolutions during sintering is simulated by defining two types of PF variables, namely, the density field, ρ, and orientation field, ηi (i = 1, 2, 3,. N, where N is the number of particles). ρ represents the normalized local mass density. To satisfy mass conservation, ρ is defined as a conserved variable. ηi identifies the individual particle and is treated as a non-conserved variable because the volume of each particle can change via grain growth. The values of ρ and ηi smoothly change from 0 to 1 in the surface or grain boundary regions. The total free energy, F, is defined as F =∫v{fbulk(ρ, ηi) +κ2ρ2|∇ρ|2 +∑iκ2η2|∇ηi|2+λ22(∇2ρ)2}dV, (1) where fbulk is the chemical free-energy density in the bulk. The second and third terms on the right-hand side in Eq. (1) are the gradient energy densities. The fourth term indicates the corner energy density, which is incorporated to avoid singularity formation by adding high-order derivative regularization [39,40]. λ is the corner energy regularization parameter. The gradient energy coefficients κρ and κη are defined as [18, 28]. κρ =̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅34(2γ̂ surf − γ̂gb)δ√and (2) κη =̅̅̅̅̅̅̅̅̅̅̅34γ̂ gbδ√, (3) where δ is the width of the diffuse interfacial region of the particle surface and grain boundary. ̂γsurf is the anisotropic surface energy, which is defined as γ̂ surf =∑iηiγsurf,i(θi,φi)∑iηi, (4) where γsurf,i(θi,φi) is the surface energy, which is dependent on the crystal grain orientation of the ith particle. θi and φi are the angles in the local coordinate system for each particle. As shown in Fig. 2, these angles determine the direction of the normal vector at an arbitrary point on the surface of the ith particle. γ̂gb is the anisotropic grain boundary energy, and is defined as [18,41]. γ̂gb =∑i∑j>1γijη2i η2j∑i∑j>1η2i η2j, (5) where γij is the grain boundary energy between the ith and jth particles. Note that the PF model has a limitation on interfacial energy, i.e., 2γ̂surf > γ̂gb. fbulk can be defined as [16]. fbulk = Aρ2(1 − ρ)2+ B[ρ2 + 6(1 − ρ)∑iη2i − 4(2 − ρ)∑iη3i+ 3(∑iη2i)2 ], (6) where A and B are parameters given as [18,28]. A =12γ̂surf − 7γ̂gbδand (7) B =γ̂gbδ. (8) The time evolution of ρ is calculated using the following equation, which is the Cahn-Hilliard equation [42] with an advection term added to analyze the rigid-body motions of the particles: ∂ρ∂t= ∇⋅(Mρ∇δFδρ − ρ∑ivi), (9) where vi denotes the advection velocity of the ith particle. Mρ is the diffusion mobility, which as defined as Mρ = Mvolh(ρ) + Mvap{1 − h(ρ) } + M̂ surfρ2(1 − ρ)2+∑i∑j∕=iMijηiηj, (10) where Mvol and Mvap are the mobilities related to the volume and vapor diffusion of the atoms, respectively. M̂surf is the anisotropic surface diffusion mobility, which depends on the crystal grain orientation. Mij is the diffusion mobility at the grain boundary between the ith and jth particles. The details of M̂surf and Mij used in this study are explained in Section 3. In this study, the anisotropy of volume diffusion is neglected because its contribution to microstructural evolution can be assumed to be significantly smaller than that of the surface and grain boundary diffusion. h(ρ) is the interpolation function, given by h(ρ) = ρ3(10 −15ρ + 6ρ2). The time evolution equation of ηi is derived by adding the advection term to the Allen-Cahn equation [43]: ∂ηi∂t= − M̂ηδFδηi− ∇⋅(ηivi), (11) where M̂η is the anisotropic mobility related to the migration of grain boundary. During sintering, vacancies are induced by the surface tension and negative curvature on surfaces where particles are in contact. The concentration of generated vacancies increases with the surface tension and curvature [1]. Because the atoms diffuse such that the concentration gradient of vacancies decreases, a neck grows between the particles in contact. In contrast, vacancies flow from the neck surface into grain boundaries, resulting in vacancy oversaturation at grain boundaries. The oversaturated vacancies, which are annihilated by dislocations, lead to Fig. 2. Definition of θi and φi angles. A. Ishii et al. Materials Today Communications 35 (2023) 1060614rigid-body motions, which drive the macroscopic densification of a sintered compact. The PF model proposed by Wang [16] describes these sintering phenomena by calculating the advection velocity of each particle. The PF model developed in this study also uses advection velocity to simulate macroscopic densification. Rigid-body motion is assumed to consist of translation and rotation. Thus, vi in Eqs. (9) and (11) is given by vi = vtri + vroi , (12) where vtri and vroi are the translational and rotational velocities of the ith particle, respectively. These are defined as [16]. vtri =mtrViFiηi and (13) vroi =mroViTi ×[r − rc,i]ηi, (14) where mtr and mro are the translation and rotation mobilities, respectively. Fi and Ti are the force and torque acting on the ith particle, respectively. Vi is the ith particle volume. r is the position vector and rc,i is the center of the ith particle. Fi is calculated as Fi =∫Vkst∑j∕=i(ρ − ρ0)〈ηiηj〉[∇ηi − ∇ηj]d3r, (15) where kst is the stiffness constant relating the magnitude of the force to the degree of vacancy oversaturation at the grain boundary. ρ0 is the equilibrium value of ρ at the grain boundary. The operation 〈ηiηj〉 is defined as 〈ηiηj〉={1(ηiηj ≥ c)0(ηiηj < c) , (16) where c is the threshold for identifying grain boundary regions based on the product of ηi and ηj. Ti is defined as Ti =∫V[r − rc,i]× dFi. (17) The time evolution equations of ρ and ηi given by Eqs. (9) and (11), respectively, are solved using the first-order Euler finite difference method for time integration and the fourth-order central finite difference method for spatial discretization. 2.2. Calculation methodology for orientation rotation To calculate the orientation rotation computationally, unit quaternions are used. The rotation calculation using unit quaternions is faster and simpler than that using Euler angles. Moreover, this method produces small computational errors without calculation defects such as a gimbal lock. The operation of rotation by angle ϕ around an arbitrary unit vector, v, is defined by the following equation using a unit quaternion, q: q = cosϕ2+ v sinϕ2, (18) where an arbitrary vector u expressed using a quaternion p with p = u is considered. The rotational operation represented by q can be applied to u using the following equation: p′ = qpq∗, (19) where the superscript “* ” denotes a conjugate quaternion, and p’ is a quaternion representing the vector obtained by rotating u. Rotation calculations for sintered particles using quaternions require the axis and angle of rotation to act on particles. Because mroTi/Vi in Eq. (14) represents the angular velocity vector, the rotation axis can be obtained from the direction of Ti and the rotation angle can be calculated from the product of mro|Ti|/Vi and the time increment of the numerical computation, Δt. Therefore, this study proposes to calculate the evolution of the crystal grain orientation based on the rotation of the rigid-body motion by defining q as follows: q = cosϕ2+Ti|Ti|sinϕ2(ϕ =mro|Ti|ViΔt). (20) Quaternions also enable the calculation of the minimum angle between two particles, Δϕij, which matches the crystal grain orientation of one particle with that of the other. By assuming that the crystal grain orientations of the two particles are represented by quaternions qi and qj, Δϕij can be calculated using the following equation [44]: Δϕij = min{2arccos(Okqi(Olqj)∗)0, 2arccos(Olqj(Okqi)∗)0}, (21) where Ok and Ol (k, l = 1, 2, …) are quaternions, called symmetry operators, that are introduced to treat equivalent crystal grain orientations in crystal structures (for example, a cubic crystal system has 24 symmetric orientations). The maximum values of k and l are equal to the number of independent symmetries. The subscript 0 indicates the scalar element of quaternions. After the misorientation angles between adjacent particles are calculated, the deviation angles from the specific misorientations can also be calculated. Therefore, this PF model can simulate microstructural evolution with coincidence site lattice (CSL) grain boundaries. The deviation angle of Δϕij from an arbitrary Σm CSL grain boundary can be calculated using the following equation [44]: Δϕij,∑m = min{2arccos(Okqmin(Olq∑m)∗)0,2arccos(Olq∑m(Okqmin)∗)0},(22) where qmin is the quaternion that provides Δϕij using Eq. (21), and qΣm is the unit quaternion, which represents the rotation axis and angle of the Σm grain boundary. 3. Conditions 3.1. Parameters To validate the developed PF model, hypothetical particles with a cubic crystal system are used for sintering simulations in this study. The parameters included in the PF model are nondimensionalized for the spacing of the finite-difference grid, Δl, and the height of the energy barrier, Δf, at the interface between the solid and vapor phases [16]. The corner energy regularization parameter is set to λ = 7. The interface width is δ = 6Δl. The mobility of volume diffusion is Mvol = 0.01 and that of vapor diffusion is Mvap = 0.001 [16]. The translational and rotational mobilities are set as mtr = 500 and mro = 1, respectively [16]. The stiffness constant is kst = 10. The equilibrium value of ρ at grain boundaries is ρ0 = 0.9816 [16]. The threshold value is c = 0.14 [16]. The values of other parameters corresponding to the surface and grain boundary properties are given in Sections 3.2 and 3.3, respectively. The nondimensionalized time increment for each computational step is Δt = 1 × 10− 4. 3.2. Anisotropy of surface energy and mobility The anisotropic surface energy of a particle with a cubic crystal structure is determined using the following equation [37]: γsurf,i(θi,φi) =k0 + k1(n2xn2y + n2yn2z + n2z n2x)+ k2(n2xn2yn2z)+ k3(n2xn2y + n2yn2z + n2z n2x)2,(23) where k0, k1, k2, and k3 are coefficients. nx, ny, and nz are defined as follows: A. Ishii et al. Materials Today Communications 35 (2023) 1060615nx = sin θi cos ϕi, (24) ny = sin θi sin ϕi, and (25) nz = cos θi. (26) In this study, the nondimensionalized tentative values of k0 = 1.0, k1 = 2.0, k2 = 0.4, and k3 = 0.1 are used to allow the formation of facets. Fig. 3 shows the Wulff plot calculated using Eq. (23) and the stated coefficients. We assume that the larger the interfacial energy, the larger the interfacial mobility. Therefore, the anisotropic diffusion mobility of the surface, M̂surf , in Eq. (10) can be defined as M̂ surf∝Msurf ⋅γ̂ surf , (27) where Msurf is the reference surface diffusion mobility, being Msurf = 4.0 [16]. 3.3. Anisotropy of grain boundary energy and mobility Because grain boundaries have five crystallographic degrees of freedom, the anisotropic grain boundary energy and mobility are determined using three types of misorientations and two types of inclinations [45]. However, for simplicity, it is assumed that the grain boundary energy and mobility are only dependent on the misorientation angle of Δϕij [44,46,47]. The grain boundary energy, γij, with an energy cusp at a Σm CSL grain boundary is defined using the Read-Shockley equation [48]: where γmingb is the minimum grain boundary energy and γgb is the grain boundary energy at high-angle grain boundaries. Δϕh is the minimum misorientation angle of the high-angle grain boundaries, which is set as Δϕh = 15◦ [44,46,47]. D1 and W1 represent the depth and width of the grain boundary energy cusp, respectively, at the Σm grain boundary. The mobility of grain boundary diffusion, Mij, can be defined using the Humpherey equation [49]: Mij = Mgb⋅{1 − H(Δϕij,Δϕij,∑m)}, (29) where Mgb is the grain boundary diffusion mobility of the high-angle grain boundaries. H(Δϕij, Δϕij,Σm) is defined as H(Δϕij,Δϕij,∑m)=⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩exp{− 5(ΔϕijΔϕh)4}(Δϕij ≤ Δϕh)D2exp{− 5(Δϕij,∑mW2)2} (Δϕij,∑m ≤ W2)0(Δϕij > Δϕh and Δϕij,∑m > W2),(30) where D2 and W2 are the depth and width of the mobility change, respectively, at the Σm grain boundary. The Σ5 grain boundary about the rotation axis of 〈001〉 is considered in this study. Fig. 4 shows the variations in γij and Mgb obtained using ΔϕΣ5 = 36.87◦ [50] and the assumed values of D1 = D2 = 0.5, W1 = 8◦, W2 = 10◦, γmingb = 0.4, γgb = 0.8, and Mgb = 0.4 [16]. These values are used in the subsequent simulations. The definition of M̂η in the following equation allows for a dependence on the misorientation to be introduced into the mobility related to grain boundary migration: M̂η = Mη −∑i∑j∕=iMη(Δϕij)ηiηj, (31) where Mη is the reference mobility set as 10 [16]. Furthermore, Mη(Δϕij) is the mobility dependent on Δϕij and is defined as Mη(Δϕij)= Mη⋅H(Δϕij,Δϕij,∑m). (32) Although the aforementioned assumed definitions of interfacial energies, mobilities, and parameters are only used to validate the developed PF model, arbitrary settings of interfacial anisotropy can be used for the model. Therefore, using the appropriate settings of the target material, the PF model can perform sintering simulations for any anisotropic material. Fig. 3. Wulff plot calculated using Eq. (23) and coefficients of k0 = 1.0, k1 = 2.0, k2 = 0.4, and k3 = 0.1. γij = γmingb +(γgb − γmingb)⋅⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ΔϕijΔϕh{1 − ln(ΔϕijΔϕh)}(Δϕij ≤ Δϕh)1 − D1[1 −Δϕij,∑mW1{1 − ln(Δϕij,∑mW1)}](Δϕij,∑m ≤ W1)1(Δϕij > Δϕh and Δϕij,∑m > W1), (28) A. Ishii et al. Materials Today Communications 35 (2023) 10606164. Results and discussions 4.1. Validation To validate the developed PF model for simulating the morphological change in sintered particles depending on the surface energy anisotropy, a 3D PF simulation is performed using a single spherical particle. Theoretically, the equilibrium shape of an isolated crystal can be determined using a Wullf plot [33]. When a sintered particle has strong anisotropic surface energy, as shown in Fig. 3, its equilibrium shape is a cube. Therefore, if the simulation results show that the spherical particle approaches a cubic shape, the validity of the PF model for surface anisotropy is confirmed. The computational domain size is 643. A dimensionless particle with a nondimensionalized radius of 20 is arranged at the center of the domain. The zero Neumann condition is applied to all boundaries of the computational domain. The calculations are performed for up to 107 computational steps. Fig. 5 shows the morphological changes in the sintered particles obtained from the PF simulation. The particle shape is indicated by the isosurface ρ = 0.5, which is colored based on the distribution of the surface energy. To decrease the total free energy of the system, the surface changes to a facet in the low-energy direction and to a convex shape in the high-energy direction. Consequently, the shape of the particle changes from spherical to cubic. Fig. 6 compares the Wullf plot, equilibrium shape predicted from the Wulff plot, and particle shape obtained after 107 steps of simulation on the ab-plane through the center of the particle. The PF simulation accurately reproduces the equilibrium shape facets. The corners show a difference between the equilibrium shape and simulation result, indicating that the simulation result has not yet reached the equilibrium state. Although the particle shape is expected to approach the equilibrium shape by continuing the calculation, the shape change rate is significantly slow after 107 steps under the present simulation conditions. Therefore, it is concluded from the results shown in Figs. 5 and 6 that the developed PF model can reasonably simulate morphological changes depending on the anisotropic surface Fig. 4. Variations in grain boundary energy, γij, and mobility, Mij, calculated using Eqs. (28)–(30), with Δϕh = 15◦, ΔϕΣ5 = 36.87◦, D1 = D2 = 0.5, W1 = 8◦, W2 = 10◦, γmingb = 0.4, γgb = 0.8, and Mgb = 0.4. Fig. 5. Snapshot of morphological changes in sintered particle due to surface energy anisotropy. Fig. 6. Comparison of the Wullf plot, equilibrium crystal shape, and simulated particle shape on the ab-plane through the center of the particle. A. Ishii et al. Materials Today Communications 35 (2023) 1060617energy. The PF model is further validated by performing 3D PF simulations under three conditions where the particles have various crystal grain orientations, and the dihedral angles at the triple junctions of grain boundaries are investigated. Fig. 7(a) shows the initial arrangement of the three particles used in simulations. At the beginning of the simulation, all dihedral angles are 120◦. Fig. 7(b)–(d) show the orientations of the particles and distributions of the grain boundary energy calculated using Eqs. (5) and (28). The particles shown in Fig. 7(b) have the same crystal grain orientation. In the initial state shown in Fig. 7(c), Particle 3 is rotated by 45◦ around the c-axis; hence, there are high-angle grain boundaries between Particles 1 and 3 and between Particles 2 and 3. The grain boundary between Particles 1 and 2, shown in Fig. 7(d), is a Σ5 CSL grain boundary because the misorientation angle between the particles is 36.87◦ around the c-axis. Fig. 8 shows the distributions of grain boundaries on the cross- section through the center of the particles obtained after 105 computation steps. The grain boundaries are indicated as regions in which the sum of the squares of ηi is 0.5. The dihedral angles obtained from the PF simulations using the initial conditions shown in Fig. 7(b), (c), and (d) are 120◦, 151◦, and 136◦, respectively. These dihedral angles are equal to the analytical results calculated using Young’s equation. These results demonstrate that the PF model provides an accurate grain boundary formation that satisfies the balance of grain boundary energies. Fig. 7. (a) Shape and arrangement of three particles, as indicated by the isosurface with ρ = 0.5. (b–d) Crystal grain orientations of particles and distributions of grain boundary energy. The grain boundary energy is displayed only in regions judged to be a grain boundary using Eq. (16) on the cross-section through the center of particles. Blue, white, and red arrows indicate the directions of the a-, b-, and c-axes, respectively. Fig. 8. Distributions of particles and grain boundaries on the cross-section through the center of particles obtained by the PF simulations with 106 computational steps. Grain boundary distributions are indicated by the sum of the square of ηi. Fig. 9. Morphological changes during sintering of two particles with (a) isotropic and anisotropic surface energies for (b) Case 1, (c) Case 2, and (d) Case 3. A. Ishii et al. Materials Today Communications 35 (2023) 10606184.2. Investigation of neck growth and densification Conventionally, neck growth and densification caused by sintering between two particles have been explained using the sintering theory [1]. However, this theory does not consider the strong surface anisotropy that causes facet formation. Because the strong anisotropy of surface energy causes a shape change in sintered particles (see Fig. 5), the anisotropy of the surface is expected to affect the sintering behavior. To investigate the effects of surface anisotropic energy and mobility on neck growth and densification, solid-state sintering PF simulations are performed using two spherical particles. The size of the 3D computational domain is 128 × 64 × 64. The zero Neumann condition is applied to all boundaries of the computational domain. The radius of particles is 20, and the distance between the initial particle centers is 40. Although there are innumerable combinations of crystal grain orientations of two adjacent particles, this study focuses on the case in which the initial crystallographic orientations are equal (that is, Δϕij = 0◦ and γij = 0.4), and PF simulations are performed for the three cases. Additionally, a simulation using particles with isotropic surface energies and mobilities is also performed for comparison with the aforementioned three cases. The isotropic surface energy is set to be constant at 1.2. The neck growth and densification rates are investigated by analyzing the variation in the cross-sectional area of the neck and the approaching distance of the particle centers from simulation results. Fig. 9 shows the initial states of the PF simulations and the simulation results obtained after 106 computational steps. When the surface energy and mobility are isotropic (Fig. 9(a)), neck growth and densification occur according to the conventional theory [1,16]. Fig. 9(b)–(d) show the three cases in which particles have anisotropic surface energy, as shown in Fig. 3. Hereafter, these three cases are denoted as Cases 1–3. In Case 1, the initial particles are in contact with each other on the surface with the lowest surface energy. Under such conditions, facets are formed on the neck. In contrast, the particles in Case 2 are in contact with each other on the highest-energy surface. In this case, the neck has a concave shape, because the neck surface is oriented with the large surface energy. In Case 3, the particles contact each other on the surfaces with high energy, which are edges of the cubic equilibrium shape. The results of the simulation for Case 3 show that the neck forms facets along the z-axis and is concave along the y-axis. Fig. 10 (a) shows the variation in the cross-sectional area of the neck as a function of the computational step. The area of the neck is evaluated on the cross-section of the center of the computational domain in the x- axis direction shown in Fig. 9. To clearly show the difference between the four results after 50 computational steps, the area surrounded by dotted lines in Fig. 10 (a) is enlarged in Fig. 10 (b). Fig. 10 (c) shows the variation in the distance between particle centers. The distance is calculated as |rc,1 − rc,2|. The hollow circles indicate the respective values at each computation step up to 10 steps, and every 20 steps thereafter. The dashed lines are the approximate line connecting the hollow circles. Black line indicates the result of the isotropic case, corresponding to the results in Fig. 9(a). Red, blue, and yellow lines indicate Cases 1, 2, and 3, corresponding to Fig. 9(b)–(d), respectively. The neck area in Case 1 is smaller than that in the isotropic case until shortly before the end of the calculation. Two factors contribute to this decline in neck growth. First, the facets are formed by the slow growth of the neck surfaces perpendicular to the y- and z-axis directions, where the surface energy is low. Second, the formed facets grow slowly because the mobility of the surface diffusion is set to be proportional to the surface energy in this study, as shown in Eq. (27). However, although the neck growth in the isotropic case slows down with sintering time, the neck continues to grow in Case 1 such that the surface of the sintered compact forms a facet, resulting in a larger neck area in Case 1 than that in the isotropic case at 106 steps. For Case 1, the reduction in the distance between the particle centers is also slowed down from the isotropic case because the negative curvature of the neck surface, which causes rigid body motion (explained in Section 2), rapidly approaches zero owing to the formation of facets. In Case 2, the diffusion mobility of the neck surface is large because of the high surface energy of the neck. However, the neck growth is slow because increasing the neck surface area prevents a decrease in total free energy. However, the decrease in the particle center distance is faster than that in the isotropic case because the concave neck has a high negative curvature. The neck in Case 3 has facets in the z-axis direction, as in Case 1, and is concave in the y-axis direction, as in Case 2. Based on the simulation results, the rate of neck growth is mostly consistent with that in Case 2, whereas the particle center distance decreases more slowly in Case 3 than in Cases 1 and 2. The aforementioned results demonstrate that even if the misorientation angle between adjacent particles is 0◦, the anisotropic surface energy and mobility change the neck growth and densification rates depending on the orientation at which the particles contact each other. We suggest that the change in the neck curvature owing to anisotropic surface properties is an important factor that affects the rates of neck growth and densification. During the actual sintering process, particles are sintered at various misorientation angles between adjacent particles. The sintering behavior depends on not only surface anisotropy but also Fig. 10. Transitions of (a and b) neck area and (b) distance between particle centers in sintering simulations. The area enclosed by the dotted line in (a) is enlarged in (b). A. Ishii et al. Materials Today Communications 35 (2023) 1060619misorientation-dependent grain boundary anisotropy and direction- dependent interface diffusion anisotropy [18]. Additionally, anisotropic properties depend on the material. Therefore, it is difficult to exhaustively predict the sintering behavior of particles with anisotropic interfaces based on sintering theory. To address such difficulties, the PF model developed in this study can provide a method for predicting the sintering behavior, because it allows simulations to be performed under various conditions. 4.3. Simulation of multiparticle sintering A simulation using 200 particles is performed to demonstrate that the PF model can reproduce a sintered structure with faceted particles. The particles are initially spherical and their radii are randomly set in the range of 15–25. The initial particle configuration and crystal orientations are randomly set. The computational domain size is 2563. The zero Neumann condition is applied to all boundaries of the computational domain. The simulation was performed using the Wisteria/BDEC-01 supercomputer system [51]. Fig. 11 shows snapshots of morphological changes and grain boundary evolution during the PF simulation. The particle shape shown in Fig. 11 (a) is described by the isosurface of ρ = 0.5. In Fig. 11 (b), the crystal grain orientation of each particle included in the region indicated by the dashed line shown in Fig. 11 (a) is described by arrows. The blue, white, and red arrows indicate the directions of the a-, b-, and c-axes, respectively. The origins of the arrows indicate the centers of mass of the particles. Fig. 11 (c) shows grain boundary distributions on the cross- section at x = 124 for the inside of the sintered compact (ρ ≥ 0.5). The sintered structure is formed by densification of particles, whereas their shape changes from spherical to cubic with orientation rotation. Grain boundary formulation, void annihilation, and grain growth occurred simultaneously in the interior of the structure. This result demonstrates the possibility that the developed PF model can simulate the formulation of actual sintered structures, such as those shown in Fig. 1, which have not yet been predicted. 5. Conclusions This study developed a new PF model for solid-state sintering capable of analyzing the microstructural evolution and macroscopic densification affected by anisotropic properties of the surface and grain boundaries, which depend on the crystal grain orientation and misorientation, respectively. Using quaternions to analyze crystal orientations with rotation due to the rigid-body motion of the sintered particles, 3D PF simulations of solid-state sintering with interface anisotropy were facilitated. The result of PF simulations using the developed model led to Fig. 11. Time evolution of (a) shape of sintered particles, (b) orientations of particles in the area enclosed by the dashed line shown in (a), and (c) microstructure on cross-section at x = 124 obtained from the PF simulation using 200 particles. Blue, white, and red arrows indicate the directions of a-, b-, and c-axes, respectively. The origins of the arrows indicate the centers of mass of particles. A. Ishii et al. Materials Today Communications 35 (2023) 10606110the following conclusions. 1) Some validations demonstrated that the developed PF model reasonably reproduces both the theoretically predicted equilibrium shape and grain boundary formation based on the strong anisotropy of the surface energy and the balance of grain boundary energies, respectively. 2) The PF simulations using two particles indicated that the change in the neck curvature due to anisotropic surface properties is a key factor in determining the rates of neck growth and densification. Moreover, the simulation results demonstrated that the developed PF model provides a means to predict neck growth and densification affected by surface and grain boundary anisotropies under various combinations of crystal grain orientations. 3) A sintering simulation using 200 particles demonstrated the potential of the developed PF model to reproduce the experimentally observed sintered structure with faceted crystals. This study used assumed values for multiple material parameters for the validation of the proposed PF model. Using the actual values for material parameters identified via state-of-the-art techniques such as data assimilation [44,52–58], the developed PF model can play a significant role in shedding light on the microstructural evolution during solid-state sintering in the real world. CRediT authorship contribution statement Akimitsu Ishii: Conceptualization, Methodology, Software, Validation, Investigation, Data Curation, Writing – Original Draft, Visualization, Formal analysis. Kondo Kyoyu: Validation, Investigation, Writing – Review & Editing. Akiyasu Yamamoto: Writing – Review & Editing, Supervision, Project administration, Funding acquisition. Akinori Yamanaka: Conceptualization, Methodology, Validation, Investigation, Resources, Writing – Review & Editing, Supervision, Project administration. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability The data that has been used is confidential. Acknowledgment This work was supported by the Strategic Basic Research Programs, Core Research for Evolutional Science and Technology (CREST), [Revolutional Materials Development] Revolutional material development by fusion of strong experiments with theory/data science, “Development of polycrystalline superconducting materials and magnets based on superconducting materials informatics” (Funding agency: JST, Japan) [Grant number: JPMJCR18J4]. References [1] R.M. German, Sintering: From Empirical Observations to Scientific Principle, Butterworth-Heinemann,, 2014. [2] M. Sagawa, S. Hirosawa, H. Yamamoto, S. Fujimura, Y. Matsuura, Nd–Fe–B permanent magnet materials, Jpn. J. Appl. Phys. 26 (1987) 785–800, https://doi. org/10.1143/JJAP.26.785. [3] K.J. Strnat, Rare earth-cobalt permanent magnets, in: E.P. Wohlfarth, K.H. J. Buschow (Eds.), Handbook of Ferromagnetic Materials, vol. 4, Elsevier, 1988, pp. 131–209, https://doi.org/10.1016/S1574-9304(05)80077-X. [4] Y. Takano, H. Takeya, H. Fujii, H. Kumakura, T. Hatano, K. Togano, H. Kito, H. Ihara, Superconducting properties of MgB2 bulk materials prepared by high- pressure sintering, Appl. Phys. Lett. 78 (2001) 2914–2916, https://doi.org/ 10.1063/1.1371239. [5] S. Tokuta, Y. Hasegawa, Y. Shimada, A. Yamamoto, Enhanced critical current density in K-doped Ba122 polycrystalline bulk superconductors via fast densification, iScience 25 (2022), 103992, https://doi.org/10.1016/j. isci.2022.103992. [6] D. Grossin, A. Montón, P. Navarrete-Segado, E. Özmen, G. Urruth, F. Maury, D. Maury, C. Frances, M. Tourbin, P. Lenormand, G. Bertrand, A review of additive manufacturing of ceramics by powder bed selective laser processing (sintering / melting): calcium phosphate, silicon carbide, zirconia, alumina, and their composites, Open Ceram. 5 (2021), 100073. [7] N. Tuncer, A. Bose, Solid-state metal additive manufacturing: a review, JOM 72 (2020) 3090–3111, https://doi.org/10.1007/s11837-020-04260-y. [8] T. Katase, Y. Ishimaru, A. Tsukamoto, H. Hiramatsu, T. Kamiya, K. Tanabe, H. Hosono, Advantageous grain boundaries in iron pnictide superconductors, Nat. Commun. 2 (2011) 409, https://doi.org/10.1038/ncomms1419. [9] T.T. Sasaki, T. Ohkubo, K. Hono, Structure and chemical compositions of the grain boundary phase in Nd-Fe-B sintered magnets, Acta Mater. 115 (2016) 269–277, https://doi.org/10.1016/j.actamat.2016.05.035. [10] M. Bah, R. Podor, R. Retoux, F. Delorme, K. Nadaud, F. Giovannelli, I. Monot- Laffez, A. Ayral, Real-time capturing of microscale events controlling the sintering of lead-free piezoelectric potassium-sodium niobate, Small 18 (2022), e2106825, https://doi.org/10.1002/smll.202106825. [11] Y. Zuo, C. Zhao, A. Robador, M. Wickham, S.H. Mannan, Quasi-in-situ observation of the grain growth and grain boundary movement in sintered Cu nanoparticle interconnects, Acta Mater. 236 (2022), 118135, https://doi.org/10.1016/j. actamat.2022.118135. [12] T. Takahashi, F. Sakamoto, J. Tatami, M. Iijima, In situ observation of evolution of internal structure of alumina during sintering by swept-source OCT, Int. J. Appl. Ceram. Technol. 19 (2022) 1171–1179, https://doi.org/10.1111/ijac.13909. [13] X.L. Phuah, J. Jian, H. Wang, X. Wang, X. Zhang, H. Wang, Ultra-high heating rate effects on the sintering of ceramic nanoparticles: an in situ TEM study, Mater. Res. Lett. 9 (2021) 373–381, https://doi.org/10.1080/21663831.2021.1927878. [14] R. Kobayashi, Modeling and numerical simulations of dendritic crystal growth, Phys. D. 63 (1993) 410–423, https://doi.org/10.1016/0167-2789(93)90120-P. [15] D. Choudhuri, L. Blake, Particle curvature effects on microstructural evolution during solid-state sintering: phenomenological insights from phase-field simulations, J. Mater. Sci. 56 (2021) 7474–7493, https://doi.org/10.1007/s10853- 021-05802-8. [16] Y.U. Wang, Computer modeling and simulation of solid-state sintering: a phase field approach, Acta Mater. 54 (2006) 953–961, https://doi.org/10.1016/j. actamat.2005.10.032. [17] S. Biswas, D. Schwen, J. Singh, V. Tomar, A study of the evolution of microstructure and consolidation kinetics during sintering using a phase field modeling based approach, Extrem. Mech. Lett. 7 (2016) 78–89, https://doi.org/ 10.1016/j.eml.2016.02.017. [18] S. Biswas, D. Schwen, H. Wang, M. Okuniewski, V. Tomar, Phase field modeling of sintering: role of grain orientation and anisotropic properties, Comput. Mater. Sci. 148 (2018) 307–319, https://doi.org/10.1016/j.commatsci.2018.02.057. [19] S. Biswas, D. Schwen, V. Tomar, Implementation of a phase field model for simulating evolution of two powder particles representing microstructural changes during sintering, J. Mater. Sci. 53 (2018) 5799–5825, https://doi.org/10.1007/ s10853-017-1846-3. [20] Z. Zhao, X. Zhang, H. Zhang, H. Tang, Y. Liang, Numerical investigation into pressure-assisted sintering using fully coupled mechano-diffusional phase-field model, Int. J. Solids Struct. 234–235 (2022), 111253, https://doi.org/10.1016/j. ijsolstr.2021.111253. [21] V. Ivannikov, F. Thomsen, T. Ebel, R. Willumeit-Römer, Capturing shrinkage and neck growth with phase field simulations of the solid state sintering, Model. Simul. Mater. Sci. Eng. 29 (2021), 075008, https://doi.org/10.1088/1361-651X/ac1f87. [22] R. Termuhlen, X. Chatzistavrou, J.D. Nicholas, H.-C. Yu, Three-dimensional phase field sintering simulations accounting for the rigid-body motion of individual grains, Comput. Mater. Sci. 186 (2021), 109963, https://doi.org/10.1016/j. commatsci.2020.109963. [23] W. Ma, Y. Shen, A phase field model for the solid-state sintering with parametric proper generalized decomposition, Powder Technol. 419 (2023), 118345, https:// doi.org/10.1016/j.powtec.2023.118345. [24] J. Hötzer, M. Seiz, M. Kellner, W. Rheinheimer, B. Nestler, Phase-field simulation of solid state sintering, Acta Mater. 164 (2019) 184–195, https://doi.org/10.1016/ j.actamat.2018.10.021. [25] M. Plapp, Unified derivation of phase-field models for alloy solidification from a grand-potential functional, Phys. Rev. 84 (2011), 031601, https://doi.org/ 10.1103/PhysRevE.84.031601. [26] I. Greenquist, M.R. Tonks, L.K. Aagesen, Y. Zhang, Development of a microstructural grand potential-based sintering model, Comput. Mater. Sci. 172 (2020), 109288, https://doi.org/10.1016/j.commatsci.2019.109288. [27] I. Greenquist, M. Tonks, M. Cooper, D. Andersson, Y. Zhang, Grand potential sintering simulations of doped UO2 accident-tolerant fuel concepts, J. Nucl. Mater. 532 (2020), 152052, https://doi.org/10.1016/j.jnucmat.2020.152052. [28] F. Abdeljawad, D.S. Bolintineanu, A. Cook, H. Brown-Shaklee, C. DiAntonio, D. Kammler, A. Roach, Sintering processes in direct ink write additive manufacturing: a mesoscopic modeling approach, Acta Mater. 169 (2019) 60–75, https://doi.org/10.1016/j.actamat.2019.01.011. [29] Y. Yang, O. Ragnvaldsen, Y. Bai, M. Yi, B.-X. Xu, 3D non-isothermal phase-field simulation of microstructure evolution during selective laser sintering, npj Comput. Mater. 5 (2019) 81, https://doi.org/10.1038/s41524-019-0219-7. A. Ishii et al. http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref1http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref1https://doi.org/10.1143/JJAP.26.785https://doi.org/10.1143/JJAP.26.785https://doi.org/10.1016/S1574-9304(05)80077-Xhttps://doi.org/10.1063/1.1371239https://doi.org/10.1063/1.1371239https://doi.org/10.1016/j.isci.2022.103992https://doi.org/10.1016/j.isci.2022.103992http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref6http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref6http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref6http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref6http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref6https://doi.org/10.1007/s11837-020-04260-yhttps://doi.org/10.1038/ncomms1419https://doi.org/10.1016/j.actamat.2016.05.035https://doi.org/10.1002/smll.202106825https://doi.org/10.1016/j.actamat.2022.118135https://doi.org/10.1016/j.actamat.2022.118135https://doi.org/10.1111/ijac.13909https://doi.org/10.1080/21663831.2021.1927878https://doi.org/10.1016/0167-2789(93)90120-Phttps://doi.org/10.1007/s10853-021-05802-8https://doi.org/10.1007/s10853-021-05802-8https://doi.org/10.1016/j.actamat.2005.10.032https://doi.org/10.1016/j.actamat.2005.10.032https://doi.org/10.1016/j.eml.2016.02.017https://doi.org/10.1016/j.eml.2016.02.017https://doi.org/10.1016/j.commatsci.2018.02.057https://doi.org/10.1007/s10853-017-1846-3https://doi.org/10.1007/s10853-017-1846-3https://doi.org/10.1016/j.ijsolstr.2021.111253https://doi.org/10.1016/j.ijsolstr.2021.111253https://doi.org/10.1088/1361-651X/ac1f87https://doi.org/10.1016/j.commatsci.2020.109963https://doi.org/10.1016/j.commatsci.2020.109963https://doi.org/10.1016/j.powtec.2023.118345https://doi.org/10.1016/j.powtec.2023.118345https://doi.org/10.1016/j.actamat.2018.10.021https://doi.org/10.1016/j.actamat.2018.10.021https://doi.org/10.1103/PhysRevE.84.031601https://doi.org/10.1103/PhysRevE.84.031601https://doi.org/10.1016/j.commatsci.2019.109288https://doi.org/10.1016/j.jnucmat.2020.152052https://doi.org/10.1016/j.actamat.2019.01.011https://doi.org/10.1038/s41524-019-0219-7Materials Today Communications 35 (2023) 10606111[30] Y. Yang, P. Kühn, M. Yi, H. Egger, B.-X. Xu, Non-isothermal phase-field modeling of heat–melt–microstructure-coupled processes during powder bed fusion, JOM 72 (2020) 1719–1733, https://doi.org/10.1007/s11837-019-03982-y. [31] X. Wang, Y. Liu, L. Li, C.O. Yenusah, Y. Xiao, L. Chen, Multi-scale phase-field modeling of layer-by-layer powder compact densification during solid-state direct metal laser sintering, Mater. Des. 203 (2021), 109615, https://doi.org/10.1016/j. matdes.2021.109615. [32] L.-X. Lu, N. Sridhar, Y.-W. Zhang, Phase field simulation of powder bed-based additive manufacturing, Acta Mater. 144 (2018) 801–809, https://doi.org/ 10.1016/j.actamat.2017.11.033. [33] T.L. Einstein, Equilibrium shape of crystals, in: T. Nishinaga (Ed.), Handbook of Crystal Growth, second ed., Elsevier, 2015, pp. 215–264. [34] A. Delanoë, S. Lay, Evolution of the WC grain shape in WC–Co alloys during sintering: effect of C content, Int. J. Refract. Met. Hard Mater. 27 (2009) 140–148, https://doi.org/10.1016/j.ijrmhm.2008.06.001. [35] S. Kim, S.-H. Han, J.-K. Park, H.-E. Kim, Variation of WC grain shape with carbon content in the WC–Co alloys during liquid-phase sintering, Scr. Mater. 48 (2003) 635–639, https://doi.org/10.1016/S1359-6462(02)00464-5. [36] H.B. Lin, Y.M. Zhang, H.B. Rong, S.W. Mai, J.N. Hu, Y.H. Liao, L.D. Xing, M.Q. Xu, X.P. Li, W.S. Li, Crystallographic facet- and size-controllable synthesis of spinel LiNi0.5Mn1.5O4 with excellent cyclic stability as cathode of high voltage lithium ion battery, J. Mater. Chem. A 2 (2014) 11987–11995, https://doi.org/10.1039/ C4TA01810A. [37] R.S. Qin, H.K.D.H. Bhadeshia, Phase-field model study of the effect of interface anisotropy on the crystal morphological evolution of cubic metals, Acta Mater. 57 (2009) 2210–2216, https://doi.org/10.1016/j.actamat.2009.01.024. [38] H. Li, Y. Du, J. Long, Z. Ye, Z. Zheng, H. Zapolsky, G. Demange, Z. Jin, Y. Peng, 3D phase field modeling of the morphology of WC grains in WC–Co alloys: the role of interface anisotropy, Comput. Mater. Sci. 196 (2021), 110526, https://doi.org/ 10.1016/j.commatsci.2021.110526. [39] S. Wise, J. Kim, J. Lowengrub, Solving the regularized, strongly anisotropic Cahn- Hilliard equation by an adaptive nonlinear multigrid method, J. Comput. Phys. 226 (2007) 414–446, https://doi.org/10.1016/j.jcp.2007.04.020. [40] A.A. Wheeler, Phase-field theory of edges in an anisotropic crystal, Proc. R. Soc. A 462 (2006) 3363–3384, https://doi.org/10.1098/rspa.2006.1721. [41] N. Moelans, B. Blanpain, P. Wollants, Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems, Phys. Rev. B 78 (2008), 024113, https://doi.org/10.1103/ PhysRevB.78.024113. [42] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, J. Chem. Phys. 28 (1958) 258–267, https://doi.org/10.1063/1.1744102. [43] S.M. Allen, J.W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Met. 27 (1979) 1085–1095, https://doi.org/10.1016/0001-6160(79)90196-2. [44] A. Yamanaka, Y. Maeda, K. Sasaki, Ensemble Kalman filter-based data assimilation for three-dimensional multi-phase-field model: estimation of anisotropic grain boundary properties, Mater. Des. 165 (2019), 107577, https://doi.org/10.1016/j. matdes.2018.107577. [45] H. Zhang, M.I. Mendelev, D.J. Srolovitz, Mobility of Σ5 tilt grain boundaries: Inclination dependence, Scr. Mater. 52 (2005) 1193–1198, https://doi.org/ 10.1016/j.scriptamat.2005.03.012. [46] E. Miyoshi, T. Takaki, Y. Shibuta, M. Ohno, Bridging molecular dynamics and phase-field methods for grain growth prediction, Comput. Mater. Sci. 152 (2018) 118–124, https://doi.org/10.1016/j.commatsci.2018.05.046. [47] T. Hirouchi, T. Tsuru, Y. Shibutani, Grain growth prediction with inclination dependence of 〈110〉 tilt grain boundary using multi-phase-field model with penalty for multiple junctions, Comput. Mater. Sci. 53 (2012) 474–482, https:// doi.org/10.1016/j.commatsci.2011.08.030. [48] W.T. Read, W. Shockley, Dislocation models of crystal grain boundaries, Phys. Rev. 78 (1950) 275–289, https://doi.org/10.1103/PhysRev.78.275. [49] F.J. Humphreys, A unified theory of recovery, recrystallization and grain growth, based on the stability and growth of cellular microstructures—I. The basic model, Acta Mater. 45 (1997) 4231–4240, https://doi.org/10.1016/S1359-6454(97) 00070-0. [50] H. Zheng, X.-G. Li, R. Tran, C. Chen, M. Horton, D. Winston, K.A. Persson, S.P. Ong, Grain boundary properties of elemental metals, Acta Mater. 186 (2020) 40–49, https://doi.org/10.1016/j.actamat.2019.12.030. [51] Wisteria/BDEC-01 supercomputer system. 〈https://www.cc.u-tokyo.ac.jp/superc omputer/wisteria/service/〉. (accessed February 13, 2023). [52] K. Sasaki, A. Yamanaka, S. Ito, H. Nagao, Data assimilation for phase-field models based on the Ensemble Kalman Filter, Comput. Mater. Sci. 141 (2018) 141–152, https://doi.org/10.1016/j.commatsci.2017.09.025. [53] A. Yamanaka, K. Takahashi, Data assimilation for three-dimensional phase-field simulation of dendritic solidification using the local ensemble transform Kalman filter, Mater. Today Commun. 25 (2020), 101331, https://doi.org/10.1016/j. mtcomm.2020.101331. [54] K. Takahashi, A. Yamanaka, Quantitative three-dimensional phase-field modeling of dendritic solidification coupled with local ensemble transform Kalman filter, Comput. Mater. Sci. 190 (2021), 110296, https://doi.org/10.1016/j. commatsci.2021.110296. [55] A. Ishii, A. Yamanaka, E. Miyoshi, Y. Okada, A. Yamamoto, Estimation of solid- state sintering and material parameters using phase-field modeling and ensemble four-dimensional variational method, Modell. Simul. Mater. Sci. Eng. 29 (2021), 065012, https://doi.org/10.1088/1361-651X/ac13cd. [56] A. Ishii, A. Yamanaka, E. Miyoshi, A. Yamamoto, Efficient estimation of material parameters using DMC-BO: application to phase-field simulation of solid-state sintering, Mater. Today Commun. 30 (2022), 103089, https://doi.org/10.1016/j. mtcomm.2021.103089. [57] A. Yamamura, S. Sakane, M. Ohno, H. Yasuda, T. Takaki, Data assimilation with phase-field lattice Boltzmann method for dendrite growth with liquid flow and solid motion, Comput. Mater. Sci. 215 (2022), 111776, https://doi.org/10.1016/j. commatsci.2022.111776. [58] Y. Matsuura, Y. Tsukada, T. Koyama, Adjoint model for estimating material parameters based on microstructure evolution during spinodal decomposition, Phys. Rev. Mater. 5 (2021), 113801, https://doi.org/10.1103/ PhysRevMaterials.5.113801. A. Ishii et al. https://doi.org/10.1007/s11837-019-03982-yhttps://doi.org/10.1016/j.matdes.2021.109615https://doi.org/10.1016/j.matdes.2021.109615https://doi.org/10.1016/j.actamat.2017.11.033https://doi.org/10.1016/j.actamat.2017.11.033http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref33http://refhub.elsevier.com/S2352-4928(23)00752-3/sbref33https://doi.org/10.1016/j.ijrmhm.2008.06.001https://doi.org/10.1016/S1359-6462(02)00464-5https://doi.org/10.1039/C4TA01810Ahttps://doi.org/10.1039/C4TA01810Ahttps://doi.org/10.1016/j.actamat.2009.01.024https://doi.org/10.1016/j.commatsci.2021.110526https://doi.org/10.1016/j.commatsci.2021.110526https://doi.org/10.1016/j.jcp.2007.04.020https://doi.org/10.1098/rspa.2006.1721https://doi.org/10.1103/PhysRevB.78.024113https://doi.org/10.1103/PhysRevB.78.024113https://doi.org/10.1063/1.1744102https://doi.org/10.1016/0001-6160(79)90196-2https://doi.org/10.1016/j.matdes.2018.107577https://doi.org/10.1016/j.matdes.2018.107577https://doi.org/10.1016/j.scriptamat.2005.03.012https://doi.org/10.1016/j.scriptamat.2005.03.012https://doi.org/10.1016/j.commatsci.2018.05.046https://doi.org/10.1016/j.commatsci.2011.08.030https://doi.org/10.1016/j.commatsci.2011.08.030https://doi.org/10.1103/PhysRev.78.275https://doi.org/10.1016/S1359-6454(97)00070-0https://doi.org/10.1016/S1359-6454(97)00070-0https://doi.org/10.1016/j.actamat.2019.12.030https://www.cc.u-tokyo.ac.jp/supercomputer/wisteria/service/https://www.cc.u-tokyo.ac.jp/supercomputer/wisteria/service/https://doi.org/10.1016/j.commatsci.2017.09.025https://doi.org/10.1016/j.mtcomm.2020.101331https://doi.org/10.1016/j.mtcomm.2020.101331https://doi.org/10.1016/j.commatsci.2021.110296https://doi.org/10.1016/j.commatsci.2021.110296https://doi.org/10.1088/1361-651X/ac13cdhttps://doi.org/10.1016/j.mtcomm.2021.103089https://doi.org/10.1016/j.mtcomm.2021.103089https://doi.org/10.1016/j.commatsci.2022.111776https://doi.org/10.1016/j.commatsci.2022.111776https://doi.org/10.1103/PhysRevMaterials.5.113801https://doi.org/10.1103/PhysRevMaterials.5.113801 Phase-field modeling of solid-state sintering with interfacial anisotropy 1 Introduction 2 Formulation 2.1 Phase-field model 2.2 Calculation methodology for orientation rotation 3 Conditions 3.1 Parameters 3.2 Anisotropy of surface energy and mobility 3.3 Anisotropy of grain boundary energy and mobility 4 Results and discussions 4.1 Validation 4.2 Investigation of neck growth and densification 4.3 Simulation of multiparticle sintering 5 Conclusions CRediT authorship contribution statement Declaration of Competing Interest Data availability Acknowledgment References