# Fileset

[片桐_Results_in_Engineering.pdf](https://mdr.nims.go.jp/filesets/eb414104-006a-4fba-83ff-78c8d53f721c/download)

## Creator

[Jun Katagiri](https://orcid.org/0000-0002-6399-1951), [Sukeharu Nomoto](https://orcid.org/0000-0001-6503-4699), [Masahiro Kusano](https://orcid.org/0000-0002-5061-0195), [Makoto Watanabe](https://orcid.org/0000-0002-5064-9583)

## Rights

[Creative Commons BY-NC Attribution-NonCommercial 4.0 International](https://creativecommons.org/licenses/by-nc/4.0/)

## Other metadata

[Effect of mushy zone constant in Voller-Prakash model on keyholing behaviour in laser powder bed fusion simulation](https://mdr.nims.go.jp/datasets/3ba222a0-5244-49a8-931b-4cf5fd8dbad7)

## Fulltext

Effect of mushy zone constant in Voller-Prakash model on keyholing behaviour in laser powder bed fusion simulationContents lists available at ScienceDirectResults in Engineeringjournal homepage: www.sciencedirect.com/journal/results-in-engineeringResearch paperEffect of mushy zone constant in Voller-Prakash model on keyholing behaviour in laser powder bed fusion simulationJun Katagiri, Sukeharu Nomoto, Masahiro Kusano, Makoto Watanabe ∗Additive Manufacturing Group, Research Center for Structural Materials, National Institute for Materials Science, Sengen 1-2-1, Tsukuba, 305-0047, Ibaraki, JapanA R T I C L E I N F O A B S T R A C T Keywords:Laser powder bed fusionMushy zone flow resistanceVoller-Prakash modelMultiphysics CFD simulationKeyhole pore formationThis paper is a study of the model constant of the Voller-Prakash model: the mushy zone constant (𝐶), which determines the flow resistance in the mushy zone during metal solidfication. A drag force-liquid fraction relation was derived, assuming that the dynamic and static pressures were equal. The 𝐶 values for the three morphologies were estimated by analytically calculating the morphological parameters of the Kozeny-Carman equation, which was the basis of the Voller-Prakash model. As a result, much smaller values of 𝐶 were obtained compared to the commonly used range of 𝐶 > 105. This was due to the drag force model included a parameter related to the grid size of the simulation. Based on previous studies of the scale of the dendrite structure and the grid size widely used in simulations, we concluded that the rectangular morphology was reasonable. To evaluate the effect of 𝐶 value on keyholing behaviour, a series of the multiphysics CFD simulations was performed varying 𝐶 , including the 𝐶for rectangular morphology and the commonly-used 𝐶 from past studies. The flow resistance force in the mushy zone increased with larger 𝐶 values, resulting in significant restriction of pore displacement. Consequently, pores formed with smaller 𝐶 values were immobilised farther away from the laser position. For larger values of 𝐶 , the flow resistance force had a significant effect, causing the mushy zone to behave like a solid. The pores that formed were trapped by the solid-like mushy zone before they could reach a stable shape and size.1. IntroductionLaser powder bed fusion (L-PBF) is a type of metal 3D additive manufacturing technology. Its advantage is that it can easily produce complex-shaped parts compared to conventional metal manufacturing processes. As L-PBF is a new processing technology, there are many unclear factors that affect the production of high-quality parts. Therefore, extensive research and development of L-PBF is being conducted worldwide. Note that there are a number of different additive manufacturing techniques. One of the most representative techniques is the fused deposition method (FDM). Researches and developments for the FDM have been extensively performed in the past studies. For example, Manola et al. have applied machine learning for optimisation of process parameters in the FDM [1]. Ramlee et al. have evaluated effect of post-printing material on mechanical property and reported that cement-reinforced printed object exhibited higher strength than the objects reinforced with the other materials [2]. Lestari et al. have investigated the impact of the FDM process parameters on printing time and weight of printed objects, and optimised the process parameters using the Taguchi method and the * Corresponding author.E-mail address: WATANABE.Makoto@nims.go.jp (M. Watanabe).response surface method [3]. Moreover, Nagaraja et al. have reviewed the utilisation of natural materials such as in the additive manufacturing [4]. Abd Aziz et al. have used the 3D printer for validation of the finite element method (FEM) simulations of mechanical properties of 3D scanned femur model [5]. Following the direct comparison of the FEM and experiment, they evaluated fixation methods for femur fractures based on the FEM simulations.The defects generated strongly depend on the laser processing conditions (laser power and scanning speed). A process map or process window, which is a diagram illustrating the defects produced by combining these two variables, is crucial for practical manufacturing [6]. It is important to determine the process window before manufacturing the product, as it varies depending on the material used. Generating the process window is a time-consuming task that involves single track experiments at different laser powers and scan speeds, track cutting and polishing, and SEM observation. To assist the experimental observation, multi-physics computational fluid dynamics (CFD) simulation of L-PBF has been extensively studied. Such simulations have advantages such as the ability to observe physical quantities that are difficult to observe https://doi.org/10.1016/j.rineng.2024.103567Received 23 January 2024; Received in revised form 8 August 2024; Accepted 26 November 2024 Results in Engineering 24 (2024) 103567 Available online 2 December 2024 2590-1230/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC license ( http://creativecommons.org/licenses/by- nc/4.0/ ). http://www.ScienceDirect.com/http://www.sciencedirect.com/journal/results-in-engineeringmailto:WATANABE.Makoto@nims.go.jphttps://doi.org/10.1016/j.rineng.2024.103567https://doi.org/10.1016/j.rineng.2024.103567http://creativecommons.org/licenses/by-nc/4.0/http://creativecommons.org/licenses/by-nc/4.0/J. Katagiri, S. Nomoto, M. Kusano et al. in experiments. Moreover, these numerical simulations, in conjunction with computer-aided design (CAD), offer a further notable advantage: they can significantly accelerate the design and prototype testing processes. For example, in the biomedical field, there have been studies that have employed FEM-based mechanical simulations to investigate materials and designs for total hip replacement techniques [7--9]. Moreover, the FEM [10] and TRIZ-based techniques [11] have been utilised in the design and development of mechanical components for surgical chairs. Note that CFD is a simulation method that is widely used in the design and prototyping processes across a broad range of industries, not only the metal 3D additive manufacturing but also mechanical [12,13] and the biomedical industries [14--16].The L-PBF is a complex multi-physics phenomenon involving heat conduction with melting and solidfication, therm-fluid dynamics of the molten metal, laser heating and powder particle dynamics. Such simulations require the coupling of several physical models. The coupling simulation of heat transfer, fluid dynamics and laser heating has been performed in previous numerous studies, e.g., [17--25]. Note that this study focuses on defect formation during metal solidfication in LPBF fabrication. The other type of simulation, e.g. thermo-elasto-plastic mechanical simulations, were also performed to evaluate the residual deformation after L-PBF fabrication [26] and solidfication crack [27]. In other cases, finite element simulation of heat transfer, taking into account melting and solidfication, is used to evaluate molten pool dimensions (e.g. molten pool width and depth) [6,28]. Thus, numerical simulation is becoming an essential tool for evaluating the phenomena occurring during L-PBF fabrication.The multiphysics CFD simulation naturally uses several sets of governing equations and a variety of physical models. Among such physical models, we focus on the flow resistance model in the solid-liquid coexistence zone, the so-called mushy zone. The mushy zone flow resistance model proposed by Voller and Prakash [29] has been used in a number of studies of multiphysics CFD simulation of L-PBF. The Voller-Prakash model is described as follows.𝑓𝑑,𝑖 = 𝐶(1 − 𝑓𝑙)2𝑓𝑙3 + 𝑞𝑢𝑖, (1)where 𝑓𝑑,𝑖, 𝑓𝑙 , 𝑞 and 𝑢𝑖 are the fluid drag force, the volume fraction of the liquid phase, a constant to prevent division by zero and the flow velocity, respectively. The parameter 𝐶 , referred to as the mushy zone constant, is assumed to be a constant determined according to the solidification morphology. The Kozeny-Carman equation [30--32], which is the basis of the Voller-Prakash model, includes morphological parameters such as the specific surface area of the porous media, the hydraulic tortuosity and the shape factor. The shape factor is a constant related to the cross-sectional shape perpendicular to the flow [30]. In Eq. (1), 𝐶combines a group of morphological parameters in the Kozeny-Carman equation into a single empirical parameter. Hence, the reduction in flow velocity due to the decrease in liquid fraction depends almost exclusively on 𝐶 . Furthermore, while Eq. (1) is the relationship between drag force and liquid fraction, it should be noted that the original KozenyCarman equation is a linear relationship between pressure gradient and liquid fraction (porosity) in a porous medium. Even with the same porosity, the morphology of porous media can be different. A representative morphological parameter included in the Kozeny-Carman equation is the specific surface area. The pressure gradient acts on the surface area of the porous media, which means that changes in morphology could lead to variations in drag force, even for the same porosity. However, we can hardly find the studies that investigated the constant 𝐶 in Eq. (1) from the perspective of morphology. Previous studies have used different values of 𝐶 , such as 𝐶 = 2.5×105 for the simulation of a latent heat thermal energy storage system [33], and 𝐶 = 106 [25], 5.56×106 [18], 105-107 [17] and 1014 [21] for the simulations of the L-PBF. In addition, previous studies have highlighted the ifluence of the mushy zone constant on melting and solidfication behaviour. Shmueli et al. have reported that 𝐶 significantly affects the area of the molten region inside the cylindrical tube for latent heat energy storage system (LHTES) application [34]. For L-PBF simulation, the authors of these studies have performed a series of multiphysics simulations varying 𝐶 and reported that 𝐶 has little effect on the dimensions of the molten pool when 𝑃= 300 W and 𝑉 = 1000 mm/s for the material of Inconel738LC [28]. Wang et al. have reported that the parameter 𝐶 can affect the collapse behaviour of the keyhole pore [18]. They have also reported that the value of 𝐶 affects the surface velocity distribution for the keyhole mode molten pool, implying that 𝐶 should affect the formation and displacement of the keyhole pores.The mushy zone constant 𝐶 is a parameter that rflects the ifluence of porous morphology; however, we have not found examples of previous studies that have dfined 𝐶 based on morphology. It is well known that dendrites form during metal solidfication. The dendrites are known to develop multiple columnar structures towards the normal direction of the solid/liquid interface. This means that the value of 𝐶 should be based on morphological parameters that rflect this columnar structure.In the light of the above, the following issues need to be considered for the mushy zone constant.• The relationship between pressure gradient and liquid fraction should be converted to a relationship between drag force and liquid fraction.• The impact of solidfication morphology should be taken into consideration.The aim of this study is to address the aforementioned issues. A relation between drag force and liquid fraction was derived by applying the widely known fluid drag equation and Bernoulli’s theorem to the Kozeny-Carman equation. The values of 𝐶 were then calculated for three morphology types of solidfication: (1) rectangular shape, (2) a cone, and (3) an assembly of cones, respectively. Flow resistance forces were determined for three different morphologies and compared with those calculated using the commonly used 𝐶 value, as a flow Reynolds number of 1 is the limit of the Kozeny-Carman equation. Furthermore, the impact of the mushy zone constant on keyholing behaviour in L-PBF was evaluated using the multiphysics CFD simulation.2. Kozeny-Carman-based mushy zone flow resistance model2.1. Drag force-liquid fraction relationThe Voller-Prakash model (Eq. (1)) is based on the Kozeny-Carman equation, which represents the relationship between the porosity of the porous media and the hydraulic conductivity (𝑘).𝑘 = 𝜌𝑔𝜇1 𝜙1 𝑆𝑠21 𝜏2𝜖3(1 − 𝜖)2, (2)where 𝑘, 𝜌, 𝑔, 𝜇, 𝜙, 𝑆𝑆 , 𝜏 and 𝜖 are the hydraulic conductivity of the porous media, the fluid density, the gravitational acceleration (9. 81 𝑚∕𝑠2), fluid viscosity, specific surface area (the ratio of the surface area to the volume of the porous media), hydraulic tortuosity, and porosity (determined by the ratio of pore volume to bulk volume of the porous media), respectively. The Kozeny-Carman equation is a commonly used equation for soil and rocks, with the fluid phase being water in porous media. This study specifically focuses on the melting and solidfication of metallic materials, assuming that the mushy zone is represented by either solid or liquid metal. When denoting the volume fractions of the solid and liquid phases as 𝑓𝑠 and 𝑓𝑙 , Eq. (3) can be obtained.𝑓𝑠 + 𝑓𝑙 = 1. (3)The relationship between the pore volume (𝑉𝑣), the porous volume (𝑉𝑠), the bulk volume (𝑉 ), and the porosity (𝜖) is as follows.𝑉𝑠 + 𝑉𝑣𝑉=𝑉𝑠𝑉+ 𝜖 = 1. (4)Results in Engineering 24 (2024) 103567 2 J. Katagiri, S. Nomoto, M. Kusano et al. From Eqns. (3) and (4), 𝑓𝑙 is used instead of 𝜖 in Eq. (2) for the purposes of this study.Darcy’s Law describes the relationship between the pressure drop due to porous flow (Δ𝑝) and the porous flow velocity (𝑈 ).𝑈 = 𝑘Δ𝑝 𝜌𝑔𝐿, (5)where 𝐿 is the length of the porous material in the direction of flow. By applying Eqns. (2) and (5), we can derive Eq. (6).Δ𝑝𝐿 = 𝜇𝜙𝑆𝑠2𝜏2(1 − 𝑓𝑙)2𝑓𝑙3 𝑈, (6)Voller and Prakash [29] formulated Eq. (1) under the assumption that the drag force is related to the pressure drop.In this study, we introduce the following assumption (Eq. (7)) to convert Eq. (6) into the relationship between drag force and liquid fraction.Δ𝑝 = 12𝜌𝑈2. (7)Eq. (7) means that the static pressure is perfectly converted into the dynamic pressure. Eq. (7) is known as Bernoulli’s law if the fluid is a non-viscous fluid (perfect fluid). The molten metal is a viscous fluid; however, Eq. (7) does not consider the dissipation due to viscosity. The justfication of Eq. (7) for the viscous fluid will be verfied by the CFD simulations in our future work.A commonly used equation for fluid drag is as follows.𝐹𝐷 = 12𝜌𝑈2𝐴𝐶𝐷, (8)where 𝐴 and 𝐶𝐷 are the cross-sectional area perpendicular to the flow direction and the drag coefficient, respectively. Substituting Eqns. (6) and (7) into Eq. (8) gives Eq. (9).𝐹𝐷 = 𝜇𝜙𝑆𝑠2𝜏2𝐿𝐴𝐶𝐷(1 − 𝑓𝑙)2𝑓𝑙3 𝑈. (9)The equation should be structured as follows: It is important that Eq. (9) includes the porous morphological parameters: 𝜏 , 𝑆𝑆 , 𝜙 and 𝐴. For this study, we calculated the specific values for the morphological parameters of simple solidfication geometries. Additionally, it should be noted that 𝐿, which refers to the grid size in numerical simulations, is also significant. Eq. (9) shows that the drag force is dependent on the scale, which distinguishes it from Eqns. (1) and (9).2.2. Rectangular morphologyFig. 1 shows an illustration of a morphology with a rectangular shape. Previous studies have observed the dendrite columnar structures. According these studies, the dendrite length is estimated to be in the range of a few hundred micrometres [35--40]. In contrast, the primary dendrite arm spacing (the distance between two dendrites) has been measured to be several micrometres [41--43]. In this study, it is assumed that dendrites begin to grow in the range of a few micrometres and extend to measured lengths of hundreds of micrometres or more. Therefore, in this section, the solidfication morphology is modelled as a rectangular shape at the onset of dendrite growth. The multiphysics CFD simulation of the L-PBF was performed using a grid resolution of several micrometres, which is significantly smaller than the several hundred micrometres dendrite structure observed in the experiment. Hence, the micrometre-grid provides only a partial surface prfile of the dendrite. Using the height (ℎ𝑠𝑜𝑙), area (𝐴𝑠𝑜𝑙) and volume (𝑉𝑠𝑜𝑙) of the solidfied metal part from Fig. 1, Eq. (10) can be used to express the specific surface area (𝑆𝑆 ).𝑆𝑠 =𝐴𝑠𝑜𝑙𝑉𝑠𝑜𝑙= 𝐿2ℎ𝑠𝑜𝑙𝐿2 = 1 ℎ𝑠𝑜𝑙. (10)Fig. 1. A rectangular-shaped solidification morphology. The drag coefficient (𝐶𝐷) is ifluenced by both the shape of the object and the Reynolds number (Re) of the fluid. The Kozeny-Carman equation is applicable for low Re flow regimes [44--46]. For the flow through porous media, the Re is typically between 1 and 10 [44]. In this study we have assumed Re = 1 to satisfy the constraints of the Kozeny-Carman equation. Note that Re = 𝑈𝐷∕𝜇, where 𝐷 is the characteristic length, for example, the average particle size is often used as 𝐷 for granular material. In this study we have chosen to use the grid size 𝐿 as for 𝐷. By setting 𝐿 to 2.5 μm, the velocity 𝑈 is about 0.14 m/s to obtain Re = 1. It is worth noting that the grid size of 2.5 μm was previously implemented in our multiphysics CFD simulations [47]. It is widely recognised that 𝐶𝐷 depends on the Re. In the low Re flow regime, specifically characterised by laminar flow, the 𝐶𝐷of a sphere is obtained by calculating 24∕Re. According to Stokes’ law, we have 𝐶𝐷 = 24 when Re = 1. According to a previous study using two-dimensional CFD simulation [48], it was found that the 𝐶𝐷 of a square is greater than that of a circular shape, with the respective 𝐶𝐷s being 11.72 and 12.9. In addition, the 𝐶𝐷 of a disc and a plate are of the order of Re = 𝑂(101) [49]. Therefore a value of 𝐶𝐷 = 10 was used in this study.The 𝜏 characterises the degree of twisting of the streamlines. In Fig. 1, 𝜏 = 1 because there is no flow twisting in the direction parallel to the solidfication interface (grey in Fig. 1). The flow is stopped by the solidfication interface in the direction perpendicular to the interface. Hence, we choose 𝜏 = 1 for Fig. 1.The 𝐴 is equal to 𝐿2 for the direction perpendicular to the solidification interface, but the flow along this direction is impeded by the interface. For the direction parallel to the interface, 𝐴 = ℎ𝑠𝑜𝑙𝐿. The 𝜇denotes the viscosity of the molten metal, which in this study has a fixed value of 0.005 Pas.The determination of 𝜙 is theoretically difficult; therefore, the first author of this study and co-authors determined the 𝜙 value by backcalculation with Eq. (2) using the hydraulic conductivity (𝑘) measured by the CFD simulations [50,51]. Unfortunately, the investigation of the 𝜙 value for metallic materials can not be found in the literature. Therefore, we set 𝜙 = 1, which means that the effect of 𝜙 is not considered in this study. To quantitatively determine the shape factor, we plan to obtain 3D dendrite structures by x-ray CT and perform CFD simulation using the 3D dendrite structure in our future work.The liquid fraction, denoted 𝑓𝑙 , can be expressed by Eq. (11) as shown in Fig. 1.𝑓𝑙 = 1 −ℎ𝑠𝑜𝑙𝐿 . (11)Substituting 𝜏 = 1, 𝜙 = 1, 𝜇 = 0.005 Pas, 𝐴 = ℎ𝑠𝑜𝑙𝐿, and Eqns. (10) and (11) into Eq. (9), we obtained the drag force-liquid fraction relationship for Fig. 1.Results in Engineering 24 (2024) 103567 3 J. Katagiri, S. Nomoto, M. Kusano et al. Fig. 2. A cone-shaped solidification morphology. 𝐹𝐷 = 0.1𝐿(1 − 𝑓𝑙)𝑓𝑙3 𝑈, (12)2.3. Conical morphologyFig. 2 displays a conical-shaped morphology. The model assumes a single dendrite tip grows within a single grid, meaning no multiple micro-cones exist in the grid. Eq. (13) calculates the 𝑆𝑆 for the geometry presented in Fig. 2.𝑆𝑠 =𝐴𝑠𝑜𝑙𝑉𝑠𝑜𝑙=√5𝜋𝐿24 𝜋𝐿212 = 3√5. (13)In the cone morphology, the flow field varies in different directions: left-right (and front-back) and top-bottom. The 𝜏 has a value of approximately one at the top of the cone in the left-right direction. In contrast, at the bottom of the cone, the streamline is half the circumference of the cone base, giving 𝜏 = 𝜋𝐿∕𝐿. Assuming that a representative 𝜏 for the left-right direction is at the centre of the cone height, then 𝜏 = (𝜋 + 1)∕2. The 𝐴 for the left-right direction equals 𝐿2∕2. For the topbottom direction of Fig. 2, the streamlines assume to be straight along the ridge of the triangle (projection view from the side of the cone); as a result, 𝜏 =√3∕2. The 𝐴 for the top-bottom direction is the circle’s area (the cone’s base), with a value of 𝜋𝐿2∕4.We presumed to 𝐶𝐷 = 10 for both the left-right and top-bottom directions although the 𝐴 is different for the left-right and top-bottom directions. In both of the left-right and top-bottom directions, 𝜙 is assumed to be 1. Eqns. (9), (14) and (15) can be derived for the left-right and top-bottom directions by incorporating the aforementioned morphological parameters.𝐹𝐷 = 4.82𝐿3(1 − 𝑓𝑙)2𝑓𝑙3 𝑈, (14)𝐹𝐷 = 2.65𝐿3(1 − 𝑓𝑙)2𝑓𝑙3 𝑈. (15)According to Eqns. (14) and (15), the difference in 𝐹𝐷 values between the left-right and top-bottom directions is insignificant. The key implication of Eqns. (14) and (15) is that the 𝐹𝐷 is contingent on 𝐿3, which becomes significantly smaller when the grid size is a few micrometres.2.4. Cone assemblyIn this section, we consider a geometry composed of multiple dendrite columns. Fig. 3 displays that the cone assembly morphology, i.e., Fig. 3. An example of solidfication morphology composed of cone assembly: 25 cones (𝑛 = 5).cones are regularly arranged on the bottom plane. 𝑛 cones are aligned along one direction of the parallel to the bottom plane; as a result, a total of 𝑛2 cones are placed in the grid. Fig. 3 is an example for 𝑛 = 5. To calculate 𝑆𝑆 for cone assembly composed of 𝑛 cones, we use Eq. (16).𝑆𝑆 =𝐴𝑠𝑜𝑙𝑉𝑠𝑜𝑙=(𝑛𝜋𝐿22 √1 + 1 4𝑛2)𝜋𝐿212 = 6𝑛√1 + 1 4𝑛2, (16)For the left-right direction, the value of 𝜏 is expressed by Eq. (17), assuming that the streamline length at the centre of the cone height represents the average of the entire geometry.𝜏 = 14((1 + 𝜋) + 1𝑛 )(for 𝑛 > 1). (17)Moreover, the 𝐴 is a collection of 𝑛 triangles; hence 𝐿2∕2.For the top-bottom direction, 𝜏 = √1 + 1∕4𝑛2 which is the ridge length of one of the triangles. The 𝐴 is 𝜋𝐿2∕4, which is the area of a group of 𝑛2 circles.For the both left-right and top-bottom directions, the values of 𝜙 and the 𝐶𝐷 are set to 1 and 10, respectively. The 𝐶𝐷 in the left-right as well as in the top-bottom directions are expressed in Eqns. (18) and (19), respectively.𝐹𝐷 = 0.90𝑛2((1 + 𝜋)𝑛+ 14𝑛 )2(1 − 1 4𝑛2)𝐿3(1 − 𝑓𝑙)2𝑓𝑙3 𝑈, (18)𝐹𝐷 = 1.41𝑛2(1 − 1 4𝑛2)2𝐿3(1 − 𝑓𝑙)2𝑓𝑙3 𝑈. (19)3. Results and discussions3.1. Comparison of drag forceTable 1 shows the constant values of the mushy zone constant used in past studies. The abbreviation LHTES in Table 1 denotes Latent Heat Thermal Energy Systems. The micrometre grid was used in the multiphysics CFD simulations in the past studies, as listed in Table 1. Despite the grid size, these studies used a mushy zone constant greater than 105. Our calculations of the 𝐹𝐷 values were carried out on a 2.5 μm grid size, with 𝑈 = 0.14 m/s to achieve Re = 1, and 𝑛 in Eqns. (18) and (19) set to 5. Fig. 4 shows the relationship between 𝐹𝐷 and 𝑓𝑙 using the Voller-Prakash model with the use of 𝐶 = 106 (equation (1)), alongside Eqns. (12), (14) and (15), in addition to Eqns. (18) and (19), when 𝑛 = 5. Results in Engineering 24 (2024) 103567 4 J. Katagiri, S. Nomoto, M. Kusano et al. Table 1Mushy zone constant (𝐶 in Eq. (1)) used in the past studies.𝐶 [kg/(𝑚3s)] Grid size Phenomenon 1.6 × 105 1/80 lattice unit Thermal cavity [29] 2.5 × 105 ≈1.81 mm LHTES [33] 105-107 4 μm L-PBF [17] 105-1010 ≈0.5 mm LHTES [34] 1 × 108 2×2×1 μm3 L-PBF [22] 5.56 × 106 4 μm L-PBF [18] 1 × 106 5 μm L-PBF [25] 1 × 1014 4 μm L-PBF [21] Fig. 4. Relationship between liquid fraction (𝑓𝑙) and drag force (𝐹𝐷) calculated by the Voller-Prakash model and proposed models in this study.The drag force calculated from the proposed models in this study, represented by Eqns. (12), (14), (15), (18) and (19), are significantly lower than that of the Voller-Prakash model. This is due to the dependence of the proposed models on the grid size (𝐿). In particular, Eqns. (14), (15), (18) and (19) have an ifluence on 𝐹𝐷 of the order of 𝐿3; hence their 𝐹𝐷 values are significantly smaller than those of Eq. (1).The large 𝐶 value has been used in most of the past studies. Li et al. have performed the multi-physics simulation with different 𝐶 values. They have reported the 𝐶 ifluenced the velocity fields inside the molten pool and the surface shape of the top of the molten pool [52]. In their study, the molten pool shape predicted from the simulation is reasonably agreed with that of the experimental observation.The authors of this study have cofirmed that 𝐶 has little effect on the molten pool shape predicted by the simulation [47]. The mushy zone is a region whose temperature is between the liquidus and solidus temperatures and can occur in two phases: the melting (heating) phase and the solidfication (cooling) phase. In the melting phase, the powder bed and substrate are rapidly melted and depressed by evaporation. During the vapor depression, the temperature of the vaporised gas rises rapidly and reaches more than 4,000 K [47]. The mushy zone occurs slightly inside the depressed metal surface. The high gas temperature results in a large spatial temperature gradient towards the interior of the metal. The spatial gradient between the solidus and liquidus temperatures will be cofined to a narrow region. The molten pool is the maximum area melted by the laser irradiation. Therefore, the mushy zone constant has little effect on the molten pool dimensions because the mushy zone is cofined to a very small area in the melting phase.In contrast, the mushy zone of the solidfication process will continue for a longer period of time than that of the melting phase, i.e., the dynamic and drastic process caused by the laser irradiation. As a result, the length of the molten pool along the direction of laser movement will be ifluenced by the mushy zone constant. Wang et al. have investigated the ifluence of the mushy zone constant on the keyhole behaviour and Fig. 5. Relationship between number of cones (𝑛2) and drag force calculated by Eqns. (18) and (19).reported that the mushy zone constant affects the velocity distribution of the keyhole type molten pool [18], which suggests that the mushy zone constant can ifluence keyhole pore formation and its dynamics.3.2. Solidfication morphology effectEquations (18) and (19) enable us to assess the effect of the number of cones (𝑛) on 𝐹𝐷. Fig. 5 illustrates the effect of 𝑛 on 𝐹𝐷 for 𝑓𝑙 values of 0.2, 0.5 and 0.8. Note that Eqns. (18) and (19) have been calculated using a value of 𝑈 = 0.14 m/s and 𝐿 = 2.5 μm. As the 𝑛 increases, the 𝐹𝐷 tends to increase monotonically; the increase in 𝐹𝐷 becomes moderate as 𝑛 increases. The 𝑆𝑆 strongly ifluences the change of the increasing trend. Eq. (9) proves that 𝜇, 𝜙, 𝐿 and 𝐶𝐷 are identical regardless of 𝑛: the two factors: 𝜏 and 𝑆𝑆 can modify 𝐹𝐷 .When 𝑛 = 2, the 𝜏 of the cone assembly for the left-right direction is about 2.20. The 𝜏 of the conical assemblies for the top-bottom direction is √1.25 when 𝑛 = 1. These 𝜏 values do not have a strong effect on 𝐹𝐷 . The previous studies have demonstrated that 𝜏 cannot take a large value for dense packing of granular materials and rock. The morphological parameters in the Kozeny-Carman equation obtained in the past studies and the 𝐶 values estimated from such morphological parameters in the past studies are summarised in Table 2. The sandstone formed by consolidation and particle crushing under high effective stress conditions in the deep subsurface has 𝜏 = 3.1 for the porosity of 0.25 [53]. Unlike 𝜏 , the 𝑆𝑆 is the same for the left-right and top-bottom directions. In addition, 𝑆𝑆 is not a dimensionless parameter, which means that the 𝑆𝑆 value varies significantly depending on the microstructural morphology and should be a highly ifluential factor for 𝐹𝐷 . Note that 𝐴 is the projected cross-sectional area, which is different from the surface area used to calculate 𝑆𝑆 . That is, the surface area is the total area of the solid surface in contact with the liquid surface. It is reasonable from a physical point of view that the greater the surface area, i.e. the more complex the morphology, the greater the drag force.3.3. Revisit on mushy zone constantAs shown in Table 1, different values of 𝐶 have been used in past studies, mostly greater than 105. As described in the introduction, the single parameter 𝐶 in Eq. (1) has been derived from a set of morphological parameters of the Kozeny-Carman equation. In this study, we have derived the novel fluid drag-liquid fraction relationship (Eq. (9)) in order to keep the morphological parameters of the Kozeny-Carman equation. Prior to the discussion, we will exclude 𝐿 and 𝐴 in Eq. (9) because they are grid-dependent parameters. The 𝐶 values were calculated using the measured values of 𝜙, 𝜏 and 𝑆𝑆 : 𝐶 = 𝜙𝜏2𝑆𝑆2.Results in Engineering 24 (2024) 103567 5 J. Katagiri, S. Nomoto, M. Kusano et al. Table 2The mushy zone constant values estimated from the past studies.Material 𝜖 𝜙 𝜏 𝑆𝑆 [m−1] 𝐶 [kg/(m3s)] Ref. Cylinder - 2.5 √2 4∕𝐷 80∕𝐷2 [30,31] Spheres - 2.5 √2 6∕𝐷 180∕𝐷2 [30,31] Al-Cu alloy ≈0.25 - - 1.35×105 1.82×1010𝜙𝜏2 [54] Zinc battery - - - 1.04×106 1.08×1012𝜙𝜏2 [55] Al-15.5%Cu alloy - - - 1.62×105 2.62×1010𝜙𝜏2 [56] Al-19.5%Cu alloy - - - 1.08×105 1.17×1010𝜙𝜏2 [56] Ni-Al-W alloy (C) 0.1 - 1.30 1.86×105 5.85×1010𝜙 [57] Ni-Al-W alloy (C) 0.2 - 1.18 4.94×105 3.40×1011𝜙 [57] Ni-Al-W alloy (C) 0.3 - 1.16 5.75×105 4.45×1011𝜙 [57] Ni-Al-W alloy (C) 0.4 - 1.25 4.74×105 3.51×1011𝜙 [57] Ni-Al-W alloy (C) 0.5 - 1.19 2.05×105 5.95×1010𝜙 [57] Ni-Al-W alloy (C) 0.6 - 1.24 7.37×104 8.35×109𝜙 [57] Ni-Al-W alloy (C) 0.6 - 1.27 3.73×104 2.24×109𝜙 [57] Ni-Al-W alloy (P) 0.25 - 1.51 2.63×105 1.58×1011𝜙 [57] Ni-Al-W alloy (P) 0.35 - 1.19 8.90×105 1.12×1012𝜙 [57] Ni-Al-W alloy (P) 0.45 - 1.32 6.62×105 7.64×1011𝜙 [57] Ni-Al-W alloy (P) 0.55 - 1.24 4.18×105 2.69×1011𝜙 [57] Ni-Al-W alloy (P) 0.65 - 1.10 1.67×105 3.37×1010𝜙 [57] Ni-Al-W alloy (P) 0.75 - 1.24 4.41×104 2.99×109𝜙 [57] Ni-Al-W alloy (P) 0.85 - 1.07 3.55×104 1.44×109𝜙 [57] Glass beads 0.33 4.60 1.31 3.69×104 1.07×1010 [50] Glass beads 0.37 5.90 1.26 4.15×104 1.61×1010 [50] Glass beads 0.40 8.60 1.24 4.22×104 2.35×1010 [50] Fine sand 0.39 6.17 1.29 5.48×104 2.99×1010 [50] Fine sand 0.40 7.08 1.29 5.27×104 3.27×1010 [50] Fine sand 0.42 9.35 1.28 5.33×104 4.35×1010 [50] Cylinders (C) 0.223 0.60 1.03 1.63×104 1.69×108 [51] Spheres 0.472 1.90 1.21 1.89×104 9.94×108 [51] Sandpack LV60A 0.38 - - 2.21×104 3.54×109 [53] Sandpack F62A 0.33 - - 1.47×104 1.29×109 [53] Sandstone 0.25 - - 2.73×104 9.09×109 [53] Carman, who proposed the Kozeny-Carman equation [30,31], suggested that 𝜏 =√2 and 𝜙 = 2.5 from the geometrical consideration for a regular array of monodisperse spheres [30]. He also proposed 𝑆𝑆as 4∕𝐷 and 6∕𝐷 for monodisperse cylinders and spheres respectively. Note that 𝐷 is the diameter of the cylinder or sphere. As a result, 𝐶 for the monodisperse cylinders and spheres are 80∕𝐷2 and 180∕𝐷2 respectively. It is important to note that 𝐶 values estimated using Carman’s morphological parameters include the factor 1∕𝐷2 . As the 𝐷 becomes smaller, the 𝐶 increase in the order of 𝐷2.For the real metallic materials, the value of 𝑆𝑆 has been determined by x-ray CT. The values of 𝑆𝑆 obtained from x-ray CT images are listed in Table 2 and range between 104 and 106 𝑚−1, although the materials have been varied in the literature. From Eq. (9), the drag force depends on the square of 𝑆𝑆 . The factor 𝜙𝜏2 is assumed to be greater than 1; as a result, the value of 𝐶 will be in the order of 1010 and 1012. Madison et al. have performed flow simulations using the 3D columnar structures obtained from x-ray CT as boundary conditions and directly evaluated 𝜏 [57]. According to their simulations, 𝜏 differs between flow parallel and perpendicular to the dendrite growth direction. The term ``(C)'' in Table 2 refers to the direction perpendicular to the dendrite growth, while the term ``(P)'' refers to the same direction as the dendrite growth. Madison et al. have found that the difference in 𝜏 due to flow direction is insignificant, implying that 𝐶 is insignificantly affected by flow direction. In the model derived in this study, there is an insignificant difference in the 𝐶 values between the left-right and top-bottom directions, which follows the findings of Madison et al.Table 2 contains data obtained from x-ray CT analysis and flow simulation using the 3D porous structure of the granular material and rock. Although the morphology of the dendrite columnar structure and that of the granular material and rock are clearly different, the values of 𝜏and 𝑆𝑆 are found to be similar, which suggest that their 𝐶 values is also similar.To the best of the authors’ knowledge, there is no direct method to determine the shape factor 𝜙; therefore, the first author of this study and co-authors have estimated the 𝜙 value by back-calculation using the hydraulic conductivity, 𝜏 , 𝑆𝑆 and 𝜖 calculated by the flow simulation [51,50]. The 𝜙 values estimated by the back-calculation are listed in the Table 2 and are of the order of 𝑂(100). In the case of glass beads and fine sand, 𝐶 is found to be of the order of 𝑂(1010). In the case of dendrite structures of alloy, Madison et al. [57] did not provide the values of 𝜙. Assuming that 𝜙 = 10 is applied to the dendrite structure, the 𝐶 values are expected to be in the order of 𝑂(1010) to 𝑂(1012). Furthermore, 𝜏2 is estimated to be at most 1.5 in the CFD simulation of Madison et al. [57]. Assuming that 𝜏=1.5 and 𝜙 is 10, 𝜙𝜏2 ≈ 22.5; as a result, the 𝐶 values for the other alloys [54,56] are of the order of 𝑂(1011). These results suggest that the large 𝐶 values used in previous studies are reasonable.However, there is actually the characteristic length for the alloy dendrite structure; therefore, the solidfication morphology should be varied depending on the characteristic length, i.e., the grid size of the simulation. Several previous studies have reported that the scale of the dendrite is estimated to be between several hundred μm to several mm, for example, [35--40]. In contrast, several micrometre grid has been used in most of previous studies for the multiphysics CFD simulation, see., Table 1. As shown in Table 1, Voller and Prakash have proposed the 𝐶 value of Eq. (1) for relatively large grid sizes. When the grid size is significantly larger than the scale of the dendrites, it is reasonable to apply the commonly used 𝐶 values: between 105 and 107. However, due to the characteristic length scale in the dendrite formation mechanism, complex morphologies, such as cone assembly are difficult to form in several micrometre grids, as in the studies listed in Table 1 and the present study. For the simulation using several micrometre grid, the rectangular morphology could be reasonable among the models proposed in this study.4. Ifluence of mushy zone constant on keyholing behaviour4.1. Simulation overviewTo validate the proposed models, the keyhole pore formation simulation was carried out following Wang et al. [18]. The simulation includes Results in Engineering 24 (2024) 103567 6 J. Katagiri, S. Nomoto, M. Kusano et al. Fig. 6. Particle size distribution of the powder. multiphysics of the L-PBF such as solid-liquid-vapour phase change and its multiphase flow, heat conduction with the vaporisation-meltingsolidification, laser heating by ray tracing. In addition, the simulation with the powder bed generated by the discrete element method (DEM) simulation. A brief summary of the simulation procedure is given in this paper; the full simulation methodology and validation through systematic comparison with the experiments can be described in detail in our previous publications [47,58].4.1.1. Discrete-element modelling of powder spreadingWe used the general purpose DEM software: LIGGGHTS (version 3.8.0) in this study. The particle size distribution shown in Fig. 6 were used in the DEM simulation. The Young’s modulus and Poisson’s ratio of the particles were set to 5×106 N/m and 0.4, respectively. Both the particle-particle and the particle-wall friction coefficients were set to 0.3. The ``f'' style boundary condition was adopted in this study. The ``f'' boundary acts as a fixed smoothed wall boundary, but the particles are removed when the particles fall outside the domain. Note that it is difficult to completely reproduce the behaviour of physical phenomena in physical simulations; therefore, certain assumptions are often incorporated. For example, a homogeneity, an isotropy, and a linear elasticity have been employed in the mechanical simulation of elastic bodies [59--62]. In the context of the DEM, the particles, which would ideally be modelled as elastic bodies, are instead assumed to be rigid particles. The two elastic particles undergo deformation when in contact with one another. The elastic particles in contact are simulated by allowing for slight overlap of the rigid particles in the DEM. In the case of a low contact force, the deformation of elastic particles in contact is small. In this case, the rigid-body assumption is reasonable. In the DEM, the overlap distance between the particles mainly depends on the time increment and the particle contact stiffness (Young’s modulus). In this study, the time increment of 2 × 10−8 sec was identfied by a parametric study, as it ensures that the overlap between the particles remains within acceptable distance.The width (𝑥), depth (𝑦), and height (𝑧) of the DEM model domain are 6.2 mm, 2 mm, and 0.75 mm. The 80,000 powder particles were randomly generated to a rectangular region above a powder coating space which was a concave and rectangular space on top of the DEM model. The width, depth, and height of the powder coating stage were 4.5 mm, 1.5 mm, and 0.08 mm, respectively. Then, the powder particles were settled under the gravity. After the gravitational deposition, a rectangular-shaped coating blade moves toward 𝑥-direction with 25 mm/s speed to spread the powder particles within the powder coating space. A snapshot during powder spreading simulation is shown in Fig. 7. 4.1.2. Governing equation of multiphysics CFD simulationThe positions and radii of the powder particles were used as the initial conditions for the multiphysics CFD simulation. We used the multiphysics CFD solver based on OpenFOAM-8, which was partly developed by TERRABYTE Co. Ltd, Tokyo, Japan.In the multiphysics CFD simulation, both gas and liquid (solid and liquid metal) are assumed to be incompressible fluids. The interface between solid or metal (phase 1) and gas (phase 2) is captured by the function 𝐹 using the Volume of Fluid (VOF) method.𝜕𝐹𝜕𝑡 +𝜕𝐹𝑢𝑗𝜕𝑥𝑗= 0, (20)where, 𝐹 is the value of the VOF function (0≤𝐹≤1). The time evolution of the interface is calculated by solving the advection equation of 𝐹 .The governing equations used in the simulations are the energy equation, the continuity equation and the conservation of momentum, which are expressed in Eqns. (21), (22) and (23) respectively.𝜕(𝜌𝐶𝑝𝑇 )𝜕𝑡 +𝜕(𝜌𝑢𝑗𝐶𝑝𝑇 )𝜕𝑥𝑗= −𝜕𝑞𝑗𝜕𝑥𝑗+ ̇𝑞𝑚 + ̇𝑞𝑣, (21)𝜕𝑢𝑖𝜕𝑥𝑖= 0, (22)𝜕(𝜌𝑢𝑖)𝜕𝑡 +𝜕(𝜌𝑢𝑖𝑢𝑗 )𝜕𝑥𝑖= −𝑔𝑗(𝑥𝑗 − 𝑟𝑗) 𝜕𝜌 𝜕𝑥𝑖+ 𝜕𝜕𝑥𝑗{𝜇(𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖)}+ 𝐹𝐶𝑆𝐹 ,𝑖 + 𝐹𝑟,𝑖 +𝑆𝑖, (23)where 𝜌, 𝐶𝑝, 𝑢𝑖, 𝑟𝑗 , 𝑞𝑗 , 𝜇, ̇𝑞𝑚 and ̇𝑞𝑣 are the density, specific heat capacity, flow velocity, reference point, heat flux, viscosity, phase change heat dissipation term and phase change heat absorption term, respectively. Eqns (21), (22) and (23) are solved using the Finite Volume Method (FVM). In the FVM algorithm, a computational domain is partitioned by a set of polyhedral volume elements: control volumes. The governing equations are discretised to algebraic equations that are conservative for each control volume. The FVM in OpenFOAM is described in detail in the cited references [63--66]. The FVM has been a standard tool for CFD in both academic and industrial engineering for many years. The FVM is distinguished by its simplicity of algorithmic structure and robustness. Furthermore, the FEM is frequently employed for CFD. However, the FEM is applicable to a variety of physical phenomena described by partial differential equations, including thermal fluid behaviour of the L-PBF [20] and the mechanical behaviour of materials [67--71]. In the manufacturing process of L-PBF, the thermal fluid behaviour is of great importance. However, warping (plastic deformation) due to the residual stress cannot be evaluated by the FVM. In order to evaluate such phenomena, it is necessary to simulate thermo-elastoplastic behaviour. The FEM is suitable for performing simulations of both thermal fluid behaviour analysis and thermo-elastoplastic analysis.The terms 𝐹CSF,𝑖, 𝐹𝑟,𝑖 and 𝑆𝑖 are the surface tension term using the continuum surface force (CSF) model, the recoil pressure and the mushy zone flow resistance force, respectively. The surface tension, recoil pressure and mushy zone flow resistance models are expressed in the Eqns. (24), (25) and (26) respectively.𝐹𝐶𝑆𝐹 ,𝑖 = 𝜎𝜅𝜕𝐹𝜕𝑥𝑖, (24)𝐹𝑟,𝑖 = −0.54𝑃0 exp((𝐿𝑣𝑇 − 𝑇𝑣𝑅𝑇𝑇𝑣))𝐧, (25)𝑆𝑖 = −𝐶(1 − 𝑓𝑙)2𝑓𝑙3 𝑢𝑖, (26)where 𝜎, 𝜅, 𝑃0, 𝐿𝑣, 𝑅, 𝑇𝑣, and 𝐧 are the surface tension coefficient, the curvature of the gas-liquid interface, the gas phase pressure, the latent heat of evaporation, the gas constant, the evaporation temperature, and a unit normal vector of the gas-liquid interface, respectively.Results in Engineering 24 (2024) 103567 7 J. Katagiri, S. Nomoto, M. Kusano et al. Fig. 7. A snapshot of powder spreading DEM simulation. The ray tracing algorithm treats the laser beam as a bundle of laser beams and tracks each beam. Thermal energy is applied to the metal surface when the beam reaches the surface and is calculated from a laser rflectivity (absorptivity) of the metal surface. The laser rays repeat the mirror rflections until the applied thermal energy becomes extremely small or the laser ray does not reach the surface.In physical simulations, it is of great importance to validate the accuracy of the simulations by assessing the fidelity with which they reproduce real-world processes. For example, in simulations of elastic materials, it is essential to determine a boundary condition that is based on an actual physical process. The resulting strains or stresses should be compared with experimental data or real-world observations [72--76]. In our previous studies, we conducted a quantitative comparison between the molten pool dimensions obtained from simulations and those obtained from experiments [47,58]. These studies have demonstrated that the molten pool dimensions obtained from the simulations were influenced by a number of factors, including vaporisation-induced recoil pressure and the laser rflectivity models used and their respective values. Note that the impact of the laser rflectivity on the molten pool dimensions is significant [47]. In the case without powder bed, it has been cofirmed that the constant rflectivity value of 0.2 accurately reproduces the molten pool dimensions observed in experiments [58]. For this reason, the constant laser rflectivity value of 0.2 has also been adopted in this study.We focused on the rectangular region within the powder-coated space, whose width and depth were 1.5 mm and 0.3 mm respectively. The upper and lower regions of the powder bed were added to form gas phase and substrate regions, respectively. The resulting CFD model domain was 1.5 mm wide, 0.3 mm deep and 0.4 mm high. The entire model domain was partitioned into 5 μm and then the inner region of 1.4 mm width, 0.2 mm depth and 0.31 mm height was partitioned into 2.5 μm grid to ensure calculation accuracy.A constant temperature of 298.15 K was applied to the top 𝑥-𝑦 plane, and an adiabatic condition was applied to the other five planes for the thermal boundary conditions. An atmospheric pressure (101,325 Pa) was adopted for the upper 𝑥-𝑦 plane, and a zero-gradient condition was adopted for the other five sides for the pressure boundary conditions. For the flow velocity boundary conditions, the pressureInletOutletVelocity boundary condition implemented in OpenFOAM was adopted for the upper 𝑥-𝑦 plane. The pressureInletOutletVelocity condition assumed a zero gradient condition when the flow is outward direction. When the inward flow, the velocity is obtained from the flow with the specfied inlet direction. The zero-gradient condition was applied to the other five sides. These boundary conditions were set based on past studies; however, they do not entirely rflect the actual conditions that would occur in a real-world setting. For Example, the effect of a vapor plume has not been incorporated in this study. The vapor plume is the metal vapor generated by laser irradiation. According to previous study, the metal vapor forms like a cloud and can reach temperatures about 5500 K [77]. Table 3Materials properties of Inconel738LC and the simulation parameters used in this study.Property or Parameter Value [Unit] Density of solid 8820 [kg/m3] Density of liquid 8220 [kg/m3] Density of gas 1.176 [kg/m3] Solidus temperature 1371.15 [K] Liquidus temperature 1620.15 [K] Evaporation temperature 3034.0 [K] Thermal conductivity 46.1 [W/(m⋅K)] Latent heat of fusion 25,000 [J/kg] Latent heat of evaporation 734,000 [J/kg] Specific heat capacity of metal 710 [J/kg⋅K] Specific heat capacity of gas 1007 [J/kg⋅K] Surface tension coefficient (𝜎) 1.2 [N/m] Liquid viscosity (𝜇) 0.005 [Pa⋅s] Gas viscosity 1.865×10−5 [Pa⋅s] Liquid-gas heat transfer coefficient 40 [W/(m2⋅K)] Emissivity of thermal radiation 0.3 [-] Laser diameter 80.0×10−6 [m] Value of 𝐶 in Eq. (26) 10−6,100,1012,1018 [kg/(m3s)] Such high-temperature metal vapor may result in secondary heating the powder bed and the substrate. Indeed, our experiments revealed that the molten pool width near the top (i.e., near the substrate) was observed to be greater than that observed at the bottom or middle, for example, [28]. In the simulation, the upper thermal boundary condition along 𝑧 direction is maintained at 298.15 K, which results in the continuous cooling of the molten metal. Consequently, the molten pool width near the substrate may be narrower than that at the middle region. To accurately reproduce the molten pool shape observed in the experiments, it may be necessary to incorporate the heating effect by the vapor plume into the boundary conditions.The laser energy distribution is a Gaussian distribution and its standard deviation is set to 1.0. The laser power and scanning speed were set to 300 W and 500 mm/s respectively. The laser was initially positioned at 𝑥 = 2.7 mm, 𝑦 = 1.0 mm, and 𝑧 = 0.1 mm, respectively, and then displaced horizontally toward the 𝑥 direction by 1.3 mm for 2.6×10−3sec. Then the laser was stopped and the spontaneous cooling simulation was performed for 0.65 sec. The computation time of the typical simulation case was about 7 days using a 96-core parallel computer. The value of 𝐶 was varied by 10−6, 100, 106, 1012 and 1018. As in the previous study [47], we used the physical properties of Inconel738LC obtained using thermodynamics software: JMatPro. The physical properties and calculation parameters are listed in Table 3. Fig. 8 shows snapshots during the simulation with 𝐶=106. The laser process parameter condition of 300 W and 500 mm/s results deep vapour depression. The tip of the vapour depression then splits and is trapped inside the substrate, which is a typical keyhole pore formation behaviour cofirmed in the previous study, e.g. [78].Results in Engineering 24 (2024) 103567 8 J. Katagiri, S. Nomoto, M. Kusano et al. Fig. 8. Snapshots during the simulation for 𝐶 = 106 kg/(m3s). 4.2. Results and discussionsAs a result of the simulations, we cofirmed that the number, size and position of the pores were different when just the 𝐶 value was varied. Therefore, we evaluated the position and size of each individual pore. We used the visualisation software Paraview (ver. 5.10.1) to generate a bounding box for each pore and measured its length along the 𝑥, 𝑦 and 𝑧 axes. We calculated the pore size (𝐿𝑝𝑜𝑟𝑒) using the following formula𝐿𝑝𝑜𝑟𝑒 =𝑙𝑏𝑥 + 𝑙𝑏𝑦 + 𝑙𝑏𝑧3 , (27)where 𝑙𝑏𝑥, 𝑙𝑏𝑦 and 𝑙𝑏𝑧 are the bounding box lengths along the 𝑥, 𝑦and 𝑧 directions respectively. Note that the bounding box is centred at (𝑙𝑏𝑥∕2, 𝑙𝑏𝑦∕2, 𝑙𝑏𝑧∕2). We have calculated the positions and sizes of all the pores. The average values and standard deviations are shown in Fig. 9.There is no correlation between the position of the 𝑦 axis and 𝐶 , nor between the position of the 𝑧 axis and 𝐶 . However, there is a correlation between the position of the 𝑥 axis (along the direction of the laser scan) and 𝐶 . Fig. 10 shows the velocity vectors in the 𝑥-𝑧 cross section at 𝑦=0.15 mm for the simulation with 𝐶=10−6 and 1018. The two solid lines in the figure represent the contours of the liquidus and solidus temperatures respectively. The area between the two contours is the mushy zone. The direction of flow in both the molten pool and the mushy zone is in the negative 𝑥 direction; hence the pore displaces in the opposite direction to the laser scan. The use of a large value of 𝐶 results in strong flow resistance on the melt, which significantly restricts pore displacement in the mushy zone. As a result, the position of the pore along the 𝑥 axis is smaller when the smaller 𝐶 value is used.The pore size becomes slightly smaller with increasing 𝐶 and is almost the same for 𝐶 ≥ 106. The standard deviations of the pore size are larger for 𝐶 ≥ 100 than for 𝐶 < 100. As shown in Fig. 10, the flow velocity distribution in the mushy zone is very small when 𝐶 is large, which means that the mushy zone can be considered almost solid. In other words, the length of the molten pool varies as 𝐶 increases. Huang et al. visualised the pore formation and its dynamics during the L-PBF experiment using synchrotron x-ray CT [79]. According to them, the pores pinched off from the tips of the keyhole grow rapidly. Huang et Fig. 9. Influence of mushy zone constant (𝐶) on position and size of pores. Results in Engineering 24 (2024) 103567 9 J. Katagiri, S. Nomoto, M. Kusano et al. Fig. 10. Velocity distributions in the 𝑥-𝑧 cross section at 𝑦=0.15 mm for the simulation for different 𝐶 values. al. suggest that this is due to pressure equalisation. According to the hypothesis of Huang et al., the pores generated in the simulation rapidly change shape. Subsequently, during the cooling process, Huang et al. proposed the hypothesis that the pore diameter shrinks due to condensation of vapour inside the pore. The shape change of the pore immediately after its formation in the molten region: temperature between 𝑇𝑙 and 𝑇𝑣. Thereafter, the shape change gradually decreases as the pores pass through the mushy zone. As 𝐶 is increased, the pores are trapped in the solid before it becomes stable in shape: spherical in most cases. The mushy zone is considered to be the solid metal because of the high flow resistance when 𝐶 > 106. Therefore, we assume that the mean and variation of pore size are almost identical when 𝐶 > 106. In other words, the dynamics after pore formation will not be reproduced by the simulation using the large 𝐶 value. By using an appropriate (plausible) value of 𝐶 , we believe we can reproduce the phenomenon of pore dynamics, i.e. the pores emerge from the molten pool after the shape change and their displacement is gradually slowed down in the mushy zone. In this study, the mushy zone constant of 10−6 was proposed. In order to verify the acceptability of this value, it is necessary to carry out experiments and visualisations under the same conditions as the simulation and to evaluate the pore dynamics.As described above, the Kozeny-Carman equation, which is the basis of the Voller-Prakash model, has been applied to very slow porous flows: Re ≈ 1. However, as shown in Fig. 10, the flow velocity at the bottom of the molten pool was high. For such a high flow velocity regime, we should consider the applicability of the Kozeny-Carman equation. The pressure gradient - porosity (liquid fraction) relationship that can be applied to the high flow velocity regime (inertial flow) has been proposed in the past studies, e.g. Ergun’s equation [80] and Di Felice’s equation [81]. The simulations using the Ergun or Di Felice relations will also be considered in our future task. In the present study, the mushy zone constants were analytically determined for the simple morphologies. To justify the use of the simple morphologies, we should visualise the real dendrite structure formed in the L-PBF process. However, the L-PBF process is a rapid process and it is difficult to visualise its 3D dendrite structures. It has been confirmed that 𝐶 is related to pore size, which means that the mushy zone constant, i.e. the solidfication morphology, may be determined by backcalculation if the pore size distribution of the L-PBF experiment can be obtained. Until now, it has been difficult to determine the mushy zone constant on the basis of evidence, but the results of the present study have demonstrated the possibility of estimation from the experimental observations. In this study, the laser processing condition was just under a single condition of 300 W and 500 mm/s. We will perform the simulation under the other conditions and also evaluate the pore dynamics in detail in our future work.This study has focused on the therm-fluid dynamics behaviour during the L-PBF fabrication process. Following the fabrication using L-PBF, post-processing such as surface finishing is often performed, which can alter the material properties. For example, previous studies have demonstrated that surface modfications due to post-processing can affect the material mechanical properties [82--85]. It should be noted that this study does not assess the effect of such post-processing on the material characteristics.5. ConclusionsThis paper investigates the mushy zone constant of the VollerPrakash model, which is used in multiphysics CFD simulations involving phase changes, such as LPBF. A drag force model that includes the morphological parameters of the Kozeny-Carman equation is derived, assuming that dynamic and static pressures are equivalent. The model is characterised by the grid size of the simulation, in addition to the morphological parameters. The model assumes three simple morphologies: rectangular, cone, and cone assembly. Their morphological parameters and the mushy zone constant are estimated analytically. The morphological parameters are compared with those obtained in previous studies. Results in Engineering 24 (2024) 103567 10 J. Katagiri, S. Nomoto, M. Kusano et al. The rectangular morphology is a suitable choice for the proposed drag model. This is because the micrometre-scale grid is commonly used in multiphysics CFD simulations of LPBFs, and the alloy dendrite structure is typically on the scale of a few hundred micrometres to a few millimetres, as observed in past experimental visualisations such as x-ray CT. Furthermore, we have discovered that the mushy zone constant for the rectangular morphology is significantly smaller than the value conventionally used in previous multiphysics CFD simulations.CRediT authorship contribution statementJun Katagiri: Writing -- review & editing, Writing -- original draft, Visualization, Validation, Methodology, Formal analysis, Data curation, Conceptualization. Sukeharu Nomoto: Writing -- review & editing, Validation, Methodology, Formal analysis. Masahiro Kusano: Writing -- review & editing, Validation, Software, Formal analysis. Makoto Watanabe: Writing -- review & editing, Validation, Supervision, Project administration, Funding acquisition, Formal analysis.Declaration of competing interestThe authors declare no coflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.AcknowledgementFinancial support from Innovative Science and Technology Initiative for Security, Grant Number JPJ004596 implemented by the Acquisition, Technology, and Logistics Agency (ATLA).Data availabilityData will be made available on request.References[1] M.S. Manola, B. Singh, M.K. Singla, J.S. Chohan, R. Kumar, R. Kumar, Y.S. Bisht, R. Kumar, M.Q. Alkahtani, S. Islam, M.I. Ammarullah, M.I. Ammarullah, M.I. Ammarullah, Investigation of melt flow index and tensile properties of dual metal reinforced polymer composites for 3d printing using machine learning approach: biomedical and engineering applications, AIP Adv. 14 (5 2024), https://doi.org/10.1063/5.0207551.[2] M.H. Ramlee, M.I. Ammarullah, N.S.M. Sukri, N.S.F. Hassan, M.H. Baharuddin, M.R.A. Kadir, Investigation on three-dimensional printed prosthetics leg sockets coated with different reinforcement materials: analysis on mechanical strength and microstructural, Sci. Rep. 14 (12 2024), https://doi.org/10.1038/s41598-024-57454-8.[3] W.D. Lestari, N. Adyono, A.K. Faizin, A. Haqiyah, K.H. Sanjaya, A. Nugroho, W. Kusmasari, M.I. Ammarullah, Optimization of 3d printed parameters for socket prosthetic manufacturing using the Taguchi method and response surface methodology, Results Eng. 21 (3 2024), https://doi.org/10.1016/j.rineng.2024.101847.[4] S. Nagaraja, P.B. Anand, K.M. Kumar, M.I. Ammarullah, Synergistic advances in natural fibre composites: a comprehensive review of the eco-friendly bio-composite development, its characterization and diverse applications, RSC Adv. 14 (2024) 17594--17611, https://doi.org/10.1039/d4ra00149d.[5] A.U.A. Aziz, M.I. Ammarullah, B.W. Ng, H.S. Gan, M.R.A. Kadir, M.H. Ramlee, Unilateral external fixator and its biomechanical effects in treating different types of femoral fracture: a finite element study with experimental validated model, Heliyon 10 (2 2024), https://doi.org/10.1016/j.heliyon.2024.e26660.[6] R. Seede, D. Shoukr, B. Zhang, A. Whitt, S. Gibbons, P. Flater, A. Elwany, R. Arroyave, I. Karaman, An ultra-high strength martensitic steel fabricated using selective laser melting additive manufacturing: densfication, microstructure, and mechanical properties, Acta Mater. 186 (2020) 199--214, https://doi.org/10.1016/j.actamat.2019.12.037.[7] M.I. Ammarullah, T. Hidayat, M. Danny, P. Lamura, J. Jamari, Relationship between deformation and running-in wear on hard-on-hard bearings from metal, ceramic, and diamond materials for total hip prosthesis, J. Tribol. 38 (2023) 69--81.[8] Z.F.M. Salaha, M.I. Ammarullah, N.N.A.A. Abdullah, A.U.A. Aziz, H.S. Gan, A.H. Abdullah, M.R.A. Kadir, M.H. Ramlee, Biomechanical effects of the porous structure of Gyroid and Voronoi hip implants: a finite element analysis using an experimentally validated model, Materials 16 (5 2023), https://doi.org/10.3390/ma16093298.[9] M.I. Ammarullah, R. Hartono, T. Supriyono, G. Santoso, S. Sugiharto, M.S. Permana, Polycrystalline diamond as a potential material for the hard-on-hard bearing of total hip prosthesis: von Mises stress analysis, Biomedicines 11 (3 2023), https://doi.org/10.3390/biomedicines11030951.[10] G. Santoso, M.I. Ammarullah, S. Sugiharto, R.M. Rachayu, A. Mughni, A.P. Bayuseno, J. Jamari, Von Mises stress analysis of surgery chair designed for laparoscopic surgeon with lifting mechanism, AIP Adv. 14 (4 2024), https://doi.org/10.1063/5.0188663.[11] G. Santoso, M.I. Ammarullah, S. Sugiharto, T. Hidayat, S. Khoeron, A.P. Bayuseno, J. Jamari, Triz-based method for developing a conceptual laparoscopic surgeon’s chair, Cogent Eng. 11 (2024), https://doi.org/10.1080/23311916.2023.2298786.[12] M. Tauviqirrahman, J. Jamari, S. Susilowati, C. Pujiastuti, B. Setiyana, A.H. Pasaribu, M.I. Ammarullah, Performance comparison of Newtonian and non-Newtonian fluid on a heterogeneous slip/no-slip journal bearing system based on cfd-fsi method, Fluids 7 (7 2022), https://doi.org/10.3390/fluids7070225.[13] E.R. Babu, N.C. Reddy, A. Babbar, A. Chandrashekar, R. Kumar, P.S. Bains, M. Alsubih, S. Islam, S.K. Joshi, A. Rizal, M.I. Ammarullah, Characteristics of pulsating heat pipe with variation of tube diameter, filling ratio, and sio2 nanoparticles: biomedical and engineering implications, Case Stud. Therm. Eng. 55 (3 2024), https://doi.org/10.1016/j.csite.2024.104065.[14] A.T. Prakoso, H. Basri, D. Adanta, I. Yani, M.I. Ammarullah, I. Akbar, F.A. Ghazali, A. Syahrom, T. Kamarul, The effect of tortuosity on permeability of porous scaffold, Biomedicines 11 (2 2023), https://doi.org/10.3390/biomedicines11020427.[15] M. Muchammad, M. Tauviqirrahman, M.I. Ammarullah, M. Iqbal, B. Setiyana, J. Jamari, Performance of textured dual mobility total hip prosthesis with a concave dimple during Muslim prayer movements, Sci. Rep. 14 (12 2024), https://doi.org/10.1038/s41598-023-50887-7.[16] M. Arun, D. Barik, S.S. Chandran, N. Govil, P. Sharma, T.M.Y. Khan, R.U. Baig, B.J. Bora, B.J. Medhi, R. Kumar, A. Rizal, M.I. Ammarullah, Twisted helical tape’s impact on heat transfer and friction in zinc oxide (zno) nanofluids for solar water heaters: biomedical insight, Case Stud. Therm. Eng. 56 (4 2024), https://doi.org/10.1016/j.csite.2024.104204.[17] M. Bayat, S. Mohanty, J.H. Hattel, Multiphysics modelling of lack-of-fusion voids formation and evolution in in718 made by multi-track/multi-layer l-pbf, Int. J. Heat Mass Transf. 139 (2019) 95--114, https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.003.[18] L. Wang, Y. Zhang, H.Y. Chia, W. Yan, Mechanism of keyhole pore formation in metal additive manufacturing, npj Comput. Mater. 8 (12 2022), https://doi.org/10.1038/s41524-022-00699-6.[19] D. Grange, A. Queva, G. Guillemot, M. Bellet, J.D. Bartout, C. Colin, Effect of processing parameters during the laser beam melting of inconel 738: comparison between simulated and experimental melt pool shape, J. Mater. Process. Technol. 289 (3 2021), https://doi.org/10.1016/j.jmatprotec.2020.116897.[20] A. Queva, G. Guillemot, C. Moriconi, C. Metton, M. Bellet, Numerical study of the impact of vaporisation on melt pool dynamics in laser powder bed fusion - application to in718 and ti–6al--4v, Addit. Manuf. 35 (10 2020), https://doi.org/10.1016/j.addma.2020.101249.[21] Z. Ren, D.Z. Zhang, G. Fu, J. Jiang, M. Zhao, Hig-fidelity modelling of selective laser melting copper alloy: laser rflection behavior and therma-fluid dynamics, Mater. Des. 207 (9 2021), https://doi.org/10.1016/j.matdes.2021.109857.[22] E.L. Li, L. Wang, A.B. Yu, Z.Y. Zhou, A three-phase model for simulation of heat transfer and melt pool behaviour in laser powder bed fusion process, Powder Technol. 381 (2021) 298--312, https://doi.org/10.1016/j.powtec.2020.11.061.[23] P. Ninpetch, P. Kowitwarangkul, S. Mahathanabodee, P. Chalermkarnnon, P. Rattanadecho, Computational investigation of thermal behavior and Molten metal flow with moving laser heat source for selective laser melting process, Case Stud. Therm. Eng. 24 (4 2021), https://doi.org/10.1016/j.csite.2021.100860.[24] T.F. Flint, J.D. Robson, G. Parivendhan, P. Cardiff, laserbeamfoam: laser ray-tracing and thermally induced state transition simulation toolkit, SoftwareX 21 (2 2023), https://doi.org/10.1016/j.softx.2022.101299.[25] W. Alphonso, M. Baier, S. Carmignato, J. Hattel, M. Bayat, On the possibility of doing reduced order, therm-fluid modelling of laser powder bed fusion (l-pbf) -- assessment of the importance of recoil pressure and surface tension, J. Manuf. Process. 94 (2023) 564--577, https://doi.org/10.1016/j.jmapro.2023.03.040, https://www.sciencedirect.com/science/article/pii/S1526612523002542.[26] M. Carraturo, B. Lane, H. Yeung, S. Kollmannsberger, A. Reali, F. Auricchio, Numerical evaluation of advanced laser control strategies ifluence on residual stresses for laser powder bed fusion systems, Integr. Mater. Manuf. Innov. 9 (2020) 435--445, https://doi.org/10.1007/s40192-020-00191-3.[27] H. Kitano, M. Kusano, M. Tsujii, A. Yumoto, M. Watanabe, Process parameter optimization framework for the selective laser melting of hastelloy x alloy considering defects and solidfication crack occurrence, Crystals 11 (2021), https://doi.org/10.3390/cryst11060578, https://www.mdpi.com/2073-4352/11/6/578.[28] J. Katagiri, M. Kusano, S. Minamoto, H. Kitano, K. Daimaru, M. Tsujii, M. Watanabe, Melt pool shape evaluation by single-track experiments and finite-element thermal analysis: balling and lack of fusion criteria for generating process window of inconel738lc, Materials 16 (2 2023), https://doi.org/10.3390/ma16041729.[29] V.R. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems, Int. J. Heat Mass Transf. 30 (1987) 1709--1719.Results in Engineering 24 (2024) 103567 11 https://doi.org/10.1063/5.0207551https://doi.org/10.1063/5.0207551https://doi.org/10.1038/s41598-024-57454-8https://doi.org/10.1038/s41598-024-57454-8https://doi.org/10.1016/j.rineng.2024.101847https://doi.org/10.1039/d4ra00149dhttps://doi.org/10.1016/j.heliyon.2024.e26660https://doi.org/10.1016/j.actamat.2019.12.037https://doi.org/10.1016/j.actamat.2019.12.037http://refhub.elsevier.com/S2590-1230(24)01810-3/bib18727875CC2F30D28B2B83451D8B6BB8s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib18727875CC2F30D28B2B83451D8B6BB8s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib18727875CC2F30D28B2B83451D8B6BB8s1https://doi.org/10.3390/ma16093298https://doi.org/10.3390/biomedicines11030951https://doi.org/10.3390/biomedicines11030951https://doi.org/10.1063/5.0188663https://doi.org/10.1063/5.0188663https://doi.org/10.1080/23311916.2023.2298786https://doi.org/10.3390/fluids7070225https://doi.org/10.1016/j.csite.2024.104065https://doi.org/10.3390/biomedicines11020427https://doi.org/10.1038/s41598-023-50887-7https://doi.org/10.1038/s41598-023-50887-7https://doi.org/10.1016/j.csite.2024.104204https://doi.org/10.1016/j.csite.2024.104204https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.003https://doi.org/10.1016/j.ijheatmasstransfer.2019.05.003https://doi.org/10.1038/s41524-022-00699-6https://doi.org/10.1038/s41524-022-00699-6https://doi.org/10.1016/j.jmatprotec.2020.116897https://doi.org/10.1016/j.addma.2020.101249https://doi.org/10.1016/j.addma.2020.101249https://doi.org/10.1016/j.matdes.2021.109857https://doi.org/10.1016/j.powtec.2020.11.061https://doi.org/10.1016/j.csite.2021.100860https://doi.org/10.1016/j.softx.2022.101299https://doi.org/10.1016/j.jmapro.2023.03.040https://www.sciencedirect.com/science/article/pii/S1526612523002542https://www.sciencedirect.com/science/article/pii/S1526612523002542https://doi.org/10.1007/s40192-020-00191-3https://doi.org/10.3390/cryst11060578https://doi.org/10.3390/cryst11060578https://www.mdpi.com/2073-4352/11/6/578https://doi.org/10.3390/ma16041729http://refhub.elsevier.com/S2590-1230(24)01810-3/bib4A16906693803A6DA0FBB463BA1A2804s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib4A16906693803A6DA0FBB463BA1A2804s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib4A16906693803A6DA0FBB463BA1A2804s1J. Katagiri, S. Nomoto, M. Kusano et al. [30] P.C. Carman, Fluid flow through granular beds, Trans. Inst. Chem. Eng. 15 (1937) 150--166, URL citeulike-article-id:5101159.[31] P.C. Carman, Permeability of saturated sands, soils and clays, J. Agric. Sci. 29 (1939) 262--273, https://doi.org/10.1017/S0021859600051789.[32] W. Carrier, Goodbye, Hazen; hello, Kozeny-Carman, J. Geotech. Geoenviron. Eng. 129 (2003) 1054--1056, https://doi.org/10.1061/(ASCE)1090-0241(2003)129:11(1054).[33] S. Tiari, S. Qiu, M. Mahdavi, Discharging process of a finned heat pipe-assisted thermal energy storage system with high temperature phase change material, Energy Convers. Manag. 118 (2016) 426--437, https://doi.org/10.1016/j.enconman.2016.04.025.[34] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: numerical investigation and comparison with experiments, Int. J. Heat Mass Transf. 53 (2010) 4082--4091, https://doi.org/10.1016/j.ijheatmasstransfer.2010.05.028.[35] L. Yu, G.L. Ding, J. Reye, S.N. Ojha, S.N. Tewari, G. Ding, S. Ojha, S. Tewari, Mushy zone morphology during directional solidfication of pb-5.8 wt pct sb alloy, Metall. Mater. Trans. A 31 (2000) 2275--2285.[36] M.A. Azeem, P.D. Lee, A.B. Phillion, S. Karagadde, P. Rockett, R.C. Atwood, L. Courtois, K.M. Rahman, D. Dye, Revealing dendritic pattern formation in ni, fe and co alloys using synchrotron tomography, Acta Mater. 128 (2017) 241--248, https://doi.org/10.1016/j.actamat.2017.02.022.[37] R.H. Mathiesen, L. Arnberg, X-ray radiography observations of columnar dendritic growth and constitutional undercooling in an Al-30wt%Cu, Acta Mater. 53 (2005) 947--956, https://doi.org/10.1016/j.actamat.2004.10.050.[38] H. Yasuda, K. Morishita, N. Nakatsuka, T. Nishimura, M. Yoshiya, A. Sugiyama, K. Uesugi, A. Takeuchi, Dendrite fragmentation induced by massive-like 𝛿--𝛾 transformation in fe–c alloys, Nat. Commun. 10 (12 2019), https://doi.org/10.1038/s41467-019-11079-y.[39] G. Reinhart, D. Grange, L. Abou-Khalil, N. Mangelinck-Noël, N.T. Niane, V. Maguin, G. Guillemot, C.A. Gandin, H. Nguyen-Thi, Impact of solute flow during directional solidfication of a ni-based alloy: in-situ and real-time x-radiography, Acta Mater. 194 (2020) 68--79, https://doi.org/10.1016/j.actamat.2020.04.003.[40] S. Karagadde, C.L.A. Leung, P.D. Lee, Progress on in situ and operando X-ray imaging of solidfication processes, Materials 14 (5 2021), https://doi.org/10.3390/ma14092374.[41] M.L.N.M. Melo, E.M.S. Rizzo, R.G. Santos, Predicting dendrite arm spacing and their effect on microporosity formation in directionally solidfied al-cu alloy, J. Mater. Sci. 40 (2005) 1599--1609.[42] Y.M. Wang, T. Voisin, J.T. McKeown, J. Ye, N.P. Calta, Z. Li, Z. Zeng, Y. Zhang, W. Chen, T.T. Roehling, R.T. Ott, M.K. Santala, P.J. Depond, M.J. Matthews, A.V. Hamza, T. Zhu, Additively manufactured hierarchical stainless steels with high strength and ductility, Nat. Mater. 17 (2018) 63--70, https://doi.org/10.1038/NMAT5021.[43] M. Yang, L. Wang, W. Yan, Phas-field modeling of grain evolutions in additive manufacturing from nucleation, growth, to coarsening, npj Comput. Mater. 7 (12 2021), https://doi.org/10.1038/s41524-021-00524-6.[44] A. Trykozko, M. Peszynska, M. Dohnalik, Modeling non-Darcy flows in realistic pore-scale proppant geometries, Comput. Geotech. 71 (2016) 352--360, https://doi.org/10.1016/j.compgeo.2015.08.011, http://www.sciencedirect.com/science/article/pii/S0266352X15001883.[45] R. Sivanesapillai, H. Steeb, A. Hartmaier, Transition of effective hydraulic properties from low to high Reynolds number flow in porous media, Geophys. Res. Lett. 41 (2014) 4920--4928, https://doi.org/10.1002/2014GL060232.[46] G.A. Narsilio, O. Buzzi, S. Fityus, T.S. Yun, D.W. Smith, Upscaling of Navier–Stokes equations in porous media: theoretical, numerical and experimental approach, Comput. Geotech. 36 (2009) 1200--1206, https://doi.org/10.1016/j.compgeo.2009.05.006, http://www.sciencedirect.com/science/article/pii/S0266352X0900086X.[47] J. Katagiri, M. Kusano, S. Nomoto, M. Watanabe, Ifluence of recoil pressure, mushy zone flow resistance and rflectivity on melt pool shape in laser powder bed fusion simulation, Case Stud. Therm. Eng. 50 (2023) 103477, https://doi.org/10.1016/j.csite.2023.103477, https://www.sciencedirect.com/science/article/pii/S2214157X23007839.[48] M.I. Yuce, D.A. Kareem, A numerical analysis of fluid flow around circular and square cylinders, J. Am. Water Works Assoc. 108 (10) (2016) E546--E554, https://doi.org/10.5942/jawwa.2016.108.0141, https://awwa.onlinelibrary.wiley.com/doi/pdf/10.5942/jawwa.2016.108.0141.[49] S.F. Hoerner, Fluid-Dynamic Drag, 2nd edition, 1965, Midland Park, N.J., first edition published under the title: Aerodynamic Drag.[50] J. Katagiri, S. Kimura, S. Noda, Significance of shape factor on permeability anisotropy of sand: representative elementary volume study for pore-scale analysis, Acta Geotech. 15 (2020) 2195--2203, https://doi.org/10.1007/s11440-020-00912-0.[51] J. Katagiri, Y. Konno, J. Yoneda, N. Tenma, Pore-scale modeling of flow in particle packs containing grain-coating and por-filling hydrates: verfication of a Kozeny–Carman-based permeability reduction model, J. Nat. Gas Sci. Eng. 45 (2017) 537--551, https://doi.org/10.1016/j.jngse.2017.06.019.[52] E. Li, Z. Zhou, L. Wang, H. Shen, R. Zou, A. Yu, Particle scale modelling of melt pool dynamics and pore formation in selective laser melting additive manufacturing, Powder Technol. 397 (2022) 117012, https://doi.org/10.1016/j.powtec.2021.11.056, https://www.sciencedirect.com/science/article/pii/S0032591021010111.[53] P. Mostaghimi, M.J. Blunt, B. Bijeljic, Computations of absolute permeability on micro-ct images, Math. Geosci. 45 (2013) 103--125, https://doi.org/10.1007/s11004-012-9431-4.[54] J.W. Gibbs, K.A. Mohan, E.B. Gulsoy, A.J. Shahani, X. Xiao, C.A. Bouman, M.D. Graef, P.W. Voorhees, The three-dimensional morphology of growing dendrites, Sci. Rep. 5 (7 2015), https://doi.org/10.1038/srep11824.[55] V. Yfit, F. Tariq, D.S. Eastwood, M. Biton, B. Wu, P.D. Lee, N.P. Brandon, Operando visualization and multi-scale tomography studies of dendrite formation and dissolution in zinc batteries, Joule 3 (2019) 485--502, https://doi.org/10.1016/j.joule.2018.11.002.[56] E. Khajeh, D.M. Maijer, Numerical determination of permeability of al-cu alloys using 3d geometry from X-ray microtomography, Mater. Sci. Technol. 26 (2010) 1469--1476, https://doi.org/10.1179/174328409X411718.[57] J. Madison, J.E. Spowart, D.J. Rowenhorst, L.K. Aagesen, K. Thornton, T.M. Pollock, Fluid flow and defect formation in the three-dimensional dendritic structure of nickel-based single crystals, Metall. Mater. Trans. A, Phys. Metall. Mater. Sci. 43 (2012) 369--380, https://doi.org/10.1007/s11661-011-0823-8.[58] J. Katagiri, S. Nomoto, M. Kusano, M. Watanabe, Particle size effect on powder packing properties and Molten pool dimensions in laser powder bed fusion simulation, J. Manuf. Mater. Process. 8 (4 2024), https://doi.org/10.3390/jmmp8020071.[59] M.I. Ammarullah, G. Santoso, S. Sugiharto, T. Supriyono, O. Kurdi, M. Tauviqirrahman, T.I. Winarni, J. Jamari, Tresca stress study of cocrmo-on-cocrmo bearings based on body mass index using 2d computational model, J. Tribol. 33 (2022) 31--38.[60] T. Hidayat, M.I. Ammarullah, E. Saputra, M.D.P. Lamura, C.N. Chethan, R. Ismail, A.P. Bayuseno, J. Jamari, A method for estimating the contact area of a dual-mobility total hip prosthesis, AIP Adv. 14 (1 2024), https://doi.org/10.1063/5.0188638.[61] T. Hidayat, R. Ismail, M. Tauviqirrahman, E. Saputra, M.I. Ammarullah, M.D.P. Lamura, A.P. Bayuseno, Jamari, Running-in behavior of dual-mobility cup during the gait cycle: a finite element analysis, Proc. Inst. Mech. Eng., H J. Eng. Med. 238 (2024) 99--111, https://doi.org/10.1177/09544119231216023, pMID: 38156402.[62] M. Tauviqirrahman, M.I. Ammarullah, J. Jamari, E. Saputra, T.I. Winarni, F.D. Kurniawan, S.A. Shiddiq, E. van der Heide, Analysis of contact pressure in a 3d model of dual-mobility hip joint prosthesis under a gait cycle, Sci. Rep. 13 (12 2023), https://doi.org/10.1038/s41598-023-30725-6.[63] H. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Pearson Education Limited, 2007, https://books.google.co.jp/books?id=RvBZ-UMpGzIC.[64] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys. 12 (6) (1998) 620--631, https://doi.org/10.1063/1.168744, https://pubs.aip.org/aip/cip/article-pdf/12/6/620/7865493/620_1_online.pdf.[65] J.H. Ferziger, M. Perić, Computational Methods for Fluid Dynamics, Springer-Verlag, Berlin, 1996.[66] T. Marić, J. Höpken, K.G. Mooney, The OpenFOAM Technology Primer, Zenodo, 2021.[67] R.U. Putra, H. Basri, A.T. Prakoso, H. Chandra, M.I. Ammarullah, I. Akbar, A. Syahrom, T. Kamarul, Level of activity changes increases the fatigue life of the porous magnesium scaffold, as observed in dynamic immersion tests, over time, Sustainability (Switzerland) 15 (1 2023), https://doi.org/10.3390/su15010823.[68] J. Jamari, M.I. Ammarullah, G. Santoso, S. Sugiharto, T. Supriyono, M.S. Permana, T.I. Winarni, E. van der Heide, Adopted walking condition for computational simulation approach on bearing of hip joint prosthesis: review over the past 30 years, Heliyon 8 (12 2022), https://doi.org/10.1016/j.heliyon.2022.e12050.[69] M.I. Ammarullah, G. Santoso, S. Sugiharto, T. Supriyono, D.B. Wibowo, O. Kurdi, M. Tauviqirrahman, J. Jamari, Minimizing risk of failure from ceramic-on-ceramic total hip prosthesis by selecting ceramic materials based on Tresca stress, Sustainability (Switzerland) 14 (10 2022), https://doi.org/10.3390/su142013413.[70] J. Jamari, M.I. Ammarullah, G. Santoso, S. Sugiharto, T. Supriyono, E. van der Heide, In silico contact pressure of metal-on-metal total hip implant with different materials subjected to gait loading, Metals 12 (8 2022), https://doi.org/10.3390/met12081241.[71] M.I. Ammarullah, I.Y. Afif, M.I. Maula, T.I. Winarni, M. Tauviqirrahman, I. Akbar, H. Basri, E. van der Heide, J. Jamari, Tresca stress simulation of metal-on-metal total hip arthroplasty during normal walking activity, Materials 14 (12 2021), https://doi.org/10.3390/ma14247554.[72] T. Hidayat, R. Ismail, M. Tauviqirrahman, E. Saputra, M.I. Ammarullah, M. Danny, P. Lamura, A.P. Bayuseno, J. Jamari, Investigation of mesh model for a finite element simulation of the dual-mobility prosthetic hip joint keywords abstract mesh model element size fea contact pressure von Mises artficial hip joint, J. Tribol. 38 (2023) 118--1401.[73] M.D.P. Lamura, T. Hidayat, M.I. Ammarullah, A.P. Bayuseno, J. Jamari, Study of contact mechanics between two brass solids in various diameter ratios and friction coefficient, Proc. Inst. Mech. Eng., Part J J. Eng. Tribol. 237 (2023) 1613--1619, https://doi.org/10.1177/14657503221144810.[74] M.D.P. Lamura, M.I. Ammarullah, T. Hidayat, M.I. Maula, J. Jamari, A.P. Bayuseno, Diameter ratio and friction coefficient effect on equivalent plastic strain (peeq) during contact between two brass solids, Cogent Eng. 10 (2023), https://doi.org/10.1080/23311916.2023.2218691.[75] M.D.P. Lamura, M.I. Ammarullah, M.I. Maula, T. Hidayat, A.P. Bayuseno, J. Jamari, The effect of load, diameter ratio, and friction coefficient on residual stress in a Results in Engineering 24 (2024) 103567 12 http://refhub.elsevier.com/S2590-1230(24)01810-3/bib5F761ED13B484DF49D6D76847938AC00s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib5F761ED13B484DF49D6D76847938AC00s1https://doi.org/10.1017/S0021859600051789https://doi.org/10.1061/(ASCE)1090-0241(2003)129:11(1054)https://doi.org/10.1061/(ASCE)1090-0241(2003)129:11(1054)https://doi.org/10.1016/j.enconman.2016.04.025https://doi.org/10.1016/j.enconman.2016.04.025https://doi.org/10.1016/j.ijheatmasstransfer.2010.05.028http://refhub.elsevier.com/S2590-1230(24)01810-3/bib89D46E8A9726FB3ABA7E5985BF9591EEs1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib89D46E8A9726FB3ABA7E5985BF9591EEs1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib89D46E8A9726FB3ABA7E5985BF9591EEs1https://doi.org/10.1016/j.actamat.2017.02.022https://doi.org/10.1016/j.actamat.2017.02.022https://doi.org/10.1016/j.actamat.2004.10.050https://doi.org/10.1038/s41467-019-11079-yhttps://doi.org/10.1038/s41467-019-11079-yhttps://doi.org/10.1016/j.actamat.2020.04.003https://doi.org/10.3390/ma14092374https://doi.org/10.3390/ma14092374http://refhub.elsevier.com/S2590-1230(24)01810-3/bibBFD7767A72AA33EF3B2285F6510BDE83s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bibBFD7767A72AA33EF3B2285F6510BDE83s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bibBFD7767A72AA33EF3B2285F6510BDE83s1https://doi.org/10.1038/NMAT5021https://doi.org/10.1038/NMAT5021https://doi.org/10.1038/s41524-021-00524-6https://doi.org/10.1016/j.compgeo.2015.08.011https://doi.org/10.1016/j.compgeo.2015.08.011http://www.sciencedirect.com/science/article/pii/S0266352X15001883http://www.sciencedirect.com/science/article/pii/S0266352X15001883https://doi.org/10.1002/2014GL060232https://doi.org/10.1016/j.compgeo.2009.05.006https://doi.org/10.1016/j.compgeo.2009.05.006http://www.sciencedirect.com/science/article/pii/S0266352X0900086Xhttps://doi.org/10.1016/j.csite.2023.103477https://doi.org/10.1016/j.csite.2023.103477https://www.sciencedirect.com/science/article/pii/S2214157X23007839https://www.sciencedirect.com/science/article/pii/S2214157X23007839https://doi.org/10.5942/jawwa.2016.108.0141https://awwa.onlinelibrary.wiley.com/doi/pdf/10.5942/jawwa.2016.108.0141https://awwa.onlinelibrary.wiley.com/doi/pdf/10.5942/jawwa.2016.108.0141http://refhub.elsevier.com/S2590-1230(24)01810-3/bib310E71BA1872F828303A0F0101C10F1Cs1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib310E71BA1872F828303A0F0101C10F1Cs1https://doi.org/10.1007/s11440-020-00912-0https://doi.org/10.1007/s11440-020-00912-0https://doi.org/10.1016/j.jngse.2017.06.019https://doi.org/10.1016/j.powtec.2021.11.056https://doi.org/10.1016/j.powtec.2021.11.056https://www.sciencedirect.com/science/article/pii/S0032591021010111https://doi.org/10.1007/s11004-012-9431-4https://doi.org/10.1007/s11004-012-9431-4https://doi.org/10.1038/srep11824https://doi.org/10.1016/j.joule.2018.11.002https://doi.org/10.1016/j.joule.2018.11.002https://doi.org/10.1179/174328409X411718https://doi.org/10.1007/s11661-011-0823-8https://doi.org/10.3390/jmmp8020071http://refhub.elsevier.com/S2590-1230(24)01810-3/bibF297E187BC1CD9F3D727302C8A066063s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bibF297E187BC1CD9F3D727302C8A066063s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bibF297E187BC1CD9F3D727302C8A066063s1https://doi.org/10.1063/5.0188638https://doi.org/10.1177/09544119231216023https://doi.org/10.1038/s41598-023-30725-6https://doi.org/10.1038/s41598-023-30725-6https://books.google.co.jp/books?id=RvBZ-UMpGzIChttps://books.google.co.jp/books?id=RvBZ-UMpGzIChttps://doi.org/10.1063/1.168744https://pubs.aip.org/aip/cip/article-pdf/12/6/620/7865493/620_1_online.pdfhttps://pubs.aip.org/aip/cip/article-pdf/12/6/620/7865493/620_1_online.pdfhttp://refhub.elsevier.com/S2590-1230(24)01810-3/bib010ED95BB57CF759C72328745801D598s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib010ED95BB57CF759C72328745801D598s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib9FBFB5ED22D3D737784FB352DEABE0BBs1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib9FBFB5ED22D3D737784FB352DEABE0BBs1https://doi.org/10.3390/su15010823https://doi.org/10.1016/j.heliyon.2022.e12050https://doi.org/10.3390/su142013413https://doi.org/10.3390/met12081241https://doi.org/10.3390/met12081241https://doi.org/10.3390/ma14247554https://doi.org/10.3390/ma14247554http://refhub.elsevier.com/S2590-1230(24)01810-3/bib41149A8B011EE8086F16FCB0505F5B33s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib41149A8B011EE8086F16FCB0505F5B33s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib41149A8B011EE8086F16FCB0505F5B33s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib41149A8B011EE8086F16FCB0505F5B33s1http://refhub.elsevier.com/S2590-1230(24)01810-3/bib41149A8B011EE8086F16FCB0505F5B33s1https://doi.org/10.1177/14657503221144810https://doi.org/10.1080/23311916.2023.2218691https://doi.org/10.1080/23311916.2023.2218691J. Katagiri, S. Nomoto, M. Kusano et al. hemispherical contact for application in biomedical industry, J. Mater. Eng. Perform. (2024), https://doi.org/10.1007/s11665-024-09330-9.[76] T. Hidayat, M.I. Ammarullah, R. Ismail, E. Saputra, M.D.P. Lamura, K.N. Chethan, A.P. Bayuseno, J. Jamari, Investigation of contact behavior on a model of the dualmobility artficial hip joint for asians in different inner liner thicknesses, World J. Orthop. 15 (2024) 321--336, https://doi.org/10.5312/wjo.v15.i4.321.[77] Y. Cai, H. Heng, F. Li, M. Wang, The ifluences of welding parameters on the metal vapor plume in fiber laser welding based on 3d reconstruction, Opt. Laser Technol. 107 (2018) 1--7, https://doi.org/10.1016/j.optlastec.2018.05.016.[78] C. Zhao, N.D. Parab, X. Li, K. Fezzaa, W. Tan, A.D. Rollett, T. Sun, Critical instability at moving keyhole tip generates porosity in laser melting, Science 370 (2020) 1080--1086, https://doi.org/10.1126/science.abd1587.[79] Y. Huang, T.G. Fleming, S.J. Clark, S. Marussi, K. Fezzaa, J. Thiyagalingam, C.L.A. Leung, P.D. Lee, Keyhole fluctuation and pore formation mechanisms during laser powder bed fusion additive manufacturing, Nat. Commun. 13 (2022) 1170, https://doi.org/10.1038/s41467-022-28694-x.[80] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89--94, https://doi.org/10.1029/JB088iS01p0B353.[81] R.D. Felice, The voidage function for fluid-particle interaction systems, Int. J. Multiph. Flow 20 (1994) 153--159, https://doi.org/10.1016/0301-9322(94)90011-6, http://www.sciencedirect.com/science/article/pii/0301932294900116.[82] R. Kodandappa, S. Nagaraja, M.M. Chowdappa, M. Krishnappa, G.S. Poornima, M.I. Ammarullah, Application of the Taguchi method and rsm for process parameter optimization in awsj machining of cfrp composite-based orthopedic implants, Open Eng. 14 (1 2024), https://doi.org/10.1515/eng-2024-0057.[83] M.U. Farooq, S. Anwar, H.A. Bhatti, M.S. Kumar, M.A. Ali, M.I. Ammarullah, Electric discharge machining of ti6al4v eli in biomedical industry: parametric analysis of surface functionalization and tribological characterization, Materials 16 (6 2023), https://doi.org/10.3390/ma16124458.[84] K. Mughal, M.P. Mughal, M.U. Farooq, S. Anwar, M.I. Ammarullah, Using nanofluids minimum quantity lubrication (nf-mql) to improve tool wear characteristics for efficient machining of cfrp/ti6al4v aeronautical structural composite, Processes 11 (5 2023), https://doi.org/10.3390/pr11051540.[85] S. Roseno, M.I. Ammarullah, S. Rohman, F. Kurniawati, T. Wahyudi, A.H.S. Wargadipura, M. Masmui, D. Budiyanto, M.D. Effendi, W. Wahyudin, E. Kalembang, H. Hernawan, S. Subari, S. Habibie, T.P.H. Simanjuntak, H. Santoso, A. Ahmad, A.L. Juwono, The effects of carbon fiber surface treatment by oxidation process for enhanced mechanical properties of carbon fiber/epoxy composites for biomedical application, AIP Adv. 14 (1 2024), https://doi.org/10.1063/5.0183153.Results in Engineering 24 (2024) 103567 13 https://doi.org/10.1007/s11665-024-09330-9https://doi.org/10.5312/wjo.v15.i4.321https://doi.org/10.1016/j.optlastec.2018.05.016https://doi.org/10.1126/science.abd1587https://doi.org/10.1038/s41467-022-28694-xhttps://doi.org/10.1038/s41467-022-28694-xhttps://doi.org/10.1029/JB088iS01p0B353https://doi.org/10.1016/0301-9322(94)90011-6http://www.sciencedirect.com/science/article/pii/0301932294900116https://doi.org/10.1515/eng-2024-0057https://doi.org/10.3390/ma16124458https://doi.org/10.3390/pr11051540https://doi.org/10.1063/5.0183153 Effect of mushy zone constant in Voller-Prakash model on keyholing behaviour in laser powder bed fusion simulation 1 Introduction 2 Kozeny-Carman-based mushy zone flow resistance model 2.1 Drag force-liquid fraction relation 2.2 Rectangular morphology 2.3 Conical morphology 2.4 Cone assembly 3 Results and discussions 3.1 Comparison of drag force 3.2 Solidification morphology effect 3.3 Revisit on mushy zone constant 4 Influence of mushy zone constant on keyholing behaviour 4.1 Simulation overview 4.1.1 Discrete-element modelling of powder spreading 4.1.2 Governing equation of multiphysics CFD simulation 4.2 Results and discussions 5 Conclusions CRediT authorship contribution statement Declaration of competing interest Acknowledgement Data availability References