Computational Modeling of Droplet-Based Printing Using Multiphase Volume of Fluid (VOF) Method: Prediction of Flow, Spread Behavior, and Printability

Abstract

The evolution of droplets during the printing process is modeled using the volume of fluid (VOF) method, which involves solving the Navier–Stokes and continuity equations for incompressible flow with multiple immiscible phases on a finite volume grid. An indicator function tracks the interfaces and calculates surface tension forces. A grid independence study confirmed the convergence and efficacy of the solutions. The computational model agreed well with experimental data, accurately capturing the impact, spreading, and recoiling of droplets on a solid surface. Additionally, the model validated the interaction of droplets with hydrophilic and hydrophobic surfaces for both constant and dynamic contact angles. Key non-dimensional numbers (Re, We, Oh) were considered to study the interplay of forces during droplet impact on a solid surface. The final print quality is influenced by droplet dynamics, governed by body forces (surface tension, gravity), contact angle, dissipative forces due to motion, and material properties. Computational studies provide insights into the overall process performance and final print quality under various process conditions and material properties.

1. Introduction

The dynamics of droplet behavior and its spreading mechanism have been extensively explored through experimental and computational methods. The droplets interact with the printing substrate in various ways, including spreading, rebounding, recoiling, and splashing, as reviewed by Yarin [1]. These interactions vary depending on the process conditions and material properties. Despite significant research on droplet spreading simulations, there has been less emphasis on studying the recoiling behavior. This study aims to fill that gap by providing a comprehensive three-dimensional (3D) analysis of droplet impact, spreading, recoiling, and interaction with different surfaces.

The deposition process is characterized by the ‘Spreading Factor’, which is defined as the ratio of the maximum spread diameter to the initial drop diameter. This factor directly reflects the printability and quality of droplet-based printing. Experimenting with droplet-based printing can be challenging, costly, and cumbersome. However, computational fluid dynamics (CFD) analyses offer a practical and efficient alternative for exploring the potential applications of droplet-based printing in additive manufacturing. Although numerical methods for accurately representing the interfaces between different phases are still evolving, they hold promise for investigating droplet behavior and predicting the outcomes of droplet-based printing.

Limited research has been conducted on the oscillatory behavior of droplets after colliding with a solid surface. Elliot and Ford [2] studied how droplets spread and recoil upon collision, observing changes in contact angles until equilibrium is reached. Similarly, Fukai et al. [3] used both experimental and theoretical models to analyze the fluid mechanics of droplet impact and spreading on a solid surface. While some authors have reported empirical and semi-empirical relations for the spreading diameter based on fluid properties and process conditions, these studies often exclude surface characteristics and dynamic contact angles [4,5,6,7,8]. Pasandideh-Fard et al. [9] and Feng [10] presented improved empirical relations incorporating advancing contact angles and energy conservation factors for predicting the final spreading diameter.

Various modeling techniques are available for numerical investigations of free surface flow with deformations, including explicit interface tracking and implicit Eulerian interface capturing methods. The implicit Eulerian methods, particularly the volume of fluid (VOF) model, are highly effective in capturing flows with deformations, including interface breakups. Most numerical studies employ the VOF model alongside the solution of Navier–Stokes equations and a dynamic contact angle model to accurately depict droplet impact and spreading. For instance, Adeniyi et al. [11] used a coupled Euler–Lagrange approach to model droplet-to-film interactions in aeroengine bearing chambers, highlighting the challenges of modeling droplet impact phenomena in finite-volume CFD techniques.

Pasandideh-Fard et al. [9] conducted simulations using VOF and finite difference Navier–Stokes solvers, concluding that dynamic contact angles yield more accurate results than constant or equilibrium contact angles. Yokoi et al. [12] used the finite volume method and VOF model, incorporating a switching function for dynamic contact angles based on the velocity of the interface normal near the wall. Gujjula [13] employed similar methods with different dynamic contact angle models, comparing their results to those of Pasandideh-Fard et al. [9] and finding that the spread diameter was underestimated. Wang et al. [14] introduced a morphology-adaptive multifield two-fluid model to simulate annular flow, proposing a new droplet entrainment model based on the shear-off entrainment mechanism at the gas–liquid interface.

The present study also examines the influence of various forces (surface tension, gravitational forces, contact angle forces, and other dissipative forces), material properties (density, viscosity, and surface tension), and processing conditions (droplet size, impact velocity, and nozzle-substrate distance) on the spreading and recoiling behavior of droplets, which dictate the final printing resolution. Building on previous work by Shah and Mohan [15], the present work extends the computational domain to three dimensions. It investigates the effects of relevant non-dimensional numbers (Reynolds, Weber, and Ohnesorge numbers) on droplet dynamics to analyze the interplay of different forces in droplet-based printing. In related work, Vontas [16] explored the practical applications of droplet impact in fuel injection systems, using enhanced VOF-based CFD simulations to model the interaction of fuel spray droplets with solid surfaces in combustion chambers.

During continuous printing, the print nozzle moves in a traverse direction depending on the print speed and deposits the subsequent droplets. The print setup and traverse speed decide the position of consequent droplets deposited on the substrate. It is essential to study the dynamics of droplets being printed, with the droplet interaction with the substrate and other droplets. The droplet–droplet interaction and the dynamics of coalescence determine the final print resolution for continuous droplet-based printing. The coalescence mechanism involves rupture, deformation, and interface merging of the droplets and needs further investigation, focusing on droplet-based printing.

Applications:

Droplet-based printing has various applications across various fields due to its precision and versatility. In biomedical engineering, it is used for bioprinting to fabricate complex 3D structures with living cells for tissue engineering and regenerative medicine, and to create drug delivery systems with controlled release properties [17,18]. In microfluidics, droplet-based printing is employed to fabricate lab-on-a-chip devices for biochemical analyses and to create monodisperse emulsions for chemical reactions and material synthesis [19,20]. In flexible electronics, it produces printed sensors for wearable devices and prints conductive inks for flexible circuits and displays [21,22]. In the aerospace and automotive industries, droplet-based printing is utilized to manufacture lightweight components, improve thermal management systems, and optimize fuel injection processes [23,24,25]. These applications highlight technology’s potential to enhance manufacturing processes, improve product performance, and enable innovative solutions across various sectors.

This paper is structured to provide a detailed exploration of droplet dynamics in the context of droplet-based printing. The main content is organized as follows: Section 2 outlines the numerical methods employed, including the governing equations and computational framework. Section 3 presents the results of the computational model, validated against experimental studies, and discusses the effects of key non-dimensional numbers (Re, We, Oh) on droplet dynamics. Section 4 delves into the fluid–fluid interaction and final print resolution prediction, examining droplets’ coalescence behavior during continuous printing. Finally, Section 5 offers a comprehensive discussion of the findings, highlighting the implications for additive manufacturing and potential areas for future research. The conclusions summarize the key insights gained from the study and propose directions for further investigation.

2. Numerical Method

2.1. Governing Equations

The Navier–Stokes equations describe the behavior of an incompressible fluid and elucidate the dynamics of a droplet impinging upon a solid surface. This is complemented by the volume fraction function, denoted as the indicator function α. The pertinent equations governing this isothermal modeling framework are delineated as follows:

Continuity equation:

(1)$\begin{array}{c}\nabla ·U=0\end{array}$

Momentum equation:

(2)$\begin{array}{c}\rho {\displaystyle \frac{DU}{Dt}}=-\nabla p+\rho g+\nabla ·\mu \left(\nabla U+\nabla U^T\right)+F_{\sigma }\end{array}$

Conservation of phase fraction:

(3)$\begin{array}{c}{\displaystyle \frac{D\alpha }{Dt}}=0\end{array}$

where U is the velocity vector, t is time, p is the lumped pressure, g is the acceleration due to gravity, and F_σ are the body (surface tension and gravity) forces

In this computational framework, the gas–liquid interface is conceptualized as a dispersed layer with a thickness commensurate with a finite volume cell. This interface is characterized by the indicator function, wherein the gradient of the volume fraction α delineates the phase discontinuity at the interface. The volume fraction, continuously distributed within the domain (0 ≤ α ≤ 1) except at the interface, is instrumental in determining phase properties such as distributed density and viscosity. In fluid modeling, the distributed density is typically calculated using the volume fraction of the phases involved. For a two-phase system, the density and viscosity at any point in the domain are a weighted average of the densities and viscosities of the individual phases, respectively, based on their corresponding volume fractions. Mathematically, this can be expressed as follows:

(4)$\begin{array}{c}\rho ={\rho }_d\alpha +{\rho }_g\left(1-\alpha \right)\end{array}$

(5)$\begin{array}{c}\mu ={\mu }_d\alpha +{\mu }_g\left(1-\alpha \right)\end{array}$

where α is ‘1’ for the liquid phase, ‘0’ for the gas/air phase, and takes an intermediate value at the interface.

As the resolution of the interface is equal to the size of the finite volume cell, it needs to be sharpened in the momentum equation. The ‘interfoam’ solver (OpenFOAM^® v2021) improves the resolution of the interface using an additional compression velocity term in the conservation equation of the phase fraction (Equation (3)) [26]. The surface tension is considered as body force in the momentum equation and modeled as follows:

(6)$\begin{array}{c}{F}_{\sigma }=\sigma kn=\sigma \left(\nabla n\right)·n=\sigma \nabla \left({\displaystyle \frac{\nabla \alpha }{\left|\nabla \alpha \right|}}\right)\end{array}$

where σ is the surface tension coefficient, k is the interfacial curvature, and ‘n’ is the interfacial normal.

The transient terms in the Navier–Stokes equations are discretized using the Crank–Nicolson scheme, a second-order implicit method known for its stability and accuracy in time integration. The time step is determined based on the specified Courant–Friedrichs–Lewy (CFL) number, ensuring numerical stability by maintaining the CFL condition [16]. Meanwhile, the spatial terms are discretized using a second-order central differencing scheme, which provides a balance between accuracy and computational efficiency by approximating derivatives at the cell centers. This method is particularly effective in minimizing numerical diffusion and maintaining sharp gradients.

The pressure gradient is solved using the conjugate gradient method, an efficient iterative solver for large, sparse linear systems typically encountered in computational fluid dynamics (CFDs). This method leverages the symmetry and positive definiteness of the pressure matrix to converge rapidly. Pressure–velocity coupling is achieved through the pressure implicit with splitting of operators (PISO) method. The PISO method is an extension of the semi-implicit method for pressure-linked equations (SIMPLE), designed to enhance the convergence rate and accuracy of the solution. In the PISO algorithm, the corrected velocities are used to update the pressure field iteratively. This iterative process ensures that the pressure and velocity fields satisfy the momentum and continuity equations, achieving a consistent and accurate solution. The PISO method addresses the challenge of not having a separate equation for the pressure term in the Navier–Stokes equations. Instead, it updates the pressure field after correcting the velocities to satisfy the continuity equation, followed by further corrections using the momentum equation. This iterative correction process ensures the pressure–velocity coupling is accurately resolved.

Furthermore, the transport and volume of fluid (VOF) equations are solved using the Gauss–Seidel algorithm, an iterative method that solves linear systems by successively updating the solution vector. This algorithm is particularly effective for solving the discretized transport equations in CFD.

The multidimensional universal limiter for explicit solution (MULES) technique is employed to achieve interface compression in the VOF method [26]. MULES is designed to maintain a sharp interface between immiscible fluids by limiting numerical diffusion and ensuring that the volume fraction remains between 0 and 1. This technique is crucial for accurately capturing the interface dynamics in multiphase flow simulations.

By incorporating these advanced numerical methods and algorithms, the solver ensures robust and accurate solutions to the Navier–Stokes equations, enabling detailed and reliable fluid dynamics simulations in various engineering applications. Several researchers have utilized the ‘interFoam’ VOF solver provided by the open-source CFD software OpenFOAM^® in their studies [26]. This solver offers numerous benefits, including enhanced representation of the interface, an efficient advection scheme, and the ability to handle large density ratios. These features make ‘interFoam’ particularly suitable for simulating multiphase flows with complex interface dynamics.

Fujimoto et al. [27] achieved good agreement between their numerical simulations and experimental data by utilizing a three-dimensional finite volume solver with the VOF method. Their study incorporated a switching dynamic contact angle approach, similar to that used by Yokoi et al. [12]. This approach allows for more accurate modeling of contact angle hysteresis, which is crucial for simulating droplet behavior on solid surfaces. The dynamic contact angle model incorporates the contact line velocity as a critical parameter, thereby accounting for the influence of flow velocity and accurately capturing the behavior of droplets on both hydrophilic and hydrophobic surfaces. Feng [10] also employed the finite volume method with the VOF model using the ‘interFoam’ solver. Their results were comparable to empirical correlations, demonstrating the solver’s capability to accurately predict fluid dynamics in multiphase systems. The study highlighted the solver’s robustness in handling complex flow phenomena, such as droplet impact and spreading.

The present work focuses on analyzing the behavior of droplets upon impact and spreading using the ‘interFoam’ solver provided by OpenFOAM^®. The main objective is to investigate droplets’ spreading and recoiling behavior, particularly in the context of additive manufacturing, and to predict the final spread diameter. Such analysis involves detailed numerical simulations to capture the intricate dynamics of droplet deformation, spreading, and recoiling, which are critical for optimizing processes in additive manufacturing.

By leveraging the advanced capabilities of the ‘interFoam’ solver, this study aims to provide deeper insights into the fluid mechanics of droplet interactions with surfaces. The findings are expected to contribute to developing more efficient and precise additive manufacturing techniques, where control over droplet behavior is essential for achieving high-quality prints.

2.2. Computational Domain

The computational domain and associated boundary conditions are depicted in Figure 1. The domain consists of a rectangular block representing the air phase, while the dark blue region indicates the water droplet. The bottom wall of the domain represents a solid substrate with ‘no slip’ boundary conditions, ensuring that the fluid velocity at the wall is zero.

Constant and dynamic contact angle boundary conditions are applied at the bottom wall to accurately capture the fluid–solid interaction. The dynamic contact angle model from the OpenFOAM^® library is utilized, which is described by the following equation:

(7)$\begin{array}{c}\theta ={\theta }_e+\left({\theta }_a-{\theta }_r\right)\tanh \left({\displaystyle \frac{U_w}{U_{\theta }}}\right)\end{array}$

where ${\theta }_e$, ${\theta }_a$, and ${\theta }_r$ represent the equilibrium, advancing, and receding contact angles, respectively, $U_w$ denotes the normal velocity of the interface near the wall. and $U_{\theta }$ is the characteristic velocity scale. Feng [10] validated with the experimental data that the model successfully predicts the contact angle evolution on various surfaces, including oscillations and other dynamic behaviors. The current study follows the empirical relation reported by Feng [10] with a simplified approach by specifying a static contact angle (${\theta }_0$) at the contact line with ${\theta }_a$ = ${\theta }_0$ + 5°, ${\theta }_r$ = ${\theta }_0$ − 5°, and U_w = 1 m/s.

The top wall and four vertical walls are treated as open to the atmosphere faces with the ‘PressureInletOutletVelocity’ boundary conditions for velocities and pressure, as specified by the OpenFOAM CFD solver [26]. The pressure settings are maintained at constant atmospheric pressure, ensuring that the pressure at the inlets/outlets remains stable and does not fluctuate. The velocity settings are zero uniform velocities, meaning there is no net flow through these boundaries. This setup allows for handling the open-to-atmosphere boundary conditions, allowing flow to freely enter and exit the domain.

The rectangular domain is of size 10 × 10 × 10 mm, which is divided into 125 × 125 × 75 finite volume cells after a convergence study is conducted for the present numerical scheme. A domain independence study is conducted to ensure that the computational results are unaffected by the dimensions of the computational domain. This study evaluated multiple domain sizes and systematically compared outcomes to verify that the selected dimensions yield accurate and stable solutions. Given that the initial droplet diameter ranged from 2.5 to 3.5 mm, we assessed various domain sizes to capture critical aspects of droplet dynamics, including maximum diameter spread and droplet velocity. The grid independence study confirmed that a domain size of 10 × 10 × 10 mm is adequate to encapsulate the fundamental dynamics of droplet behavior, ensuring reliable and precise computational results. A comprehensive grid independence study was conducted to ensure the accuracy and reliability of the simulation results, involving the evaluation of computational performance and solution convergence across multiple grid resolutions. A detailed summary is provided in Appendix A. The governing equations and required numerical methods are well described in Section 2.1.

3. Results

The computational method is initially validated by comparing the results with previously reported experimental studies by Kim and Chun [28], which investigated the interaction of water droplets with polycarbonate surfaces. Their study provided a comprehensive analysis of droplet dynamics, including spreading and recoiling behaviors, using high-speed video systems to capture the intricate details of the interactions.

We have conducted a validation study and performed a parametric analysis for a more comprehensive understanding of droplet dynamics by incorporating additional relevant non-dimensional numbers, such as the Reynolds (Re), Weber (We), and Ohnesorge (Oh) numbers, and their effects on droplet dynamics. These dimensionless parameters are crucial for characterizing the fluid dynamics of droplet impact and spreading, providing valuable insights into the underlying mechanism. By incorporating these non-dimensional numbers, the study aims to provide a more detailed and accurate prediction of droplet dynamics, particularly in additive manufacturing. The analysis focuses on droplets’ spreading and recoiling behavior, aiming to predict the final spread diameter and understand the influence of various fluid properties and impact conditions.

In Section 3.1, we present the validation studies for the droplet impact on a hydrophilic and hydrophobic surface, implicating the efficacy of the model to capture the complex physics of the droplet impact and spreading.

3.1. Flow and Spreading Behavior Method Validation: Impact of a Water Droplet on a Hydrophilic Surface

3.1.1. Method Validation: Impact of a Water Droplet on a Hydrophilic Surface

The numerical results were validated by comparing them with the experimental study reported by Kim and Chun [28]. In their study, a water droplet with a diameter of 3.6 mm was made to impact a polycarbonate surface at an initial velocity of 0.77 m/s. The experimental setup consisted of a solid surface onto which liquid droplets were projected using an inkjet printer at varying velocities. The recoiling behavior of the droplets upon impact was recorded using high-speed cameras and subsequently analyzed using advanced image processing techniques. Detailed descriptions of the experimental procedures can be found in the work published by Kim and Chun [28].

The equilibrium contact angle between the water and polycarbonate surface was measured to be θ₀ = 87.4°. The dynamic contact angles, including the advancing and receding angles, were specified as θ₀ + 5° and θ₀ − 5°. These contact angles are critical parameters for accurately modeling the fluid–solid interaction and predicting the droplet dynamics.

The results of the droplet shape evolution over time were compared with the experimental observations, as illustrated in Figure 2. This comparison demonstrated good agreement between the numerical simulations and the experimental data, validating the accuracy of the computational model. The results highlight the importance of accurately capturing the contact angle dynamics to predict the spreading and recoiling behavior of droplets upon impact.

To further validate the computational results, Figure 3 presents the superimposed images of the simulated droplet shape extracted at the same time intervals as the experimental images. This overlay allows for a direct visual comparison, highlighting discrepancies or agreements between the simulation and experimental results. The validation study demonstrated a strong agreement between the present numerical simulations and the experimental data reported by Kim and Chun [28]. The superimposed droplet shapes from the simulations closely matched the experimental images at various time steps, indicating the accuracy and fidelity of the computational model. These findings validate the effectiveness of the ‘interFoam’ solver in capturing the complex dynamics of droplet impact and spreading, providing a robust foundation for further investigations in this domain.

As illustrated in Figure 2 and Figure 3, the simulation accurately captures both the spreading and recoiling phases during the shape evolution of a water droplet impinging on a polycarbonate surface. The spreading phase is characterized by a rapid expansion, reaching the maximum spread diameter at 11.6 ms, and is marked by the formation of a rim at the furthest point from the droplet’s center. Conversely, the recoiling phase is slower and exhibits an oscillatory behavior in the spread diameter before stabilizing at equilibrium at 39.8 ms.

This study effectively captures the droplet’s detailed shape evolution and dynamic behavior throughout the impact process. The qualitative analysis of the shape evolution is complemented by a quantitative analysis of the droplet’s spread diameter over time. For the quantitative analysis, we use normalized base and normalized time as defined below:

Normalized base: $D_b^{\mathrm{*}}={\displaystyle \frac{D_b}{D_0}}$

Normalized time: $t^{\mathrm{*}}={\displaystyle \frac{t\sqrt{\sigma }}{\sqrt{\rho D_o^3}}}$

The normalized base is plotted against the normalized time to further elucidate the dynamics. This plot provides insights into the orders of magnitude changes in the spreading diameter over time, highlighting the intricate interplay between inertial, viscous, and surface tension forces during the droplet’s impact and subsequent evolution. Figure 4 below shows the variation in the normalized base with normalized time for constant and dynamic contact angle models for comparison with the experimental results reported by Kim and Chun [28].

The constant contact angle model provides a reasonable prediction of the initial spreading behavior of the droplet; however, it tends to overpredict the recoiling behavior. In contrast, the dynamic contact angle model better captures the initial spreading and recoiling behaviors. Additionally, it more accurately represents the post-recoiling oscillations compared to the constant contact angle model.

The computational model, while slightly overpredicting the equilibrium base diameter, as observed from the qualitative shape evolution analysis, effectively captures the overall dynamics of the droplet. The predicted dynamics include the intricate details of the spreading and recoiling phases. The dynamic contact angle model’s ability to better capture the transient behaviors and oscillations highlights its robustness and accuracy in simulating droplet dynamics.

Overall, the computational results provide a comprehensive and effective prediction of the spreading and recoiling behavior of a water droplet impinging on a polycarbonate surface. This validation underscores the model’s capability to accurately simulate complex fluid–solid interactions, which is crucial for applications in various engineering and scientific studies.

3.1.2. Method Validation: Impact of a Water Droplet on a Hydrophobic Surface

After examining the droplet dynamics on a hydrophilic surface in the previous section, this study extends the analysis to a hydrophobic surface. The computational model is tested for a scenario where a water droplet impacts a hydrophobic surface. Azimi et al. [29] investigated the nanostructured modification of material surfaces and the induced hydrophobicity of rare-earth oxide ceramics. They fabricated ceramic nanostructures, consisting of a thin layer on cubical silicon micro-posts, which exhibited significant hydrophobicity with contact angles of 160 ± 2°, as illustrated in Figure 5.

In the computational domain, a water droplet with a diameter of 2.5 mm is made to fall on the hydrophobic surface with an initial velocity of 1.6 m/s. The bottom wall is modeled as a hydrophobic surface with a contact angle of 160°, while the other boundary conditions remain consistent with those described in Section 2.1. Due to the hydrophobic nature of the nanostructured surface, the water droplet bounces off upon impact. The computational model accurately captures this behavior, as shown in Figure 5.

The study demonstrates the model’s capability to simulate the complex interactions between droplets and hydrophobic surfaces. The results provide valuable insights into the influence of surface properties on droplet dynamics, which are crucial for applications in surface engineering, coating technologies, and fluid dynamics research.

The comparison of the base diameter variation with time between the experimental data and the simulations is illustrated in Figure 6. The results demonstrate a high degree of correlation, with the simulations accurately capturing the temporal evolution of the base diameter, thereby validating the computational model’s effectiveness in predicting droplet dynamics on hydrophobic surfaces.

Based on the experimental data reported by Azimi et al. [29], the base diameter of the droplet increases post-impact, reaching a maximum of 5.71 mm at 2.48 ms before decreasing as the droplet bounces off the surface at 20.8 ms. The experiments did not document or report the subsequent rebound and oscillatory behavior stages. However, the computational model replicates the observed trend, with the dynamic contact angle model providing more accurate results than the constant contact angle model. The constant contact angle model tends to overpredict the base diameter. In contrast, the dynamic contact angle model offers better predictions due to its more sophisticated approximation of the contact line velocity at the solid–liquid interface.

Overall, the multiphase CFD solver effectively simulates a water droplet’s flow and spreading behavior impinging on a polycarbonate surface. These findings are crucial for understanding the variation in base diameter across different fluid–solid interactions and for conducting parametric studies to optimize process parameters. The subsequent section will explore the effects of non-dimensional numbers, such as Reynolds, Weber, and Ohnesorge numbers, to investigate the interplay of various forces associated with this phenomenon.

3.2. Effect of Non-Dimensional Reynolds (Re), Weber (We), and Ohnesorge (Oh) Numbers

The numerical study has examined droplet dynamics in relation to spreading and recoiling behavior and fluid–solid interactions (contact angle). The computational models have demonstrated good agreement with experimental results. This section investigates and discusses the effect of relevant non-dimensional numbers on droplet dynamics. These dimensionless numbers are crucial in fluid dynamics as they elucidate the interplay of different forces associated with flow phenomena.

Key non-dimensional numbers pertinent to droplet dynamics include the Reynolds (Re), Weber (We), and Ohnesorge (Oh) numbers. These numbers are defined as ratios of different forces and are dimensionless. The Reynolds number (Re) is the ratio of inertial to viscous forces and is essential for understanding flow profiles and fall velocity in droplet dynamics. The Weber number (We) represents the ratio of inertial to surface tension forces and characterizes the droplet shape influenced by inertia (impact velocity). The Ohnesorge number (Oh) is the ratio of viscous forces to inertial and surface tension forces. The definitions of Re, We, and Oh numbers are provided below:

$Re={\displaystyle \frac{inertia}{viscous force}}={\displaystyle \frac{\rho uD}{\mu }}$

$We={\displaystyle \frac{inertia}{surface tension force}}={\displaystyle \frac{\rho u^{2}D}{\sigma }}$

$Oh={\displaystyle \frac{viscous force}{\sqrt{inertial.surface tension force}}}={\displaystyle \frac{\mu }{\rho \sigma D}}$

where ρ is the fluid density, U is the characteristic velocity, L is the characteristic length, μ is the dynamic viscosity, and σ is the surface tension. These dimensionless numbers play a pivotal role in characterizing the droplet dynamics and understanding the influence of various forces on the behavior of droplets upon impact.

The case study presented in Section 3.1 analyzes droplet dynamics for different Reynolds numbers (Re) ranging from 1400 to 3000. The focus is primarily on the initial spreading and recoiling behavior, influenced by the interplay of inertial and viscous forces. As depicted in Figure 7, the initial spreading is governed by the Reynolds number, with the spreading increasing as the Reynolds number rises. The peak spread diameter increases with higher Reynolds numbers, indicating that the base diameter spread is directly proportional to the Reynolds number. This behavior is attributed to the increased inertial forces and the more significant impact of the collision force associated with higher Reynolds numbers. However, it is important to note that further increases in the Reynolds number can lead to splash generation, which is undesirable for printing applications.

Interestingly, there is minimal variation in the equilibrium base diameter across different Reynolds numbers, with convergence to a similar equilibrium base diameter despite the varying peaks. These results are consistent with the literature and experimental data, confirming the Reynolds number’s significant role in dictating the droplet’s initial spread. The recoiling and equilibrium stages do not exhibit significant variations with changes in the Reynolds number.

We subsequently study the effect of the Weber number (We) and Ohnesorge number (Oh) on droplet dynamics to understand the interplay of inertial, surface tension, and viscous forces, as shown in Figure 8.

The Weber number (We), which ranges from 40 to 200, is examined for its impact on the droplet’s initial spreading and recoiling phases. The Weber number, the ratio of inertial to surface tension forces, is examined primarily for its impact on the initial spreading and recoiling phases. The Weber number primarily influences the initial spreading phase. As shown in Figure 8a, the normalized base diameter is plotted against the normalized time for different Weber numbers. The peak spread diameter increases as the Weber number rises due to the increased inertial effect. However, the recoiling phase remains relatively unaffected, with recoiling time almost constant across different Weber numbers. Further investigation is required to fully understand the role of the Weber number in determining recoiling behavior. The Weber number is particularly significant in dictating droplet generation at microscales, where the interplay of inertial and surface tension forces is critical.

The Ohnesorge number (Oh) is used to predict the nature of the recoiling phase and analyze the spreading and recoiling behavior of droplet impingement on a solid surface. The Ohnesorge number (Oh) is defined as the ratio of viscous forces to inertial and surface tension forces, and the Ohnesorge number quantifies the relative importance of these forces for a given inertia. Figure 8b compares the base diameter evolution for different Ohnesorge numbers, with all cases having the same inertia, to distinctly compare the effects of viscous and surface tension forces.

Figure 8b shows that the spreading phase exhibits minimal variation across different Ohnesorge numbers, indicating that the Ohnesorge number does not significantly influence the initial spreading when the inertial force is constant. However, the recoiling phase is shorter for lower Ohnesorge numbers as the viscous and surface tension forces directly scale in this phase. The Ohnesorge number determines the resistance force to the recoiling behavior in the droplet impingement phenomenon.

One potential explanation for the counterintuitive behavior observed is the relationship between the Ohnesorge number and surface tension. A higher Oh value indicates stronger viscous forces relative to surface tension forces. In droplet spreading, the dominant forces are viscous and surface tension forces. A higher Ohnesorge number implies that viscous forces play a more significant role, potentially weakening the influence of surface tension in restraining the droplet’s spreading. Additionally, a higher Oh value may correspond to a decrease in the capillary length scale, further diminishing the importance of surface tension relative to viscous forces and promoting spreading.

It is crucial to note that the interaction between the Ohnesorge number and surface tension is complex and can be influenced by factors such as the Weber number. The combined influence of the Ohnesorge and Weber numbers on droplet behavior warrants comprehensive analysis to accurately elucidate the underlying physics. Further investigation and detailed analysis are necessary to understand why a higher Oh value, indicative of more dominant viscous forces, enhances droplet spreading rather than inhibiting it. Exploring the interplay between the Ohnesorge number, Weber number, and surface tension can elucidate the intricate dynamics governing droplet behavior in various scenarios.

In summary, these non-dimensional numbers are significant in determining the interplay of forces associated with droplet dynamics, primarily inertia, surface tension, and viscous forces. The Reynolds and Weber numbers directly affect the spreading phase due to the dominance of inertial forces, while the Ohnesorge number governs the recoiling phase.

3.3. Fluid–Fluid Interaction (Droplet Coalescence) and Prediction of Final Print Resolution in Droplet-Based Printing

During continuous printing, the print nozzle traverses in a direction dictated by the print speed, depositing subsequent droplets accordingly. The print setup and traverse speed determine the precise positioning of the droplets on the substrate. It is crucial to investigate the dynamics of droplets during printing, particularly their interactions with the substrate and other droplets. The interactions between droplets and the dynamics of coalescence are pivotal in determining the final print resolution in continuous droplet-based printing.

The coalescence mechanism, which encompasses rupture, deformation, and interface merging of droplets, warrants further investigation within the context of droplet-based printing. Extensive studies have been conducted on the dynamics of coalescence, focusing on droplet–droplet interactions [30,31,32], interactions with solid surfaces [33,34,35,36], and interactions with liquid films [32,35,37,38]. Various computational approaches have been employed to capture the interfaces during coalescence, including the solution of Navier–Stokes equations and computational fluid dynamics (CFDs) techniques [39,40,41], the lattice Boltzmann method, and molecular dynamics (MDs) simulations [42,43,44,45].

The present work extends these computational approaches by utilizing the volume of fluid (VOF) method to solve flow equations and capture the merging of interfaces driven by surface tension forces. The methodology involves the addition of a droplet, with its center location determined by the print speed or traverse speed. Initially, a sessile droplet is deposited and allowed to impinge upon the substrate, followed by the deposition of a satellite droplet, whose position is also governed by the print speed.

A schematic diagram illustrating the positions of the droplets and the center-to-center distance of the subsequent droplet is presented in Figure 9 below.

3.3.1. Flow and Spreading Behavior of Two Coalescing Droplets

Extending the current studies, the computational model rigorously simulates the dynamics and coalescence of two water droplets impinging upon a polycarbonate substrate to gain insights into the flow and spreading behavior. The computational domain and boundary conditions employed are consistent with those listed in Section 3.

The traverse speed and the corresponding lateral migration of the print nozzle along the X-axis govern the positional determination of the subsequent droplet. Observations indicate that the impinging droplets exhibit characteristic spreading and recoiling phases analogous to those observed in the case of a solitary uniform droplet. However, the offset distance between the two droplets notably influences the overall dynamics.

Figure 10 provides a comprehensive overview of the overall spreading and coalescence mechanisms within the two-dimensional computational domain, highlighting the intricate interplay of forces and the resultant flow patterns and morphological changes.

Although the droplet dynamics adhere to the aforementioned spreading and recoiling phases, the overall behavior exhibits notable asymmetry due to the addition of a secondary droplet. The subsequent impingement of the satellite droplet significantly influences the initial spreading of the sessile droplet. This collision introduces a higher degree of dynamism to the flow, with initial observations indicating air entrapment at the interface. However, as the flow progresses, the entrapped air dissipates, culminating in the formation of a crown which represents the apex of phase accumulation.

The crown formation predominantly occurs at the first droplet’s impingement site, indicating a localized accumulation of phases at the collision point. This phase is followed by the formation of an unequal rim at the conclusion of the spreading phase. The recoiling behavior is characterized by the inward movement of this unequal rim towards the center. Notably, the recoiling phase is prolonged compared to a solitary uniform droplet, attributable to the flow’s increased complexity and dynamic nature.

The system exhibits multiple oscillations during the subsequent spreading and recoiling phases, resulting in a longer final deposition time than a single droplet. The present study highlights the intricate spreading and recoiling phases and the asymmetrical flow behavior of coalescing droplets. These findings underscore the necessity for extending the analysis to a three-dimensional framework to fully capture the complexities of the phenomenon.

3.3.2. Prediction of the Flow Dynamics with Relevance to Droplet-Based Printing

The computational model is extended to a three-dimensional domain, maintaining the same boundary conditions as those utilized in the previous section. A more viscous fluid (similar to printing ink) is considered to study the flow behavior and predict the final print resolution, emulating the inkjet printing process. A droplet with a diameter of 3.6 mm is considered to impinge on a polycarbonate surface with an initial velocity of 0.68 m/s. The fluid properties include a viscosity of 4.98 × 10⁻⁵ Pa.s and a surface tension of 0.0798 N/m. The equilibrium contact angle between the printing ink and the substrate is set at 70° to mimic a hydrophilic surface.

The flow behavior is analyzed through the evolution of the droplets during coalescence, as depicted in Figure 11. The three-dimensional analysis reveals that the droplet follows a single uniform droplet’s spreading and recoiling behavior. However, the higher viscosity of the fluid resists the spreading and recoiling phases, underscoring the significant role of viscosity in droplet dynamics. The flow behavior is further examined using velocity streamlines, which elucidate the material flow patterns during droplet impact and coalescence.

As shown in Figure 11, the initial droplet begins to spread upon impact, and subsequently, the satellite droplet impinges on the sessile droplet, altering the overall dynamics. The coalescence process demonstrates the merging of the two droplets, with the interfaces disappearing as mixing occurs. The spreading phase is prominent, while the recoiling phase is less dominant due to the higher viscosity of the fluid. The final print resolution is determined by the coalescence and merging of the interfaces of these two droplets.

The velocity streamlines provide insights into the flow patterns of the print material throughout the spreading and recoiling phases. Initially, the velocity streamlines are directed outwards for the sessile droplet and downwards for the satellite droplet, as observed at 8 ms in Figure 11. Upon the satellite droplet’s impact on the sessile droplet, the flow pattern remains consistent instantaneously. However, as the sessile droplet continues to spread, the flow patterns for the satellite droplet exhibit both outward and downward directions. This behavior is followed by a strong mixing of the velocity streamlines at the interface where the two droplets merge. This information is crucial for determining the material flow in different regions, aiding in the prediction of the surface characteristics of the final printed part. The flow pattern reverses during the recoiling phase, although it is less dominant than the spreading phase. The final print resolution and surface finish result from this dynamic flow pattern of the printing material.

Furthermore, the effect of print speed on the print resolution is analyzed at various print speeds. The print speed directly influences the center-to-center distance between two droplets. Different print speeds ranging from 0.5 m/s to 5 m/s were modeled, representing varying center-to-center distances between the droplets, as illustrated in Figure 12. This analysis is crucial for predicting the final resolution and surface characteristics at different printing speeds.

Three distinct cases are examined, ranging from drop-on-drop deposition to higher print speeds. The comparison of the final resolution for these cases is presented in Figure 12. This comparative analysis provides valuable insights into how print speed affects the coalescence dynamics, droplet interaction, and, ultimately, the surface finish of the printed part.

This study aims to optimize the printing process to achieve desired resolutions and surface qualities by understanding the relationship between print speed and droplet spacing. The findings highlight the importance of precise control over print speed to ensure high-quality outcomes in droplet-based printing applications.

The print resolution increases with the increase in print speed, which has significant implications for manufacturing time and surface finish. Higher printing speeds enable faster manufacturing processes, thereby reducing production times. However, this acceleration often comes at the expense of print quality. Extreme print speeds can lead to undesirable outcomes such as puddle formation, poor surface finish, and wastage of print material.

Therefore, it is essential to strike a balance between printing speed and achievable print resolution. The parametric study conducted in this research helps identify the optimal printing speeds that yield the desired print resolutions while maintaining acceptable surface quality. By analyzing different print speeds and their effects on droplet coalescence and interaction, this study provides valuable insights into optimizing the printing process for various applications.

The findings underscore the importance of carefully selecting print speeds to ensure high-quality outcomes in droplet-based printing. This balance is crucial for achieving efficient manufacturing processes without compromising the integrity and precision of the printed structures.

4. Discussion

The numerical model extends the qualitative analysis to a three-dimensional study. This quantitative study agrees well with the experimental studies conducted by Kim and Chum [28] for hydrophilic surfaces and by Azimi et al. [29] for hydrophobic surfaces. Notably, the computational model captures the drop evolution with impingement on nanostructured hydrophobic surfaces, as reported by Azimi et al. [29]. Additionally, an in-depth analysis of dimensionless numbers (Re, We, and Oh) quantifies the role of different forces associated with printing. This non-dimensional number-based analysis sheds light on the mechanics of droplet dynamics and the printability of materials specific to additive manufacturing applications.

The numerical model and modeling framework agree with the experimental data by following the shape evolution of a droplet on a polycarbonate surface and capturing the mechanisms in droplet-based printing. The present work illustrated different stages of spreading, recoiling, lamella formation, and wetting for droplet dynamics in printing, as reported in other experimental and numerical studies.

A three-dimensional model provided an effective modeling framework and improved the quantitative analysis of droplet-based printing in additive manufacturing. The numerical model extended the qualitative analysis reported by Shah and Mohan [15] to a three-dimensional study, showing good agreement with the experimental studies. The computational model even captured the various stages of the interaction of droplets with the polycarbonate surface in agreement with the experimental study. The accuracy of the VOF method heavily relies on the resolution of the computational grid. When dealing with small droplets or complex geometries, capturing the precise shape and position of the droplet interface becomes challenging. Insufficient grid resolution may result in inaccuracies, especially near the contact line between the droplet and the solid substrate.

Furthermore, the study predicted the final printing resolution depending on the process parameters and material properties. In addition, a dimensionless analysis of the effect of characteristic numbers (Re, We, and Oh) revealed the role of different forces associated with printing. The inertial forces mainly govern the initial spreading phenomenon, while the recoiling is characterized by the interplay of inertia, viscous, and surface tension forces. This study gives an insight into the mechanics of droplet dynamics, in general, and the printability of materials specific to additive manufacturing applications.

The numerical method captures the flow and spreading behavior of two coalescing droplets and shows the asymmetric flow pattern with the crown formation and uneven rim around the center. The flow behavior and the velocity streamlines reveal the intermixing of the material during the printing process. Higher printing speeds lead to less manufacturing time, but with poor resolution and surface finish. The parametric studies are conducted to achieve the desired printing resolution with optimal printing speeds. This information is much needed and helps decipher the conditions for optimizing the printing process. Furthermore, the study helps predict the final printing resolution depending on the process parameters and material properties.

Future studies could consider understanding the spreading behavior with multiple droplets, incorporating the material phase change and solidification, and the surface variation effects when the materials are printed on the subsequent layers. The inclusion of the phase-change model into the present work can extend this study to analyze the photopolymerization of polymers in 3D printing and its effect on print surfaces and overall print quality. The effect of temperatures can be studied further, and their implications on the final printing quality can be examined. Many industrial materials like polymers, bio-inks, and other commercial printing inks can also be analyzed with parametric studies for better print quality and fidelity.

The current modeling framework with complex computational models, phenomenological models for surface characteristics effects, and higher computing times is required for such analyses.

5. Conclusions

The present investigations on multiphase modeling of flow and spreading behavior in droplet-based printing yield several key insights. The numerical model and modeling framework demonstrate strong agreement with experimental data, accurately capturing the shape evolution of droplets on both hydrophilic and hydrophobic surfaces. The model effectively captures the mechanisms involved in droplet-based printing, including different stages of spreading, recoiling, lamella formation, and wetting, as reported in other experimental and numerical studies.

The extension to a three-dimensional model provides an effective framework for improving the quantitative analysis of droplet-based printing in additive manufacturing. The computational model successfully captures the interaction of droplets with nanostructured hydrophobic surfaces, providing essential information for optimizing the 3D printing process.

The study predicts the final printing resolution based on process parameters and material properties. A dimensionless analysis of characteristic numbers (Reynolds number, Weber number, and Ohnesorge number) reveals the roles of different forces associated with printing. Inertial forces primarily govern the initial spreading phenomenon, while the interplay of inertial, viscous, and surface tension forces characterizes the recoiling phase. This analysis provides valuable insights into the mechanics of droplet dynamics and the printability of materials specific to additive manufacturing applications.

The study of droplet coalescence in the context of 3D printing reveals asymmetric flow behavior around the center. Flow and velocity streamlines indicate the mixing and dynamics of interface merging and droplet coalescence. Higher printing speeds result in shorter manufacturing times but poorer print quality. The effect of print speeds on final resolution and surface finish is analyzed, emphasizing the need to optimize the printing process for different speeds.

While our current study provides valuable insights into droplet coalescence, a comprehensive quantitative analysis of the surface tension-driven coalescence phenomenon is necessary. The studies can be extended to perform process optimization by conducting a sensitivity analysis for various process parameters, including traverse speed, droplet size, and impinging velocities. This future scope will enhance our understanding of droplet coalescence and its implications for print quality and resolution. Future studies should investigate the effects of temperature on final printing quality. Incorporating a phase-change model into the present work could extend the analysis of the photopolymerization of polymers in 3D printing and its impact on print surfaces and overall print quality. Many industrial materials, such as polymers, bio-inks, and other commercial printing inks, can be analyzed through parametric studies to achieve better print quality and fidelity. This study can also be extended to bio-ink materials and biological cells to predict and characterize printability and print quality in bioprinting.

Future work should incorporate non-Newtonian models to simulate compound droplets and their applications in droplet-based printing. It could also focus on understanding the spreading behavior with multiple droplets, incorporating material phase change and solidification, and examining surface variation effects when materials are printed on subsequent layers.

Author Contributions

Conceptualization, R.V.M.; Methodology, R.S. and R.V.M.; Software, R.S.; Formal analysis, R.S.; Investigation, R.S.; Resources, R.V.M.; Writing—original draft, R.S.; Writing—review & editing, R.V.M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request. The data are not publicly available due to Researcher Intent.

Acknowledgments

The authors acknowledge the support provided for this work at the Joint School of Nanoscience and Nanoengineering, a member of the National Nanotechnology Coordinated Infrastructure (NNCI), which is supported by the National Science Foundation. We thank the University of North Carolina at Chapel Hill and the Research Computing group for providing computational resources and support that have contributed to these research results. Computing resource support from the University of North Carolina at Charlotte is also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

VOF	Volume of Fluid
CFD	Computational Fluid Dynamics
CFL	Courant–Friedrichs–Lewy
PISO	Pressure Implicit with Splitting of Operators
SIMPLE	Semi-Implicit Method for Pressure-Linked Equations
MULES	Multidimensional Universal Limiter for Explicit Solution
Re	Reynolds number
We	Weber number
Oh	Ohnesorge number

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1 Schematic showing the computational domain and associated boundary conditions for a droplet impinging on a solid surface.

Figure 2 Comparison of experimental and numerical results for a water droplet impinging on a polycarbonate surface. Experimental results are shown in black-and-white images [28], while the numerical results are blue, representing the evolution of the water droplet shape. Experimental results are reproduced with permission from [28]. Copyright 2001 AIP Publishing LLC.

Figure 3 Comparison of shape evolution of the droplet: The surface of the droplet captured in the numerical study is compared with experimental results. Experimental results are shown in gray area images [28], while the numerical results are blue contours representing the water droplet shape evolution. Experimental results are reproduced with permission from [28]. Copyright 2001 AIP Publishing LLC.

Figure 4 Comparison of experimental studies [25] and numerical results for a water droplet impinging on a polycarbonate surface. Quantities with superscript * denote the normalized quantities, as defined earlier.

Figure 5 Comparison of experiments and the computational results. (Top) High-speed camera images of a water droplet impinging on the nanograss-covered micro-posts sputtered with ceria (impact velocity ~1.6 m s⁻¹); scale bar, 2.5 mm [29]. (Bottom) Shape evolution of the water droplet emulating the experimental study using a multiphase CFD solver. Experimental results are reproduced with permission from Azimi et al. [29]. Copyright 2013 Nature Publishing Group.

Figure 6 Comparison of experimental data reported by Azimi et al. [29] for a water droplet impinging on a hydrophobic surface and the computational results using a multiphase CFD solver for the same case study.

Figure 7 Droplet evolution (base diameter) comparison for different Reynolds numbers. Quantities with superscript * denote the normalized quantities, as defined earlier.

Figure 8 Droplet evolution (base diameter) comparison for different (a) Weber and (b) Ohnesorge numbers. Quantities with superscript * denote the normalized quantities, as defined earlier.

Figure 9 Sketching of the in-print conditions. (a) The sessile droplet position with print speed v_t. (b) Center-to-center distance for two consecutive impinged droplets.

Figure 10 Selected frames from the simulation depicting salient features for the coalescence of two droplets.

Figure 11 The flow behavior of two droplets during impact and coalescence, along with corresponding velocity streamlines.

Figure 12 Comparison of droplet coalescence and shape evolution for different print speeds (0.5 m/s, 3 m/s, and 5 m/s from left to right), snapshots taken at 2.4 ms, 8 ms, 15.2 ms, 18.2 ms, and 40 ms (from top to bottom).

Appendix A. Grid Independence Study Results

Convergence Analysis: We conducted a grid independence study using multiple grid resolutions to ensure that the results are independent of the mesh size. The study involved testing grid sizes of 50 × 50 × 35, 80 × 80 × 45, 100 × 100 × 60, 125 × 125 × 75, and 150 × 150 × 90 cells. The results showed that the key parameters, such as droplet spread diameter and velocity profiles, converged with increasing grid resolution. Specifically, the differences in these parameters between the 125 × 125 × 75 and 150 × 150 × 90 grids were less than 1%, indicating that further refinement does not significantly affect the simulation results.

Computational Time Analysis: The parallel computations were executed on the Longleaf high-performance computing (HPC) cluster at the University of North Carolina at Chapel Hill, NC, USA. Each simulation utilizes 8 nodes with 36 processing cores to efficiently distribute the computational workload. The above-mentioned grid sizes were employed to analyze the computational time, providing insights into the scalability and performance of the simulations across different mesh resolutions. The 100 × 100 × 60 grid required approximately two hours of computational time, while the 125 × 125 × 75 grid required 3.5 hours, and the 150 × 150 × 90 grid required 5 hours. The 125 × 125 × 75 grid was chosen as the optimal resolution, balancing accuracy and computational efficiency.

Table A1

Grid Independence Study.

Mesh Resolution	Cells	Computational Time (Approx Hours)	Key Parameter (Droplet Spread Diameter in mm)
Coarse	50 × 50 × 35	1.0	6.532
Medium	80 × 80 × 45	1.5	6.010
Fine	100 × 100 × 60	2.0	5.791
Finer	125 × 125 × 75	3.5	5.767
Finest	150 × 150 × 90	5.0	5.720

References

1. Yarin, A.L. DROP IMPACT DYNAMICS: Splashing, Spreading, Receding, Bouncing. Annu. Rev. Fluid Mech.; 2006; 38, pp. 159-192. [DOI: https://dx.doi.org/10.1146/annurev.fluid.38.050304.092144]

2. Elliot, T.A.; Ford, D.M. Dynamic contact angles. J. Chem. Soc.; 1972; 68, pp. 1814-1823. [DOI: https://dx.doi.org/10.1016/0021-9797(67)90183-X]

3. Fukai, J.; Shiiba, Y.; Yamamoto, T.; Miyatake, O.; Poulikakos, D.; Megaridis, C.M.; Zhao, Z. Wetting effects on the spreading of a liquid droplet colliding with a flat surface: Experiment and modeling. Phys. Fluids; 1995; 7, pp. 236-247. [DOI: https://dx.doi.org/10.1063/1.868622]

4. Asai, A.; Shioya, M.; Hirasawa, S.; Okazaki, T. Impact of an Ink Drop on Paper. J. Imaging Sci. Technol.; 1993; 37, 205.

5. Bayer, I.S.; Megaridis, C.M. Contact angle dynamics in droplets impacting on flat surfaces with different wetting characteristics. J. Fluid Mech.; 2006; 558, pp. 415-449. [DOI: https://dx.doi.org/10.1017/S0022112006000231]

6. Scheller, B.; Bousfield, W. Newtonian Drop Impact with a Solid Surface. AIChE J.; 1995; 41, pp. 1357-1367. [DOI: https://dx.doi.org/10.1002/aic.690410602]

7. Roisman, I. Inertia Dominated Drop Collisions. II. An Analytical Solution of the Navier–Stokes Equations for a Spreading Viscous Film. Phys. Fluids; 2009; 21, 052104. [DOI: https://dx.doi.org/10.1063/1.3129283]

8. Toivakka, M. Numerical Investigation of Droplet Impact Spreading in Spray Coating of Paper. Proceedings of the 2003 TAPPI 8th Advanced Coating Fundamentals Symposium; Chicago, IL, USA, 8–10 May 2003.

9. Pasandideh-Fard, M.; Qiao, Y.M.; Chandra, S.; Mostaghimi, J. Capillary Effects During Droplet Impact on a Solid Surface. Phys. Fluids; 1996; 8, pp. 650-659. [DOI: https://dx.doi.org/10.1063/1.868850]

10. Feng, J.Q. A Computational Study of High-Speed Microdroplet Impact onto a Smooth Solid Surface. J. Appl. Fluid Mech.; 2017; 10, pp. 243-256. [DOI: https://dx.doi.org/10.18869/acadpub.jafm.73.238.26440]

11. Adeniyi, A.A.; Morvan, H.P.; Simmons, K.A. A coupled Euler-Lagrange CFD modeling of droplets-to-film. Aeronaut. J.; 2017; 121, pp. 1897-1918. [DOI: https://dx.doi.org/10.1017/aer.2017.107]

12. Yokoi, K.; Vadillo, D.; Hinch, J.; Hutchings, I. Numerical Studies of the Influence of the Dynamic Contact Angle on a Droplet Impacting on a Dry Surface. Phys. Fluids; 2009; 21, 072102. [DOI: https://dx.doi.org/10.1063/1.3158468]

13. Gujjula, S. Numerical Investigation of Droplet Spread: Effect of Contact Angle Models. Master’s Thesis; University of Cincinnati: Cincinnati, OH, USA, 2015.

14. Wang, L.S.; Krull, B.; Lucas, D.; Meller, R.; Schlegel, F.; Tekavcic, M.; Xu, J.Y. Simulation of droplet entrainment in annular flow with a morphology adaptive multifield two-fluid model. Phys. Fluids; 2023; 35, 103312. [DOI: https://dx.doi.org/10.1063/5.0169288]

15. Shah, R.G.; Mohan, R.V. Multiphase Modeling of Droplet-Based 3D Printing: Predicting Printability, Resolution and Shape Fidelity in Additive Manufacturing Processes. Proceedings of the ASME 2023 International Mechanical Engineering Congress and Exposition; New Orleans, LA, USA, 29 October–2 November 2023.

16. Vontas, K. Enhanced VOF-based CFD simulations of droplet impact. Bachelor’s Thesis; Piraeus University of Applied Sciences, University of Brighton: Brighton, UK, 2016.

17. Teng, T.; Joralmon, D.; Anyigbo, T.; Li, X. Acoustic Levitation-Assisted Contactless Printing of Microdroplets for Biomedical Applications. J. Manuf. Sci. Eng.; 2024; 146, 011002. [DOI: https://dx.doi.org/10.1115/1.4062971]

18. Wu, D.; Xu, C. Predictive Modeling of Droplet Formation Processes in Inkjet-Based Bioprinting. Trans. ASME.; 2018; 140, 101007. [DOI: https://dx.doi.org/10.1115/1.4040619]

19. Aladese, A.; Jeong, H. Recent Developments in 3D Printing of Droplet-Based Microfluidics. BioChip J.; 2021; 15, pp. 313-333. [DOI: https://dx.doi.org/10.1007/s13206-021-00032-1]

20. Britel, A.; Tomagra, G.; Apra, P.; Varzi, V.; Sturari, S.; Amine, N.; Olivero, P.; Picollo, F. 3D printing in microfluidics: Experimental optimization of droplet size and generation time through flow focusing, phase, and geometry variation. RSC Adv.; 2024; 14, pp. 7770-7778. [DOI: https://dx.doi.org/10.1039/D4RA00752B]

21. Jiang, J.; Chen, X.; Mei, Z.; Chen, H.; Chen, J.; Wang, X.; Li, S.; Zhang, R.; Zheng, G.; Li, W. Review of droplet printing technologies for flexible electronic devices: Materials, control, and applications. Micromachines; 2024; 15, 333. [DOI: https://dx.doi.org/10.3390/mi15030333]

22. Abdolmaleki, H.; Kidmose, P.; Agarwala, S. Droplet-based techniques for printing of functional inks for flexible physical sensors. Adv. Mater.; 2021; 33, 2006792. [DOI: https://dx.doi.org/10.1002/adma.202006792]

23. Zhang, H.; Cui, L.; Lee, P.-H.; Kim, Y.; Moon, S.K.; Choi, J.P. Data-Driven autonomous printing process optimization and real-time abnormality identification in aerosol jet-deposited droplet morphology. Virtual Phys. Prototyp.; 2024; 19, e2429530. [DOI: https://dx.doi.org/10.1080/17452759.2024.2429530]

24. Tofan, T.; Borodinas, S.; Kačianauskas, R.; Jasevičius, R. Modeling 3D droplet movement using a drop-on-demand inkjet printhead model. Processes; 2022; 10, 1467. [DOI: https://dx.doi.org/10.3390/pr10081467]

25. Jia, Y.; Cheng, X.; Song, W.; Zhou, Y.; Zhao, H. A droplet state prediction method for inkjet printing under small sample conditions based on the two-stage TrAdaBoost.R2 algorithm. AIP Adv.; 2025; 15, 015120. [DOI: https://dx.doi.org/10.1063/5.0246942]

26. The OpenFOAM Foundation. OpenFOAM v2021 User Guide. Available online: https://cfd.direct/openfoam/user-guide (accessed on 10 January 2023).

27. Fujimoto, H.; Shiotani, Y.; Tong, A.Y.; Hama, T.; Takuda, H. Three-Dimensional Numerical Analysis of the Deformation Behavior of Droplets Impinging onto a Solid Substrate. Int. J. Multiph. Flow; 2007; 33, pp. 317-332. [DOI: https://dx.doi.org/10.1016/j.ijmultiphaseflow.2006.06.013]

28. Kim, H.-Y.; Chun, J.-H. The recoiling of liquid droplets upon collision with solid surfaces. Phys. Fluids; 2001; 13, pp. 643-659. [DOI: https://dx.doi.org/10.1063/1.1344183]

29. Azimi, G.; Dhiman, R.; Kwon, H.-M.; Paxson, A.T.; Varanasi, K.K. Hydrophobicity of rare-earth oxide ceramics. Nat. Mater.; 2013; 12, pp. 315-320. [DOI: https://dx.doi.org/10.1038/nmat3545] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23333998]

30. Hsu, A.S.; Roy, A.; Leal, L.G. Drop-size effects on coalescence of two equal-sized drops in a head-on collision. J. Rheol.; 2008; 52, pp. 1291-1310. [DOI: https://dx.doi.org/10.1122/1.2980013]

31. Acevedo-Malavé, A.; García-Sucre, M. Coalescence collision of liquid drops I: Off-center collisions of equal-size drops. AIP Adv.; 2011; 1, 032117. [DOI: https://dx.doi.org/10.1063/1.3624553]

32. Tang, C.; Zhang, P.; Law, C.K. Bouncing, coalescence, and separation in head-on collision of unequal-size droplets. Phys. Fluids; 2012; 24, 022101. [DOI: https://dx.doi.org/10.1063/1.3679165]

33. Josserand, C.; Thoroddsen, S.T. Drop impact on a solid surface. Annu. Rev. Fluid Mech.; 2016; 48, pp. 365-391. [DOI: https://dx.doi.org/10.1146/annurev-fluid-122414-034401]

34. Antonini, C.; Amirfazli, A.; Marengo, M. Drop impact and wettability: From hydrophilic to superhydrophobic surfaces. Phys. Fluids; 2012; 24, 102104. [DOI: https://dx.doi.org/10.1063/1.4757122]

35. Rein, M. The transitional regime between coalescing and splashing drops. J. Fluid Mech.; 1996; 306, pp. 145-165. [DOI: https://dx.doi.org/10.1017/S0022112096001267]

36. Thoroddsen, S.T.; Etoh, T.G.; Takehara, K. Air entrapment under an impacting drop. J. Fluid Mech.; 2003; 478, pp. 125-134. [DOI: https://dx.doi.org/10.1017/S0022112002003427]

37. Thoroddsen, S.T.; Takehara, K. The coalescence cascade of a drop. Phys. Fluids; 2000; 12, pp. 1265-1267. [DOI: https://dx.doi.org/10.1063/1.870380]

38. Mohamed-Kassim, Z.; Longmire, E.K. Drop coalescence through a liquid/liquid interface. Phys. Fluids; 2004; 16, pp. 2170-2181. [DOI: https://dx.doi.org/10.1063/1.1735686]

39. Chen, Y.; Shen, C.; Peterson, G.P. Hydrodynamics and morphologies of droplet coalescence. Ind. Eng. Chem. Res.; 2015; 54, pp. 9257-9262. [DOI: https://dx.doi.org/10.1021/acs.iecr.5b01459]

40. Shen, C.; Chen, Y.; Yu, C.; Liu, X. Numerical study on the liquid-liquid interface evolution during droplet coalescence. Microgravity Sci. Technol.; 2020; 32, pp. 737-748. [DOI: https://dx.doi.org/10.1007/s12217-020-09805-z]

41. Zimmermann, P.; Mawbey, A.; Zeiner, T. Calculation of droplet coalescence in binary liquid-liquid systems: An incompressible Cahn–Hilliard/Navier–Stokes approach using the non-random two-liquid model. J. Chem. Eng.; 2020; 65, pp. 1083-1094. [DOI: https://dx.doi.org/10.1021/acs.jced.9b00493]

42. Baroudi, L.; Kawaji, M.; Lee, T. Effects of initial conditions on the simulation of inertial coalescence of two drops. Comput. Math. Appl.; 2014; 67, pp. 282-289. [DOI: https://dx.doi.org/10.1016/j.camwa.2013.05.002]

43. Lim, S.J.; Choi, M.C.; Weon, B.M.; Gim, B. Lattice Boltzmann simulations for water coalescence. Appl. Phys. Lett.; 2017; 111, 101602. [DOI: https://dx.doi.org/10.1063/1.4985264]

44. Gross, M.; Steinbach, I.; Raabe, D.; Varnik, F. Viscous coalescence of droplets: A lattice Boltzmann study. Phys. Fluids.; 2013; 25, 052101. [DOI: https://dx.doi.org/10.1063/1.4803178]

45. Perumanath, S.; Borg, M.K.; Chubynsky, M.V.; Sprittles, J.E.; Reese, J.M. Droplet coalescence is initiated by thermal motion. Phys. Rev. Lett.; 2019; 122, 104501. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.104501]