This study aims to provide understanding of the influence of external factors, such as gravity, during sintering of three dimensional (3D)-printed parts in which the initial relative density and cohesion between the powder particles are lower compared with those present in the green parts produced by traditional powder technologies. A developed model is used to predict shrinkage and shape distortion of 3D-printed powder components at high sintering temperatures.
Three cylindrical shape connector geometries are designed, including horizontal and vertical tubes of different sizes. Several samples are manufactured by binder jetting to validate the model, and numerical results are compared with the measurements of the sintered shape.
Simulations are consistent with empirical data, proving that the continuum theory of sintering can effectively predict sintering deformation in additively manufactured products.
This work includes the assessment of the accuracy and limits of a multiphysics continuum mechanics–based sintering model in predicting gravity-induced distortions in complex-shaped additively manufactured components.
1. Introduction
The term “additive manufacturing” (AM) technology groups different processes that are described in the ISO and BS ASTM 52900: 2015, 52900: 2015 (2021) as “process of joining materials to make parts from three dimensional model data, usually layer upon layer, as opposed to subtractive manufacturing and formative manufacturing methodologies.” Among the seven categories classified, one of the processes that has received great interest from both the industrial and research points of view is binder jetting (BJ) (Dahmen et al., 2020; Nastac and Klein, 2017; Lee et al., 2020; Mostafaei et al., 2021; Ziaee and Crane, 2019; Shakor et al., 2022; Chen et al., 2022).
BJ was developed in 1993 at the Massachusetts Institute of Technology (Sachs et al., 1993). In BJ, a liquid binder is selectively ejected over a powder spread layer by layer to induce the joining between the particles. This low-energy process offers significant advantages, including the ability to print a wide range of material classes, such as metals, ceramics, and composites. Additionally, it achieves a relatively high build rate compared to other AM technologies (Martin et al., 2021; Cooke et al., 2020; Crane, 2020). Because the process is conducted at room temperature, it avoids all those unwanted effects typically present in processes like laser and electron beam melting, where the melting and rapid solidification can induce distortions, cracks and reactions (Nguyen et al., 2002; Sachs et al., 1990; Ziaee et al., 2017; Moat et al., 2011; Yu et al., 2019). The printed component composed of particles bound by a binder is typically called a “green” part, which, to achieve the final mechanical, physical and dimensional properties, requires the postprocessing steps: debinding and sintering. During conventional sintering, the densification of the part is achieved, and as widely discussed in the literature, it depends on the powder morphology, initial density distribution, printing and sintering conditions, as summarized in Bordia et al. (2017) work. During sintering and densification of the parts produced by binder jetting, it is difficult to control and predict the shrinkage and distortion of printed parts. It is widely shown in the literature that shrinkage in BJ parts shows anisotropic dimensional changes, as shown by the works of Cabo Rios et al. (2022a, 2022b, 2022c) and Jamalkhani et al. (2022) analyzing parts made of stainless steel 316 L, and Emanuelli et al. (2023) for 17–4 PH stainless steel. Therefore, considering its influence on the final part is critical and has to be taken into account in the design of the components manufactured by binder jetting technology, as shown by the work of Zago et al. (2021a, 2021b).
This anisotropic densification is mainly conducted by the “layered” pores/particles’ structure obtained during printing and by the particle rearrangement during spreading, which influence the accuracy of the printed part (Maximenko et al., 2021; Maidaniuk et al., 2022) as well the shrinkage along the different directions (Barthel et al., 2021; Lee et al., 2022).
Another phenomenon that is less studied but is important for the final outcome of real printed parts is the role played by gravity to induce distortions during sintering (German et al., 2023; Torresani et al., 2022; Borujeni et al., 2022; Paudel et al., 2023). This phenomenon is difficult to capture by observing the typical specimens (e.g. small cubes and cylinders) which are used to study this technology; these distortions are visible when the specimens have overhanging elements such as beams (Torresani et al., 2022; Borujeni et al., 2022; Paudel et al., 2023). Therefore, to achieve precise manufacturing, these deformations developed during sintering need to be taken into account in the initial design prior to the printing process.
In literature, many trial-and-error studies are present, where experimental data obtained at different printing and sintering conditions are extrapolated (Wang and Zhao, 2017; Bai et al., 2017 Huber et al.,2021; Lecis et al., 2021); these data aim to find the conditions to minimize or compensate for the deformation. Other works (Zago et al., 2021a, 2021b, 2022a, 2022b) derived theoretical models to predict the dimensional accuracy of the samples with respect to the building direction. This approach can be expensive and time-consuming because many trial-and-error and compensation iterations may be needed. Therefore, numerical simulations can be a cost- and time-effective method to redesign the parts compensating for the distortions.
The introduction of Finite Element Method (FEM) between the end of 1950s and the beginning of 1960s (Turner et al., 1956; Zienkiewicz and Hollister, 1965; Greene et al., 1969; Finlayson, 1974; Eduardo and De Arantes, 1968) induced, starting from 1970s, the development of the field defined as continuum mechanics of powder pressing to describe the compaction of porous bodies stemming from the theory of plasticity of porous bodies (Kuhn, 1971; Shima and Oyane, 1976; Green, 1972). Successively, since the 1980s, with the growth of computational power, the continuum mechanics approach started to be used to model the sintering process, where the diffusion-flow mechanism involved during the densification correlates to the macroscopic viscous flow of the sintering body, either crystalline or amorphous (Olevsky, 1998; Hsueh et al., 1986; Abouaf et al., 1988; Riedel and Blug, 2000). Models based on these laws have been largely used (Olevsky, 1998; Boccaccini and Olevsky, 1997; Olevsky et al., 2012 and Olevsky et al., 2013; Alvarado-Contreras et al., 2013; Alvarado-Contreras et al., 2014; Giuntini et al., 2016; Manière et al., 2017; Lee et al., 2022; Torresani et al., 2022; Sahli et al., 2015; Van Nguyen et al., 2016; Van Nguyen et al., 2017; Braginsky et al., 2005; Kraft and Riedel, 2004), thanks to their versatility and capability to describe part’s shrinkage and evolution during the sintering.
In the present work, the influence of gravity on the densification, dimensional changes and macrostructural evolution developed during sintering of binder jetted complex pipe connectors parts made by 316 L stainless steel have been investigated. The model framework based on the continuum mechanics of sintering developed to study densification of binder jetted simple cubic samples proposed by Cabo Rios et al. (2022a, 2022b, 2022c) and Cabo Rios et al. (2023) has been implemented in FEM software and used to simulate the evolution of the complex-shaped components subjected to the gravitational force during sintering. The model results have been compared with those obtained after sintering, to assess the model’s capability to describe the final dimensions and distortion with adequate accuracy of complex-shaped parts.
2. Materials and methods
2.1 Parts printing and sintering processes
In this study, 316 L stainless steel powder was used. The particle size distribution is 27 μm (D90), 15.3 μm (D50) and 8.1 μm (D10), as declared by Digital Metals powder provider. Three connector geometries were designed having different dimensions but consisting of three main features: a basement and a vertical tube intersecting a horizontal tube. The three geometries from here on are called G1, G2 and G3. Figure 1 displays the frontal and side view of the geometry 1 (G1), indicating the most relevant dimensions of the components. Since the other two geometries have the same features, the equivalent relevant dimensions were measured.
Designed tee connector geometry, indicating the relevant dimensions of the component
Designed tee connector geometry, indicating the relevant dimensions of the component
Tables 1 and 2 report the linear and diameter nominal dimensions at green state of the three connector geometries, referring to Figure 1. The geometries were scaled to compensate for the expected dimensional shrinkages during sintering. Different scaling factors were applied along X, Y and Z printing directions, according to the suggestions of machine provider.
Scaled linear dimensions of CAD geometries in millimeter
| Geometry | LX [mm] | LY [mm] | LZ [mm] | W [mm] |
|---|---|---|---|---|
| G1 | 83.393 | 31.079 | 63.706 | 40.505 |
| G2 | 83.393 | 31.079 | 61.302 | 40.505 |
| G3 | 59.567 | 7.172 | 33.656 | 16.679 |
| Geometry | LX [mm] | LY [mm] | LZ [mm] | W [mm] |
|---|---|---|---|---|
| G1 | 83.393 | 31.079 | 63.706 | 40.505 |
| G2 | 83.393 | 31.079 | 61.302 | 40.505 |
| G3 | 59.567 | 7.172 | 33.656 | 16.679 |
Scaled “diameters” of CAD geometries in millimeter
| Geometry | ||||||||
|---|---|---|---|---|---|---|---|---|
| G1 | 31.079 | 31.252 | 27.493 | 27.646 | 30.975 | 31.079 | 27.401 | 27.493 |
| G2 | 31.079 | 31.252 | 27.493 | 27.646 | 11.913 | 11.953 | 9.531 | 9.563 |
| G3 | 7.172 | 7.212 | 4.781 | 4.808 | 7.148 | 7.172 | 4.765 | 4.781 |
| Geometry | ||||||||
|---|---|---|---|---|---|---|---|---|
| G1 | 31.079 | 31.252 | 27.493 | 27.646 | 30.975 | 31.079 | 27.401 | 27.493 |
| G2 | 31.079 | 31.252 | 27.493 | 27.646 | 11.913 | 11.953 | 9.531 | 9.563 |
| G3 | 7.172 | 7.212 | 4.781 | 4.808 | 7.148 | 7.172 | 4.765 | 4.781 |
The anisotropic scaling factor determined an elliptical section of the horizontal and vertical tubes, consisting of two different axis lengths, which are identified with the respective scaling axis in Table 2 (e.g., ellipse z-axis of internal horizontal cylinder ).
In this study X, Y and Z directions correspond to the binder injection, powder spreading and building direction, respectively. For each geometry, four replicates were printed by a Digital Metal (P1701) machine. Samples were located in the printing box having the horizontal axis of the tube aligned with binder injection direction (X). The axis of the vertical tube was aligned with the building direction (Z). Printing parameters are reported in Table 3.
Printing parameters
| Layer thickness | Printing speed | Powder applicator speed | Bed temperature | Resolution |
|---|---|---|---|---|
| 42 μm | 200 mm/s | 30 mm/s | 80 °C | 1,200 DPI |
| Layer | Printing | Powder | Bed | Resolution |
|---|---|---|---|---|
| 42 μm | 200 mm/s | 30 mm/s | 80 °C | 1,200 DPI |
The binder saturation level was controlled by the Dark Body parameter, which is set equal to 3, corresponding to the saturation of the 69% of the pixels (Dahmen et al., 2021). Because measuring green density of complex shapes is challenging, 10-mm cubical samples were printed. For these geometries, an average green density of 4.38 g/cm3 was calculated using the mass-to-volume ratio.
After the printing process, the building box was heat treated to evaporate the aqueous fraction content of the binder and to cure the polymeric fraction. Subsequently, de-powdering was carefully performed to extract the green products. After the dimensional and geometrical characterization of green products, all samples were thermically debinded in air atmosphere at 640°C, and successively sintered in a batch furnace under hydrogen. An isothermal sintering temperature of 1,377°C was reached during the sintering cycle which was held for 210 min.
2.2 Parts characterization and measurements
Each sample was measured by a coordinate measuring machine (CMM), before and after sintering, aiming at reconstructing the whole geometry. A global DEA image 07-07-07 CMM machine was used, equipped with an SP600 Renishaw head mounting a stylus of 6 mm of diameter to measure connectors 1 and 2. A stylus of 1 mm diameter was used for geometry 3. The machine was calibrated according to ISO 10360-2, guaranteeing a maximum permissible error of 1.5 + L/333 μm.
An automatic measurement routine was developed. In the beginning, the Datum Reference Frame (DRF) was defined. The axis of the external cylinder of the horizontal tube was used to establish the x-axis of the part. The axis of the external cylinder of the vertical tube was used to align the z-axis. The intersection of the two axes determines the center of the DRF. After the definition of the DRF, the routine was programmed to reconstruct the part geometry. Planar surfaces were inspected by a point-by-point strategy. Using the measured points, planes were reconstructed by a best-fit least squares method. Using the distance between the reconstructed planes, linear dimensions were estimated. Subsequently, horizontal and vertical tubes were reconstructed. A parametric procedure was programmed to measure the external cylinder at three levels equally spaced along the distance H as shown in Figure 2. Both left and right sides of horizontal cylinder, with respect to the axis of the vertical one, were measured at three levels. The same strategy was adopted for internal cylinder to measure circle sections at four levels. By this parametric approach, it is possible, in principle, to measure the same cylinder section both at green and sintered state.
Illustration of the measurement levels along the vertical and horizontal cylinders
Illustration of the measurement levels along the vertical and horizontal cylinders
In total, 24 points were acquired at each level for reconstructing each single closed curve. Points were acquired around the circumference at equidistant angular position. Each ellipse’s circumference was reconstructed by a best-fit least squares algorithm as shown in Figure 3. At sintered state, the measured points of distorted cylinder were fitted by equation (1):
where y0 and z0 are the coordinates of the intersection of the ellipse’s axes, a and b are the major and minor semi-axis of the ellipse, ζ is the ellipse rotation angle around the x-axis. In this way, shape distortion was experimentally derived. During this work, the major and minor semi-axes were identified by the closest parallel direction of the component as can be seen in Figure 3(a).
Assuming that circle profile had evolved into an elliptical shape, measured points were fitted to derive the minor and major axes at three levels: 3H/4, 2H/4 and H/4 [see Figure 3 (a)]. In parallel, the sintering simulation allows monitoring the diameter variations [see Figure 3(b)] during sintering as a function of the time and temperature. Consequently, the experimental CMM measurements were compared against the simulation results at the end of sintering (t = 755 min).
The final density of each sintered part was measured through Archimedes method.
3. Sintering modeling
In the present macroscale model, through the coupling between thermal, mechanical (stress-strain) and densification was possible to model the influence of the sintering parameters (e.g. thermal cycle) and external conditions (e.g. gravity and friction) to determine the final densification and deformation of the studied components. The model, which is based on the continuum theory of sintering (Olevsky, 1998; Cabo Rios et al., 2023), was incorporated in COMSOLTM finite elements software and its framework and boundary conditions are described below.
3.1 Main constitutive equation
The main constitutive equation for the free sintering process is defined as:
in which σij is the stress tensor (Pa) and represents the externally applied stresses, is the first strain rate tensor invariant (s−1), which corresponds to the volumetric shrinkage, η0 is the shear viscosity of the fully dense material constituting the powder (Pa〮s), and δij is the Kronecker delta. The normalized shear and bulk viscosity modulus, φ and ψ, respectively, and the effective sintering stress PL (Pa), can be written as function of porosity θ according to the modified Skorokhod–Olevsky model (Cabo Rios et al., 2023):
In the last equation, the effective sintering stress or Laplace pressure is proportional to the material’s surface tension α (J〮m−2) because this represents the driving force for the sintering, and is inversely proportional to the average grain radius, rg (m).
The strain rate tensor (s−1) is given in terms of matrix velocity gradient as:
Assuming isotropy condition, the strain rate tensor can be rewritten as:
The first invariant of the strain rate tensor corresponds to the rate of volume change in the porous material, and can be correlated to the porosity evolution with the continuity equation:
As an external load, the gravity force has been introduced; therefore, the equilibrium equation can be written as:
where fi is the force per unit of mass, given in terms of density (kg m−3) of the part and gravity acceleration vector (m s−2).
3.2 Grain growth and delta-ferrite kinetics
To evolution of the grain growth during the sintering process can be described using the kinetic equation proposed by (Olevsky et al., 2012; Cabo Rios et al., 2023):
In this nonlinear equation, the dependence of the grain growth kinetic on the temperature is described by the pre-exponential term k0 (m3 s−1) and activation energy Qg (kJ mol−1), and the dependence from the porosity evolution on the presence of the critical porosity θc. This last parameter represents the limit when the pore depinning from the grain boundary is relevant in influencing both grain growth and densification kinetics. These parameters were experimentally determined in a previous work (Cabo Rios et al., 2023) and are shown in Table 4.
Parameters for equation (9) identified for binder jetted samples made by 316 L (Cabo Rios et al., 2023)
| k0 (m3 s−1) | Qg (kJ mol−1) | θc |
|---|---|---|
| 2.97〮10−22 | 164.8 | 5.20% |
| k0 (m3 s−1) | Qg (kJ mol−1) | θc |
|---|---|---|
| 2.97〮10−22 | 164.8 | 5.20% |
A previously proposed Arrhenius-like piecewise function was used for the material viscosity η0 parameter, where a transition temperature is defined for the evolution between two temperature domains (Cabo Rios et al., 2022a, 2022b, 2022c):
The transition temperature from the experimentally derived material viscosity agreed with the delta-ferrite transition temperature from thermodynamic calculations. Furthermore, the model showed good agreement with the densification behavior characterized by dilatometry experiments below and above the transformation temperature, where the increase of shrinkage rate was accurately described by the model.
3.3 Boundary conditions
Proper motion and equilibrium conditions are applied to the outer boundaries in the x-y-z reference system. Through the interaction between the lower surfaces of the components with a support were imposed the conditions:
where σr is the surface load opposing any motion during the densification or slumping, which is corelated to the normal surface load σz, induced by the gravitational acceleration, through the coefficient of Coulomb friction μ. In this preliminary study the influence of the friction has been considered negligible (μ = 0).
3.4 Experimental determination of sintering model parameters
The sintering model parameters were determined through dilatometric test as shown in previous works (Cabo Rios et al., 2022a, 2022b, 2022c, 2023). Through the dilatometric test the effect of gravity can been neglected, which is reasonable due to the small dimensions of the specimens. Therefore, equation (2) can be rewritten for a pressure-less isotropic case:
Combining equation (14) with the continuity equation expressed by equation (8) allows to derive the relationship, which describes the porosity evolution during sintering:
This equation describes the porosity variation as function of porosity θ, effective stress PL, normalized bulk viscosity ψ and material shear viscosity η0. As mentioned above, PL and ψ are function of the porosity. The powder material shear viscosity η0 can be described by an Arrhenius-like function to represent its temperature-dependent nature, and includes the influence of the delta-ferrite phase formation on the powder material behavior by introducing the piece-wise function from equation (11), where the pre-exponential factor A0 (Pa s K−1) and activation energy Q (kJ mol−1) are materials constant and should be determined. These constants were previously obtained from similar 316 L BJ samples printed using the same BJ and powder system and sintered at different heating rates (from 2 to 15°C/min) (Cabo Rios et al., 2022a, 2022b, 2022c). Dilatometry results at different heating rates showed negligible influence of the δ-ferrite transformation kinetics on the material viscosity. Thus, the powder material shear viscosity is now defined for the complete range of the temperature (up to 1,370°C) used for the sintering of the component studied in the present work. The constants used are detailed in Table 5.
Parameters for equation (9) identified for binder jetted samples made by 316 L (Cabo Rios et al., 2023)
| i | Ai/α [s m−1K−1] | Qi [K J mol−1] | TT [K] |
|---|---|---|---|
| 1 | 8.993 | 201.7 | 1,583.62 |
| 2 | 5.335e-32 | 1,178.7 |
| i | Ai/α [s m−1K−1] | Qi [K J mol−1] | TT [K] |
|---|---|---|---|
| 1 | 8.993 | 201.7 | 1,583.62 |
| 2 | 5.335e-32 | 1,178.7 |
4. Results and discussion
4.1 Sintering simulation results and qualitative comparison with sintered components
The finite element mesh chosen for this study consists of ten-node tetrahedron elements with the following characteristics: maximum element size of 9.64 mm, minimum element size of 1.6 mm, maximum element growth rate of 1.6, curvature factor of 0.7 and resolution in narrow regions of 0.4. The meshing strategy for the convergence of the simulation has followed the method from Grippi et al. (2024) with ten-nodes tetrahedron. The finite element simulation computes the sintering for 755 min with 5-min time steps.
As mentioned previously, gravity plays an important role in shape deviation, and the simulation must consider gravity as an external load. A preliminary study was conducted to compare the same part with and without gravity. The visual comparison, shown in Figure 4, clearly demonstrates that if the FEM model can compute densification, it does not display any slumping or shape distortion of the studied part without the inclusion of gravity in the simulation.
Shape distortion of the geometry 1 (G1) vertical and horizontal tube: (a and c) numerical result without computation of the gravity and (b and d) numerical result under gravity
Shape distortion of the geometry 1 (G1) vertical and horizontal tube: (a and c) numerical result without computation of the gravity and (b and d) numerical result under gravity
The final shapes of the three different geometries, as by experiment and simulation consecutively, can be observed in Figures 5(G1), 6(G2) and 7(G3). In these figures, on the simulation side, the sintered shape is colored expressing the distribution of the relative density at the end of sintering. Also, the shape evolution of external and internal profiles from simulations is reported at different depth positions of horizontal and vertical tube (following method detailed in Figure 3). Monitoring the dimensions of the different ellipse’s axis, the shape deformation can be clearly highlighted.
Shape distortion of the geometry 1 (G1) vertical and horizontal tube: (a and c) photo, (b and d) numerical result, and (e) evolution during sintering simulation of the external and internal axis dimensions
Shape distortion of the geometry 1 (G1) vertical and horizontal tube: (a and c) photo, (b and d) numerical result, and (e) evolution during sintering simulation of the external and internal axis dimensions
Shape distortion of the geometry 2 (G2) vertical and horizontal tube: (a and c) photo, (b and d) numerical result and (e) evolution during sintering simulation of the external and internal axis dimensions
Shape distortion of the geometry 2 (G2) vertical and horizontal tube: (a and c) photo, (b and d) numerical result and (e) evolution during sintering simulation of the external and internal axis dimensions
Shape distortion of the geometry 3 (G3) vertical and horizontal tube: (a and c) photo, (b and d) numerical result and (e) evolution during sintering simulation of the external and internal axis dimensions
Shape distortion of the geometry 3 (G3) vertical and horizontal tube: (a and c) photo, (b and d) numerical result and (e) evolution during sintering simulation of the external and internal axis dimensions
As shown in Figure 5, due to its thin wall thickness, the first geometry has significant sintering distortions of both horizontal and vertical cylinder. Counterintuitively, the upper cylinder shows more distortion than the horizontal one as it will be thoroughly analyzed in subsection 4.2. The simulation helps to understand that the upper vertical cylinder distortion is enhanced by the lower horizontal one that applies additional tension, in the XZ plane, by its collapsing. The experimental measurements of the density for the sintered part through Archimedes’ method provided the average value of 7.89 (± 0.02) g cm−3, which corresponds to a relative density of 99.25 (±0.3) %, this value is in agreement with the modeling result shown in Figure 4(b).
The second geometry (G2) shows most of its shape distortions on the horizontal cylinder. The vertical cylinder does not show important cylindrical shape distortions. However, the concentrated weight of this upper part emphasizes the distortion of the lowest cylinder. This behavior is opposite to the deformation observed in G1 [Figure 5(e)] and shows the larger slumping of the horizontal cylinder closer to the vertical cylinder due to the concentrated weight of the vertical cylinder with lower diameter.
The relative density measurement for this part provided an average value of 7.91 (± 0.06) g cm−3, which corresponds to a relative density of 99.44 (±0.8) %, also in this case the measured values are in agreement with the model results shown in Figure 5(b).
Note that on the sintered components G1 and G2 displayed in Figures 5(a) and 6(a), the cross-sectional profile of horizontal tube exhibits a slight rotation around the cylinder axis, resulting in an additional shape distortion (Zago et al., 2022a). This phenomenon is larger in extent in G1 and can be ascribed to an inhomogeneous green density distribution along the part, which deviation leads to different densification and shrinkage on sintering. Another explanation considers a slight misalignment between the direction of the gravity load and the axis of the vertical tube, giving rise to a tangential force that induces the observed distortion. However, the simulation results do not show this effect due to the symmetry of the geometries studied, which hinder the possible small deviations during the sintering experiments detailed above. Further studies can be performed to analyze the effect of density inhomogeneities and small deviations of the vertical axis direction on the cylindrical shape deformation during sintering.
In Figure 7, the third geometry (G3) shows a negligible slumping of the cylindrical shapes caused by gravity. The vertical cylinder does not show significant shape distortion during the sintering simulation, but the initial minor differences between the ellipse’s axis dimensions applied by the geometry scaling factors were not reduced. This is due to the use of isotropic sintering model while the applied factors are anisotropic. However, the anisotropic shrinkage of these components is very small compared to typical BJ samples studied in literature such as in the works of Cabo Rios et al. (2022a, 2022b, 2022c), Zago et al. (2021a, 2021b) and Lee et al. (2022).
The measured density for this component is 7.95 (± 0.01) g cm−3, which corresponds to a relative density of 99.94 (±0.1) %, also in this case the experimental and model results [Figure 6(b)] are in agreement.
Apart from the details of the cylindrical shape deformation, the geometries experience high levels of shrinkages during the sintering simulations. The shrinkage evolution of the component’s global dimensions (LX, LY, LZ and W) is detailed in Figure 8. In general, the shrinkage of dimensions LX and W showed a direct correlation with the densification behavior, without any shrinkage occurring after full densification is reached. The shrinkage of the vertical dimension LZ increases due to densification from temperatures ∼ 900°C and keeps increasing during the holding time due to the slumping of the geometry caused by the gravitation forces. This effect is larger for the geometry G2 when compared to G1, as previously discussed in this section, while it is negligible for G3. On the contrary, the dimension LY shows a shrinkage during the densification process, which tends to decrease close to the highest temperature. Then, densification shrinkage stops, and LY increases during the holding time, as observed by the decrease on the dimensional shrinkage. This is caused by the lateral slumping of the base caused by low viscosity at the dwell temperature and the gravitational forces, which may be increased by the absence of frictional forces in the simulation.
Shrinkage evolution of the component’s global dimensions during the sintering simulations
Shrinkage evolution of the component’s global dimensions during the sintering simulations
4.2 Quantitative comparison between coordinate measuring machine experimental measurements with simulation results
A quantitative analysis of sintering shrinkages and shape distortions is presented here according to the procedure described in subsection 2.2. Figure 9 gathers the relative (%) and absolute (mm) differences between the experimental and simulation results obtained by the method depicted in Figure 3.
Nominal and percentage difference between experimental and simulation results of the axis of the ellipse shape profile at depth levels H/4, 2H/4 and 3H/4 of external (a) and internal (b) surface of horizontal tube; and external (c) and internal (d) surface of vertical tube
Nominal and percentage difference between experimental and simulation results of the axis of the ellipse shape profile at depth levels H/4, 2H/4 and 3H/4 of external (a) and internal (b) surface of horizontal tube; and external (c) and internal (d) surface of vertical tube
In general, the simulation results reveal a very good prediction accuracy of the model with all absolute deviations from experimental values below 1 mm. The largest deviation along the horizontal and vertical cylinders axis measurements are −0.96 (4.9%) mm and −0.39 (5.27%) mm, respectively. It is also noticeable the dominant presence of negative values in Figure 9, meaning that the dimensions measured in the simulation results are larger than the experiments in general. However, these deviations may come from differences in the densification shrinkage and/or cylindrical shape deformation during sintering. In particular, G3 has minimal shape deformation during sintering and all the deviations from experimental results are negative (see Figure 9), suggesting that the initial green density used was slightly overestimated. The different deviations at various depth levels in G1 and G2 geometries are more likely related to the prediction of cylindrical shape distortions rather than the densification shrinkage. Therefore, to isolate the cylinder’s shape distortions, the axis ratios were calculated and plotted in Figure 10.
Ratio between the ellipses axis dimensions from (a) simulations and (b) CMM measurements of each geometry (G1, G2 and G3)
Ratio between the ellipses axis dimensions from (a) simulations and (b) CMM measurements of each geometry (G1, G2 and G3)
The vertical cylinder axis ratios in G1 show the same increasing trend along the cylinder in the experiment and simulation results [Figure 10(a)] with largest deformation at 3H/4. However, all the values are larger and the increasing trend stronger on the sintered component (CMM). This reveals the systematically larger shape distortion of the sintered component compared to the simulation result. Also, G2 and G3 show larger values for the CMM measurements, but the distortions are minimal.
In Figure 10(b), G1 and G2 horizontal cylinder ratios derived from CMM measurements show a consistent decreasing trend from H/4 to 3H/4, revealing the largest deformation close to the vertical cylinder. This trend is stronger on the G2 component, due to the concentrated weight of the vertical cylinder. The axis ratios derived from simulation results [Figure 10(a)] reveals the same trend for the G2 geometry, but with lower values and difference between the H/4 and 3H/4 depth levels. This again indicates a larger shape distortion on the sintered components. On the contrary, the simulation showed the opposite trend along the G1 horizontal cylinder [Figure 10(a)], where largest distortions were found further from the connection between cylinders (at 3H/4). The reason of this discrepancy is under investigation.
As shown in Figure 9, the comparison between axis ratios, for both external and internal profiles, reveals similar predicted and experimental results. The best prediction is made for the less distorted geometry (Geometry 3), but simulation remains within acceptable divergence range also for the most distorted one (Geometry 1). In summary, the simulation results show consistent lower ratios and weaker trends for most components and depth levels. This lower simulation shape distortion could be caused by the overestimation of the material shear viscosity η0. Further refinement of the material constant could be done to improve the simulation performance.
Figure 11 shows the differences between experimental and simulation on linear dimensions at the final time step, and Figure 12 compares the experimental and simulated shrinkages in the three directions.
Nominal and percentage difference between experimental and simulation results of the linear global dimensions from the different geometries (G1, G2 and G3)
Nominal and percentage difference between experimental and simulation results of the linear global dimensions from the different geometries (G1, G2 and G3)
Sintering shrinkages measured for the linear global dimensions of the different samples (G1, G2 and G3), comparison between experimental and simulation shrinkages
Sintering shrinkages measured for the linear global dimensions of the different samples (G1, G2 and G3), comparison between experimental and simulation shrinkages
As shown in Figures 11 and 12, the highest accuracy for the prediction of shrinkages is observed along the X direction, where the LX dimension largest deviation is 0.3% (maximum difference 0.21 mm on the G2 component). The lower accuracy in dimensional prediction of the LY and LZ dimensions is affected by the slumping of the components due to gravity. LY shrinkages in G1 and G2 are lower for the experimental measurements (∼11%) compared to simulation results (∼15%), related to larger horizontal slumping (along Y) of the sintered components. This may be due to the overestimation of the model material shear viscosity and the absence of friction forces in the preliminary simulation studies presented in this work. Theoretically, the shrinkages of LZ dimension are consistently larger for all the components because of the sintering anisotropy. However, Figure 11 shows that LZ shrinkages are considerably larger for the G1 and G2 components compared to G3 due to the shape deformation of the horizontal tubes. The largest deviation between the LZ shrinkages was observed for the G2 component with an overestimation of 2.7% from the simulation results when compared to the experiment. The larger LZ shrinkages of the experimental measurements may be caused by a combination of shrinkage anisotropy (not considered in the isotropic model used) and the underestimation of the material shear viscosity in the model.
The measurements and model values for the various geometries and elements examined in this study are provided in the Appendix.
5. Conclusions
This paper aims to assess the accuracy of the proposed continuum mechanics-based model in predicting sintering deformation of complex-shaped metal binder jetted products, particularly the shape distortion induced by gravity load. Here are the fundamental conclusions:
Sintering simulations are consistent with experimental results. The model effectively estimates gravity-induced deformation of the stainless-steel material behavior at high sintering temperatures. In general, simulations capture the deformation interaction of the different vertical and horizontal cylinder’s geometrical features.
The largest deviation between experimental and simulation results consists of −0.96 mm (4.9%) and −0.39 mm (5.27%) on the horizontal and vertical tubes, respectively. Generally, simulations predict lower shrinkage than experimental results. A more accurate characterization of initial green density and modeling the temperature distribution and convection in the sintering furnace could potentially provide a better match between empirical and numerical results.
For the horizontal tube (cylinder axis perpendicular to gravity load), the prediction accuracy of the external and internal profiles tends to decrease moving to inner depth levels of G1 and G2 geometry, with the minor axis showing lower precision. Also, the axis ratio from CMM measurements is consistently larger than the simulation results ratio. This larger experimental distortion could be explained by an underestimation of the material softening (higher shear viscosity) utilized in the model, which results in lower slumping, especially close to the vertical tube due to its weight effect.
For the vertical tube (cylinder axis parallel to gravity load), simulations can precisely predict the severe deformation that occurred in G1. In this geometry, accuracy tends to increase moving to inner depth levels, due to the increasing difference between the axis ratio of the experimental and simulation results. This again may be related to the underestimation of the model material softening as previously described. No significant distortion occurred on the vertical tubes of G2 and G3 axes.
In conclusion, the work demonstrates that the sintering model used can effectively predict shape deformation occurring in complex sintered BJ products. Hence, the current model demonstrates its potential as a predictive tool for designing BJ parts, forecasting the deformations caused by gravity during the sintering process, and eliminating the need for time- and resource-consuming experimental trial and error. The experimental characterization of green density and the actual sintering temperature is challenging but essential for a good model prediction. Moreover, the work did not consider that dimensional and geometrical errors can be introduced during the printing operation. Friction coefficients were neglected; moreover, densification was assumed isotropic instead of anisotropic. Further study will be directed to implement such aspects to increase the accuracy of model predictions.
Declaration of competing interest: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References
Appendix
The experimental CMM measurements of the sintered specimens and simulation results, extracted at the last time step (t = 755 min), corresponding to the horizontal and vertical cylindrical sections are detailed in Tables 1 and 2, respectively. These results were acquired following the methods depicted in Figure 3, which focuses on the elliptical distortion of the hollow cylindrical profiles. The results for the external linear dimensions of the sintered geometries from CMM measurements and simulations are depicted in Table A3.
Experimental linear dimensions of the sintered geometries (G1, G2 and G3) derived from the CMM measurements
| Green (as-printed) [mm] | Sintered [mm] | |||||||
|---|---|---|---|---|---|---|---|---|
| Geometry | LX | LY | LZ | W | LX | LY | LZ | W |
| SIM | ||||||||
| G1 | 83.39 | 31.08 | 63.71 | 40.51 | 68.92 | 26.46 | 50.95 | 33.43 |
| G2 | 83.39 | 31.08 | 61.30 | 41.51 | 68.74 | 26.63 | 48.17 | 33.93 |
| G3 | 59.57 | 7.17 | 33.66 | 16.68 | 49.14 | 5.95 | 27.62 | 13.77 |
| EXP | ||||||||
| G1 | 83.43 | 31.06 | 63.75 | 40.59 | 69.11 | 27.54 | 50.29 | 33.84 |
| G2 | 83.44 | 31.10 | 61.30 | 40.60 | 68.95 | 27.36 | 46.90 | 33.75 |
| G3 | 59.59 | 7.15 | 33.60 | 16.88 | 49.14 | 5.86 | 27.06 | 13.83 |
| Green (as-printed) [mm] | Sintered [mm] | |||||||
|---|---|---|---|---|---|---|---|---|
| Geometry | LX | LY | LZ | W | LX | LY | LZ | W |
| SIM | ||||||||
| G1 | 83.39 | 31.08 | 63.71 | 40.51 | 68.92 | 26.46 | 50.95 | 33.43 |
| G2 | 83.39 | 31.08 | 61.30 | 41.51 | 68.74 | 26.63 | 48.17 | 33.93 |
| G3 | 59.57 | 7.17 | 33.66 | 16.68 | 49.14 | 5.95 | 27.62 | 13.77 |
| EXP | ||||||||
| G1 | 83.43 | 31.06 | 63.75 | 40.59 | 69.11 | 27.54 | 50.29 | 33.84 |
| G2 | 83.44 | 31.10 | 61.30 | 40.60 | 68.95 | 27.36 | 46.90 | 33.75 |
| G3 | 59.59 | 7.15 | 33.60 | 16.88 | 49.14 | 5.86 | 27.06 | 13.83 |
Axis dimensions of the horizontal cylinder of sintered geometries (G1, G2 and G3) derived from the CMM measurements and simulations
| Geometry | Profile | Axis | Horizontal cylinder | |||||
|---|---|---|---|---|---|---|---|---|
| EXP (mean ± SD) [mm] | SIM (Value at t = 755 min) [mm] | |||||||
| H/4 | 2H/4 | 3H/4 | H4 | 2H/4 | 3H/4 | |||
| G1 | Ext | Y | 27.66 ± 0.16 | 27.26 ± 0.17 | 27.58 ± 0.17 | 27.46 | 27.46 | 27.49 |
| Z | 23.21 ± 0.21 | 23.32 ± 0.21 | 23.35 ± 0.22 | 23.73 | 23.57 | 23.43 | ||
| Int | Y | 24.42 ± 0.12 | 24.40 ± 0.11 | 24.37 ± 0.09 | 24.51 | 24.52 | 24.54 | |
| Z | 20.04 ± 0.22 | 20.16 ± 0.21 | 20.20 ± 0.21 | 20.76 | 20.60 | 20.46 | ||
| G2 | Ext | Y | 28.07 ± 0.05 | 27.85 ± 0.07 | 27.63 ± 0.08 | 27.62 | 27.55 | 27.46 |
| Z | 22.75 ± 0.06 | 23.04 ± 0.05 | 23.24 ± 0.03 | 23.45 | 23.49 | 23.51 | ||
| Int | Y | 24.78 ± 0.05 | 24.63 ± 0.07 | 24.47 ± 0.07 | 24.67 | 24.667 | 24.51 | |
| Z | 19.53 ± 0.06 | 19.86 ± 0.05 | 20.07 ± 0.04 | 20.48 | 20.48 | 20.55 | ||
| G3 | Ext | Y | 5.93 ± 0.03 | 5.92 ± 0.03 | 5.94 ± 0.04 | 5.95 | 5.95 | 5.96 |
| Z | 5.82 ± 0.02 | 5.82 ± 0.01 | 5.82 ± 0.02 | 5.89 | 5.89 | 5.88 | ||
| Int | Y | 3.91 ± 0.03 | 3.90 ± 0.04 | 3.89 ± 0.08 | 3.98 | 3.98 | 3.98 | |
| Z | 3.80 ± 0.01 | 3.81 ± 0.01 | 3.79 ± 0.04 | 3.92 | 3.91 | 3.91 | ||
| Geometry | Profile | Axis | Horizontal cylinder | |||||
|---|---|---|---|---|---|---|---|---|
| EXP (mean ± SD) [mm] | SIM (Value at t = 755 min) [mm] | |||||||
| H/4 | 2H/4 | 3H/4 | H4 | 2H/4 | 3H/4 | |||
| G1 | Ext | Y | 27.66 ± 0.16 | 27.26 ± 0.17 | 27.58 ± 0.17 | 27.46 | 27.46 | 27.49 |
| Z | 23.21 ± 0.21 | 23.32 ± 0.21 | 23.35 ± 0.22 | 23.73 | 23.57 | 23.43 | ||
| Int | Y | 24.42 ± 0.12 | 24.40 ± 0.11 | 24.37 ± 0.09 | 24.51 | 24.52 | 24.54 | |
| Z | 20.04 ± 0.22 | 20.16 ± 0.21 | 20.20 ± 0.21 | 20.76 | 20.60 | 20.46 | ||
| G2 | Ext | Y | 28.07 ± 0.05 | 27.85 ± 0.07 | 27.63 ± 0.08 | 27.62 | 27.55 | 27.46 |
| Z | 22.75 ± 0.06 | 23.04 ± 0.05 | 23.24 ± 0.03 | 23.45 | 23.49 | 23.51 | ||
| Int | Y | 24.78 ± 0.05 | 24.63 ± 0.07 | 24.47 ± 0.07 | 24.67 | 24.667 | 24.51 | |
| Z | 19.53 ± 0.06 | 19.86 ± 0.05 | 20.07 ± 0.04 | 20.48 | 20.48 | 20.55 | ||
| G3 | Ext | Y | 5.93 ± 0.03 | 5.92 ± 0.03 | 5.94 ± 0.04 | 5.95 | 5.95 | 5.96 |
| Z | 5.82 ± 0.02 | 5.82 ± 0.01 | 5.82 ± 0.02 | 5.89 | 5.89 | 5.88 | ||
| Int | Y | 3.91 ± 0.03 | 3.90 ± 0.04 | 3.89 ± 0.08 | 3.98 | 3.98 | 3.98 | |
| Z | 3.80 ± 0.01 | 3.81 ± 0.01 | 3.79 ± 0.04 | 3.92 | 3.91 | 3.91 | ||
Axis dimensions of the vertical cylinder of sintered geometries (G1, G2 and G3) derived from the CMM measurements and simulations
| Geometry | Profile | Axis | Vertical cylinder | |||||
|---|---|---|---|---|---|---|---|---|
| EXP (mean ± SD) [mm] | SIM (Value at t = 755 min) [mm] | |||||||
| H/4 | 2H/4 | 3H/4 | H4 | 2H/4 | 3H/4 | |||
| G1 | Ext | X | 27.50 ± 0.15 | 27.93 ± 0.07 | 28.37 ± 0.07 | 27.51 | 27.73 | 27.96 |
| Y | 23.42 ± 0.15 | 22.95 ± 0.15 | 22.46 ± 0.15 | 23.48 | 23.18 | 22.84 | ||
| Int | X | 24.29 ± 0.10 | 24.60 ± 0.12 | 24.94 ± 0.15 | 24.57 | 24.79 | 25.02 | |
| Y | 20.12 ± 0.18 | 19.72 ± 0.14 | 19.18 ± 0.14 | 20.53 | 20.23 | 19.89 | ||
| G2 | Ext | X | 9.98 ± 0.05 | 9.95 ± 0.03 | 9.92 ± 0.03 | 9.91 | 9.85 | 9.81 |
| Y | 9.64 ± 0.07 | 9.70 ± 0.07 | 9.72 ± 0.07 | 9.75 | 9.80 | 9.82 | ||
| Int | X | 7.57 ± 0.49 | 7.51 ± 0.50 | 7.45 ± 0.50 | 7.94 | 7.88 | 7.85 | |
| Y | 7.63 ± 0.05 | 7.69 ± 0.05 | 7.72 ± 0.05 | 7.78 | 7.83 | 7.85 | ||
| G3 | Ext | X | 5.87 ± 0.04 | 5.89 ± 0.03 | 5.92 ± 0.05 | 5.92 | 5.92 | 5.92 |
| Y | 5.80 ± 0.01 | 5.82 ± 0.02 | 5.85 ± 0.05 | 5.90 | 5.90 | 5.90 | ||
| Int | X | 3.82 ± 0.02 | 3.82 ± 0.03 | 3.80 ± 0.04 | 3.95 | 3.95 | 3.94 | |
| Y | 3.86 ± 0.01 | 3.87 ± 0.03 | 3.85 ± 0.03 | 3.93 | 3.93 | 3.93 | ||
| Geometry | Profile | Axis | Vertical cylinder | |||||
|---|---|---|---|---|---|---|---|---|
| EXP (mean ± SD) [mm] | SIM (Value at t = 755 min) [mm] | |||||||
| H/4 | 2H/4 | 3H/4 | H4 | 2H/4 | 3H/4 | |||
| G1 | Ext | X | 27.50 ± 0.15 | 27.93 ± 0.07 | 28.37 ± 0.07 | 27.51 | 27.73 | 27.96 |
| Y | 23.42 ± 0.15 | 22.95 ± 0.15 | 22.46 ± 0.15 | 23.48 | 23.18 | 22.84 | ||
| Int | X | 24.29 ± 0.10 | 24.60 ± 0.12 | 24.94 ± 0.15 | 24.57 | 24.79 | 25.02 | |
| Y | 20.12 ± 0.18 | 19.72 ± 0.14 | 19.18 ± 0.14 | 20.53 | 20.23 | 19.89 | ||
| G2 | Ext | X | 9.98 ± 0.05 | 9.95 ± 0.03 | 9.92 ± 0.03 | 9.91 | 9.85 | 9.81 |
| Y | 9.64 ± 0.07 | 9.70 ± 0.07 | 9.72 ± 0.07 | 9.75 | 9.80 | 9.82 | ||
| Int | X | 7.57 ± 0.49 | 7.51 ± 0.50 | 7.45 ± 0.50 | 7.94 | 7.88 | 7.85 | |
| Y | 7.63 ± 0.05 | 7.69 ± 0.05 | 7.72 ± 0.05 | 7.78 | 7.83 | 7.85 | ||
| G3 | Ext | X | 5.87 ± 0.04 | 5.89 ± 0.03 | 5.92 ± 0.05 | 5.92 | 5.92 | 5.92 |
| Y | 5.80 ± 0.01 | 5.82 ± 0.02 | 5.85 ± 0.05 | 5.90 | 5.90 | 5.90 | ||
| Int | X | 3.82 ± 0.02 | 3.82 ± 0.03 | 3.80 ± 0.04 | 3.95 | 3.95 | 3.94 | |
| Y | 3.86 ± 0.01 | 3.87 ± 0.03 | 3.85 ± 0.03 | 3.93 | 3.93 | 3.93 | ||

