## Abstract

The fundamental operation in binder jet three-dimensional printing is the deposition of liquid binder into a powder layer to selectively bond particles together. Upon droplet impact, the binder spreads into the powder bed forming a bound network of wetted particles called a primitive. A computational fluid dynamics framework is proposed to directly simulate the capillary and hydrodynamic effects of the interfacial flow that is responsible for primitive formation. The computational model uses the volume-of-fluid method for capturing dynamic binder-air interfaces, and the immersed boundary method is adopted to include particle geometries on numerical Cartesian grids. Three-phase contact angles are prescribed through an interface extension algorithm. Binder droplet impact on powder beds of varying contact angle are simulated. Furthermore, the numerical model is used to simulate liquid bridges connecting binary and ternary particle systems, and the resulting capillary and hydrodynamic forces are validated by comparison with published experimental and theoretical model results.

## 1 Introduction

Additive manufacturing (AM) describes a broad class of technologies that are based on the successive addition or consolidation of thin material layers to generate three-dimensional objects in accordance to computer-aided-design (CAD) models. AM allows for the fabrication of highly complex geometries, which can be prohibitively expensive, inefficient, or even impossible to produce with traditional manufacturing. Furthermore, most AM machines are indifferent to the geometry of the part that it has been tasked to build; therefore, designs can be adapted or altered without major retooling or machine modifications. In contrast with conventional manufacturing methods based on milling, casting, or forging, the general “bottom-up” paradigm of AM relaxes many manufacturing imposed restrictions on topology and shape optimization. With this added flexibility, part designs can incorporate unique lattice structures and exotic geometries that permit significant mass reduction while maintaining strength and performance. In the aerospace industry, where weight is a primary design constraint, AM is being utilized for the creation of strong, lightweight, high-performance aircraft and spacecraft components. A variety of other engineering applications, such as those in the biomedical, defense, and energy sectors, are also leveraging the unique advantages offered by AM for current or future part production.

Binder jet 3D printing (BJ3DP) is a powder bed AM method that relies on the jetting of liquid binder into a powder layer to selectively bind particles together to form solid, three-dimensional objects. The method is carried out in a series of repeating steps where the desired shape is built up one cross-sectional layer at a time. The general procedure starts with the delivery of a powder feedstock on top of a build platform. Next, a cylindrical roller spreads the newly delivered powder heap into a smooth, thin layer. After spreading, a printhead scans over the powder bed and jets binder droplets into regions of desired solidification. Once the jetting stage is complete, the build platform lowers by a small, user-specified layer thickness, and the cycle repeats many times over until the complete geometry is formed. Typically, the parts must then be sintered at high temperature to achieve densification and strengthening.

Compared to other AM methods, BJ3DP is amenable to a wide range of materials, including metal superalloys and advanced ceramics [1]. Additionally, it boasts excellent build rates and cost-effectiveness. Despite these and further benefits of the process, which are described in the exhaustive BJ3DP review by Mostafaei et al. [2], the final end-use parts generally exhibit inferior mechanical properties and performance characteristics to those produced by other powder bed AM approaches, such as powder bed fusion (PBF) either by selective laser melting or electron beam melting [3]. As a result, process improvements are required before BJ3DP finds wide-scale adoption for manufacturing high-performance and mission-critical components.

The interfacial fluid-particle interaction that occurs between the powder, binder, and surrounding air is a predominant factor governing the final quality of BJ3DP parts. Upon impact of the jetted binder, momentum transferred from the droplet to the powder bed can result in ballistic particle ejection, which may give rise to increased porosity and reduced dimensional accuracy [4,5]. As binder is imbibed into the bed, the liquid forms an interconnected network of mutually wetted particles in the form of liquid bridges. Attractive forces by these liquid bridges can occur due to surface tension and capillary pressure induced by the curvature of the fluid interface. The resulting particle agglomeration bound together by these cohesive forces is referred to as a *primitive* and is the fundamental building element of the part [6]. In some cases, the capillary forces are large enough to overcome the gravitational and frictional forces holding the powder bed in place. This results in particle rearrangement, whereby the wetted particles are accelerated inward toward each other forming a roughly spherical shaped primitive resting on the powder bed surface. The primitive shape and binder distribution within it have a direct effect on part resolution, dimensional accuracy, and green part (pre-sintered) strength.

Liquid-mediated adhesion of capillary bridges between two surfaces is an important tribological problem in applications such as microelectromechanical (MEMS) systems, atomic force microscopy, and magnetic disk drives [7]. The same fundamental problem applies to BJ3DP primitive formation where liquid bridges connecting the powder particles produce attractive capillary forces. The binder saturation level, which is the relative volume of binder to the pore space in a predefined envelope within the powder bed, affects the configuration of the liquid network inside the primitive. The pendular, funicular, capillary, and droplet primitive regimes shown in Fig. 1 are dependent on the saturation level prescribed as an input process parameter to the machine. A crucial factor governing the final quality of BJ3DP parts is the selection of an optimal combination of process parameters in each stage of the method for a specific powder material, geometry, and function [8]. Saturation is an important parameter in the jetting stage because a sufficient amount of binder must be deposited to develop adequate inter-particle binding; however, excessive saturation will cause the liquid to flow out of the desired part boundary resulting in an increase in minimum feature size (i.e., decreased resolution) and loss of dimensional accuracy [9]. Assuming the capillary regime in Fig. 1, Miyanaji et al. [10] devised a combined empirical and theoretical model to determine the optimal saturation input for two powder feedstocks. The model succeeded in predicting the optimal saturation level for single-layer prints with Ti-6Al-4V powder, but it failed for 420 stainless steel (SS) powder and in multi-layer prints with either Ti-6Al-4V or 420 SS. The authors ascribe the model’s failure to its inability to completely incorporate solid particle geometry; however, we also theorize that significant error could be attributed to neglecting particle rearrangement caused by capillary forces. Several direct experimental approaches were taken to determine the effect of saturation, as well as other BJ3DP process parameters, on the resulting surface roughness, dimensional accuracy, and strength of the final parts [11–16].

To date, most of the advancements in BJ3DP process parameter optimization have been made in the aforementioned empirical studies or through basic trial and error experimentation. These experiments are often arduous and can be expensive if a significant amount of consumable material is required. Moreover, experiments can be difficult to setup and various physical aspects limit the amount of information that can be gathered from them. Computational simulation based on the numerical solution of the governing differential equations may provide new insight into the complicated physics of BJ3DP, which could help guide process modifications and improvements. This work is primarily focused on developing a computational model for simulating interfacial fluid flow in BJ3DP to ultimately be used in numerical experiments for determining optimal combinations of printing parameters.

The majority of computational modeling as it relates to BJ3DP has been deployed for investigating particle dynamics in the powder spreading stage with the discrete element method (DEM) [17–21]. Direct numerical simulation of fluid-particle interaction in BJ3DP using a coupled computational fluid dynamics (CFD) and DEM approach has currently not been published. Tan [22] performed simulations of binder droplet impact on fixed powder beds but did not consider the fluid-induced forces acting on the particles and therefore did not incorporate particle motion. This work seeks to develop a CFD algorithm sufficient to simulate the hydrodynamic and capillary flows occurring during primitive formation in any liquid bridge regime, i.e., we are not restricted to pendular liquid bridges of well-defined shape. The CFD algorithm numerically solves the Navier–Stokes equations of incompressible flow and employs a volume-of-fluid (VOF) strategy for capturing the fluid interface dynamics [23]. The immersed boundary (IB) method is included to model the geometry of the solid powder particles within the flow field. The numerical framework is applied to several test cases involving multi-particle liquid bridge systems and compared with published experimental results to gauge its ability to model capillary and hydrodynamic forces in BJ3DP.

## 2 Governing Equations

**u**and pressure

*p*fields of a fluid in a given spatial domain are governed by the time-dependent Navier–Stokes (N–S) equations:

*a*)

*b*)

*ρ*is density, $\tau $ is the viscous stress tensor, and

**I**is the identity matrix. The terms

**f**

_{g}, $f\sigma $, and

**f**

_{ib}are body forces induced by gravity, surface tension, and the presence of an immersed boundary, respectively. The fluid is considered to be Newtonian with the viscous stress tensor defined as

*μ*is the dynamic viscosity. The N–S equations are a set of coupled, nonlinear, partial differential equations (PDEs) that enforce mass and momentum conservation, Eqs. (1a) and (1b), respectively. For incompressible flows, Eq. (1a) yields the following condition:

*Dρ*/

*Dt*= 0); however, the fluid properties, namely

*ρ*and

*μ*, are not constant throughout the domain as they would be for a single-phase. To incorporate these quantities into Eq. (1), the VOF method developed at the Los Alamos National Laboratory by Hirt and Nichols [23] is used. In the VOF approach, a scalar field

*c*is introduced that indicates which fluid (fluid 1 or fluid 2) is present at any point in the domain, such that

**x**is the position vector of the given point. Numerical schemes for solving partial differential equations typically require that the domain be discretized into many small volumetric cells (or elements). The volume fraction

*C*, which we often refer to as the VOF function, of fluid 1 in a cell Ω is then

*C*= 1 for a cell completely full of fluid 1,

*C*= 0 for a cell completely absent of fluid 1 (i.e., full of fluid 2) and 0 <

*C*< 1 for a cell containing a mixture of the two phases. For immiscible fluids, a mixed cell implicitly indicates the presence of an interface. The density and viscosity fields can then simply be found by

*a*)

*b*)

*σ*is the surface tension coefficient,

*κ*is the interfacial curvature, $n^\Gamma $ is the interface unit normal vector pointing outward from fluid 1, and $\delta \Gamma $ is a Dirac delta function that is zero everywhere except for at the interface. The subscript Γ indicates the fluid–fluid interface.

We use the IB method to model the geometry of solid powder particles embedded within the flow domain [31]. Generally, an IB method is any approach that incorporates fictitious body forces in Eq. (1b) to implicitly apply no-slip velocity boundary conditions (BCs) on solid surfaces; thus, forcing the flow to behave as if there was a solid object present. Today, there exists many variations of Peskin’s [31] original method (see Ref. [17] for a review of popular IB methodologies). The IB strategy we adopt here is discussed further in Sec. 3.4.

*cl*) must also be included. Therefore, the total fluid force on the particle is

*s*and

*A*

_{s}is the surface area. The first term accounts for the force caused by the fluid pressure and viscous stresses, and the second term is the singular contact line force. The fluid-induced torque can be determined in a similar manner; however, we do not consider it in this work (see Ref. [32] for more details on the torque computation).

## 3 Numerical Methodology

Analytical solutions to the N–S equations have only been found for a small number of idealized cases that permit extensive simplifications of the PDEs; therefore, numerical solutions are usually required for practical engineering and scientific applications. We use the finite volume method (FVM) to discretize the governing equations on three-dimensional Cartesian grids. Additionally, we incorporate adaptive mesh refinement (AMR) with block-structured composite grids in order to localize grid refinement only in regions where it is needed, which leads to reduced computational expense. For simplicity, the following presentation of our numerical algorithm is in the context of a single Cartesian grid. For more information on the particular AMR strategy used, see Ref. [33].

The variables are arranged in accordance to the staggered marker-and-cell (MAC) layout [34], which is shown in Fig. 2. The scalar values, such as *p* and *C*, are stored at the center of each computational cell, and the *x*, *y*, and *z* velocity vector components, *u*, *v*, and *w*, respectively, are stored at the centroids of the corresponding cell faces. On structured Cartesian grids, all cell-centered nodes can be identified by a coordinate triplet (*i*, *j*, *k*), and face-centered nodes are offset by half-indices in each direction. Staggering the nodes in this way prevents velocity-pressure decoupling that results in the classic checkerboard oscillations [35]. Furthermore, FVM is based on tracking fluxed quantities through cell boundaries; therefore, positioning the fluxing velocities at the cell faces simplifies the FVM discretization.

The unsteady governing equations are integrated in time such that the velocity, pressure, and volume fraction fields are solved at discrete time-steps. We denote the current time-step as *n*; thus, *n* + 1 is the next time-step where the solution is updated. In the case that *n* = 0, appropriate initial conditions must be supplied to begin time marching.

### 3.1 Volume-of-Fluid Interface Advection.

*C*. As described by Rider and Kothe [36], the modified equation includes a divergence correction term resulting in

_{i,j,k}, which we define as the computational cell shown in Fig. 2. Applying Gauss’s divergence theorem yields:

*f*is a cell face, $n^f$ is the outward pointing unit normal vector to its respective face, $\u2212V\Omega i,j,k$ is the cell volume, and

*A*

_{f}is the area of the face, e.g.,

*A*

_{f}=

*h*

_{y}

*h*

_{z}for cell faces whose normal vector is parallel to the

*x*-axis. The first term on the right-hand side is the net volume fraction flux out of the cell and the second term is the discrete divergence correction.

*a*)

*a*)

*b*)

*t*is the discrete time-step and Φ

^{x}, Φ

^{y}, and Φ

^{z}are the one-dimensional volume fraction fluxes at each face. Note that for brevity, we drop the indicial subscripts of cell face nodes that coincide with cell centers. Instead of using the actual volume fraction in the divergence correction term, Weymouth and Yue [37] found that mass conservation could be improved by replacing it with the following:

A well-known difficulty of the VOF method is maintaining a narrow transition region as the volume fraction varies from zero to unity along its gradient, i.e., the interface must remain sharp as it is advected with the flow. Most established numerical schemes applied to Eq. (10) are apt to either excessively diffuse the interface (lower-order methods) or introduce nonphysical oscillations (higher-order methods) [38]. Ultimately, the source of this problem is Eq. (5), where we are attempting to represent the interface step function *c* as its cell volume average, effectively neglecting the position of the interface at length scales smaller than the cell size. Since the original introduction of the VOF approach, there has been a multitude of strategies devised to mitigate this issue, some of the prevalent methods are reviewed by Tryggvason et al. [38].

To preserve the interface sharpness, we implement the tangent of hyperbola interface capturing (THINC) technique [39] for the computation of cell face fluxes in Eq. (12). The method is based on fitting one-dimensional hyperbolic tangent functions to the VOF field in each coordinate direction allowing for analytical flux equations by integration of the smooth sigmoid functions. The multidimensional interface orientation is accounted for using the slope-weighting extension of the original THINC method (THINC-SW) given by Xiao et al. [40]. The THINC-SW technique adequately maintains sharp interfaces and is a simple and robust approach for both two- and three-dimensional flow fields. For further details on the THINC algorithm, see Refs. [41,42].

### 3.2 Incompressible Flow Solver

#### 3.2.1 Time Integration.

*n*+ 1. The first step is to discretize Eq. (1b) in time while temporarily neglecting the pressure and the immersed boundary terms:

**u*** is a provisional velocity field that is not divergence-free. The nonlinear advection term, viscous term, and gravitational body force are evaluated explicitly, while the surface tension force is considered implicit since the updated VOF function is used for the curvature in Eq. (8)—surface tension is described in greater detail below. The second step of the algorithm is the projection step, where the provisional velocity field is corrected to be divergence-free by means of the pressure gradient assessed implicitly:

The solution to the PPE is the pressure that enforces the incompressibility constraint on the velocity in Eq. (15). In general, solving the elliptic PPE is computationally expensive and is usually the primary bottleneck of the entire algorithm. For practical grid resolutions, the linear system resulting from the discrete PPE is too large for direct matrix inversion to be feasible; therefore, iterative solvers are standard for projection methods and in CFD more broadly. Convergence is accelerated with the use of a geometric multigrid solver designed for Poisson equations having variable and discontinuous coefficients [44,45].

#### 3.2.2 Spatial Discretization.

**A**and the diffusion term

**D**are given as

*a*)

*b*)

In this case, Eq. (17) is a vector equation where the velocity components lie on the cell faces rather than at the cell center. Therefore, the control volume for each vector component does not coincide with the cell shown in Fig. 2, but rather it is shifted by half the cell size such that it becomes centered on the respective staggered node. For example, the control volume Ω_{i+1/2,j,k} is centered on *u*_{i+1/2,j,k}, so that the net flux out of this staggered control volume is used to calculate *u*_{i+1/2,j,k}*. The derivation of the discrete PPE is carried out in the same way as shown by Tryggvason et al. [38], along with more details on the discretizations presented here.

Because of the staggered control volumes, the fluxes in Eq. (18) require velocity values at the cell centers and edges, as well as density and viscosity values at cell faces. Linear interpolation is used for all values except for the advected momentum density (*ρ***u**)^{n} in Eq. (18a). Away from the interface, this quantity is interpolated at the needed location with the upwind van Leer scheme to maintain stability and second-order accuracy [46]. For any cell that contains an interface, simple upwind interpolation is used for the advected momentum density. Additionally, it is noted that the velocity gradients in the viscous stress tensor are replaced with centered finite differences.

Due to the fact that the VOF advection equation given by Eq. (7) is just a different form of the mass continuity equation (1a), care must be taken to ensure the momentum fluxes are consistent with the VOF (i.e., mass) fluxes in Eq. (12) [47,48]. Inconsistent mass and momentum coupling can lead to severe interface distortions, especially when the density ratio between the two phases is large [49], e.g., *ρ*_{1}/*ρ*_{2} = *O*(10^{3}) for water and air. In this case, we have developed our own approach for maintaining mass-momentum coupling similar to the method introduced by Fuster et al. [50], which will be presented in future work.

The velocity BCs at the domain walls are either no-slip, free-slip, outflow, or periodic. Homogeneous Neumann BCs for pressure are supplied for the solution of the PPE, except when the coinciding velocity BC is prescribed as an outflow condition. For which case, a Dirichlet BC for pressure that is set to an arbitrary reference pressure is used.

### 3.3 Surface Tension.

*h*subscript indicates that the gradient operator is approximated with finite differences. The unit interface normal is also estimated with the VOF gradient by

*κ*is by numerically calculating the divergence of the interface normals:

*parasitic*velocities. These spurious currents are nonphysical and can deteriorate the interface, especially at small length scales where much of the dominating physics in BJ3DP occur. Although it remains an active research topic, various alternatives to Eq. (22) have been shown to significantly reduce parasitic velocities such that surface tension driven flows can be simulated with the CSF method [51–55]. Presently, we follow the height function approach given by Lopez et al. [56] to find

*κ*, which helps minimize the deleterious effects of spurious velocities.

### 3.4 Immersed Boundary Method.

*ϕ*is a scalar field expressing the fractional volume of solid that is present in each cell, just as the VOF function does for the reference fluid phase, and

**U**

_{ib}is the velocity of the particular IB in question. We note that each component of

**f**

_{ib}is located on its respective face-centered node, but the scalar

*ϕ*value is cell-centered. Therefore, the actual value of

*ϕ*in Eq. (23) is found by averaging the volume fractions of the neighboring cells sharing the common face. The force is used to correct the velocity field to account for the embedded flow obstacles by

### 3.5 Three-Phase Contact Line Dynamics.

Although fluid–fluid interface dynamics are accounted for with the combination of the VOF and CSF methods, additional capability is required to handle the three-phase contact line such that capillary effects are adequately modeled. The method employed here is based on the interface extension technique first developed by Sussman [61] for the level set method, but has since been applied to several VOF algorithms with little additional effort [58,59,62–64]. The balance of interfacial tensions at the trijunction of the solid–gas (SG), solid–liquid (SL), and liquid–gas (LG) interfaces suggests a contact angle *θ* that the LG interface must make with the SL interface to maintain equilibrium, as described by Young’s equation: *σ*_{SL} − *σ*_{SG} + *σ*_{LG} cos *θ* = 0. Note that our previous use of the symbol *σ* implies *σ* = *σ*_{LG}. This means that *θ* is a property of the solid–liquid–gas system and therefore should be an input into the computational framework.

The numerical approach for prescribing *θ* at three-phase contact lines is explained by first considering the example VOF field given in Fig. 3(a), which represents a fluid–fluid interface intersecting an IB without a prescribed *θ*. The assumed *θ* that is shown is simply the result of truncating a circular VOF field with a circular solid volume fraction field. Two issues become immediately evident: (1) the apparent *θ* must somehow be adjusted to fulfill the actual *θ* of the material system being modeled and (2) the VOF function does not have physical meaning inside the solid because there is no available volume for either fluid to exist in. This results in the VOF function erroneously indicating that gas (or fluid 2) is present inside the solid since *C* = 0 here. Because of this, an artificial LG interface occurs that coincides with the IB surface. This would be representative of a completely hydrophobic system (*θ* = 180 deg), and if the simulation proceeded as is, surface tension forces attributed to the artificial LG interface would be generated normal to the SL interface causing the droplet to quickly retract and “jump” off of the solid.

One way of addressing both of these concerns is to extend the interface several cells into the solid such that the correct *θ* is satisfied at the contact line (or contact point in 2D), and the artificial LG interface is displaced inside of the solid where the velocities that it generates are corrected by the force in Eq. (24). Fujita et al. [65] proved that integrating the surface tension term over the extended *virtual* interface produces an identical total surface tension force as would result from integrating over the contact line provided that the virtual surface remains completely within the IB and crosses the IB interface at the required *θ*. This means that an interface extension method implicitly handles the capillary effect due to three-phase contact lines and does not require further capability.

*τ*= min(

*h*

_{x},

*h*

_{y},

*h*

_{z}). An expression for the artificial extension velocity

**u**

_{e}was given by Sussman [61], which reads

**u**

_{e}are cell-centered. In principal, the same VOF advection solver introduced in Sec. 3.1 can solve Eq. (30), only now the extension velocity is used in place of the actual flow velocity. In practice, however, we find that the nonzero divergence of the extension velocity field results in a crude extension of

*C*into the solid. A simple semi-Lagrangian approach, as described by Sun and Sakai [58], was found to be better suited for this task. The semi-Lagrangian scheme updates the volume fraction in each cell using:

**x**

_{cc}is the position vector of the cell center where the current

*C*value is being updated. The right-hand side of Eq. (29) is evaluated with trilinear interpolation from the surrounding nodes at the location of an upwind point cast backwards one pseudo-time-step along the material point’s trajectory. The extension is performed for a given number of iterations

*n*

_{e}, which will be equal to the number of cells inside the solid that the interface is extended by. Unless otherwise specified, we take

*n*

_{e}= 5.

The result of applying the extension technique to the VOF field in Fig. 3(a) with *θ* = 30 deg is shown in Fig. 3(b). At the contact point, the interface has been extended along a path that approximates the prescribed contact angle. Away from the contact point, the VOF function is extended inward along the solid surface normal. The extension procedure does result in a slight change of the contact line position in which a small amount of mass may be added or subtracted. To account for mass conservation errors resulting from the extension, the simple correction approach taken by Sun and Sakai [58] is used to correct volume excesses or deficits.

**u**

_{c}is an artificial velocity field that serves to compress the smeared interface. The quantity

*C*(1 −

*C*) localizes the compression in regions closest to the 0.5 contour of the

*C*field. The compression velocity is set to

*K*being a constant that controls the compression strength. In most cases, we find that

*K*= 2 provides the best results. The extension algorithm then becomes

*a*)

*b*)

*C** in Eq. (32a) and is then compressed in Eq. (32b). We note that the compression velocity found in Eq. (31) uses

*C** for the interface normal calculation. It was found that the compression may lead to small overshoots and undershoots of the VOF bounds (e.g.,

*C*> 1 or

*C*< 0). For these cases, the volume fractions are simply set to zero or unity. The result of the compression procedure applied to the example in Fig. 3(b) is shown in Fig. 3(c). It is clear that the extended interface is now sharper, and the contact angle is better defined.

The extension algorithm for applying three-phase contact angles requires less numerical overhead than other approaches since the contact line does not need to be explicitly located. Furthermore, the finite difference stencils used in the surface tension and curvature calculations do not need to be modified in the vicinity of the IB, which would be required if the interface were not extended.

We recognize that dynamic contact angles and contact angle hysteresis can have a significant impact on the overall flow dynamics for certain material combinations. Presently, we neglect velocity dependent contact angles; however, the VOF-IB approach with contact line extension is amenable to adding such capability [62,63].

### 3.6 Fluid-Induced Forces.

*a*)

*b*)

_{ib}is the control volume containing the IB solid. Substituting the alternative surface tension force given by Eq. (33b) into Eq. (1b), then integrating it over Ω

_{ib}and relating it to Eq. (35) allows for a numerically convenient equation for the total fluid force on the solid given as

## 4 Results

### 4.1 Binder Droplet Impact on a Fixed Powder Bed.

In this section, the simulation results for a binder droplet impacting a fixed polydispersed powder bed are shown. Two cases are presented with identical input parameters except for the value of *θ* applied at the particle-binder-air interface. Each test case is prescribed with three-phase contact angles of *θ* = 30 deg and *θ* = 150 deg, representing hydrophilic and (super)hydrophobic powder beds, respectively.

The problem setup is similar to what was done by Tan [22], except here we use a powder bed with a polydispersed particle size and a random packing arrangement of approximately 1000 spherical particles. The particle diameters were generated from a uniform distribution between 10 and 30 *μ*m. The packing density was calculated as 0.402. The input parameters that were used in both test cases are given in Table 1. In Fig. 4, the initial droplet and simulation domain are shown. The grid spacing at the binder-air interface is set so that *R*_{b}/*h* ≈ 38, where *R*_{b} is the initial binder droplet radius and *h* is the uniform cell size. It is once again noted that AMR is leveraged in our simulations; therefore, the grid refinement adapts to the interface evolution and is coarser in areas where it is not needed. The nondimensional Bond number for this problem is $Bo=\rho lgD2/\sigma \u22481\xd710\u22123$ indicating that gravity has very little influence compared to surface tension; therefore, we neglect it here. We do the same for all of the simulations hereafter as well.

The droplet evolution during impact is shown for both test cases in Fig. 5. In Figs. 5(a) and 5(b), we see the same initial droplets just before impacting the bed at *t* = 0 *μ*s, and then subsequent impact and penetration into the powder bed as time progresses. Top views for both impacts are given in Figs. 5(c) and 5(d), where binder spreading as a result of inertial and capillary forces may be observed. In the range of *t* = 0 to *t* = 5 *μ*s, both droplets penetrate into their respective powder beds at approximately the same depth and achieve a similar spreading diameter. After *t* = 5 *μ*s, the droplets begin to behave much differently as a result of the contrasting particle contact angles we prescribed for each test. In the hydrophilic case, the droplet continues to spread radially in the *x*–*y* plane and penetrate in the −*z* direction until around *t* = 10 *μ*s, when the migration slows considerably. This is indicative of the transition from the inertial to capillary regimes described by Tan [22] and Das et al. [68], where the initial high kinetic energy that is driving the flow is dissipated and the continued imbibition of binder into the bed is then governed solely by capillarity. In the case of the hydrophobic powder, the binder droplet begins to exhibit “finger” instabilities around *t* = 10 *μ*s, which are typically encountered for droplet impacts onto flat hydrophobic surfaces with significant roughness [69]. Because of the high contact angle, the liquid–air interfaces within the pores assume a convex shape resulting in surface tension forces being generated that oppose the flow into the pores. Between *t* = 5 and *t* = 10 *μ*s, the droplet begins to retract radially and liquid fingers in the powder bed pores are drawn back upward to the surface. At *t* = 20 *μ*s, the binder on the hydrophilic powder has formed an interconnected network of wetted particles on and beneath the powder bed surface; however, on the superhydrophobic powder, no liquid remains in the pores and the deformed droplet is concentrated completely on top of the powder bed surface. In fact, the droplet is in the process of rebounding off the hydrophobic powder at *t* = 20 *μ*s. In Fig. 6, we show a longer time range of the hydrophobic impact simulation up to *t* = 55 *μ*s, where a complete droplet rebound in the opposite direction of the initial velocity is observed.

### 4.2 Validation of Hydrostatic Liquid Bridge Forces.

As described in Sec. 1, the capillary forces acting on the wetted particles result in the formation of BJ3DP primitives. In order to choose optimal machine process parameters for a particular binder-powder system, it is important to understand the mechanisms by which the primitives are formed. In this section, we investigate the model’s ability to compute liquid-mediated cohesive forces on particles connected by liquid bridges. Simple pendular bridges form between two particles, but the more complicated liquid bridges that connect three or more particles are in the funicular regime (Fig. 1). Since both regimes are often encountered in BJ3DP, we simulate liquid bridge formation on binary and ternary particle systems. We readily admit that the use of CFD is seemingly overkill and not the most efficient tool for the study of hydrostatic liquid bridges. For such cases, the full system can be completely described by the minimization of surface tension and gravitational energies subject to a constant volume constraint and boundary conditions for the contact line angle. A popular numerical tool for energy minimization is SURFACE EVOLVER [70], which can be used to compute equilibrium interface configurations and resulting forces much faster than numerically solving the N–S equation. However, if the fluid or particle dynamics are of interest, then energy minimization is not sufficient in describing the system. Of course, a time-dependent CFD solution should still result in a minimal surface energy configuration if it is allowed to run long enough for the interface to reach a steady-state. With this in mind, we believe directly simulating a liquid bridge that attains its equilibrium shape is an important test for the robustness of the framework.

#### 4.2.1 Binary Particle System (pendular).

The capillary forces of liquid bridges are described as being made up of two components, the Laplace pressure force and the singular contact line force. The Laplace pressure force occurs due to the pressure difference across a curved interface, which, at equilibrium, is described by the famous Young-Laplace equation, *p*_{i} − *p*_{o} = *σκ*, where *p*_{i} and *p*_{o} are the pressures inside and outside of the liquid bridge, respectively. This force can be attractive or repulsive depending on the interface curvature. The other component, which is always attractive, is due to surface tension acting tangentially to the interface along the contact line. A schematic of a binary particle liquid bridge system is shown in Fig. 7(a), and a graphical explanation of the two components of the total capillary force is given in Fig. 7(b).

Pitois et al. [71] constructed an experimental apparatus to measure the static liquid bridge force between two equal-sized particles as a function of the spacing between the particle surfaces *S* and the liquid bridge volume $\u2212Vlb$. In this section, we directly simulate these experiments and compare our results. The apparatus was designed to basically reproduce the two-particle configuration shown in Fig. 7(a). In the experiments, both spheres are a polished ruby material and each has a radius of *R* = 0.004 mm. The bottom sphere is connected to a rod that can move vertically in the (+/− )*z*-direction by a motorized micrometer screw. After a precise volume of liquid is deposited into the gap between the spheres, the liquid is stretched as the bottom sphere moves at a constant velocity in the −*z* direction. During this process, a scale continuously measures the resulting force on the top sphere. To assume a quasi-static state, the bottom sphere is set to move at only −0.01 *μ*m/s.

We begin the simulation by initializing a “blob” of liquid between the particles at a prescribed spacing. The blob is a sphere of fluid that is truncated by the particles such that after the truncation it has the correct volume of liquid used in the experiments. When the simulation begins, the truncated spherical shape of the liquid is not consistent with the equilibrium shape for the specified contact angle; therefore, velocity will be generated that begins to move the liquid to its equilibrium configuration. The simulation is allowed to run until the fluid ceases to move and the force converges to a steady-state. Unlike the quasi-static experiment, the simulations are completely static. Eleven separate cases are performed at various spacing distances. In principle, we could have exactly replicated the quasi-static experiments by moving the bottom sphere at −0.01 *μ*m/s; however, it takes around 22 h to move a distance of 8.0 × 10^{−4} m (0.2 · *R*) at that velocity, so a simulation (with an execution time much longer than the actual simulated time) would not be feasible. All of the input parameters to the model are the same as those used in the experiments, which can be found in Table 2. We note that the density and viscosity of the fluids do not affect the equilibrium capillary force (although they do affect the dynamics before reaching steady-state); thus, the values for these variables are the same as those given in Table 1. Since this is an axisymmetric problem, we only simulate one-quarter of the liquid bridge with appropriate symmetry (free-slip) boundary conditions for computational efficiency. The grid spacing at the AMR level that tracks the fluid interface is such that *R*/*h* ≈ 64. The capillary force on the particles is computed with Eq. (36).

The simulation results for the steady-state liquid bridge shapes are shown in Fig. 8, and plots of the normalized capillary force with the normalized grid spacing are shown in Fig. 9. We note that the force was measured from the top sphere, but confirm that equal values of opposite sign are obtained from the bottom sphere. In Fig. 9, good agreement is found between the simulations and experimental results, especially for normalized surface spacing between 0.01 and 0.1. The largest deviation with the experiment exists at zero spacing, i.e., when the spheres are in contact. This is likely due to the fact that, in the model, there will be computational cells that contain solid volume fraction attributed to both spheres. In such cases, the surface normal calculation that is based on the gradients of solid volume fraction suffers as the gradient no longer represents the solid interface normal. Although rigid spheres will theoretically only contact each other at a singular point, the model represents the spheres by their volume fraction within finite cells; therefore, several cells near the theoretical contact point will contain volume fraction from both spheres. Beyond the normalized spacing of 0.1, the forces are slightly overestimated in the simulation. However, we were able to closely capture the distance at which the bridge becomes unstable and ruptures causing a sharp decrease in force. In the experiment, rupture occurs just under *S*/*R* = 0.25, and in the simulation, we find rupture occurring between the last two data points of *S*/*R* = 0.225 and *S*/*R* = 0.25.

#### 4.2.2 Ternary Particle System (funicular).

In this section, we deploy the modeling framework to investigate more complicated funicular liquid bridges between three equal-sized particles. Wang et al. [72] performed experiments to determine how the capillary force from a static liquid bridge connecting three particles varies with spacing distance, liquid volume, contact angle, and capillary pressure. These experimental results were also compared with corresponding SURFACE EVOLVER solutions. Here, we execute direct simulations of three-particle liquid bridges at various spacing distances and liquid volumes until equilibrium is reached. The resulting capillary forces are compared with the experimental and SURFACE EVOLVER results given by Wang et al. [72].

As shown in Fig. 10, the simulation is setup to exactly replicate the experimental apparatus described by Wang et al. [72] (see Table 3 for simulation parameters). As in the previous section, the experiment is considered quasi-static; however, the spheres in the simulations are completely static for the entirety of each test. After steady-state is reached for a particular spacing distance, the bottom particle is regenerated to adjust *S* for the next simulation. Liquid volumes of $\u2212Vlb=3.0\xd710\u221210m3$ (0.3 *μ*L) and $\u2212Vlb=6.0\xd710\u221210m3$ (0.6 *μ*L) are tested. Unlike binary liquid bridges, funicular bridges are not axisymmetric; therefore, the entire particle system is within the simulated domain for this test. The grid spacing at the AMR level that tracks the fluid interface is such that *R*/*h* ≈ 22.

In their experiments, Wang et al. [72] observed significant variations in contact angle on each sphere with varying spacing distance. The contact angles were optically measured and fit with third-order polynomials. The input contact angles for each simulation are obtained from these polynomials, which range from 10–80 deg. We point the reader to Wang et al. [72] for the specific polynomial coefficients.

The simulated steady-state funicular bridge shapes for the constant liquid bridge volume of $\u2212Vlb=6.0\mu $L are shown in Fig. 11. The corresponding normalized capillary forces measured on the bottom particle versus normalized spacing are given for both volumes in Fig. 12. Comparing Figs. 12(a) and 12(b), we observe that a larger liquid bridge volume results in reduced capillary force at small spacing distances and increased capillary force at larger spacing distances. The simulation results are in good agreement with the SURFACE EVOLVER solution and in reasonable agreement with the experimentally measured values. The termination point of each SURFACE EVOLVER plot is the distance for which the liquid bridge ruptures. It is clear that the sudden drop in force shown in both Figs. 12(a) and 12(b) corresponds closely with the rupture distance predicted by SURFACE EVOLVER; however, the experimental results indicate an appreciable increase in rupture distance compared to our simulations and SURFACE EVOLVER. The explanation given by Wang et al. [72] for the departure of the SURFACE EVOLVER forces from the experimental ones was that there were likely errors in the optical method used for measuring the contact angles at different spacing distances; therefore, the contact angles that were input into SURFACE EVOLVER were different from the actual physical values. For this reason, it makes sense that our simulation results are more closely aligned with the SURFACE EVOLVER solution than the experiments. Numerical experiments indicate that the slight deviation between the present simulation and the SURFACE EVOLVER solutions is reduced as the computational grid is further refined.

### 4.3 Validation of Hydrodynamic Liquid Bridge Forces.

The last simulation we present in this work is that of a dynamic liquid bridge between two spherical particles separating at a constant velocity. During separation, viscous hydrodynamic forces in the thin liquid film attempt to impede the motion of the particles; thus, augmenting the attractive capillary force and attributing to the strength of the particle agglomeration or primitive. For this section, we compare the forces obtained from the direct numerical simulation of a stretching liquid bridge with the theoretical model of Pitois et al. [71,73] for the combined capillary and viscous forces acting on the particles.

*H*(

*r*) is the film thickness given by

*r*indicates a point on the radial axis. With appropriate boundary conditions (see Ref. [74]), Eq. (38) can be integrated twice with respect to

*r*to yield an approximation for the viscous lubrication force as:

*S*/d

*t*. The above integration was taken from

*r*= 0 to

*r*= ∞, so it only considers an infinite film of a single fluid. A correction was derived by Matthewson [75] assuming a finite, cylindrical liquid bridge. The modified viscous force then becomes

*a*is the radius of the wetted contact area.

*H*and $\u2212Vlb$ by

To compare the increase in force due to the dynamic stretching with the static liquid bridge force, the same parameters for the binary system in Sec. 4.2.1 and Table 2 are applied here; however, in this simulation, each sphere is completely contained inside the domain, i.e., we are not applying an axisymmetric approximation for this case. The finest AMR grid size that evolves with the interface is such that *R*/*h* ≈ 20. In the first part of the simulation, the liquid is allowed to attain its equilibrium configuration at *S* = 0 before either sphere is moved. When equilibrium is achieved, the velocities are then quickly ramped up to their constant values, which are prescribed to the bottom and top spheres as $Vzbottom=\u22120.05$ m/s and $Vztop=0.05$ m/s.

Snapshots of the dynamic stretching simulation are shown in Fig. 13, where necking and subsequent rupture are observed at $S/R\u22480.4$. In Fig. 14, the dynamic normalized forces as a function of normalized spacing that are obtained from the simulation and the theoretical model described by Eq. (44) are given. Unsurprisingly, it is observed that the viscous force is most prevalent at small spacing distances, but then quickly dies out as the spacing increases, leaving only the capillary force to withstand further stretching. For more viscous liquids or faster separation velocities, the lubrication force will have a larger effect on the system for greater spacing distances. The theoretical model agrees closely with the simulation at small spacing, but then consistently underestimates the simulated force as time increases, suggesting that the main deviation lies in the capillary force model. Pitois et al. [73] also found Eq. (42) to underestimate the capillary force at larger separation distances.

## 5 Conclusions

A three-dimensional interfacial flow solver has been presented to simulate the capillary and hydrodynamic effects central to the BJ3DP process. The CFD framework uses the VOF method to capture the dynamics of the binder-air interface, and the IB method is applied for the inclusion of solid powder particles inside the fluid domain. Capillary effects due to three-phase contact lines are made possible by a VOF extension algorithm. Several numerical experiments are performed to demonstrate the model’s efficacy of simulating various fluid dynamical phenomena in BJ3DP. First, binder droplet impact on fixed hydrophilic and hydrophobic powder beds is simulated. Inertial and capillary imbibition of the binder into the hydrophilic bed is observed. In the hydrophobic case, the growth of finger instabilities and complete droplet rebound is captured. Next, the CFD framework is tasked with simulating static liquid bridges in two- and three-particle systems. The forces obtained in these simulations are validated by direct comparison with published experimental data and results from surface energy minimization software. Finally, a dynamic liquid bridge system is simulated and compared with a theoretical force model based on the addition of capillary and viscous lubrication forces.

With the CFD solver experimentally validated, future work will include the addition of particle motion with DEM so that full fluid-particle interactions in BJ3DP may be simulated. We point out that all of the simulations presented here are idealized cases of the fundamental physical phenomena in BJ3DP. However, they are necessary to ensure that the framework adequately captures the dominant physics of the process before it is deployed to make new discoveries and advancements in BJ3DP.

## Acknowledgment

This work was supported by a NASA Space Technology Research Fellowship. The authors would like to thank Craig Smith (NASA Glenn Research Center) for his support and valuable input on this project.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*c*=discontinuous fluid phase indicator function

**f**=body force vector in N–S

*g*=gravitational acceleration

*h*=computational cell size

*n*=time-step number

*p*=pressure

*r*=radius of curvature or coordinate of radial axis

*s*=solid surface

*t*=time

**u**=fluid velocity vector

**x**=position vector

**A**=advection vector term in N-S

*A*=surface area

*C*=fluid phase volume fraction

**D**=diffusion (viscous) vector term in N-S

**F**=total fluid-induced force vector on solid

*R*=particle radius

*S*=particle surface separation distance

**U**=solid body velocity vector

*i*,*j*,*k*=nodal indices

- $n^,t^$ =
unit normal and tangent vectors

*u*,*v*,*w*=fluid velocity vector components

*α*=half-filling angle

- Γ =
fluid-fluid interface

*δ*=dirac delta function

*θ*=equilibrium three-phase contact angle

*κ*=curvature

*μ*=dynamic viscosity

- Π =
capillary pressure tensor

*ρ*=density

*σ*=surface tension coefficient

- $\tau $ =
viscous stress tensor

*τ*=pseudo-time

- Φ =
VOF flux

*ϕ*=solid volume fraction

- Ω =
control volume

- $\u2207$ =
gradient operator

- $\u2212V$ =
volume

## References

*x*-ray Imaging

*arXiv preprint arXiv:1811.12327*.