# Fileset

[jmmp-09-00383.pdf](https://mdr.nims.go.jp/filesets/4896e5b4-6697-45eb-9a6e-44bc6fa19911/download)

## Creator

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

## Rights

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

## Other metadata

[Bayesian Optimization-Based Parameter Identification for Discrete Element Method Simulation of Consolidation and Its Application to Powder Spreading Analysis](https://mdr.nims.go.jp/datasets/8159e097-c2fe-48d8-b210-c0b2bed296dc)

## Fulltext

Bayesian Optimization-Based Parameter Identification for Discrete Element Method Simulation of Consolidation and Its Application to Powder Spreading AnalysisAcademic Editor: Shuting LeiReceived: 24 October 2025Revised: 18 November 2025Accepted: 20 November 2025Published: 21 November 2025Citation: Katagiri, J.; Kusano, M.;Watanabe, M. BayesianOptimization-Based ParameterIdentification for Discrete ElementMethod Simulation of Consolidationand Its Application to PowderSpreading Analysis. J. Manuf. Mater.Process. 2025, 9, 383. https://doi.org/10.3390/jmmp9120383Copyright: © 2025 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license(https://creativecommons.org/licenses/by/4.0/).ArticleBayesian Optimization-Based Parameter Identification forDiscrete Element Method Simulation of Consolidation and ItsApplication to Powder Spreading AnalysisJun Katagiri 1,* , Masahiro Kusano 2 and Makoto Watanabe 21 Degradation Diagnostics Research Team, Integrated Research Center for Resilient Infrastructure, NationalInstitute of Advanced Industrial Science and Technology, 1-1-1 Higashi Tsukuba, Ibaraki 305-8567, Japan2 Additive Manufacturing Group, Research Center for Structural Materials, National Institute for MaterialsScience, 1-2-1 Sengen, Tsukuba 305-0047, Japan; watanabe.makoto@nims.go.jp (M.W.)* Correspondence: j-katagiri@aist.go.jpAbstractThis study develops a Bayesian optimization framework to calibrate two discrete elementmethod (DEM) parameters—the cohesion-related surface energy coefficient (k) and therolling resistance coefficient (µr)—based on experimental void ratio data obtained frompowder consolidation tests. The optimized parameter set reproduces the void ratio ob-tained from the consolidation experiment, demonstrating efficient and physically plausiblecalibration under confined loading. However, when these parameters are applied to pow-der spreading simulations, the resulting powder beds become excessively cohesive, leadingto poor layer uniformity. This discrepancy is attributed to (i) the mismatch in mechanicalscaling (σ/E) between the experimental and simulated conditions and (ii) the shift indominant particle-scale mechanics from normal-force-controlled consolidation to shear-dominated spreading. The results indicate that DEM parameter calibration for powderbed-based additive manufacturing should incorporate shear-related experimental metricsand scaling considerations rather than rely solely on consolidation-based fitting.Keywords: powder bed fusion laser melting (PBF-LM); Bayesian optimization; discreteelement method (DEM); parameter identification; consolidation; powder spreading1. IntroductionThe deterioration of transportation infrastructure such as bridges has become a crit-ical issue in many countries worldwide, including Japan. Repairing such infrastructurerequires restoration measures that are commensurate with the extent and nature of thedegradation. In recent years, additive manufacturing (AM) technologies have attractedsignificant attention as innovative approaches for the maintenance and repair of bridges.This study focuses on laser powder bed fusion (PBF-LM), one of the AM processes. Inthe PBF-LM process, thin layers of metal powder are sequentially spread and selectivelymelted and solidified by a moving laser beam, thereby enabling the fabrication of arbi-trarily shaped three-dimensional objects. In such manufacturing processes, it is naturallyanticipated that both the condition of the spread powder layer, often referred to as thepowder bed, and the laser processing parameters, including laser scanning speed and laserpower, significantly influence the quality of the fabricated parts.Conventionally, laser technology has been extensively investigated in the field ofmetal laser welding, and correspondingly, numerous studies have been conducted onJ. Manuf. Mater. Process. 2025, 9, 383 https://doi.org/10.3390/jmmp9120383https://doi.org/10.3390/jmmp9120383https://doi.org/10.3390/jmmp9120383https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/https://www.mdpi.com/journal/jmmphttps://www.mdpi.comhttps://orcid.org/0000-0002-6399-1951https://orcid.org/0000-0002-5061-0195https://orcid.org/0000-0002-5064-9583https://doi.org/10.3390/jmmp9120383https://www.mdpi.com/article/10.3390/jmmp9120383?type=check_update&version=1J. Manuf. Mater. Process. 2025, 9, 383 2 of 16the influence of the laser processing parameters on build quality in PBF-LM. It has beenreported that the dimensions of the molten pool formed by the laser—specifically its depth,width, length, and their ratios—are correlated with defects occurring during the fabricationprocess [1,2]. Consequently, research efforts have been directed towards identifying theprocess window, i.e., the range of laser processing parameters under which various defectsdo not occur [2–6].In contrast, comparatively fewer investigations, such as [7], have addressed the effectsof the powder bed condition on build quality relative to those concerning laser parameters.One feasible reason for this is that the field of powder engineering has been historicallyoriented towards the management of large quantities of powder as bulk material. Thissuggests that applications involving the thin and widespread distribution of powder havehistorically been less prevalent. In view of this, we aim to evaluate the influence of powdercharacteristics such as particle size, morphology, and adhesion on the powder bed quality.Metallic powders with an average particle diameter of several tens of mm have beenoften employed in PBF-LM manufacturing. At the micrometer scale, adhesion forcesarising from electrostatic interactions, van der Waals forces, and humidity become non-negligible compared to the gravitational force acting on the powder particles [8]. Toaccurately simulate the behavior of the metal powders, it is essential to incorporate adhesioneffects within the particle contact force models of the discrete element method (DEM) [9].Furthermore, it is imperative to ascertain the model parameters that govern the adhesionforce in order to accurately reproduce experimental observations.Conventionally, the process of parameter selection has involved multiple simulationruns with varied parameters, followed by a comparison with experiments. For example,past studies have calibrated DEM parameters to match the experimental angle of repose(e.g., [10–12]), while Lupo et al. calibrated the rolling friction coefficient against powderspreading measurements of packing density [13]. In past studies, parameters to be iden-tified are chosen a priori and identified via manual parametric exploration and expertassessment. However, such a manual procedure frequently entails substantial effort andtime for parameter identification. From a statistical viewpoint, the parameters may besystematically varied and simulations conducted in an exhaustive manner. Nevertheless,DEM simulations tend to be computationally expensive, and consequently, it is unfeasibleto run a large number of cases. The objective of this study is to utilize Bayesian optimizationto achieve parameter identification with a minimal number of simulation runs.The paper is organized as follows. The second section provides a comprehensiveoverview of the materials utilized in the study, along with a detailed description of theexperimental procedures and a thorough exposition of the DEM analysis methodology,and the parameter identification approach, which combines Bayesian optimization withDEM simulations. In Section 3, the identified parameter values are discussed. Also, weconducted the powder spreading simulations using the identified parameters, comparedthe results with previous studies, and verified the validity of the identified parameters.2. Materials and MethodsIn this study, the squared error between the void ratio obtained from the consolidationexperiment and that derived from the DEM simulation of consolidation was employed asthe objective function. The optimal parameter set of k and µr was determined by minimizingthis objective function through Bayesian optimization. The identified parameters weresubsequently applied to powder spreading simulations to examine their transferability.A schematic illustration of the entire procedure is presented in Figure 1, and detaileddescriptions of each step are provided in the following subsections.J. Manuf. Mater. Process. 2025, 9, 383 3 of 16Figure 1. A schematic of the flow of the consolidation experiment, the DEM simulation of consolida-tion, Bayesian optimization, and its application to powder spreading simulation.2.1. Materials UsedWe used gas-atomized Hastelloy X (HX) powder (AMPERPRINT 0228.074, HöganäsAB, Höganäs, Skåne Län, Sweden) with a density of 8220 kg/m3. The particle morphol-ogy is nearly spherical. The particle diameters were mostly distributed in the range of10–40 mm, with an average particle diameter of approximately 29 µm [14].2.2. Consolidation ExperimentA consolidation experiment was conducted using an FT4 powder rheometer (FreemanTechnology Ltd., Tewkesbury, Gloucestershire, UK) [15–17]. The cylindrical vessel wasinitially filled with the powder via a funnel. Subsequently, the powder was compressed bythe top wall under a constant applied stress of 9 kPa, as shown in Figure 2. Subsequent tothe process of consolidation, the void ratio (e), defined as the ratio of the post-consolidationvoid volume to the total volume of the particle assembly, was calculated. Smaller values ofe indicate a denser packing state.Figure 2. A schematic of the consolidation experiment.J. Manuf. Mater. Process. 2025, 9, 383 4 of 162.3. Discrete Element Method OverviewThe open-source DEM software LIGGGHTS-PUBLIC (ver. 3.8.0) [18] was employedin this study. The simulation parameters are shown in Table 1. In order to reproducethe cohesion of the powders, a simplified Johnson–Kendall–Roberts (SJKR)-type adhesionmodel (the sjkr2 model) and a rolling resistance model (cdt model) were also adopted [18].The JKR model and the cdt model are expressed by Equations (1) and (2), respectively.F = 4kπδnrirjri+rj, (1)T rf = µr kn δnωr,shear∣∣∣∣ωr,shear∣∣∣∣ rirjri+rj, (2)Table 1. The DEM model parameters and parameters to be set for the simulation.Parameter ValueContact law Hertz–Mindlin (no history)Density 8220 kg/m3Young’s modulus 5 × 106 PaPoisson’s ratio 0.3Restitution coefficient 0.4Particle–particle friction coefficient 0.3Particle–wall friction coefficient 0.3Time increment 1 × 10−6 sGravity 9.81 m/s2Neighbor list radius 40 µmIntegrator nve/sphereBoundary condition x: periodic, y: periodic, z: fixedWall consolidation stress 0.9 kPaServo wall maximum velocity 0.1 m/sWall servo control parameters (kp) 0.1k is the surface energy density. ri and rj are the radii of i-th and j-th particles, respectively. dn is the overlap depthbetween i-th and j-th particles. mr, kn, and wr ,shear are the rolling friction coefficient, normal spring constant, andcontact tangential component of the relative rotation vector between i-th and j-th particles, respectively. The JKRmodel requires the input of parameter k, whilst the rotational resistance model requires mr. In the case of a hybridmodel of the two models, it is necessary to set two model parameters.2.4. Consolidation SimulationA simulation was conducted to reproduce a certain volume of the element subjectedto the consolidation test. Figure 3 shows an example snapshot of the final state of theconsolidation simulation. A computational domain of 2 mm in the x and y directions wasdefined, with periodic boundary conditions. The initial height was set to 5 mm, and rigidwalls were placed at z = 0 mm and z = 5 mm under fixed boundary conditions. In total,500,000 particles inserted from the upper part of the computational domain were settledunder gravity. Subsequently, the powder layer was compressed by the servo controllingthe upper rigid wall, and then we calculated e.For the HX considered in this study, Young’s modulus is about 200 GPa [19]. In otherwords, a consolidation pressure of 9 kPa is maintained by a particle assembly with a stiffnessof approximately 2 × 1011 Pa. However, employing a stiffness of 2 × 1011 Pa in DEMsimulations results in an extremely small time step, leading to unfeasible computation times.Considering the above discussion, we employed an empirically determined consolidationpressure of 0.9 kPa.J. Manuf. Mater. Process. 2025, 9, 383 5 of 16Figure 3. An example snapshot of the consolidation DEM simulation.2.5. Bayesian OptimizationPython library GpyOpt 1.2.6 (the required libraries are Gpy 1.13.2, NumPy 1.26.4,Python 3.9) [20] was utilized. Upper and lower bounds were specified for two parameters,and preliminary simulations were conducted for four cases corresponding to the combina-tions of these bounds, as well as one case using the midpoint values of both parameters,resulting in a total of five initial cases. It should be noted that the optimization was initiatedwith five initial samples. Although the use of prior information is not strictly necessary, ithas been empirically demonstrated that the incorporation of such information generallyincreases the possibility of identifying an optimal solution. To distribute data points widelyin the design space, researchers typically employ grid search or Latin hypercube sampling(LHS). Given the potential for these strategies to inflate the number of simulations, usinga coarse-level grid (i.e., a few levels per variable) is frequently a reasonable compromise.The parameter search was then executed using the results of the five initial cases as priorinformation. The parameter values in Bayesian optimization are summarized in Table 2.The acquisition weight (wa), often referred to as jitter, governs the extent to which priorinformation is considered, with larger values tending to explore regions further from theprior points and smaller values tending to search in the vicinity of the prior points. Thevalue of wa was manually adjusted for each optimization run. Then, the workflow shownin Figure 1 was conducted to identify the parameter set.It is important to note that the experimental void ratio target was obtained froma single trial. As anticipated, experimental measurements show random variability. InGPyOpt, measurement noise can be accommodated by setting “exact_feval = False”, whichallows the Gaussian process likelihood to include a noise term. When multiple experimentalreplicates are available, one may supply the replicate observations directly or use theiraverage together with an estimate of the noise variance (e.g., from the sample standarddeviation) so that the optimizer treats the objective as noisy and robustness can be assessed.J. Manuf. Mater. Process. 2025, 9, 383 6 of 16Table 2. Bayesian optimization parameters used.Parameter ValueAcquisition(acquisition_type)Expected improvement(EI)Model type(model_type)Gaussian process(GP)Allowance of duplicate proposals(de_duplication) TrueInclusion of noise in objective function(exact_feval) FalseAcquisition weight: wa(acquisition_weight) 0.01–1.0Objective function (eexp − eDEM)2Optimization type MinimizationInitial cases(initial_design_num_data, X, Y) 5 (0–1 to 0–5 in Table 3)Range of k 102–105 J/m3Range of µr 0.05–0.7Number of iterations 100Table 3. Results of preliminary analysis and Bayesian optimization.Case k [J/m3] µr [-] Void Ratio0–1 102 0.05 0.5750–2 105 0.7 1.4590–3 102 0.7 0.5660–4 105 0.05 1.1200–5 5 × 102 0.4 0.5681 5570 0.13 0.5692 70,500 0.37 1.0643 9100 0.24 0.5654 36,900 0.63 0.6585 42,800 0.31 0.7216 55,400 0.66 0.8897 17,500 0.68 0.5668 56,400 0.49 0.8979 69,800 0.56 1.08010 55,400 0.56 0.88611 42,800 0.55 0.71612 56,400 0.46 0.897best 47,500 0.5 0.7762.6. Powder Spreading SimulationA powder spreading simulation was conducted using the identified parameters. Thecomputational geometry of the powder bed formation simulation is shown in Figure 4.Initially, 50,000 powder particles were deposited into a designated powder storage regionby gravity. Thereafter, a rectangular coater blade, located behind the storage region, wastranslated at a constant velocity of 25 mm/s to disperse the powder into the powder bedformation region located in front of the storage region. The depth of the bed formationregion was set to 80 mm. Moreover, the computational domain was specified as a rectan-gular volume extending beyond the geometrical boundaries of the model, and particlesreaching the domain boundaries were deleted. This boundary condition is the f boundarystyle implemented in LIGGGHTS [18].J. Manuf. Mater. Process. 2025, 9, 383 7 of 16Figure 4. An example snapshot of the powder spreading simulation.3. Results3.1. Optimization ResultsBefore conducting Bayesian optimization, preliminary analyses of consolidation wereperformed for the cases listed in case numbers 0–1 and 0–5 in Table 3. These cases corre-spond to combinations of the minimum and maximum values of the two parameter sets, aswell as the results obtained using intermediate values of the two parameters. The Table 3also provides the void ratios that were obtained from the DEM simulations. The void ratioexhibited a range of values between 0.575 and 1.459, depending on the parameter values.The experimental value was found to be 0.783, which is believed to be achievable withinthe specified parameter ranges. Consequently, the search ranges for the parameters wereset to k = 102 to 105 J/m3 and mr = 0.05 to 0.7.Using the five preliminary analysis results as prior information, optimization calcula-tions were performed, and the results are summarized in Table 3. The Table 3 presents theresults obtained from 12 trials. It is important to note that in these 12 trials, wa was variedbetween 0.01 and 1, depending on the results. Although no parameter set was identifiedthat could exactly reproduce the experimental void ratio of 0.783 across the 12 trials, cases5, 6, 8, 10, 11, and 12 achieved values that were relatively close to the experimental one.Hence, it was determined that when k and mr were designated as 47,500 J/m3 and 0.5,respectively, the DEM simulation yielded a void ratio of 0.776.Figure 5 presents a color map of simulated void ratio esim (k,µr). The horizontal andvertical axes correspond to cohesion-related parameter k and rolling-friction coefficient µr,respectively. The color field is obtained by interpolating the simulated values in Table 3using the SciPy library (Python). White circles mark the actual simulation points listed inTable 3, and the red star denotes the previously identified optimal condition. It has beenobserved that, under the framework of Bayesian optimization, the search process exhibits apropensity to explore the vicinity of the red star, even in the initial stages of optimization.Despite the fact that points in proximity to the red star were repeatedly predicted, itproved challenging to achieve convergence toward values that more closely aligned withthe experimental results, i.e., the optimal value. In Bayesian optimization, the predictionmodel is iteratively updated based on the accumulated results. Increasing the number oftrials does not necessarily yield a better parameter set than those already obtained. Thisbehavior differs from optimization methods that rely on gradient-based approaches. In thegradient-based approaches, the next candidate point is obtained by following the gradientfrom the initial value so that repeated trials progressively approach the optimum value. InJ. Manuf. Mater. Process. 2025, 9, 383 8 of 16contrast, Bayesian optimization employs Gaussian process regression to explore regionsof the solution space that are not included in prior information. This suggests that theparameter sets that can reproduce the experimental void ratio can be proposed, even witha small number of trials. In the present study, a near-optimal value was obtained as earlyas the fifth trial in Table 3. On the other hand, depending on the wa value, subsequentcandidates are directed towards the sparsely sampled region of the solution space. Hence,an increase in the number of trials does not guarantee the approach to the optimumvalue. This characteristic is both an advantage of Bayesian optimization and a factor thatcomplicates convergence assessment.Figure 5. Simulated void ratio esim (k,µr) over the (k,µr) domain. The color map is interpolated fromthe simulation results in Table 3; white circles indicate the computed data points (cases 0–1 to 12),and the red star marks the best condition identified in Table 3.In recent years, significant attention has been focused on concepts such as human-in-the-loop (HITL) optimization, e.g., [21–23]. The term HITL is used to denote methodsin which human judgment is incorporated into the optimization process. In the presentoptimization procedure, the value of wa was adjusted by the user based on the simulationresults, which can be considered a form of HITL. Nevertheless, the evaluation of the resultswas based on a mechanistic criterion: the minimization of the squared error relative to theexperimental values. Consequently, although a human might have estimated the optimumby the fifth or sixth trial, additional trials were conducted. For problems where the solutionis a continuous function, gradient-based optimization methods guarantee convergencetowards the optimum with repeated trials, which allows fully automated optimization. Incontrast, for many engineering problems where the form of the solution is unknown, theactive incorporation of human judgment in the selection of optimization parameters andthe evaluation of objective values can be crucial to achieving successful optimization.It is important to note that the number of parameters to be identified was only two inthe present study. In many engineering optimization problems, the number of parametersmay often exceed ten. In such cases, it becomes critically important to cover the minimumand maximum ranges of the solution space with as few preliminary trials as possible. Suchtechniques, including orthogonal arrays and Latin hypercube sampling, can be employedJ. Manuf. Mater. Process. 2025, 9, 383 9 of 16to sample the parameter sets. However, prior to Bayesian optimization, a minimum ofseveral tens of preliminary simulations are likely to be necessary.In the original optimization trial, we varied wa (the EI exploration parameter); hencethe effect of wa on convergence is not clear. Consequently, a parametric study was conductedwith respect to wa. To maintain a reasonable computational cost, the DEM evaluations inthe workflow of Figure 1 were substituted with a surrogate model that was trained on thedataset presented in Table 3. The surrogate model maps the two inputs (k, µr) to the voidratio (output). Six regression algorithms were assessed in this study: Gaussian ProcessRegression (GPR), Random Forest (RF), Gradient Boosting (GB), Support Vector Regression(SVR), Ridge Regression, and k-Nearest Neighbors (kNN). The RF surrogate was selectedbecause it yielded the highest coefficient of determination (R2) on hold-out validation.A single DEM run requires about 8 h, whereas the surrogate returns a prediction inless than 1 s. This speed enhancement enables a systematic exploration over wa. Note thatthe number of simulation results were limited, which implies that higher-fidelity surrogaterequires additional cases.Figure 6 shows the evolution of the minimum objective function value during theoptimization. The execution of Bayesian optimization with wa = {0.01, 0.05, 0.1, 0.5, 1.0} heldfixed has yielded the same result: k = 48,779.3 J/m3 and mr = 0.60. These results suggestthat the impact of wa on the optimization path and the attained optimum is negligibleunder the present conditions.Figure 6. The relationship between the number of iterations and best objective value obtained fromthe Bayesian optimization process using the machine learning surrogate model.3.2. Powder Spreading Simulation Using the Parameters EstimatedFigure 7 presents snapshots of the simulations performed with the parameter setsidentified by Bayesian optimization and the simulation lacking the adhesion model. Inthe simulation with adhesion, a strong inter-particle adhesive force is apparent, resultingin the collective movement of powder particles as a uniform mass. Moreover, the strongadhesion of the powder bed results in its formation being limited to a small portion of thedomain. The powder bed formed under these conditions appears to differ to that observedin actual PBF-LM manufacturing. This finding suggests that the parameter set calibrated toreproduce the consolidation may overestimate the adhesive properties of the particles.J. Manuf. Mater. Process. 2025, 9, 383 10 of 16 Figure 7. Snapshots of powder spreading simulations with and without adhesion.Accordingly, simulations were conducted with k set to 23,750 and 15,833 J/m3; theseresults, the result with k = 47,500 J/m3, and the result without cohesion are illustrated inFigure 8. It is evident that as the magnitude of k decreases, the inter-particle adhesionweakens, which results in the formation of a broader powder bed region. Accordingto previous experimental studies [24,25], the powder beds do not become densely andhomogeneously packed but instead appear to exhibit a sparse structure. In the simulationwith k = 15,833 J/m3, a sparse structure consistent with previous experimental observationsis formed.J. Manuf. Mater. Process. 2025, 9, 383 11 of 16Figure 8. Snapshots of final step of powder spreading simulations without cohesion (a), withk = 15,833 J/m3 (b), 23,750 J/m3 (c), and 47,500 J/m3 (d).4. DiscussionsIn light of the aforementioned observations, it can be concluded that the parameterset determined to reproduce the consolidation is not able to reproduce the powder bedformation. This discrepancy is hypothesized to arise from two factors: (1) the differencein scaling and (2) difference in particle-scale dynamics between the consolidation andpowder spreading.4.1. ScalingGuided by Buckingham’s P theorem, we first specify the response of interest—here,the void ratio, e, which is dimensionless—and identify the dominant physical quantitiesgoverning uniaxial consolidation under quasi-static loading. With the representative length(particle diameter D) held fixed and other conditions kept constant, the principal scalesare the applied consolidation stress, σ, and the material stiffness (Young’s modulus E). Aconvenient dimensionless group that characterizes consolidation is thereforeΠ1 =σE. (3)Our scaling requirement is to match P1 for the experiment and simulation. Denotingthe experimental conditions by subscript “0” and DEM conditions by “1”, we imposeσ0E0=σ1E1. (4)In the experiments, σ0 = 9 kPa and E0 = 2 × 1011 Pa. To mitigate the time step constraintin the DEM simulations, we set E1 = 5 × 106 Pa. Enforcing Equation (4) yields the equivalentconsolidation stress:σ1 = σ0E1E0= 0.225 Pa. (5)Thus, by matching σ/E, we preserve the leading non-dimensional control parameterrelevant to e while permitting a smaller E (and hence a larger stable time step) in DEM. Wenote that if additional mechanisms (e.g., adhesion or rate effects) are significant, furtherP groups may be required; in the present consolidation setting, P1 = σ/E provides thefirst-order scaling used to set the simulation load.For the three cases of (1) k = 48,000 J/m3 and mr = 0.5, (2) k = 48,000 J/m3 and mr = 0,and (3) without cohesion, consolidation simulations at s = 0.225, 0.45, 9, 90, 900, and 9000 Pawere conducted, and the corresponding void ratios are shown in Figure 9. As illustratedJ. Manuf. Mater. Process. 2025, 9, 383 12 of 16by the Figure 9, when 9 kPa consolidation stress is applied, the simulated void ratio islower than in the experiment, even with the adhesion model incorporated. This can beattributed to the fact that, for the given particle stiffness, the applied pressure is relativelyhigh. This allows easier compression and results in a lower void ratio. It is hypothesizedthat the reduction in the void ratio is the result of two factors: firstly, the rearrangement ofthe particle assembly, and secondly, an increase in particle overlap.Figure 9. Results of the consolidation simulations varying the normal stress: relationship between thenormal stress and the void ratios for the cases of (1) k = 48,000 J/m3 and mr = 0.5, (2) k = 48,000 J/m3and mr = 0, and (3) without cohesion.The Figure 9 illustrates the void ratio (eexp) derived from the consolidation experimentconducted at s = 9 kPa. At σ = 0.225 Pa, where the scaling is satisfied using Young’smodulus from the DEM analysis, the void ratio is slightly larger than eexp for case (1).Of the three cases under consideration, case (1) is the closest to the cohesion parametersdetermined by Bayesian optimization under conditions that do not satisfy the scaling law.In case (1), the void ratio at σ = 0.225 Pa deviates significantly from the experimental value.This finding indicates that the DEM simulation, which did not satisfy the scaling law, mayhave resulted in excessive cohesion parameters.4.2. Difference in Particle-Scale DynamicsIn the case of consolidation, the application of s = 9 kPa corresponds to the resistanceforce exerted by the particle assembly. In this case, the lateral boundaries are cylindrical, andthe top and bottom are constrained by disks, which prevent shear deformations that wouldgenerate significant sliding along the tangential directions between particles, as shown inFigure 10. In other words, the resistance against external pressure is primarily provided byinter-particle normal forces. The JKR-type model applied in this study provides cohesionalong the contact normal direction between particles. In contrast, the rotational resistancemodel insignificantly affects the contact normal direction but suppresses rotations due totorque, thereby contributing to tangential inter-particle forces.J. Manuf. Mater. Process. 2025, 9, 383 13 of 16Figure 10. Particle-scale schematics of consolidation and direct shear in the powder spreading processes.In the consolidation scenario, once the particle arrangement is partially fixed, thenormal inter-particle forces become the dominant mechanism. Thus, an increase in kenables the reproduction of the consolidation phenomenon. In contrast, as shown inFigure 10, direct shear occurs between the coater blade and the powder particles duringpowder bed formation simulations. Within this shear zone, the separation at inter-particlecontacts and sliding at contact points become the dominant phenomena. Since the powderabove the shear band zone will disappear due to the coater movement, the powder particleswithin the bed formation region are originally within the shear zone.In the angle of repose-based calibration by Lupo et al. [13], the authors conducteda comprehensive analysis of powder behavior during powder spreading. It has beenreported that, at the proximity of the coater, shear stresses exerted between the coaterblade and the particle assembly induce irregular particle motion, thereby promoting theformation of localized voids in front of the coater. In the present study, we likewise inferthat inter-particle shear is dominant immediately before and after powder bed formation,providing evidence that supports the observations made by Lupo et al.The findings indicate that distinct optimal parameters may be required for consolidationand shear-dominated powder spreading. Prior studies typically calibrate against a singleexperimental modality—such as a consolidation metric or an angle-of-repose test—ratherthan jointly across consolidation and shear. To the best of our knowledge, few studies haveidentified DEM parameters simultaneously for both regimes on the same powder within aunified framework. This gap matters because powders dissipate energy efficiently throughinter-particle friction and sliding; consequently, the governing micromechanics—and theassociated parameter sensitivities—differ between normal-force-dominated consolidationand tangential-sliding-dominated spreading. Consistent with this, the parameter set cali-brated for consolidation does not accurately reproduce spreading, underscoring both thenonlinearity of powder behavior and the need for phenomenon-aware calibration.Guided by these insights, our planned work proceeds in two stages. First, we will iden-tify parameters under consolidation while enforcing stress–stiffness scaling (e.g., matchedσ/E) to ensure proper non-dimensional scaling. Second, if that set still mispredicts spread-ing, we will pivot to shear-relevant observables for powder bed formation (e.g., bulkdensity, area coverage, shear metrics) and reassess both scaling and the parameter valuesfor the shear regime. This strategy advances a multi-mechanism fitting paradigm in whichparameters are validated across the distinct physical mechanisms governing consolidationand spreading.J. Manuf. Mater. Process. 2025, 9, 383 14 of 165. ConclusionsIn the present study, a Bayesian optimization-based parameter identification methodwas developed to minimize the squared error between the powder bed void ratio obtainedfrom consolidation experiments and that predicted by DEM simulations. This methodologywas applied to metallic powders used in PBF-LM, and the identified parameter set wassubsequently employed for powder spreading simulations.• A sequential optimization process was conducted utilizing Bayesian optimization.Although the optimal parameter set was not directly identified within twelve trials,several parameter sets close to the optimum were obtained as early as the fifth trial.• Proposing the next candidate parameter sets based on the user’s judgment rather thanrelying solely on Gaussian process regression resulted in the successful attainment ofa parameter set consistent with the experimental values.• The parameter set identified for reproducing consolidation, in which inter-particlenormal forces were dominant, was found to be incapable of reproducing the powderbed formation.• This inconsistency may be caused by two main aspects: (1) the scaling law (the ratioof Young’s modulus to consolidation stress) between the experiments and DEM doesnot match and (2) there is a difference between the particle-scale dynamics betweenconsolidation and direct shear in powder spreading.In the present study, the identification of parameters was conducted via Bayesianoptimization for the consolidation process. Utilizing the same experimental apparatus, it isalso possible to evaluate the shear strength, or internal friction angle, of a powder layerunder confined pressure. Future work will focus on determining parameters based onshear strength and direct experimental observations of powder bed formation.In this study, we employed Bayesian optimization as the sole optimizer. Althoughnumerous gradient-based and metaheuristic methods are available, the strong nonlinearityand discontinuities of DEM simulations make reliable and robust gradient evaluationdifficult, limiting the applicability of robust gradient-based schemes. In practice, derivative-free metaheuristics—such as genetic algorithms, particle swarm optimization, and antcolony optimization—have also been widely used. A systematic assessment of how thechoice of optimizer affects convergence behavior and solution quality will be an importanttopic for future work.Author Contributions: Conceptualization, J.K., M.K. and M.W.; methodology, M.K. and J.K.; software,M.K. and J.K.; validation, M.K. and J.K.; formal analysis, J.K., M.K. and M.W.; investigation, J.K.;resources, J.K. and M.K.; data curation, J.K.; writing—original draft preparation, J.K. and M.K.;writing—review and editing, J.K., M.K. and M.W.; visualization, J.K. and M.K.; project administration,M.W.; funding acquisition, M.W. All authors have read and agreed to the published version ofthe manuscript.Funding: Financial support came from the Innovative Science and Technology Initiative for Security,Grant Number JPJ004596, implemented by the Acquisition, Technology, and Logistics Agency (ATLA).Data Availability Statement: The original contributions presented in this study are included in thearticle. Further inquiries can be directed to the corresponding author.Acknowledgments: We acknowledge the financial support from Innovative Science and TechnologyInitiative for Security, Grant Number JPJ004596 implemented by the Acquisition, Technology, andLogistics Agency (ATLA).Conflicts of Interest: The authors declare no conflicts of interest.J. Manuf. Mater. Process. 2025, 9, 383 15 of 16AbbreviationsThe following abbreviations are used in this manuscript:PBF-LM Powder bed fusion laser meltingBO Bayesian optimizationDEM Discrete element methodHX Hastelloy XSJKR Simplified Johnson–Kendall–RobertsHITL Human in the loopReferences1. DebRoy, T.; Wei, H.L.; Zuback, J.S.; Mukherjee, T.; Elmer, J.W.; Milewski, J.O.; Beese, A.M.; Wilson-Heid, A.; De, A.; Zhang, W.Additive manufacturing of metallic components—Process, structure and properties. Prog. Mater. Sci. 2018, 92, 112–224. [CrossRef]2. Seede, R.; Shoukr, D.; Zhang, B.; Whitt, A.; Gibbons, S.; Flater, P.; Elwany, A.; Arroyave, R.; Karaman, I. An ultra-high strengthmartensitic steel fabricated using selective laser melting additive manufacturing: Densification, microstructure, and mechanicalproperties. Acta Mater. 2020, 186, 199–214. [CrossRef]3. Fürstenau, J.-P.; Wessels, H.; Weißenfels, C.; Wriggers, P. Generating virtual process maps of slm using powder-scale sphsimulations. Comput. Part. Mech. 2020, 7, 655–677. [CrossRef]4. Grange, D.; Queva, A.; Guillemot, G.; Bellet, M.; Bartout, J.D.; Colin, C. Effect of processing parameters during the laser beammelting of inconel 738: Comparison between simulated and experimental melt pool shape. J. Mater. Process. Technol. 2021, 289,116897. [CrossRef]5. Alphonso, W.; Baier, M.; Carmignato, S.; Hattel, J.; Bayat, M. On the possibility of doing reduced order, thermo-fluid modelling oflaser pow- der bed fusion (l-pbf)—Assessment of the importance of recoil pressure and surface tension. J. Manuf. Process. 2023, 94,564–577. [CrossRef]6. Katagiri, J.; Kusano, M.; Minamoto, S.; Kitano, H.; Daimaru, K.; Tsujii, M.; Watanabe, M. Melt pool shape evaluation by single-track experiments and finite-element thermal analysis: Balling and lack of fusion criteria for generating process window ofinconel738lc. Materials 2023, 16, 1729. [CrossRef]7. Katagiri, J.; Nomoto, S.; Kusano, M.; Watanabe, M. Particle size effect on powder packing properties and molten pool dimensionsin laser powder bed fusion simulation. J. Manuf. Mater. Process. 2024, 8, 71. [CrossRef]8. Masuda, H. Adhesion of powder particles. Denshi Shashin Gakkaishi (Electrophotogr.) 1997, 36, 169–174. [CrossRef]9. Cundall, P.A.; Strack, O.D. A discrete element model for granular assemblies. Geotechnique 1979, 29, 47–65. [CrossRef]10. Li, T.; Peng, Y.; Zhu, Z.; Zou, S.; Yin, Z. Discrete Element Method Simulations of the Inter-Particle Contact Parameters for theMono-Sized Iron Ore Particles. Materials 2017, 10, 520. [CrossRef]11. Yim, S.; Bian, H.; Aoyagi, K.; Yamanaka, K.; Chiba, A. Spreading behavior of Ti–48Al–2Cr–2Nb powders in powder bed fusionadditive manufacturing process: Experimental and discrete element method study. Addit. Manuf. 2022, 49, 102489. [CrossRef]12. Geer, S.; Bernhardt-Barry, M.L.; Garboczi, E.J.; Whiting, J.; Donmez, A. A more efficient method for calibrating discrete elementmethod parameters for simulations of metallic powder used in additive manufacturing. Granul. Matter 2018, 20, 77. [CrossRef]13. Lupo, M.; Ajabshir, S.Z.; Sofia, D.; Barletta, D.; Poletto, M. Discrete element method model calibration and validation for thespreading step of the powder bed fusion process to predict the quality of the layer surface. Particuology 2024, 94, 261–273.[CrossRef]14. Kitano, H.; Kusano, M.; Tsujii, M.; Yumoto, A.; Watanabe, M. Process parameter optimization framework for the selective lasermelting of hastelloy x alloy considering defects and solidification crack occurrence. Crystals 2021, 11, 578. [CrossRef]15. Freeman, R. Measuring the flow properties of consolidated, conditioned and aerated powders—A comparative study using apowder rheometer and a rotational shear cell. Powder Technol. 2007, 174, 25–33. [CrossRef]16. Liu, Y.; Guo, X.; Lu, H.; Gong, X. An investigation of the effect of particle size on the flow behavior of pulverized coal. ProcediaEng. 2015, 102, 698–713. [CrossRef]17. Koynov, S.; Duda, K.; Santos, P.A.D.L.; Goldfarb, D.J. Comparative evaluation of peschl and ft4 full-bed rotational shear cells forpowder flow characterization. Powder Technol. 2025, 456, 120810. [CrossRef]18. DCS Computing GmbH. Liggghts-Public Documentation, Version 3.x. 2025. Available online: https://www.cfdem.com/media/DEM/docu/Manual.html (accessed on 15 August 2025).19. Haynes International. Hastelloy-X. 2025. Available online: https://haynesintl.com/en/alloys/alloy-portfolio/high-temperature-alloys/hastelloy-x/#physical-properties (accessed on 15 August 2025).20. The GpyOpt Authors. Gpyopt: A Bayesian Optimization Framework in Python. 2016. Available online:http://github.com/SheffieldML/GPyOpt (accessed on 15 August 2025).https://doi.org/10.1016/j.pmatsci.2017.10.001https://doi.org/10.1016/j.actamat.2019.12.037https://doi.org/10.1007/s40571-019-00296-3https://doi.org/10.1016/j.jmatprotec.2020.116897https://doi.org/10.1016/j.jmapro.2023.03.040https://doi.org/10.3390/ma16041729https://doi.org/10.3390/jmmp8020071https://doi.org/10.11370/isjepj.36.169https://doi.org/10.1680/geot.1979.29.1.47https://doi.org/10.3390/ma10050520https://doi.org/10.1016/j.addma.2021.102489https://doi.org/10.1007/s10035-018-0848-4https://doi.org/10.1016/j.partic.2024.08.010https://doi.org/10.3390/cryst11060578https://doi.org/10.1016/j.powtec.2006.10.016https://doi.org/10.1016/j.proeng.2015.01.170https://doi.org/10.1016/j.powtec.2025.120810https://www.cfdem.com/media/DEM/docu/Manual.htmlhttps://www.cfdem.com/media/DEM/docu/Manual.htmlhttps://haynesintl.com/en/alloys/alloy-portfolio/high-temperature-alloys/hastelloy-x/#physical-propertieshttps://haynesintl.com/en/alloys/alloy-portfolio/high-temperature-alloys/hastelloy-x/#physical-propertieshttp://github.com/SheffieldML/GPyOptJ. Manuf. Mater. Process. 2025, 9, 383 16 of 1621. Kim, M.; Ding, Y.; Malcolm, P.; Speeckaert, J.; Siviy, C.J.; Walsh, C.J.; Kuindersma, S. Human-in-the-loop bayesian optimization ofwearable device parameters. PLoS ONE 2017, 12, e0184054. [CrossRef]22. Sousa, J.; Sousa, A.; Brueckner, F.; Reis, L.P.; Reis, A. Hu-man-in-the-loop multi-objective bayesian optimization for directedenergy deposition with in-situ monitoring. Robot. Comput.-Integr. Manuf. 2025, 92, 102892. [CrossRef]23. Kim, D.B.; Bajestani, M.S.; Lee, J.Y.; Shin, S.-J.; Kim, G.-Y.; Sajadieh, S.M.M.; Noh, S. Human-in-the-loop in smart manufacturing(h-sm): A review and perspective. J. Manuf. Syst. 2025, 82, 178–199. [CrossRef]24. Kikuchi, K.; Tanifuji, Y.; Zhou, W.; Nomura, N.; Kawasaki, A. Experimental characterization and computational simulation ofpowder bed for powder bed fusion additive manufacturing. Mater. Trans. 2022, 63, 931–938. [CrossRef]25. Okugawa, M.; Isono, Y.; Koizumi, Y.; Nakano, T. Raking process for powder bed fusion of ti–6al–4v alloy powder analyzed bydiscrete element method. Mater. Trans. 2023, 64, 37–43. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.https://doi.org/10.1371/journal.pone.0184054https://doi.org/10.1016/j.rcim.2024.102892https://doi.org/10.1016/j.jmsy.2025.05.020https://doi.org/10.2320/matertrans.MT-Y2021005https://doi.org/10.2320/matertrans.MT-MLA2022010 Introduction  Materials and Methods  Materials Used  Consolidation Experiment  Discrete Element Method Overview  Consolidation Simulation  Bayesian Optimization  Powder Spreading Simulation  Results  Optimization Results  Powder Spreading Simulation Using the Parameters Estimated  Discussions  Scaling  Difference in Particle-Scale Dynamics  Conclusions  References