Computer simulation of pressure drop in a powder bed stabilized by a gas flow for 3D printing process in a zero gravity environment

Abstract

The process known as powder bed fusion exhibits a lack of stabilization in the initial layers. One approach to address this challenge involves introducing a gas flow across the particulate medium. This solution enables more targeted applications of the technique, especially in zero-gravity environments. The objective of this study is to validate a methodology that utilizes computational fluid dynamics (CFD) and the discrete element method (DEM) to apply a gas flow within a porous medium under gravity-free conditions. The validity of the proposed solution was assessed using Ergun’s equation, resulting in an error of 6.65%, and the Kozeny-Carman equation, resulting in an error of 10.53%. The drag force exerted on the particles in the simulations surpassed the gravitational force (1g). Consequently, the application investigated in this study represents an effective alternative for employing 3D printing in the absence of gravitational forces.

Keywords:
computational fluid dynamic; discrete element method; additive manufacturing; zero gravity

INTRODUCTION

In recent years, 3D printing technologies have gained widespread popularity worldwide, offering various methods for depositing and applying construction materials. These technologies emerged in the 1980s and 1990s^{1), (2} and have evolved from simple applications like toys to more advanced uses, including the aerospace industry³. Consequently, this technology has become a significant subject of study. Among its technological variants, each possesses distinct advantages and drawbacks relative to the others⁴. Notably, one extensively discussed and well-developed technique is powder bed fusion (PBF), wherein a laser beam is employed to fuse powdered particles. Within this category fall techniques like selective laser sintering/melting (SLS/SLM)⁵. However, a challenge associated with SLS is the necessity of a support base to stabilize the component during fabrication, leading to increased material consumption and additional time for support removal. Efforts to enhance support construction have been proposed^{6), (7}, but such approaches still demand additional time and material resources. Conversely, a methodology proposed by Zocca et al.⁸ involves applying gas through the powdered layer, ensuring greater stability and cohesion among the particles. Additionally, they suggest that this technology could be employed in zero-gravity environments, as the absence of gravitational force is counterbalanced by the drag force imparted on the particles by the gas stream. The lack of gravitational force affects particle compaction, potentially reducing the structural integrity of the manufactured part.

The SLS process stands out as one of the most rapidly advancing additive manufacturing (AM) technologies. Operating with a continuous or pulsed laser beam acting as a heat source, SLS selectively scans the powder bed, sintering particles layer by layer according to predetermined geometries and sizes. Diverse raw materials can serve as sources for part fabrication within this technique, encompassing waxes, cermet (metal-ceramic composites), ceramics, metallic-polymeric powders, metals, alloys, and steels. Notably, when ceramics are utilized, components with minimal post-sintering shrinkage can be produced^{9), (10}. An advantage of the SLS process is its capacity to work at very small scales, even in the micrometer range². This leads to a challenge related to the investigation of particle-particle and particle-fluid interactions, a difficulty that can be addressed through a combined approach involving: a) computational fluid dynamics (CFD), which analyzes fluid interactions and fluid-solid interactions; and b) the discrete element method (DEM), primarily applied to studying particle interactions.

The primary objective of this study is to validate a methodology that employs CFD and DEM to apply a gas flow within a micrometric-scale porous medium under conditions without gravity. This aims to characterize the hydrodynamic properties of the medium, enhancing understanding and increasing efficiency. Such insights can aid in making informed decisions when developing models for 3D printing via a powder bed within environments devoid of gravity.

MODELING

In this section, the models employed for the DEM simulation, geometry generation, and CFD simulations, alongside other adopted considerations are present. The study was subdivided into ten simulations, each incrementing the layer by 100 µm, which aligned with the height of the 3D printing layers utilized⁸. The uniform script was followed for each simulation: i) DEM simulation to determine particle coordinates; ii) generation of the STL file based on the coordinates; and iii) importation, geometry creation, and CFD simulation. The particles were modeled as regular spheres with a consistent diameter.

In scenarios involving interactions between particles as well as particles and walls, Newton’s laws for linear motion and Euler’s equations for rotational motion were employed for each particle. A time interval on the order of 10^-5 or 10^-6 s was utilized in this method, enabling the calculation of the new particle state¹¹. The general form of the momentum equation for a particle is:

$m_p\frac{d{v}_p}{d_t}=F_s+F_b$ (A)

where F_s is the superficial particle force, F_b is the body force, m_p is the mass of particle p, and v_p is the velocity of particle p. They are decomposed into:

$F_s=F_d+F_p+F_{vm}$ (B)

$F_b=F_g+F_u$ (C)

where F_d is the drag force, F_p is the pressure gradient force, F_vm is the virtual mass force, F_g is the gravitational force, and F_u is the pre-set body force.

DEM simulation was employed to establish the particle configuration within the powder bed. The selected material was aluminum oxide, chosen due to its ability to achieve high packing density¹². It has also been employed in experimental prints of this kind⁸, with a diameter of 38 µm, possessing a density (ρ) of 3950 kg/m³, Poisson’s ratio (ν_p) of 0.22, a coefficient of restitution (e) of 0.75, a coefficient of static friction (µ_s) of 0.5, and a coefficient of rolling friction (µ_n) of 0.1, employing the Hertz-Mindlin contact model¹³, which operates as a friction-based contact model¹⁴. The DEM time-step was set to 1.0x10^-6, and a gravitational acceleration of -9.81 m/s² was applied along the y-axis. To simulate particle compaction within the control volume, a random entry mode for particles was implemented, allowing them to interact with gravitational force. The stopping criterion used was the virtual absence of kinetic energy among the particles. The dimensions of the enclosure were x, y, z=150, 1000, 150 µm, respectively. Boundaries perpendicular to the flow direction were configured as symmetric interfaces, simulating an infinite domain. This choice eliminated wall influences as advised by Choi et al.¹⁵. Determining the actual dimensions of the entire powder bed is challenging due to extensive computational demands¹⁶. Xu and Jiang¹⁶ point out that numerically simulating the Navier-Stokes equations for fluid flow in a porous medium with all particles is immensely intricate due to the extensive calculations involved. Particle insertion within the enclosure (Fig. 1) was executed randomly, and subsequently, the particles were allowed to settle under the influence of gravity. A singular simulation sufficed, serving as the foundation for further division into CFD simulations by layers.

Figure 1
Image showing the control volume on y=1000 µm for DEM simulations.

Upon concluding the DEM simulations, the central point coordinates and particle diameter (uniform in this case) were obtained. These data were then employed in a Python script to generate geometries for subsequent CFD analyses. Given that each CFD simulation corresponded to a specific layer thickness (100, 200, ..., 1000 µm), only particles situated entirely below the layer’s height were imported. This criterion ensured that particle center points were situated beneath the layer’s height minus the sphere’s radius (y_layer - radius). The CFD analysis focused on calculating pressure drop and drag forces within the gas flow interacting with the particles. The particles’ arrangement was dictated by the outcomes of the DEM simulations. For each layer, porosity was determined by considering the cell volume (mesh cells) within the control volume. This calculated porosity value was subsequently incorporated into the theoretical Ergun and Kozeny-Carman equations to validate the method’s efficacy. It’s noteworthy that in the CFD simulations, particles were treated as stationary entities. This approach aligned with the DEM-CFD one-way coupling model employed by Choi et al.¹⁵.

The current study employed a laminar model with a segregated flow approach based on the SIMPLE algorithm, utilizing second-order discretization for convection terms. Nitrogen gas was utilized¹², and the system was unaffected by gravitational effects. The gas flow traversed the particle region in the negative y-axis direction. Pressure drop was computed between the layer height and the y=0 plane. Symmetry plane conditions were imposed on the walls. The input and output domains were configured as a velocity inlet with v_y=-0.118 m/s and a pressure outlet with p₀=0, respectively. To mitigate the influence of inlet and outlet boundaries, entry and exit regions were incorporated within the calculation zone, as indicated by Choi et al.¹⁵. These regions measured 380 µm (equivalent to 10 diameters) in length, as depicted in Fig. 2. Particle surfaces were treated as walls subject to a no-slip condition. Prism layers were applied to this region to facilitate boundary layer calculations. The convergence criteria involved attaining a steady-state solution, and continuing calculations until residuals reached values less than 10^-7 for the continuity and momentum equations in the x, y, and z directions. For the mesh, a polyhedral model with a prismatic mesh configuration was used. Additionally, a mesh convergence analysis was conducted for the first layer case, and the outcomes of this analysis were applied across all simulations.

Figure 2
Image showing the geometric model and some boundary conditions for mesh generation.

It is crucial to recognize that, despite addressing a micrometer-scale problem, the continuum hypothesis for fluids remains applicable¹⁷. In this context, the flow was presumed to be both stationary and incompressible. Based on these assumptions, the velocity fields of the fluid can be established by solving the continuity equation and the Navier-Stokes equations. For an incompressible fluid, these equations take the form:

$\frac{\partial U_i}{\partial x_i}=0$ (D)

$\frac{\partial U_i}{\partial t}+U_j\frac{\partial U_i}{\partial x_i}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+V\frac{{\partial }^2{U}_i}{\partial x_j^2}$ (E)

where U_i is the ith component of the fluid velocity, ρ is the density, v is the kinematic viscosity, and p is the fluid pressure. The particle’s drag force can be decomposed into two components: one arising from viscous stress and the other stemming from pressure effects¹⁸. In the context of CFD, the force exerted on a body in the flow direction is monitored. Consequently, the total drag force can be calculated using the following formula¹⁴:

$f={\Sigma }_f\left(f_f^{pressure}+f_f^{shear}\right).n_f$ (F)

where $f_f^{pressure}$ and $f_f^{shear}$ are the pressure and shear force vectors on the surface face, respectively, f and n_f are the user-specified direction vector that indicates the direction to compute the force.

$f_f^{pressure}=\left(p_f-p_{ref}\right).a_f$ (G)

$f_f^{shear}=T_f,a_f$ (H)

where p_f is the face static pressure, a_f is the face area vector, p_ref is the reference pressure, and T_f is the stress tensor at face f. To assess whether the resulting drag force can effectively counteract the force of gravity and stabilize the powder bed, it becomes imperative to theoretically quantify the drag force. In instances of low Reynolds number flows, where viscous forces dominate (laminar flow with Re<1), the drag force can be calculated using the Stokes equation:

$F_d=6\pi .R.u$ (I)

where F_d is the drag force, µ is the dynamic viscosity, R is the radius, and u is the gas velocity. The gravitational force is calculated according to:

$F_g=\rho .V.g$ (J)

where F_g is the gravitational force, ρ is the density of the particle, V is the particle volume (V=4π.R³/3), and g is the gravitational acceleration.

In the post-processing phase, an analysis was conducted on a laminar and incompressible flow, treating the powder bed as a porous medium. To validate this approach, Ergun’s law and the Kozeny-Carman equation were employed. Typically, Ergun’s complete equation is utilized to correlate pressure drop in compacted beds, and it does not consider wall effects¹⁹. This equation comprises two components: the first term on the right side of the equation quantifies frictional losses, while the second term accounts for inertial losses²⁰. It is mathematically defined as follows:

$\frac{\Delta p}{L}=150\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3}·\frac{\mu .U_m}{D_p^2}+1.75\frac{\left(1-\epsilon \right)}{{\epsilon }^3}·\frac{{\rho }_m.U_m^2}{D_p}$ (K)

where ΔP is the drop pressure, L is the layer height, ε is the porosity, µ is the dynamic viscosity, ρ_m is the fluid density, U_m is the fluid superficial velocity, and D_p is the particle diameter. The Kozeny-Carman equation presents a classical method for correlating permeability with the attributes of a densely packed powder. It can be regarded as a more comprehensive extension of the Ergun equation where the inertial effect is disregarded. This equation holds true for scenarios involving a random and homogeneous arrangement of particles, particularly within the context of laminar flow conditions:

$\frac{\Delta p}{L}=\frac{180\mu }{{\varnothing }^2.D^2}·\frac{{\left(1-\epsilon \right)}^2}{{\epsilon }^3}·\frac{Q}{A}$ (L)

where Q is the flow rate, A is the cross-sectional area, and ϕ is the sphericity of the particle (for spherical particles ϕ=1).

RESULTS AND DISCUSSION

The arrangement of the resulting bed is depicted in Fig. 3a. As the side walls were configured as periodic boundaries, particles possessed the capability to extend beyond these boundaries, a prerequisite to emulate an infinite domain. The behavior of the particle count within the control volume over time, up to stabilization, is illustrated in the graph found in Fig. 3b. In the process of constructing the layers (y=100, 200, ..., 1000 µm), the Y-coordinate of each particle was used as a parameter. It was assumed that the sum of the particle’s center coordinate and its radius value must be smaller than the dimension of the given layer.

Figure 3
Result of DEM simulation of packing bed at very low kinetic energy (a) and graph of total particles over time (b).

By utilizing a Python script within Blender, an open-source software, a geometry was generated and saved in the standard triangulation format (.stl), commonly employed in CAD programs. This geometry was then exported and utilized to establish the mesh for CFD simulations. In order to strike a balance between achieving mesh convergence and minimizing computational costs, a mesh convergence study was conducted. The variable under control was the pressure drop in simulations involving the first layer (y=100 µm). The boundary conditions and models employed were consistent throughout the study. The meshing techniques employed included the polyhedral mesher, prism layer mesher, and surface remesher. In the region of the particles, a cubic volumetric control was implemented for localized mesh refinement within the zone of interest. The principal attributes of the selected mesh are outlined in Table I. These values encompassed the base size, considered an absolute measure, while the others were configured in relation to it. This configuration included the cube size within the particle region, aimed at improving mesh quality and enhancing the capture of flow details. When a smaller size control was not defined, the target surface size corresponded to the face size. This adjustment proved essential for attaining greater refinement and effectively capturing intricate features within the flow²¹. The outcome of the mesh convergence study is illustrated in Fig. 4, with emphasis on the chosen mesh. The selection process was guided by identifying the point where the curve exhibited a substantial deflection, indicating computational limitations. Refinement was achieved through the reduction of the relative mesh size value within the control volume. At the highlighted point, a minor variation was observed in the pressure difference results. However, there existed a significant variation in the number of cells, which would inevitably lead to elevated computational costs.

Figure 4
Pressure drop versus the number of elements (mesh test behavior).

Table I
Mesh sizes.

Size	Value
Base size	30 µm
Target surface size	50% of base size
Minimum surface size	5% of base size
Prism layer total thickness	10% of base size
Relative size (volumetric control)	6% of base size

Fig. 5 presents three stages of mesh refinement for the simulation involving a single layer (y=100 µm). In the initial stage, a higher mesh density was evident within the particle region in contrast to the region characterized by free-flowing fluid. This discrepancy arose due to the implementation of the cubic control volume, particularly within the critical flow zone mentioned earlier. In the subsequent stage, the prismatic mesh region encompassed the particles. This prismatic mesh region was responsible for boundary layer calculations, on surfaces of the particles, where a non-slip condition was imposed. Following pre-processing, simulations were executed using Eqs. D and E, while post-processing involved fluid flow analysis. Additionally, the drag force was evaluated both theoretically (Eq. I) and numerically (Eq. F), enabling a comparison to the gravitational force acting on the particle (Eq. J). In the initial phase, the obtained results were matched against the theoretical predictions from the Ergun and Kozeny-Carman laws, which were calculated using the porosities derived from computer simulations. Notably, the initial two layers exhibited higher porosity due to the influence of the base on which they are deposited. This base acts as a barrier for particles while permitting fluid passage. As outlined by Chen et al.²², regions near walls displayed elevated porosity due to a more structured packing arrangement. However, this effect diminished moving away from the wall due to the random disposition of particles. As distance from the y=0 layer increased, porosity further decreased. If more layers were added, the anticipated trend would be a gradual uniformity in porosity.

Figure 5
Mesh refinement around a particle.

Fig. 6 portrays the relationship between pressure drop (Δp) and the height of the powder bed layer (h) across the simulations. The pressure drop was computed as the disparity in pressure between the onset and end of the powder bed layer. This pressure drop emerged from the friction between the fluid and particles, resulting in heat dissipation²³. The recorded values were juxtaposed against the theoretical predictions of the Ergun equation (Eq. K) and the Kozeny-Carman equation (Eq. L). Upon observation, it became evident the discrepancy between the pressure drop obtained via CFD and the one predicted by Ergun’s law remained minimal. The simulated values tended to be overestimated, resembling the behavior observed in slow flows as reported by Reddy and Joshi²⁴. However, notable errors emerged in the initial layers, potentially attributed to their higher porosity. The average error between the theoretical value (Ergun’s law) and CFD stood at 6.65%. In comparison, the error between values derived from the Kozeny-Carman equation exhibited a more uniform distribution, albeit higher than Ergun’s law, averaging around 11%.

Figure 6
Pressure drop variation against the height of layers.

Displayed in Fig. 7a is the velocity field pattern within the median plane (z=75 µm) for the even layers. Due to the laminar nature of the flow (Re_flow=1.13), fluid movement exhibited well-organized behavior from the inlet to the point of interaction with particles. Upon collision with particles, a disturbance occurred; however, this disturbance was insufficient to induce a transition to turbulent flow. Despite certain fluid portions undergoing velocity fluctuations, specific points within the flow exhibited velocities 20 times that of the inlet velocity. Observing Fig. 7a, it became apparent that the boundary layer remained attached downstream of the particles, a characteristic hallmark of laminar flow. A similar behavior was observed in the work conducted by Reddy and Joshi²⁴. This pattern indicated that the flow velocity was elevated within the interstitial spaces between particles. This finding aligned with the research conducted by Chen et al.²², which highlighted that the highest velocities were also observed within these ‘voids’. This pattern adhered to Bernoulli’s principle.

Figure 7
Velocity field (a) and pressure variation (b) of layers 2, 4, 6, 8, and 10 in the medium plane (z=75 µm).

The pressure field depicted in Fig. 7b revealed the relationship between the number of layers and the pressure behavior within the system. As additional layers were added, the pressure within the system upstream of the particulate phase increased. In the scenario featuring 10 layers, the pressure attained a maximum value of 497.66 Pa within the region preceding the particle region. When compared to the theoretical equations, the results exhibited an average error of approximately 7% when measured against Ergun’s law. In contrast, when compared to the Kozeny-Carman equation, a more uniform behavior was observed, albeit with a higher error of around 10.53%. Within the Kozeny-Carman equation, a correction was incorporated through a constant. Typically, this constant was chosen as 180, as noted by Choi et al.¹⁵ in the most prevalent form of the equation. However, Rong et al.²⁵ suggest that this constant can vary between 150.7 and 180, contingent on the Reynolds number (Re). By inversely assessing this constant and applying the Δp value obtained from simulations, a convergence pattern emerged towards a value of 159.4. This value aligned with the assertion made by Rong et al.²⁵ and appropriately fit the values ascertained from the simulations.

As per the drag equation, the sole variable influencing it was velocity. As illustrated earlier, in pores, the velocity surpassed the inflow velocity. This difference in velocity was responsible for the mean drag force value in the particle simulations exceeding the theoretically calculated value. An important aspect to consider is that the velocity did not experience substantial variation across most of the velocity fields for the studied case as layers increased. Consequently, the drag force exhibited minor variation with the expansion of layers. Nonetheless, the key observation was that the magnitude of the drag force induced by the gas flow on a powder bed surpassed the gravitational force acting on the particles in an environment devoid of gas flow and under standard gravity (g=9.81 m/s²). Fig. 8 displays the calculated drag force values on the particles for each layer. For each layer, the average drag force (F_d) is represented (Eq. F). The drag force computed using the Stokes equation was labeled as F_dStokes (Eq. I). To provide a basis for comparison, the gravitational force value (F_g) was also included (Eq. J). The results obtained in this study indicated the lack of necessity for gravitational force in this type of system. This revelation opens new possibilities for the technology under investigation, even in environments with gravity, such as the potential for upside-down printing, for instance.

Figure 8
Theoretical gravitational force (F_g) compared with the simulated drag force (F_d).

CONCLUSIONS

The analysis of pressure and velocity fields was conducted in correlation with the increasing number of layers. Notably, a linear relationship was observed for the pressure difference across the particulate phase with the addition of layers. The findings yielded an average error of 6.65% for the Ergun equation and 10.53% for the Kozeny-Carman equation. This study effectively provided simulations whose results aligned closely with the theoretical values derived from established equations within this research domain. Additionally, a Kozeny-Carman constant value of 159.4 was determined. The drag force exerted on particles within an environment devoid of gravitational force surpassed the gravitational force experienced by the particles in a 1g environment (g=9.81 m/s²). As a result, the application investigated emerges as a feasible solution for utilizing 3D printing in space, eliminating concerns about material structural integrity due to particle non-compaction stemming from gravitational absence. The employed methodology proved efficient in studying gas flows through particulate media, particularly applicable to 3D printing within microgravity environments. Demonstrating minimal error, the methodology enabled accurate prediction of gas hydrodynamic behavior and pressure field characteristics within the particulate medium.

ACKNOWLEDGMENT

The authors gratefully acknowledge the CAPES for financial support.

REFERENCES

• ¹
  A. Zocca, P. Colombo, C.M. Gomes, J. Günster, J. Am. Ceram. Soc. 98 (2015) 1983.
• ²
  M. Ansari, E. Jabari, E. Toyserkani, J. Mater. Process. Technol. 294 (2021) 117117.
• ³
  B. Blakey-Milner, P. Gradl, G. Snedden, M. Brooks, J. Pitot, E. Lopez, M. Leary, F. Berto, A. du Plessis, Mater. Des. 209 (2021).
• ⁴
  F. Jose, M. Rivera, A. Jose, R. Arciniegas, Int. J. Adv. Manuf. Technol. 109 (2020) 17.
• ⁵
  Y.A. Gueche, N.M. Sanchez-Ballester, S. Cailleaux, B. Bataille, I. Soulairol, Pharmaceutics 13 (2021) 1212.
• ⁶
  J. Jhabvala, E. Boillat, C. André, R. Glardon, Int. J. Adv. Manuf. Technol. 59 (2012) 137.
• ⁷
  A. Hussein, L. Hao, C. Yan, R. Everson, P. Young, J. Mater. Process. Technol. 213 (2013) 1019.
• ⁸
  A. Zocca, C.M. Gomes, T. Mühler, J. Günster, Adv. Mech. Eng. 6 (2014) 491581.
• ⁹
  S. Kumar, JOM 55 (2003) 43.
• ¹⁰
  J. Zhang, W. Zheng, J.-M. Wu, K.-B. Yu, C.-S. Ye, Y.-S. Shi, Ceram. Int. 48 (2022) 1173.
• ¹¹
  M. Sen, D. Barrasso, R. Singh, R. Ramachandran, Processes 2 (2014) 89.
• ¹²
  A. Zocca, J. Lüchtenborg, T. Mühler, J. Wilbig, G. Mohr, T. Villatte, F. Léonard, G. Nolze, M. Sparenberg, J. Melcher, K. Hilgenberg, J. Günster, Adv. Mater. Technol. 4 (2019) 1900506.
• ¹³
  M. Combarros, H.J. Feise, H. Zetzener, A. Kwade, Particuology 12 (2014) 25.
• ¹⁴
  CD-adapco, “User Guide STAR-CCM+v.11.04” (2016).
• ¹⁵
  D. Choi, S. Park, J. Han, M.Y. Ahn, Y. Lee, Y.H. Park, S. Cho, D. Sohn, Fusion Eng. Des. 143 (2019) 24.
• ¹⁶
  R.N. Xu, P.X. Jiang, Int. J. Heat Fluid Flow 29 (2008) 1447.
• ¹⁷
  G. Boccardo, D.L. Marchisio, R. Sethi, J. Colloid Interface Sci. 417 (2014) 227.
• ¹⁸
  C. Chan-Braun, M. García-Villalba, M. Uhlmann, Phys. Fluids 25 (2013) 75103.
• ¹⁹
  K. Vollmari, T. Oschmann, S. Wirtz, H. Kruggel-Emden, Powder Technol. 271 (2015) 109.
• ²⁰
  S. Ergun, Chem. Eng. Prog. 48 (1952) 89.
• ²¹
  S. Flaischlen, G.D. Wehinger, Chem Engineering 3 (2019) 1.
• ²²
  Y. Chen, L. Chen, S. Liu, G. Luo, Fusion Eng. Des. 114 (2017) 84.
• ²³
  L.G. Gibilaro, Fluidization dynamics, Elsevier (2001).
• ²⁴
  R.K. Reddy, J.B. Joshi, Chem. Eng. Res. Des. 86 (2008) 444.
• ²⁵
  L.W. Rong, K.J. Dong, A.B. Yu, Chem. Eng. Sci. 99 (2013) 44.

Publication Dates

• Publication in this collection
  22 Apr 2024
• Date of issue
  Oct-Dec 2023

History

• Received
  26 June 2023
• Reviewed
  14 Sept 2023
• Accepted
  20 Oct 2023