Numerical simulations of wall contact angle effects on droplet size during step emulsification

Abstract

Terrace-based microfluidic devices are currently used to prepare highly monodisperse micro-droplets. Droplets are generated due to the spontaneous pressure drop induced by the Laplace pressure, and so the flow rate of a dispersed phase has little effect on droplet size. As a result, control over the droplet is limited once a step emulsification device has been fabricated. In this work, a terrace model was established to study the effect of the wall contact angle on droplet size based on computational fluid dynamics simulations. The results for contact angles from 140° to 180° show that a lower contact angle induces wall-wetting, increasing the droplet size. The Laplace pressure equations for droplet generation were determined based on combining pressure change curves with theoretical analyses, to provide a theoretical basis for controlling and handling droplets generated through step emulsification.

---

Introduction

The development of low-cost methods for the high-throughput generation of controllable micron/sub-micron droplets is an important topic in microfluidic technology, as well as an urgent problem with regard to practical industrial production. The breakup mechanisms of droplets can be divided into two types: shear induced or spontaneous formation.^1–3 The widely used cross-junction/T-junction,⁴ flow-focusing⁵ and co-flow⁶ techniques are all shear-induced. In such cases, a flowing continuous phase is required to detach the droplets, necessitating accurate control of the flow rates of both dispersed phase and continuous phase. In contrast, the step emulsification process developed by Kawakatsu⁷ in 1997 relies on the interfacial tension between two phases, such that flow of the continuous phase is not necessary for the generation of droplets.^8,9 In this method, droplets break up due to a Laplace pressure difference resulting from a dimension change outside the terrace.¹⁰ This fracture mechanism reduces the effect of flow rate on the droplet size, making it more likely that micro-droplets with uniform sizes will be obtained.¹¹ Furthermore, terrace-based emulsification allows efficient use of the area of a microfluidic device to achieve the dense packing of nozzles and thus a high throughput of droplets.^12–14

Li et al.¹⁵ analysed the quasi-steady interface shape of the fluid before conversion to the step emulsification oscillation mode, and determined the relationship between the capillary number and the channel section ratio in the critical transition state. Dutka et al.¹⁶ added a contraction neck and by-pass at the end of the inlet channel to the terrace device. This modification allows the continuous phase flow to return through the bypass to cut off droplets, thus reducing the effect of the dispersed phase rate on the droplet volume.

Because of the theoretical complexity of this field of study and the limitations of experimental work, many researchers have started to use computational fluid dynamics (CFD) to study step emulsification. In 2004, Kobayashi et al.¹⁷ first used the finite volume method to simulate droplet generation in a symmetric straight flow-through microchannel, and obtained important insights into the movement of the water–oil interface and flow pressure variations associated with different channel sizes. In 2011, the same team studied the effects of the micro-channel size and the dispersed phase flow rate on the formation of droplets and also analysed the regulation of interfacial tension, flow rate and pressure.¹⁸ Chakraborty et al.¹⁹ designed an axisymmetric step emulsification apparatus that simulated the transition mechanism in two dimensions, and demonstrated that the droplet size is independent of the capillary number and almost independent of viscosity. However, to date, there have been fewer studies regarding the effect of the wall contact angle on droplet size during step emulsification.

The size parameter cannot be changed once the terrace device has been fabricated by using either soft-lithography²⁰ or 3D printing.²¹ Furthermore, the dispersed phase flow has little effect on the droplet size.²² Thus, the most convenient means of adjusting droplet size is to vary the contact angle between the dispersed phase and the channel wall. For this reason, it is important to study the effect of the contact angle during step emulsification, as this is an important aspect in the design of micro-droplet preparation equipment, and also assists in the choice of suitable raw materials for industrial processes. In experimental work, the interfacial tension and contact angle of the fluid are often changed by adding various surfactants. Although there have been in-depth studies of the effects of surface tension, the role of surfactants in droplet formation is complex and ambiguous.^23,24 Surfactants not only change the interfacial tension but can also modify other physical properties of the fluid. As an example, the addition of varying amounts of butanol can be used to adjust the surface tension, which in turn affects the viscosity of the continuous phase. Such changes will have an effect on the study of single variables. Moreover, the addition of surfactants may be strictly limited sometimes, especially in certain medical²⁵ or biological^26,27 applications.

The use of numerical calculations can mitigate the limitations imposed by experimental work. Thus, in the present work, a basic model of step emulsification was established and the validity of the model was verified by comparing experimental data with simulations. Subsequently, the effect of the wall contact angle on droplet size was assessed while fixing all other parameters. This work also investigated the effects of the terrace width on droplets generated in a step emulsification apparatus based on simulations.

Simulations and experimental work

In this study, water (density: 998.2 kg m⁻³ and viscosity: 1.003 × 10⁻³ Pa s) was used as the dispersed phase and silicone oil (density: 930 kg m⁻³ and viscosity: 9.3 × 10⁻³ Pa s) was the continuous phase. Because the Reynolds number in this scenario is small (Re < 23), the flow of the two phases can be considered laminar.

Computational domain and boundary conditions

The models used in the simulations were defined using the Design Modeller software package. To minimise computational time, one-fourth of the total domain was employed as the computational domain, with symmetry in the x = 0 and y = 0 planes, as shown in Fig. 1. Note that the channel width in this model is equal to the height, although this is not stated in the diagram. The ICEM software package was used to generate grids. Due to the abrupt change in volume from quasi-two-dimensional to three-dimensional, a grid transition area was included at the terrace outlet to reduce the possibility of divergence of the results, as well as to minimise the calculation burden. Table 1 shows the step parameters (10 μm in height and 50 μm in length with various widths) that were used for the simulations in this work, as well as the grid partitioning. The minimum cell size in the model was 1 μm.

Fig. 1 A three-dimensional schematic diagram of the step model used in the current simulation.

Table 1 Parameters of the step emulsification models used in this study

Groups of different step width	Channel length l1 (μm)	Channel height h (μm)	Terrace length l2 (μm)	Terrace width, W (μm)	Reservoir size, W × L × H (μm)	Number of cells
a Group 4 is the default size for contact angle simulation, as shown in Fig. 4.
1	30	10	50	70	70 × 100 × 110	120 479
2				90	90 × 100 × 110	153 009
3				110	110 × 100 × 110	185 539
4^a				150	150 × 100 × 110	228 664
5				170	170 × 100 × 110	233 559
6				190	190 × 100 × 110	288 964
7				210	210 × 100 × 110	318 644
8				230	230 × 100 × 110	338 925

---

Fig. 2 presents the grid partition and boundary settings used in the model for step emulsification. Here, the inlet is the region where the dispersed phase fluid enters, the wall is the boundary between the dispersed phase fluid and the terrace unit and the boundary at which the continuous phase is situated is defined as the outlet. Moreover, two symmetrical faces were employed to achieve a symmetrical model and to display the simulation results. The dispersed phase velocity at the channel inlet was set to 0.2 m s⁻¹ and the reference outlet had a pressure boundary of 0 Pa. A PISO scheme was employed for pressure–velocity coupling and a PRESTO! scheme was used to determine the pressure difference, together with a second-order upwind, which was used to compute the momentum spread. The initial time step was 10⁻⁶ s together with a convergence standard of 1 × 10⁻⁵.

Fig. 2 The grid partition and boundary settings in the model. (I) Inlet, the entrance for the dispersed phase. (II) Wall, the interface between the fluid and the terrace device. (III) Symmetry. (IV) Outlet, the reservoir filled with continuous phase.

Comparison of experimental and simulation results

A volume of fluid (VOF) model in which a volume fraction distinguishes the dispersed and continuous phases was used to simulate the droplet formation process. The step device was composed of two cover glasses connected using double-sided tape. A channel (80 μm height × 80 μm width × 1000 μm length) was connected to a terrace (80 × 230 × 15 000 μm) that ended in a much deeper reservoir filled with the continuous phase. The two cover glasses were treated with the hydrophobic agent Aquapel (PPG Industries) to ensure that water-in-oil emulsions were obtained. The terrace in the simulations had dimensions of 80 × 230 × 800 μm. Thus, the width was decreased slightly so as to reduce the computational load while not affecting the simulation results. Under certain conditions, the terrace width had little effect on the droplet size, as discussed below.

The simulation and experimental results are compared in Fig. 3, which demonstrates good agreement between the two sets of results. The simulation diagram shows that the formation of emulsion droplets at the terrace primarily involves three processes: disk expansion in the terrace, expansion and formation of the neck outside the terrace, and detachment. During disk expansion, the dispersed phase flows through the channel inlet and gradually expand in the terrace until reaching the terrace outlet, while part of the continuous phase inside the terrace is forced to flow out of the terrace. As the dispersed phase continues to expand outside the terrace, some of the continuous phase flows rapidly back into the terrace, squeezing the water phase to induce the necking phenomenon, which is defined as the necking stage. The detachment process is eventually triggered by the intensification of necking. After drop generation, the dispersed phase return to the terrace and a new cycle of droplet generation begins. The pressure changes and mechanisms associated with breakup at various stages are discussed below. The experimental results also indicate that step emulsification can generate droplets without additional actuation of the continuous phase, which is analogous to the CFD simulation results.

Fig. 3 A comparison of CFD simulations with experimental results. (I) Beginning of the disk expansion in the terrace. (II) The tip of the dispersed phase reaches the terrace outlet. (III) Necking occurs. (IV) Droplet detachment.

Results and discussion

Simulation results for different wall contact angles

The default model used to simulate the effect of the contact angle had a 30 μm channel length, l₁, 10 μm channel height, h, 50 μm step length, l₂, 150 μm terrace width, W, and a reservoir size of W × L × H = 150 × 100 × 110 μm. These dimensions are equivalent to those of entry 4 in Table 1. The sizes of the resulting droplets were determined using Nano Measurer. The diameter of the first droplet is termed d₀, because this droplet will coalesce with the second one in the simulated situation. Fig. 4 shows the results of simulations for different contact angles. In an ideal scenario, the surface contact angle between the dispersed phase and the terrace is 180° such that there is no wetting of the wall by the fluid, which is advantageous to droplet formation. As the angle is gradually reduced, the dispersed phase wets the wall and the Laplace pressure difference is decreased, such that droplet fracture becomes less likely. In addition, the droplets readily adhere to the terrace entrance after breakup and rapidly coalesce.

Fig. 4 Effects of the static contact angle on droplet diameter and diagrams showing terrace wetting: (I) initialization conditions for the simulation, (II) θ ≤ 135, (III) θ = 140°, (IV) θ = 145°, (V) θ = 150°, (VI) θ = 160°, (VII) θ = 180°.

It is obvious that increasing the wall contact angle reduces the droplet size. There is also a noticeable wetting phenomenon when the contact angle is less than 150°. Under these conditions, the dispersed phase flows out of the terrace boundary from the side at wall contact angles of less than 140°. To eliminate the effect of the terrace width on the droplet simulation results, eight groups of terrace widths were applied in conjunction with the same simulation process, as shown in Fig. 5. When the contact angle is small, the droplet diameter evidently decreases with increases in terrace width, whereas this value is almost independent of the terrace width at contact angles greater than 150°. The decreasing trend of the droplet diameters in response to changes in the contact angle also slows.

Fig. 5 Several groups of variation in droplet diameter with terrace widths at different contact angles.

Pressure analysis along a symmetric axis

To further study the effect of contact angle on the formation of droplets by step emulsification, we selected representative data at angles of 140° and 180° and monitored pressures during droplet generation. Fig. 6 plots the internal pressure of the dispersed phase along the symmetric axis (z axis). The contrast between the two figures leads to the conclusion that the internal pressure of the dispersed phase below 140° is lower and the droplet expansion stage represents a longer time span. The amount of dispersed phase flowing into the droplet prior to detachment is thus greater, leading to larger droplets. During droplet formation, the internal pressure of the dispersed phase is approximately equal to the Laplace pressure difference between water and oil in association with the step emulsification mechanism. The four timelines here correspond to the liquid droplet emulsification shown in Fig. 3, and their Laplace pressure expressions can be described as follows: curve I represents the dispersed phase expansion stage in the terrace, while curve II indicates that the dispersed phase flows to the terrace exit. According to the Young–Laplace equation, the pressures in a circular channel in this stage are:andwhere ΔP_c is the pressure difference inside the channel, ΔP_t is the pressure difference inside the terrace, σ is the interfacial tension between the water and oil phases, R_c is the circular channel radius (herein, the square channel radius is approximated by h/2), r_t is the radius of the liquid droplet extending into the terrace, R_t is the radius of the droplet disk in the terrace and θ is the wall contact angle. The specific parameters are shown in Fig. 7, where the thick dashed line indicates the terrace exit.

Fig. 6 The pressure along the symmetric axis indicated by the dashed line. (a) Pressure at a contact angle of 140° at (I) 0.08 ms, (II) 2.16 ms, (III) 3.28 ms, and (IV) 4.0 ms. (b) Pressure at a contact angle of 180° at (I) 0.1 ms, (II) 1.25 ms, (III) 1.85 ms, and (IV) 2.4 ms.

Fig. 7 Parameters used in CFD droplet simulation for step emulsification. (I) and (II) are the view from the x–z plane, and (III) is the view from the y–z plane.

The dispersed phase begins to expand in three dimensions after flowing out of the terraces and the Laplace pressure in the to-be-formed droplet is:

In this case, the Laplace pressure will decrease with increases in R_d, while the input pressure is constant. The difference gradient drives the dispersed phase inside the terrace to rapidly flow into the to-be-formed droplet to induce a neck, which is indicated by curve III. Due to the formation of the neck, P_t gradually decreases to a minimum value. Because the radius, R_t, can expand in the terrace to an essentially infinite value, the minimum pressure, P_tmin, primarily depends on r_t, which is defined as:

After reaching a minimum, P_t does not decrease with increases in the droplet radius, R_d. Thus, a large pressure gradient forms between P_tmin and P_d, which triggers a significant dispersed phase flux from the terrace. Simultaneously, the liquid flow into the terrace is determined by the difference, defined as ΔP, between the applied pressure, P_c, and P_t, which has a constant value. At this point, the liquid flow into the terrace no longer changes. As a result, droplet breakup will occur as soon as the liquid flux flow outside the terrace is greater than the flux inside the terrace, which is described by curve IV. The pressure at the platform is the Laplace pressure difference for the resulting liquid droplet.

Pressure analysis at the breaking point

Here, we selected a single point to illustrate the pressure change during the process of droplet formation. To ensure that the two groups of droplet fracture positions were consistent, the same point was selected for comparison. Fig. 8 presents the pressure variation curves at point (1,0,117) within the two droplet formation periods. The first droplet remains at the terrace entrance and coalesces with the second one, such that T₂ is much shorter than T₁. Here, T₂ is simply used to illustrate the periodicity of droplet generation. The two pressure increase points in the dotted circle correspond to droplet fracture. The lowest point in Fig. 8(a) shows the moment of droplet fracture. In this moment, the dispersed phase, which is flowing from the inlet, returns rapidly, leading to a significant negative pressure. However, some simulation cases did not write the data of instantaneous break, such that there was no negative pressure. The comparison shows that, at a contact angle of 140°, the droplet stays longer inside the terrace and the P_tmin required for necking is less, such that a longer time is required to achieve this value. This occurs primarily because r_t increases after the contact angle becomes smaller, decreasing the Laplace pressure difference inside the terrace. Under these conditions, a lower P_d is required to induce droplet rupture, allowing more time for the droplet to grow.

Fig. 8 Variations in pressure at one point. (a) Pressure at a contact angle of 140° at T₁ = 3.38 ms. (b) Pressure at a contact angle of 180° at T₂ = 1.86 ms. The point (1,0,117) was chosen to obtain a consistent droplet fracture position.

Summary

A computational fluid dynamics model intended to simulate step emulsification was established. A comparison of simulation results with experimental data confirmed the validity of the model to some extent. These simulations reproduced the movement of the water–oil interface during droplet formation as well as terrace wetting, due to reductions in the wall contact angle that increased the droplet size. Variations in pressure along the central axis and at a fixed point were monitored and the droplet fracture mechanism at different contact angles was analysed based on hydromechanics theory. The model established in this study can also be used to simulate the effects of other factors on step emulsification, such as terrace parameters, dispersed phase flow and fluid physical parameters. Furthermore, this simulation method provides a theoretical basis for terrace design and solution preparation to obtain specific droplet sizes.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 21676244) and the National Basic Research Program of China (Grant No. 2015CB057301). We thank Michael D. Judge, MSc, from Liwen Bianji, Edanz Editing China (www.liwenbianji.cn/ac), for editing the English text of a draft of this manuscript.

References

1. T. Nisisako, T. Torii and T. Higuchi, Lab Chip, 2002, 2, 24–26.
2. D. R. Link, S. L. Anna, D. A. Weitz and H. A. Stone, Phys. Rev. Lett., 2004, 92, 54503.
3. I. Kobayashi, M. Nakajima, K. Chun, Y. Kikuchi and H. Fujita, AIChE J., 2002, 48, 1639–1644.
4. T. Thorsen, R. W. Roberts, F. H. Arnold and S. R. Quake, Phys. Rev. Lett., 2001, 86, 4163–4166.
5. S. L. Anna, N. Bontoux and H. A. Stone, Appl. Phys. Lett., 2003, 82, 364–366.
6. A. Shams Khorrami and P. Rezai, Biomicrofluidics, 2018, 12, 34113.
7. T. Kawakatsu, Y. Kikuchi and M. Nakajima, J. Am. Oil Chem. Soc., 1997, 74, 317–321.
8. K. C. Van Dijke, C. Karin, P. G. H. Schroën and R. M. Boom, Langmuir, 2008, 24, 10107–10115.
9. X. Huang, M. Eggersdorfer, J. Wu, C. Zhao, Z. Xu, D. Chen and D. A. Weitz, RSC Adv., 2017, 7, 14932–14938.
10. S. Sugiura, M. Nakajima, S. Iwamoto and M. Seki, Langmuir, 2001, 17, 5562–5566.
11. A. Ofner, D. G. Moore, P. A. Rühs, P. Schwendimann, M. Eggersdorfer, E. Amstad, D. A. Weitz and A. R. Studart, Macromol. Chem. Phys., 2016, 218, 1600472.
12. I. Kobayashi, K. Uemura and M. Nakajima, Colloids Surf., A, 2007, 296, 285–289.
13. F. Malloggi, N. Pannacci, R. Attia, F. Monti, P. Mary, H. Willaime, P. Tabeling, B. Cabane and P. Poncet, Langmuir, 2010, 26, 2369–2373.
14. E. Amstad, M. Chemama, M. Eggersdorfer, L. R. Arriaga, M. P. Brenner and D. A. Weitz, Lab Chip, 2016, 16, 4163–4172.
15. Z. Li, A. M. Leshansky, S. Metais, L. M. Pismen and P. Tabeling, Lab Chip, 2015, 15, 3095.
16. F. Dutka, A. S. Opalski and P. Garstecki, Lab Chip, 2016, 16, 2044–2049.
17. I. Kobayashi, S. Mukataka and M. Nakajima, Langmuir, 2004, 20, 9868–9877.
18. I. Kobayashi, G. T. Vladisavljević, K. Uemura and M. Nakajima, Chem. Eng. Sci., 2011, 66, 5556–5565.
19. I. Chakraborty, J. Ricouvier, P. Yazhgur, P. Tabeling and A. M. Leshansky, Lab Chip, 2017, 17, 3609–3620.
20. W. C. Bauer, M. Fischlechner, C. Abell and W. T. S. Huck, Lab Chip, 2010, 10, 1814.
21. K. G. Lee, K. J. Park, S. Seok, S. Shin, D. H. Kim, J. Y. Park, Y. S. Heo, S. J. Lee and T. J. Lee, RSC Adv., 2014, 4, 32876–32880.
22. R. Dangla, S. C. Kayi and C. N. Baroud, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 853–858.
23. J. H. Xu, P. F. Dong, H. Zhao, C. P. Tostado and G. S. Luo, Langmuir, 2012, 28, 9250–9258.
24. J. Baret, Lab Chip, 2012, 12, 422–433.
25. M. C. Lee, H. Seonwoo, P. Garg, K. J. Jang, S. Pandey, H. B. Kim, S. B. Park, J. B. Ku, J. H. Kim and K. T. Lim, RSC Adv., 2018, 8, 6452–6459.
26. Y. Feng and Y. Lee, RSC Adv., 2017, 7, 49455–49462.
27. E. N. A. Latip, L. Coudron, M. B. McDonnell, I. D. Johnston, D. K. McCluskey, R. Day and M. C. Tracey, RSC Adv., 2017, 7, 49633–49648.

---