Computation of Hydrodynamic and Capillary Phenomena in Binder Jet Three-Dimensional Printing

Wagner, Joshua J.; Fred Higgs, C.

doi:10.1115/1.4050942

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].

Fig. 1

View large Download slide

The (a) pendular, (b) funicular, (c) capillary, and (d) droplet regimes of the liquid distribution in BJ3DP primitives dependent on binder saturation level

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

The Eulerian velocity u and pressure p fields of a fluid in a given spatial domain are governed by the time-dependent Navier–Stokes (N–S) equations:

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = 0

(1a)

\frac{\partial ρ u}{\partial t} + \nabla \cdot (ρ u \otimes u) = \nabla \cdot (τ - p I) + f_{g} + f_{σ} + f_{i b}

(1b)

where ρ is density,

τ

is the viscous stress tensor, and I is the identity matrix. The terms f_g,

f_{σ}

⁠, 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

τ = μ (\nabla u + \nabla u^{T})

(2)

where μ 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. and , respectively. For incompressible flows, Eq. yields the following condition:

\nabla \cdot u = 0

(3)

requiring that the velocity field remain divergence-free throughout the domain and at all time instances. Another result of the flow being incompressible is that the energy equation becomes decoupled from the mass and momentum equations, which can be solved separately using the velocities obtained from the solution of Eqs. and . In this work, we consider only isothermal flow as we are primarily concerned with simulating hydrodynamic and capillary effects; however, we acknowledge that many AM problems necessitate the inclusion of the energy equation as temperature variations and heat transfer must be accounted for, e.g., melt pool dynamics in PBF [24–26], liquid-phase sintering [27], and phase-change occurring in colloidal suspensions [28,29].

While the N–S equations described above are formulated for a single fluid, they may also be used to analyze interfacial, multifluid flows, e.g., the simultaneous flow of binder and surrounding gas in BJ3DP. In such cases, the flow is still considered to be incompressible (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. , 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

c (x, t) = {\begin{matrix} 1 if x \in fluid 1 \\ 0 if x \in fluid 2 \end{matrix}

(4)

where 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 = \frac{1}{- V_{Ω}} \int_{Ω} c (x, t) d Ω

(5)

where

- V_{Ω}

is the volume of the cell. It follows that 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

ρ = ρ_{1} C + ρ_{2} (1 - C)

(6a)

μ = μ_{1} C + μ_{2} (1 - C)

(6b)

where the subscripts indicate the respective fluid phases. Obviously, the volume fraction configuration evolves with time as the interface moves with the flow. Substituting Eq. into Eq. leads to a hyperbolic advection equation that governs the interface kinematics:

\frac{\partial C}{\partial t} + \nabla \cdot (C u) = 0

(7)

For each time-step, the acquired velocity from the solution of the N-S equations is used in Eq. to update the volume fraction field as it is advected with the flow.

The surface tension term is included as a body force using the continuum surface force (CSF) approach of Brackbill et al. [30]:

f_{σ} = σ κ {\hat{n}}_{Γ} δ_{Γ}

(8)

where σ is the surface tension coefficient, κ is the interfacial curvature,

{\hat{n}}_{Γ}

is the interface unit normal vector pointing outward from fluid 1, and

δ_{Γ}

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.

With the solution of Eqs. and at incremental time-steps, the force and torque acting on each particle by the fluid can be computed through integration of the pressure and viscous shear stress over the particle surface. For interfacial flows, a force attributed to the surface tension acting along the three-phase contact line (cl) must also be included. Therefore, the total fluid force on the particle is

F = \int_{s} (τ - p I) \cdot {\hat{n}}_{s} d A_{s} + \oint_{c l} σ {\hat{t}}_{Γ} d l

(9)

where

{\hat{n}}_{s}

is the unit normal vector of the solid surface 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.

Fig. 2

View large Download slide

A computational grid cell with a staggered nodal arrangement. The location indices of each node are also given.

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.

The time integration procedure for each step begins by solving a modified form of Eq. to update C. As described by Rider and Kothe [36], the modified equation includes a divergence correction term resulting in

\frac{\partial C}{\partial t} + \nabla \cdot (C u) = C \nabla \cdot u

(10)

Mathematically, Eqs. and are identical for incompressible flows; however, the correction term on the right-hand side accounts for nonzero discrete divergence introduced by the numerical method. The modified VOF advection equation is discretized with FVM by integrating it over a control volume Ω_i,j,k, which we define as the computational cell shown in Fig. 2. Applying Gauss’s divergence theorem yields:

- V_{Ω_{i, j, k}} \frac{\partial C}{\partial t} = - \sum_{f = 0}^{5} C (u \cdot {\hat{n}}_{f}) A_{f} + C \sum_{f = 0}^{5} (u \cdot {\hat{n}}_{f}) A_{f}

(11)

where f is a cell face,

{\hat{n}}_{f}

is the outward pointing unit normal vector to its respective face,

- V_{Ω_{i, j, k}}

is the cell volume, and A_f is the area of the face, e.g., A_f = h_yh_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.

We use the operator-split algorithm of Weymouth and Yue [37] to carry out the time integration for Eq. . Operator-splitting allows for the multidimensional advection equation to be solved with independent, one-dimensional sweeps along each coordinate direction. The operator-split algorithm is as follows:

\begin{aligned} C_{i, j, k}^{*} & = C_{i, j, k}^{n} - \frac{Δ t}{- V_{i, j, k}} [(Φ_{i + 1 / 2}^{x} - Φ_{i - 1 / 2}^{x}) \\ - C_{i, j, k}^{W Y} h_{y} h_{z} (u_{i + 1 / 2} - u_{i - 1 / 2})] \end{aligned}

(12a)

\begin{aligned} C_{i, j, k}^{* *} & = C_{i, j, k}^{*} - \frac{Δ t}{- V_{i, j, k}} [(Φ_{j + 1 / 2}^{y} - Φ_{j - 1 / 2}^{y}) \\ - C_{i, j, k}^{W Y} h_{x} h_{z} (v_{j + 1 / 2} - v_{j - 1 / 2})] \end{aligned}

(12a)

\begin{aligned} C_{i, j, k}^{n + 1} & = C_{i, j, k}^{* *} - \frac{Δ t}{- V_{i, j, k}} [(Φ_{k + 1 / 2}^{z} - Φ_{k - 1 / 2}^{z}) \\ - C_{i, j, k}^{W Y} h_{x} h_{y} (w_{k + 1 / 2} - w_{k - 1 / 2})] \end{aligned}

(12b)

where the * and ** superscripts indicate the intermediate sweeps in the first two dimensions, Δ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:

C_{i, j, k}^{W Y} = {\begin{matrix} 1 & if C_{i, j, k}^{n} > 0.5 \\ 0 & else \end{matrix}

(13)

The importance of the divergence correction term is evident from Eqs. since the velocities used in the one-dimensional sweeps are clearly not divergence-free. Without the correction term to counteract the nonzero divergence, spurious volume generation or depletion occurs with each dimensional sweep.

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.

After the VOF advection equation is solved, the spatially varying density and viscosity fields are updated with Eqs. and . A two-step Chorin projection [43] is employed to solve Eq. and enforce the incompressibility constraint given by Eq. on the new velocity field at n + 1. The first step is to discretize Eq. in time while temporarily neglecting the pressure and the immersed boundary terms:

\frac{ρ^{n + 1} u^{*} - ρ^{n} u^{n}}{Δ t} = - \nabla \cdot (ρ u \otimes u)^{n} + \nabla \cdot τ^{n} + f_{g}^{n} + f_{σ}^{n + 1}

(14)

where 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. —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:

u_{f l}^{n + 1} = u^{*} - \frac{Δ t}{ρ^{n + 1}} \nabla p^{n + 1}

(15)

where

u_{f l}^{n + 1}

is the corrected fluid velocity irrespective of the immersed boundary. Due to the constraint that

\nabla \cdot u_{f l}^{n + 1} = 0

⁠, taking the divergence of Eq. results in the following pressure Poisson equation (PPE):

\nabla \cdot \frac{1}{ρ^{n + 1}} \nabla p^{n + 1} = \frac{\nabla \cdot u^{*}}{Δ t}

(16)

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.

In the same way that we derived Eq. , Eq. is discretized in space with FVM by integrating it over a control volume and applying the divergence theorem to obtain the following:

ρ^{n + 1} u^{*} = ρ^{n} u^{n} - \frac{Δ t}{- V_{c v}} (A^{n} - D^{n}) + Δ t (f_{g}^{n} + f_{σ}^{n + 1})

(17)

where the flux differencing form of the advective term A and the diffusion term D are given as

A^{n} = \sum_{f = 0}^{5} (ρ u)^{n} (u^{n} \cdot {\hat{n}}_{f}) A_{f}

(18a)

D^{n} = \sum_{f = 0}^{5} τ^{n} \cdot {\hat{n}}_{f} A_{f}

(18b)

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)ⁿ 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., ρ₁/ρ₂ = O(10³) 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.

The interface delta function in Eq. is approximated using the gradient of the volume fraction:

δ_{Γ} = ‖ \nabla_{h} C ‖

(19)

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

{\hat{n}}_{Γ} = - \frac{\nabla_{h} C}{‖ \nabla_{h} C ‖}

(20)

Substituting Eqs. and into Eq. gives the numerical approximation to the surface tension term:

f_{σ} = - σ κ \nabla_{h} C

(21)

The simplest way of obtaining the interfacial curvature κ is by numerically calculating the divergence of the interface normals:

κ = - \nabla_{h} \cdot {\hat{n}}_{\tilde{C}}

(22)

where

{\hat{n}}_{\tilde{C}}

is the unit normal obtained from a smoothed VOF function

\tilde{C}

⁠. Unfortunately, this approach is known to produce numerical artifacts in the form of vortices around the interface, commonly referred to as 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. 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.

Up to this point, the IB body force in the momentum equation has been ignored. The divergence-free velocity obtained from Eq. does not “feel” the presence of any solid geometries besides the domain wall boundaries resulting in the generation of artificial fluid velocity inside the solid. A correction is made by forcing the solid velocity onto the nodes that are contained within the solid or on the IB surface. Because we consider only regular Cartesian grids that do not conform to the IB surface geometry (unless it is a flat surface parallel to one of the coordinate axes), the IB forces must be interpolated from the surface at nearby nodes. Several interpolation strategies were investigated by Fadlun et al. [57], which included direct injection of the solid velocity at the node nearest to the IB (stair-step approximation), volume fraction weighting, and linear interpolation along a line extending normal to the IB surface. Due to its simplicity and adherence to the VOF philosophy, we use the solid volume fraction in each cell as a weight on the amount of solid velocity that contributes to that cell. In addition to Fadlun et al. [57], this technique was also successfully implemented by Sun and Sakai [58,59] and Washino et al. [60]. With this method, the IB body force is given by

f_{i b}^{n + 1} = ρ \frac{ϕ (U_{i b} - u_{f l}^{n + 1})}{Δ t}

(23)

where ϕ 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. 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

u^{n + 1} = u_{f l}^{n + 1} + \frac{Δ t}{ρ} f_{i b}^{n + 1}

(24)

In the case that the total fluid force and torque on the object is not of interest, the IB force term does not actually have to be computed outright. This is because substituting Eq. into Eq. yields

u^{n + 1} = U_{i b} ϕ + u_{f l}^{n + 1} (1 - ϕ)

(25)

which allows for the final velocity field to be updated directly. On the other hand, the IB force is particularly useful when integrating the fluid stresses in Eq. , as we will show in Sec. 3.6.

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.

Fig. 3

View large Download slide

Examples of the interface extension algorithm for applying a contact angle at an IB surface (heavy black line). The C field is shown (a) before the extension, (b) after the semi-Lagrangian extension, and (c) after artificial compression. The 0.5 contour of C is also pictured (light black line). The prescribed contact angle is 30 deg.

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.

The extension is carried out by solving the following special advection equation:

\frac{\partial C}{\partial τ} + \nabla \cdot (C u_{e}) = 0

(26)

where a pseudo-time-step Δτ = 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} = {\begin{matrix} - \frac{{\hat{n}}_{s} - \cot (π - θ) {\hat{n}}_{2}}{‖ {\hat{n}}_{s} - \cot (π - θ) {\hat{n}}_{2} ‖} & if b < 0 \\ - \frac{{\hat{n}}_{s} + \cot (π - θ) {\hat{n}}_{2}}{‖ {\hat{n}}_{s} + \cot (π - θ) {\hat{n}}_{2} ‖} & if b > 0 \\ - {\hat{n}}_{s} & else \end{matrix}

(27)

with:

n_{1} = \frac{{\hat{n}}_{Γ} \times {\hat{n}}_{s}}{‖ {\hat{n}}_{Γ} \times {\hat{n}}_{s} ‖}, n_{2} = \frac{{\hat{n}}_{1} \times {\hat{n}}_{s}}{‖ {\hat{n}}_{1} \times {\hat{n}}_{s} ‖}, and b = {\hat{n}}_{Γ} \cdot {\hat{n}}_{2}

(28)

We note that the values u_e are cell-centered. In principal, the same VOF advection solver introduced in Sec. 3.1 can solve Eq. , 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:

C^{τ + Δ τ} = C^{τ} (x_{c c} - u_{e} Δ τ)

(29)

where x_cc is the position vector of the cell center where the current C value is being updated. The right-hand side of Eq. 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.

The volume fraction interpolation in the semi-Lagrangian method results in significant smearing of the extrapolated interface, which reduces the accuracy of the applied contact angle. As a result, we have devised a simple scheme based on the VOF advection algorithm of Weller [66] to artificially compress the interface so that a sharper contact line is produced. This is done by adding a compression term to the right-hand side of Eq. such that

\frac{\partial C}{\partial τ} + \nabla \cdot (C u_{e}) = - \nabla \cdot (u_{c} C (1 - C))

(30)

where 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

u_{c} = - K {\hat{n}}_{Γ}

(31)

with 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

C^{*} = C^{τ} (x_{c c} - u_{e} Δ τ)

(32a)

C^{τ + Δ τ} = C^{*} - Δ τ \nabla_{h} \cdot (u_{c} C^{*} (1 - C^{*})) for ϕ > 0

(32b)

Here, the smeared solution obtained from the semi-lagrangian method is set to be a provisional volume fraction field C* in Eq. and is then compressed in Eq. . We note that the compression velocity found in Eq. 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.

The total fluid force acting on a solid object in the flow field was given in Eq. . Following Sun and Sakai [59], the contact line term is replaced with an equivalent surface integral containing the capillary pressure tensor

Π

derived by Lafaurie et al. [67]:

Π = - σ (I - {\hat{n}}_{Γ} \otimes {\hat{n}}_{Γ}) δ_{Γ}

(33a)

f_{σ} = - \nabla \cdot Π

(33b)

which results in

F = \int_{s} (τ - p I) \cdot {\hat{n}}_{s} d A_{s} - \int_{s} Π \cdot {\hat{n}}_{s} d A_{s}

(34)

In the IB method, the solid object is represented on a Cartesian grid where the nodes do not necessarily align with the solid geometry. This makes numerical integration over surfaces difficult, and it is therefore desirable to transform the surface integrals to volume integrals with Gauss’s divergence theorem. Following the transformation, a summation over a control volume containing the solid can be performed to find the total force [64]. Applying the divergence theorem to Eq. gives

F = \int_{Ω_{i b}} \nabla \cdot (τ - p I) d Ω - \int_{Ω_{i b}} \nabla \cdot Π d Ω

(35)

where Ω_ib is the control volume containing the IB solid. Substituting the alternative surface tension force given by Eq. into Eq. , then integrating it over Ω_ib and relating it to Eq. allows for a numerically convenient equation for the total fluid force on the solid given as

F = - \int_{Ω_{i b}} f_{i b} d Ω - \int_{Ω_{i b}} f_{g} d Ω + \int_{Ω_{i b}} \frac{D ρ u}{D t} d Ω

(36)

The first term is the integration of the IB forcing which contains the hydrodynamic and capillary contributions, the second term is the buoyancy force, and the third term is the so-called added mass term that is attributed to the force required to move the surrounding fluid out of the way of the solid as it is accelerating. The benefit of representing the contact line force with the divergence of the capillary tensor is evidenced by the fact that it is contained in the IB force and requires no additional work to compute. Finally, the integrals in Eq. are approximated using

F = - \sum_{Ω_{i b}} (f_{i b}^{n + 1} + f_{g}^{n} + \frac{(ρ u)^{n + 1} - (ρ u)^{n}}{Δ t} + A^{n})

(37)

where a discrete summation is taken over all the cells containing nonzero solid volume fraction associated with the particular IB in consideration. The torque calculation is a straightforward extension of the same general methodology.

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 $B o = ρ_{l} g D^{2} / σ \approx 1 \times 10^{- 3}$ 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.

Fig. 4

View large Download slide

Domain setup for the simulation of a binder droplet impacting a powder bed

Table 1

Simulation parameters for the binder droplet impact problem

Parameter	Value	Unit
Droplet diameter, D	60	μm
Impact velocity, w_o	7	m/s
Binder density, ρ_b	1000	kg/m³
Air density, ρ_a	1.25	kg/m³
Binder viscosity, μ_b	4.31 × 10⁻³	Pa · s
Air viscosity, μ_a	1.0 × 10⁻⁵	Pa · s
Binder-air surface tension, σ	3.0 × 10⁻²	N/m

Parameter	Value	Unit
Droplet diameter, D	60	μm
Impact velocity, w_o	7	m/s
Binder density, ρ_b	1000	kg/m³
Air density, ρ_a	1.25	kg/m³
Binder viscosity, μ_b	4.31 × 10⁻³	Pa · s
Air viscosity, μ_a	1.0 × 10⁻⁵	Pa · s
Binder-air surface tension, σ	3.0 × 10⁻²	N/m

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.

Fig. 5

View large Download slide

Binder droplet impact on a fixed powder bed at different time instances. The left two columns are the front views of each simulation with translucent particles so that the binder penetration depth can be seen: (a) θ = 30 deg and (b) θ = 150 deg. The right two columns show the top view of the same simulations: (c) θ = 30 deg and (d) θ = 150 deg. The binder-air interface is colored by the velocity magnitude.

Fig. 6

View large Download slide

A binder droplet rebounding off a powder bed where the contact angle of each particle is prescribed to θ = 150 deg. Note: the first five snapshots are the same results shown in Fig. 5(b).

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).

Fig. 7

View large Download slide

A diagram of (a) the pendular liquid bridge between two particles showing relevant variables, and (b) the Laplace pressure (lp) and contact line (cl) components of the total capillary (cap) force

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 $- V_{l b}$ ⁠. 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⁻⁴ 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).

Table 2

Simulation parameters for the static binary liquid bridge problem

Parameter	Value	Unit
Particle radius, R	4 × 10⁻³	m
Surface tension, σ	0.021	N/m
Liquid bridge volume, $- V_{l b}$	1.5	μL
Contact angle, θ	10	°

Table 3

Simulation parameters for the static ternary liquid bridge problem

Parameter	Value	Unit
Particle radius, R	9.3 × 10⁻⁴	m
Surface tension, σ	0.073	N/m
Liquid bridge volume (a), $- V_{l b}^{a}$	3.0	μL
Liquid bridge volume (b), $- V_{l b}^{b}$	6.0	μL

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.

Fig. 8

View large Download slide

Simulated equilibrium liquid bridge shapes for a binary particle system at various normalized spacing distances. Note: hatching indicates the solid surface.

Fig. 9

View large Download slide

Normalized capillary force with normalized spacing distance as measured from the top sphere of the binary system. The experimental data were obtained from Pitois et al. [71].

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 $- V_{l b} = 3.0 \times 10^{- 10} m^{3}$ (0.3 μL) and $- V_{l b} = 6.0 \times 10^{- 10} m^{3}$ (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.

Fig. 10

View large Download slide

A diagram of a ternary liquid bridge. The dashed circle indicates the position of the bottom particle when S = 0.

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 $- V_{l b} = 6.0 μ$ 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.

Fig. 11

View large Download slide

Simulated equilibrium liquid bridge shapes for a ternary particle system at various normalized spacing distances. Only the case of $- V_{l b} = 6.0 μ L$ is shown here.

Fig. 12

View large Download slide

Normalized capillary force with normalized spacing distance as measured from the bottom sphere of the ternary particle system with liquid volumes of (a) 3.0 μL and (b) 6.0 μL. The experimental and SURFACE EVOLVER data were obtained from Wang et al. [72].

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.

In lubrication theory, the hydrodynamics of fluid films are often described by the Reynolds equation:

\frac{d}{d r} (r H^{3} (r) \frac{d p (r)}{d r}) = 12 μ r \frac{d S}{d t}

(38)

where H(r) is the film thickness given by

H (r) = S + r^{2} / R

(39)

Equations and are for an axisymmetric fluid film in cylindrical coordinates, where r indicates a point on the radial axis. With appropriate boundary conditions (see Ref. [74]), Eq. can be integrated twice with respect to r to yield an approximation for the viscous lubrication force as:

F_{visc} = - \frac{3}{2} π μ R^{2} \frac{1}{S} \frac{d S}{d t}

(40)

where the negative sign indicates that the force always acts opposite to the separation or approach velocity, which is simply defined as dS/dt. 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

F_{visc} = - \frac{3}{2} π μ R^{2} {(1 - \frac{S}{H (a)})}^{2} \frac{1}{S} \frac{d S}{d t}

(41)

where a is the radius of the wetted contact area.

The theoretical model for the capillary force described by Pitois et al. [71] is

F_{cap} = 2 π R σ \cos θ \cdot (1 - \frac{1}{\sqrt{1 + \frac{2 - V_{l b}}{π R S^{2}}}})

(42)

A key part of this derivation assumed the liquid bridge was cylindrical in order to relate H and

- V_{l b}

by

- V_{l b} = \frac{π R}{2} (H^{2} (a) - S^{2})

(43)

Adding Eqs. and and incorporating Eq. , the total force on a particle due to the combined capillary and viscous forces can therefore be approximated by

\begin{aligned} F_{tot} & = F_{cap} + F_{visc} \\ = 2 π R σ \cos θ \cdot X_{- V} + \frac{3}{2} π μ R^{2} \frac{1}{S} \frac{d S}{d t} \cdot X_{- V}^{2} \end{aligned}

(44)

with

X_{- V} = 1 - {(1 + \frac{2 - V_{l b}}{π R S^{2}})}^{- 1 / 2}

(45)

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 $V_{z}^{b o t t o m} = - 0.05$ m/s and $V_{z}^{t o p} = 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 \approx 0.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.

Fig. 13

View large Download slide

Simulation snapshots of a liquid bridge stretching as two particles separate at constant velocity

Fig. 14

View large Download slide

Normalized hydrodynamic force with normalized spacing distances as measured from the top sphere of a stretching two-particle system. The total hydrodynamic force from Eq. (44) is provided along with its individual viscous and capillary components given in Eqs. (41) and (42), respectively.

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
$\hat{n}, \hat{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
$τ$ =: viscous stress tensor
τ =: pseudo-time
Φ =: VOF flux
ϕ =: solid volume fraction
Ω =: control volume
$\nabla$ =: gradient operator
$- V$ =: volume

References

1.

Ziaee

,

M.

, and

Crane

,

N. B.

,

2019

, “

Binder Jetting: A Review of Process, Materials, and Methods

,”

Addit. Manuf.

,

28

, pp.

781

–

801

.

2.

Mostafaei

,

A.

,

Elliott

,

A. M.

,

Barnes

,

J. E.

,

Li

,

F.

,

Tan

,

W.

,

Cramer

,

C. L.

,

Nandwana

,

P.

, and

Chmielus

,

M.

,

2020

, “

Binder Jet 3D Printing–Process Parameters, Materials, Properties, and Challenges

,”

Prog. Mater. Sci.

, p.

100707

.

3.

Gokuldoss

,

P. K.

,

Kolla

,

S.

, and

Eckert

,

J.

,

2017

, “

Additive Manufacturing Processes: Selective Laser Melting, Electron Beam Melting and Binder Jetting Election Guidelines

,”

Materials

,

10

(

6

), p.

672

.

4.

Parab

,

N. D.

,

Barnes

,

J. E.

,

Zhao

,

C.

,

Cunningham

,

R. W.

,

Fezzaa

,

K.

,

Rollett

,

A. D.

, and

Sun

,

T.

,

2019

, “

Real Time Observation of Binder Jetting Printing Process Using High-speed x-ray Imaging

,”

Sci. Rep.

,

9

(

1

), pp.

1

–

10

.

5.

Wagner

,

J.

,

Shu

,

H.

,

Kilambi

,

R.

, and

Higgs III

,

C. F.

,

2019

, “

Experimental Investigation of Fluid-particle Interaction in Binder Jet 3D Printing

,”

Solid Freeform Fabrication Symposium–Additive Manufacturing

,

Austin, TX

, pp.

134

–

147

.

6.

Sachs

,

E.

,

Cima

,

M.

,

Cornie

,

J.

,

Brancazio

,

D.

,

Bredt

,

J.

,

Curodeau

,

A.

,

Fan

,

T.

,

Khanuja

,

S.

,

Lauder

,

A.

,

Lee

,

J.

, and

Michaels

,

S.

,

1993

, “

Three-Dimensional Printing: The Physics and Implications of Additive Manufacturing

,”

CIRP. Ann.

,

42

(

1

), pp.

257

–

260

.

7.

Streator

,

J. L.

,

2009

, “

A Model of Liquid-Mediated Adhesion With a 2D Rough Surface

,”

Tribol. Int.

,

42

(

10

), pp.

1439

–

1447

.

8.

Lores

,

A.

,

Azurmendi

,

N.

,

Agote

,

I.

, and

Zuza

,

E.

,

2019

, “

A Review on Recent Developments in Binder Jetting Metal Additive Manufacturing: Materials and Process Characteristics

,”

Powder. Metall.

,

62

(

5

), pp.

267

–

296

.

9.

Bredt

,

J. F.

,

1997

, “

Binder Stability and Powder/Binder Interaction in Three-Dimensional Printing

,”

Ph.D. thesis

,

Cambridge MA

.

10.

Miyanaji

,

H.

,

Zhang

,

S.

, and

Yang

,

L.

,

2018

, “

A New Physics-Based Model for Equilibrium Saturation Determination in Binder Jetting Additive Manufacturing Process

,”

Int. J. Mach. Tools. Manuf.

,

124

, pp.

1

–

11

.

11.

Chen

,

H.

, and

Zhao

,

Y. F.

,

2016

, “

Process Parameters Optimization for Improving Surface Quality and Manufacturing Accuracy of Binder Jetting Additive Manufacturing Process

,”

Rapid Prototy. J.

,

22

(

3

), pp.

527

–

538

.

12.

Shrestha

,

S.

, and

Manogharan

,

G.

,

2017

, “

Optimization of Binder Jetting Using Taguchi Method

,”

Jom

,

69

(

3

), pp.

491

–

497

.

13.

Lu

,

K.

, and

Reynolds

,

W. T.

,

2008

, “

3dp Process for Fine Mesh Structure Printing

,”

Powder. Technol.

,

187

(

1

), pp.

11

–

18

.

14.

Vaezi

,

M.

, and

Chua

,

C. K.

,

2011

, “

Effects of Layer Thickness and Binder Saturation Level Parameters on 3d Printing Process

,”

Int. J. Adv. Manuf. Technol.

,

53

(

1

), pp.

275

–

284

.

15.

Gaytan

,

S.

,

Cadena

,

M. A.

,

Karim

,

H.

,

Delfin

,

D.

,

Lin

,

Y.

,

Espalin

,

D.

,

Macdonald

,

E.

, and

Wicker

,

R.

,

2015

, “

Fabrication of Barium Titanate by Binder Jetting Additive Manufacturing Technology

,”

Ceram. Int.

,

41

(

5

), pp.

6610

–

6619

.

16.

Kafara

,

M.

,

Kemnitzer

,

J.

,

Westermann

,

H.-H.

, and

Steinhilper

,

R.

,

2018

, “

Influence of Binder Quantity on Dimensional Accuracy and Resilience in 3D-Printing

,”

Procedia Manufact.

,

21

, pp.

638

–

646

.

17.

Haeri

,

S.

, and

Shrimpton

,

J.

,

2012

, “

On the Application of Immersed Boundary, Fictitious Domain and Body-Conformal Mesh Methods to Many Particle Multiphase Flows

,”

Int. J. Multiphase. Flow.

,

40

, pp.

38

–

55

.

18.

Nan

,

W.

,

Pasha

,

M.

,

Bonakdar

,

T.

,

Lopez

,

A.

,

Zafar

,

U.

,

Nadimi

,

S.

, and

Ghadiri

,

M.

,

2018

, “

Jamming During Particle Spreading in Additive Manufacturing

,”

Powder. Technol.

,

338

, pp.

253

–

262

.

19.

Desai

,

P. S.

, and

Higgs

,

C. F.

,

2019

, “

Spreading Process Maps for Powder-Bed Additive Manufacturing Derived From Physics Model-based Machine Learning

,”

Metals

,

9

(

11

), p.

1176

.

20.

Zhang

,

W.

,

Mehta

,

A.

,

Desai

,

P. S.

, and

Higgs

,

C.

,

2017

, “

Machine Learning Enabled Powder Spreading Process Map for Metal Additive Manufacturing (AM)

,”

Solid Freeform Fabrication Symposium–Additive Manufacturing

,

Austin, TX

, pp.

1235

–

1249

.

21.

Desai

,

P. S.

,

Mehta

,

A.

,

Dougherty

,

P. S.

, and

Higgs III

,

C. F.

,

2019

, “

A Rheometry Based Calibration of a First-order Dem Model to Generate Virtual Avatars of Metal Additive Manufacturing (am) Powders

,”

Powder. Technol.

,

342

, pp.

441

–

456

.

22.

Tan

,

H.

,

2016

, “

Three-dimensional Simulation of Micrometer-Sized Droplet Impact and Penetration Into the Powder Bed

,”

Chem. Eng. Sci.

,

153

, pp.

93

–

107

.

23.

Hirt

,

C. W.

, and

Nichols

,

B. D.

,

1981

, “

Volume of Fluid (vof) Method for the Dynamics of Free Boundaries

,”

J. Comput. Phys.

,

39

(

1

), pp.

201

–

225

.

24.

Khairallah

,

S. A.

,

Anderson

,

A. T.

,

Rubenchik

,

A.

, and

King

,

W. E.

,

2016

, “

Laser Powder-Bed Fusion Additive Manufacturing: Physics of Complex Melt Flow and Formation Mechanisms of Pores, Spatter, and Denudation Zones

,”

Acta. Mater.

,

108

, pp.

36

–

45

.

25.

Bidare

,

P.

,

Bitharas

,

I.

,

Ward

,

R.

,

Attallah

,

M.

, and

Moore

,

A. J.

,

2018

, “

Fluid and Particle Dynamics in Laser Powder Bed Fusion

,”

Acta. Mater.

,

142

, pp.

107

–

120

.

26.

Tan

,

J.

,

Tang

,

C.

, and

Wong

,

C.

,

2018

, “

Study and Modeling of Melt Pool Evolution in Selective Laser Melting Process of Ss316l

,”

MRS Commun.

,

8

(

3

), pp.

1178

–

1183

.

27.

Villanueva

,

W.

,

Grönhagen

,

K.

,

Amberg

,

G.

, and

Ågren

,

J.

,

2009

, “

Multicomponent and Multiphase Simulation of Liquid-Phase Sintering

,”

Comput. Mater. Sci.

,

47

(

2

), pp.

512

–

520

.

28.

Fujita

,

M.

,

Koike

,

O.

, and

Yamaguchi

,

Y.

,

2015

, “

Direct Simulation of Drying Colloidal Suspension on Substrate Using Immersed Free Surface Model

,”

J. Comput. Phys.

,

281

, pp.

421

–

448

.

29.

Schlottke

,

J.

, and

Weigand

,

B.

,

2008

, “

Direct Numerical Simulation of Evaporating Droplets

,”

J. Comput. Phys.

,

227

(

10

), pp.

5215

–

5237

.

30.

Brackbill

,

J. U.

,

Kothe

,

D. B.

, and

Zemach

,

C.

,

1992

, “

A Continuum Method for Modeling Surface Tension

,”

J. Comput. Phys.

,

100

(

2

), pp.

335

–

354

.

31.

Peskin

,

C. S.

,

1977

, “

Numerical Analysis of Blood Flow in the Heart

,”

J. Comput. Phys.

,

25

(

3

), pp.

220

–

252

.

32.

Nangia

,

N.

,

Johansen

,

H.

,

Patankar

,

N. A.

, and

Bhalla

,

A. P. S.

,

2017

, “

A Moving Control Volume Approach to Computing Hydrodynamic Forces and Torques on Immersed Bodies

,”

J. Comput. Phys.

,

347

, pp.

437

–

462

.

33.

MacNeice

,

P.

,

Olson

,

K. M.

,

Mobarry

,

C.

,

De Fainchtein

,

R.

, and

Packer

,

C.

,

2000

, “

Paramesh: A Parallel Adaptive Mesh Refinement Community Toolkit

,”

Comput. Phys. Commun.

,

126

(

3

), pp.

330

–

354

.

34.

Harlow

,

F. H.

, and

Welch

,

J. E.

,

1965

, “

Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid With Free Surface

,”

Phys. Fluids

,

8

(

12

), pp.

2182

–

2189

.

35.

Anderson

,

J. D.

,

2010

,

Computational Fluid Dynamics: The Basics with Applications

,

McGraw-Hill

,

New York

.

36.

Rider

,

W. J.

, and

Kothe

,

D. B.

,

1998

, “

Reconstructing Volume Tracking

,”

J. Comput. Phys.

,

141

(

2

), pp.

112

–

152

.

37.

Weymouth

,

G. D.

, and

Yue

,

D. K. -P.

,

2010

, “

Conservative Volume-of-Fluid Method for Free-Surface Simulations on Cartesian-Grids

,”

J. Comput. Phys.

,

229

(

8

), pp.

2853

–

2865

.

38.

Tryggvason

,

G.

,

Scardovelli

,

R.

, and

Zaleski

,

S.

,

2011

,

Direct Numerical Simulations of Gas–Liquid Multiphase Flows

,

Cambridge University Press

,

New York

.

39.

Xiao

,

F.

,

Honma

,

Y.

, and

Kono

,

T.

,

2005

, “

A Simple Algebraic Interface Capturing Scheme Using Hyperbolic Tangent Function

,”

Int. J. Numer. Methods Fluids

,

48

(

9

), pp.

1023

–

1040

.

40.

Xiao

,

F.

,

Ii

,

S.

, and

Chen

,

C.

,

2011

, “

Revisit to the Thinc Scheme: A Simple Algebraic Vof Algorithm

,”

J. Comput. Phys.

,

230

(

19

), pp.

7086

–

7092

.

41.

Yokoi

,

K.

,

2007

, “

Efficient Implementation of Thinc Scheme: A Simple and Practical Smoothed Vof Algorithm

,”

J. Comput. Phys.

,

226

(

2

), pp.

1985

–

2002

.

42.

Ii

,

S.

,

Sugiyama

,

K.

,

Takeuchi

,

S.

,

Takagi

,

S.

,

Matsumoto

,

Y.

, and

Xiao

,

F.

,

2012

, “

An Interface Capturing Method with a Continuous Function: The Thinc Method with Multi-Dimensional Reconstruction

,”

J. Comput. Phys.

,

231

(

5

), pp.

2328

–

2358

.

43.

Chorin

,

A. J.

,

1968

, “

Numerical Solution of the Navier-Stokes Equations

,”

Math. Comput.

,

22

(

104

), pp.

745

–

762

.

44.

Wesseling

,

P.

,

1992

,

Introduction To Multigrid Methods

,

John Wiley and Sons

,

New York

.

45.

Trottenberg

,

U.

,

Oosterlee

,

C. W.

, and

Schuller

,

A.

,

2000

,

Multigrid

,

Elsevier

,

San Diego, CA

.

46.

Kothe

,

D.

, and

Mjolsness

,

R.

,

1992

, “

Ripple-a New Model for Incompressible Flows With Free Surfaces

,”

AIAA. J.

,

30

(

11

), pp.

2694

–

2700

.

47.

Bussmann

,

M.

,

Kothe

,

D. B.

, and

Sicilian

,

J. M.

,

2002

, “

Modeling High Density Ratio Incompressible Interfacial Flows

,”

Fluids Engineering Division Summer Meeting

, Vol.

36150

,

Montreal, Quebec, Canada

, pp.

707

–

713

.

48.

Rudman

,

M.

,

1998

, “

A Volume-Tracking Method for Incompressible Multifluid Flows With Large Density Variations

,”

Int. J. Numer. Methods Fluids

,

28

(

2

), pp.

357

–

378

.

49.

Pathak

,

A.

, and

Raessi

,

M.

,

2016

, “

A 3D Fully Eulerian, VOF-Based Solver to Study the Interaction Between Two Fluids and Moving Rigid Bodies Using the Fictitious Domain Method

,”

J. Comput. Phys.

,

311

, pp.

87

–

113

.

50.

Fuster

,

D.

,

Arrufat

,

T.

,

Crialesi-Esposito

,

M.

,

Ling

,

Y.

,

Malan

,

L.

,

Pal

,

S.

,

Scardovelli

,

R.

,

Tryggvason

,

G.

, and

Zaleski

,

S.

,

2018

, “

A Momentum-Conserving, Consistent, Volume-of-Fluid Method for Incompressible Flow on Staggered Grids

,” arXiv preprint arXiv:1811.12327.

51.

Popinet

,

S.

,

2009

, “

An Accurate Adaptive Solver for Surface-Tension-Driven Interfacial Flows

,”

J. Comput. Phys.

,

228

(

16

), pp.

5838

–

5866

.

52.

Renardy

,

Y.

, and

Renardy

,

M.

,

2002

, “

Prost: a Parabolic Reconstruction of Surface Tension for the Volume-of-Fluid Method

,”

J. Comput. Phys.

,

183

(

2

), pp.

400

–

421

.

53.

Liovic

,

P.

,

Francois

,

M.

,

Rudman

,

M.

, and

Manasseh

,

R.

,

2010

, “

Efficient Simulation of Surface Tension-Dominated Flows Through Enhanced Interface Geometry Interrogation

,”

J. Comput. Phys.

,

229

(

19

), pp.

7520

–

7544

.

54.

Cummins

,

S. J.

,

Francois

,

M. M.

, and

Kothe

,

D. B.

,

2005

, “

Estimating Curvature From Volume Fractions

,”

Comput. Struct.

,

83

(

6–7

), pp.

425

–

434

.

55.

Francois

,

M. M.

,

Cummins

,

S. J.

,

Dendy

,

E. D.

,

Kothe

,

D. B.

,

Sicilian

,

J. M.

, and

Williams

,

M. W.

,

2006

, “

A Balanced-force Algorithm for Continuous and Sharp Interfacial Surface Tension Models Within a Volume Tracking Framework

,”

J. Comput. Phys.

,

213

(

1

), pp.

141

–

173

.

56.

Lopez

,

J.

,

Zanzi

,

C.

,

Gomez

,

P.

,

Zamora

,

R.

,

Faura

,

F.

, and

Hernandez

,

J.

,

2009

, “

An Improved Height Function Technique for Computing Interface Curvature From Volume Fractions

,”

Comput. Methods. Appl. Mech. Eng.

,

198

(

33–36

), pp.

2555

–

2564

.

57.

Fadlun

,

E.

,

Verzicco

,

R.

,

Orlandi

,

P.

, and

Mohd-Yusof

,

J.

,

2000

, “

Combined Immersed-Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations

,”

J. Comput. Phys.

,

161

(

1

), pp.

35

–

60

.

58.

Sun

,

X.

, and

Sakai

,

M.

,

2016

, “

Numerical Simulation of Two-Phase Flows in Complex Geometries by Using the Volume-of-Fluid/Immersed-Boundary Method

,”

Chem. Eng. Sci.

,

139

, pp.

221

–

240

.

59.

Sun

,

X.

, and

Sakai

,

M.

,

2016

, “

Direct Numerical Simulation of Gas-Solid-Liquid Flows With Capillary Effects: An Application to Liquid Bridge Forces Between Spherical Particles

,”

Phys. Rev. E

,

94

(

6

), p.

063301

.

60.

Washino

,

K.

,

Tan

,

H.

,

Salman

,

A.

, and

Hounslow

,

M.

,

2011

, “

Direct Numerical Simulation of Solid–Liquid–Gas Three-Phase Flow: Fluid–Solid Interaction

,”

Powder. Technol.

,

206

(

1–2

), pp.

161

–

169

.

61.

Sussman

,

M.

,

2001

, “

An Adaptive Mesh Algorithm for Free Surface Flows in General Geometries

,”

Adaptive Method Lines

, pp.

207

–

231

.

62.

Patel

,

H.

,

Das

,

S.

,

Kuipers

,

J.

,

Padding

,

J.

, and

Peters

,

E.

,

2017

, “

A Coupled Volume of Fluid and Immersed Boundary Method for Simulating 3d Multiphase Flows With Contact Line Dynamics in Complex Geometries

,”

Chem. Eng. Sci.

,

166

, pp.

28

–

41

.

63.

Yokoi

,

K.

,

Vadillo

,

D.

,

Hinch

,

J.

, and

Hutchings

,

I.

,

2009

, “

Numerical Studies of the Influence of the Dynamic Contact Angle on a Droplet Impacting on a Dry Surface

,”

Phys. Fluids.

,

21

(

7

), p.

072102

.

64.

O’Brien

,

A.

, and

Bussmann

,

M.

,

2020

, “

A Moving Immersed Boundary Method for Simulating Particle Interactions At Fluid-Fluid Interfaces

,”

J. Comput. Phys.

,

402

, p.

109089

.

65.

Fujita

,

M.

,

Koike

,

O.

, and

Yamaguchi

,

Y.

,

2013

, “

Computation of Capillary Interactions Among Many Particles At Free Surface

,”

Appl. Phys. Express.

,

6

(

3

), p.

036501

.

66.

Weller

,

H. G.

,

2008

, “

A New Approach to VOF-Based Interface Capturing Methods for Incompressible and Compressible Flow

,”

OpenCFD Ltd. Report TR/HGW

,

4

, p.

35

.

67.

Lafaurie

,

B.

,

Nardone

,

C.

,

Scardovelli

,

R.

,

Zaleski

,

S.

, and

Zanetti

,

G.

,

1994

, “

Modelling Merging and Fragmentation in Multiphase Flows with Surfer

,”

J. Comput. Phys.

,

113

(

1

), pp.

134

–

147

.

68.

Das

,

S.

,

Patel

,

H.

,

Milacic

,

E.

,

Deen

,

N.

, and

Kuipers

,

J.

,

2018

, “

Droplet Spreading and Capillary Imbibition in a Porous Medium: A Coupled IB-VOF Method Based Numerical Study

,”

Phys. Fluids.

,

30

(

1

), p.

012112

.

69.

Bussmann

,

M.

,

Chandra

,

S.

, and

Mostaghimi

,

J.

,

2000

, “

Modeling the Splash of a Droplet Impacting a Solid Surface

,”

Phys. Fluids.

,

12

(

12

), pp.

3121

–

3132

.

70.

Brakke

,

K. A.

,

1992

, “

The Surface Evolver

,”

Exp. Math.

,

1

(

2

), pp.

141

–

165

.

71.

Pitois

,

O.

,

Moucheront

,

P.

, and

Chateau

,

X.

,

2000

, “

Liquid Bridge Between Two Moving Spheres: An Experimental Study of Viscosity Effects

,”

J. Colloid. Interface. Sci.

,

231

(

1

), pp.

26

–

31

.

72.

Wang

,

J.-P.

,

Gallo

,

E.

,

François

,

B.

,

Gabrieli

,

F.

, and

Lambert

,

P.

,

2017

, “

Capillary Force and Rupture of Funicular Liquid Bridges Between Three Spherical Bodies

,”

Powder. Technol.

,

305

, pp.

89

–

98

.

73.

Pitois

,

O.

,

Moucheront

,

P.

, and

Chateau

,

X.

,

2001

, “

Rupture Energy of a Pendular Liquid Bridge

,”

Eur. Phys. J. B-Condensed Matter Complex Syst.

,

23

(

1

), pp.

79

–

86

.

74.

Washino

,

K.

,

Chan

,

E. L.

,

Matsumoto

,

T.

,

Hashino

,

S.

,

Tsuji

,

T.

, and

Tanaka

,

T.

,

2017

, “

Normal Viscous Force of Pendular Liquid Bridge between Two Relatively Moving Particles

,”

J. Colloid. Interface. Sci.

,

494

, pp.

255

–

265

.

75.

Matthewson

,

M.

,

1988

, “

Adhesion of Spheres by Thin Liquid Films

,”

Philos. Mag. A

,

57

(

2

), pp.

207

–

216

.

Computation of Hydrodynamic and Capillary Phenomena in Binder Jet Three-Dimensional Printing

Abstract

1 Introduction

2 Governing Equations

3 Numerical Methodology

3.1 Volume-of-Fluid Interface Advection.

3.2 Incompressible Flow Solver

3.2.1 Time Integration.

3.2.2 Spatial Discretization.

3.3 Surface Tension.

3.4 Immersed Boundary Method.

3.5 Three-Phase Contact Line Dynamics.

3.6 Fluid-Induced Forces.

4 Results

4.1 Binder Droplet Impact on a Fixed Powder Bed.

4.2 Validation of Hydrostatic Liquid Bridge Forces.

4.2.1 Binary Particle System (pendular).

4.2.2 Ternary Particle System (funicular).

4.3 Validation of Hydrodynamic Liquid Bridge Forces.

5 Conclusions

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

References

Get Email Alerts

Cited By

Journals

Conference Proceedings

eBooks

Resources

Opportunities

Computation of Hydrodynamic and Capillary Phenomena in Binder Jet Three-Dimensional Printing

Abstract

1 Introduction

2 Governing Equations

3 Numerical Methodology

3.1 Volume-of-Fluid Interface Advection.

3.2 Incompressible Flow Solver

3.2.1 Time Integration.

3.2.2 Spatial Discretization.

3.3 Surface Tension.

3.4 Immersed Boundary Method.

3.5 Three-Phase Contact Line Dynamics.

3.6 Fluid-Induced Forces.

4 Results

4.1 Binder Droplet Impact on a Fixed Powder Bed.

4.2 Validation of Hydrostatic Liquid Bridge Forces.

4.2.1 Binary Particle System (pendular).

4.2.2 Ternary Particle System (funicular).

4.3 Validation of Hydrodynamic Liquid Bridge Forces.

5 Conclusions

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only